function DatapixxGammaCorrection(varargin)
% DatapixxGammaCorrection([screenNumber])
%
% Doing the gamma correction for a CRT monitor, using Datapixx
%
% Sept 3, 2010 Jenny Read, heavily based on code from Peter April of VPixx

nbits = 16; % number of bits you have - 16 for DATAPixx, 8 normally
nreps = 200; % number of repetitions to do
timetoleavelab  = 10; % time to allow user to turn off lights, get out of lab and close doors etc.
savefile = 'testfinevals.mat';

gamma = 2.267728989547310; %

t0=GetSecs;
try	
	AssertOpenGL;
	
	% Configure PsychToolbox imaging pipeline to use 32-bit floating point numbers.
	% Our pipeline will lumalso implement an inverse gamma mapping to correct for display gamma.
	PsychImaging('PrepareConfiguration');
	PsychImaging('AddTask', 'General', 'FloatingPoint32Bit');
	PsychImaging('AddTask', 'General', 'EnableDataPixxM16Output');
	PsychImaging('AddTask', 'FinalFormatting', 'DisplayColorCorrection', 'SimpleGamma');
	
	% Open our window.
	% If input argument is supplied, that is the screennumber to use:
	if nargin>0
		screenNumber=varargin{1};
	else
		% Assume that the DATAPixx is connected to the highest number screen.
		screenNumber=max(Screen('Screens'));
	end
	
	%NB You may need this if you have not used the Pipeline before:
	%BitsPlusImagingPipelineTest(screenNumber)length(relpixvals)
	
	oldVerbosity = Screen('Preference', 'Verbosity', 1);   % Don't log the GL stuff
	[win, winRect] = PsychImaging('OpenWindow', screenNumber);
	Screen('Preference', 'Verbosity', oldVerbosity);
	winWidth = RectWidth(winRect);
	winHeight = RectHeight(winRect);
	
    % Specify the window's inverse gamma value to be applied in the imaging pipeline
    PsychColorCorrection('SetEncodingGamma', win, 1/gamma);
	
       
    % Open a serial port to the photometer: - NB dev/ttyUSB1 is  the USB port on
	% the front of the PC
	port = serial('/dev/ttyUSB1', 'BaudRate', 4800,'databits',7,'parity','even','stopbits',2,'terminator','CR/LF','FlowControl','hardware');
    fopen(port);
	whitelevel = 2^nbits-1;
	
    pixtest = (round(whitelevel/2) + [0 1 2 3 4 5])/whitelevel;
    
    ntest = length(pixtest);
	unitfield = ones(winHeight,winWidth);
	
	WaitSecs('UntilTime', t0+timetoleavelab);
	
    for jrep=1:nreps
        
		% Do it in random order
        for j=randperm(ntest)
            fprintf('measuring luminance for value %5.3f',pixtest(j));
            lumvals(j,jrep) = DisplayGanzfeld(pixtest(j),unitfield,win,port);
            fprintf(', got %5.3f cd/m^2.\n',lumvals(j));
            save(savefile)
		end

    end
    
    meanvals = mean(lumvals,2);
    SEMvals = std(lumvals,[],2)/nreps;
	
    figure
    subplot(1,2,1)
    plot(pixtest,lumvals,'bx')
    hold on
    errorbar(pixtest,meanvals,SEMvals,'ro')
    xlabel('pixel value (0=black, 1=white)')
    ylabel('measured luminance (cd/m^2)')
    axis tight
    xl=xlim;
    xlim([xl(1)*0.99999  xl(2)*1.00001])

    a1=subplot(1,2,2);
    errorbar(pixtest,meanvals,SEMvals,'ro')
    xlabel('pixel value (0=black, 1=white)')
    axis tight
    xl=xlim;
    xlim([xl(1)*0.99999  xl(2)*1.00001])
    box off
    % Do a second axes to mark unnorm pix vals:
    a2=axes;
    set(a2,'pos',get(a1,'pos'),'xlim',get(a1,'xlim'),'ylim',get(a1,'ylim'),'color','none','xaxisloc','top');
    xlim(xlim*whitelevel)
    val2label = pixtest(1:2:end)*whitelevel;
    for kkk=1:length(val2label)
        label{kkk} = num2str(val2label(kkk),'%d');
    end
    set(a2,'xtick',val2label,'xticklabel',label)
    box off
    
    
    
    figure
    subplot(1,2,1)
    plot(pixtest*whitelevel,lumvals,'bx')
    hold on
    errorbar(pixtest*whitelevel,meanvals,SEMvals,'ro')
    xlabel('pixel value (0=black, 65535=white)')
    ylabel('measured luminance (cd/m^2)')
    axis tight
    for kkk=1:length(pixtest)
        label{kkk} = num2str(pixtest(kkk)*whitelevel,'%d');
    end
    set(gca,'xtick',pixtest*whitelevel,'xticklabel',label)
    title('individual measurements')
    a1=subplot(1,2,2);
    errorbar(pixtest*whitelevel,meanvals,SEMvals,'ro')
    xlabel('pixel value (0=black, 65535=white)')
    axis tight    
    set(gca,'xtick',pixtest*whitelevel,'xticklabel',label)
    title('means +/- SEM of 200 repeats')
    % Restore any system gamma tables we modified
	Screen('CloseAll');nreps
	
	% Close port to photometer:
	fclose(port);delete(port);clear port
	
catch
	lasterr
 	RestoreCluts;
	fclose(port);delete(port);nrepsclear port
	sca
	clear all	

end

return;

function err = rmserr(params,relpixvals,lumvals)
predlumvals = powerlaw(relpixvals,params);
err = sqrt(mean( (predlumvals-lumvals).^2));
function y = powerlaw(x,params)
power = params(1) ;
scale = params(2);
offset = params(3);
y = offset + (scale*x).^power;

function lumval = DisplayGanzfeld(relpixval,unitfield,win,port)
% Relpixval is a number between 0 (black) and 1 (white)

% Fill the whole screen with uniform luminance:
ganzfeld = relpixval*unitfield;
% Use floatprecision=2:
tex = Screen('MakeTexture', win, ganzfeld, [], [], 2);

% Draw the floating point texture
% Specify filter mode = 0 (nearest neighbour), so that GL doesn't interpolate pixel values.
Screen('DrawTexture', win, tex, [], [], [], 0);
Screen('Flip',win,GetSecs);
Screen('Close',tex);

% Read photometer:
fprintf(port,'MES');
[line, nlin] = fscanf(port);
j=min(findstr(line,' '));
val = str2num(line(j+1:end));
if isempty(val)
	val = 0;
end
lumval = val;


