% VisualGammaDemo.m
%  
% VisualGammaDemo quickly calibrates a video monitor's gamma function,
% dots/inch, and frame rate without requiring any extra equipment. The
% purpose is to create a minimal set of calibration data for the
% monitor to allow its use in most experiments and demonstrations,
% with reasonable, if not ultimate, accuracy in synthesis of stimuli.
% 
% VisualGammaDemo collects its results in a structure called "cal". 
% It optionally saves the data in a calibration file stored in the
% CalData subfolder, with a file name, e.g. screen0.mat, indicating by
% the screen number.  If there already was calibration information
% present, you have the option of replacing it entirely or updating
% the parts that VisualGamma knows about.  Replacing is dangerous
% because you will lose any information that VisualGamma doesn't know
% about.
% 
% The gamma function is measured by asking the viewer to adjust a
% luminance to visually null a grating's contrast to zero. The
% technique is precise (contrast threshold for the test pattern is
% 0.003) and accurate, since the light and dark bars of the grating
% will be physically identical in observer's retinal image at the
% null. This is what Brindley calls a Class A match. The grating is
% vertical. The even bars of the grating are uniform and adjustable.
% The odd bars, not adjustable, are a mixture of two colors (e.g.
% black and white) produced by having alternate scan lines of each
% color. The viewing distance is made sufficiently large that the
% spatial frequency of the scan line alternation cannot pass through
% the observer's pupil, so the retinal image of the bar is a
% homogeneous optical mixture. Matching a 50/50 mixture of black and
% white bisects the interval. Repeated bisections yield the gamma
% function describing (except for a linear transformation) the
% luminance as a function of voltage applied to the monitor.
% 
% The gamma function y=f(x) is fitted with a rectified power law: 
%   if x<=x0 then y=0
%   if x>x0 then y=((x-x0)/(1-x0))^gamma
% This fit has two parameters: the threshold x0 and the exponent
% gamma. Since this function usually provides an excellent fit, and
% only has two degrees of freedom, the 8 measurements are more than
% enough to accurately characterize the entire function. Knowing this
% gamma function, we can linearize our display, e.g. produce truly
% sinusoidal luminance waveforms, but we won't know their absolute
% luminance or contrast.
% 
% The value y is related to the display luminance L by the linear
% transformation y=(L-LMin)/(LMax-LMin). The range of y is [0,1],
% corresponding to the luminance range [LMin,LMax]. For many
% experiments, LMax is relatively unimportant, e.g. it would be enough
% to know, by appearance, that it's low to mid photopic, e.g. 30 to
% 300 cd/m^2. However, we usually do need to know absolute contrast.
% That requires that we measure the ratio LMin/LMax.
% 
% I tried asking the observer to adjust the contrast of a
% high-luminance grating to match the contrast of a low-luminance
% grating. Alas, this is a subjective, Class B match, and people are
% famously good at discounting mean luminance. Matching contrast
% across luminance turns out not to be a good task.
% 
% A better approach would be to ask the observer to supply a neutral
% density filter. We could calibrate the filter and the LMin/LMax
% ratio, by doing two luminance matches.
% 
% PLANS:
% 
% Here's our current model for monitor calibration:
%   1. 3 dacs may be linearly combined by a video attenuator,  
%      characterized by 3 gains
%   2. gamma function and LMin and LMax for each of 3 guns
%   3. 3 phosphors each have an emission spectrum. (Emissions
%      are additive.)
% 
% At present VisualGammaDemo only handles item 2, and in fact only
% measures a combined gamma, driving all three guns equally. It would
% be easy to have VisualGammaDemo measure the 3 guns individually.
% However, if the 3 guns are physically identical then their gamma
% functions will be the same, because light output is proportional to
% beam current for all phosphors. So measuring the 3 guns separately
% is only needed if the guns differ. I suspect the differences are
% very small.
% 
% For many purposes, items 1 and 3 can be estimated from published
% data, without any new measurements. The video attenuator resistors
% have +/-1% tolerance, but the gains of the 3 dacs aren't guaranteed
% to match to quite as high a tolerance. Similarly, for many purposes,
% one might reasonably assume standard emission spectra, based on
% manufacturer specs.
% 
% A subtle point is that the gamma function for a small uniform patch
% depends somewhat on the mean luminance of the rest of the display.
% Vision experiments typically maintain a constant mean luminance on
% every frame, but that mean luminance may vary between experiments.
% My past practice has been to do the calibration and experiment at
% the same mean luminance. However, since the effect is small, it
% seems likely that one could measure the gamma function at several
% mean luminances and produce a parametric description for how the
% gamma function's parameters depend on mean luminance. Then a single
% calibration could characterize performance at all luminances.
% However, it will mean that the gamma function will be a function of
% two variables, not one, and thus may not be amenable to tabulation.
% Or, at least, we will need multiple tables and a method of
% interpolation.
% 
% A closely related point is that it might be a good idea to make some
% basic measurements designed to characterize the worst limitations of
% the display, as suggested in my paper on pixel independence (Pelli,
% 1997). Even though these numbers wouldn't be routinely used,
% measuring them would draw everyone's attention to the limitations,
% so that they would be borne in mind by users in designing their
% experiments, and would be mentioned more often to monitor
% manufacturers, who might then improve their designs. What I have in
% mind are all the effects that deviate from the model of independent
% pixels, transformed by a point nonlinearity (the gamma function),
% and blurred by the point spread function. There are two effects in
% particular that seem worth measuring:
%   1. the effects of slew rate limitations on mean luminance of fine
%   vertical gratings,
%   2. the effect of the luminance of the rest of the screen on the 
%   luminance of a patch.
% The slew rate limitation probably could be modeled as a horizontal
% blur that occurs BEFORE the gamma function, but that's too hard to
% think about. The monitor manufacturer can eliminate this problem by
% increasing the slew rate, i.e. providing a faster video amplifier.
% There are 3 distinct causes for the luminance effect: incomplete DC
% restoration (let's keep asking for dc coupling), imperfect high
% voltage regulation (on cheaper monitors), and a new monitor
% "feature" that stabilizes the mean luminance (Tom Robson told me
% about this last one).
% 
% Pelli, D.G. (1997) "Pixel independence: Measuring spatial
% interactions on a CRT display." Spatial Vision 10(4):443-446
% web http://vision.nyu.edu/VideoToolbox/PixelIndependence.html
% 
% It would be nice to print a fine grating on transparency that could
% be used as a spatial-frequency filter, to block the mixture
% frequency at a convenient short viewing distance.
% 
% BUYING A LONG KEYBOARD CABLE:
% 
% For long viewing distances it's nice to be able to move the keyboard
% and mouse far from the computer. Until 1998, Apple keyboards and
% mice were connected to the computer by an ADB [Apple Desktop Bus]
% cable. Since then they use a USB cable. (In September, 2000, Apple
% is switching to an optical mouse, which hopefully will allow long
% working distances.) S-VHS cables exactly match the ADB ports and are
% available in many lengths (e.g. 20 feet) from local electronics
% companies, and Markertek, in Saugerties, NY at 800-522-2025,
% fax 914-246-1757.
% web http://vision.nyu.edu/Tips/HowTo.html#remoteKeyboard
% 
% See MonitorGamma, MonitorGammaError, LoadCalFile, SaveCalFile, ClutMovieDemo, RenderDemo.
%
% Denis Pelli 6/1/97

% 5/19/96 dgp Wrote it.
% 5/21/96 dgp Use the "input" command to get viewing distance, (thanks david!)
% 5/26/96 dgp Always use finest grating, which restricts us to bisection, but
%			  supplement original 2 colors by colors matching mixtures.
% 5/27/96 dgp Fit gamma by rectified power law.
% 5/28/96 dgp Added contrast matching to estimate LMinOverLMax.
% 5/28/96 dgp Updated to use new GetMouse, that flushes mouse events
% 6/15/96 dgp Added note about long keyboard cables, above
% 3/20/97 dgp Updated
% 4/2/97  dgp Updated. Use LoadCalFile instead of LoadCal.
% 4/26/97 dhb Cosmetic. Changed name.  Add warning before overwrite of cal file.
% 6/1/97	dgp	Use FIGURE(GCF) to make sure we get to see it.
%							Disabled contrast matching, since we're no good at it.
%							Polished Figure 1.
% 8/16/97 dgp Added "clear text" to prevent any such variable from displacing TEXT function.
% 4/10/98 dhb Updated for new calibration fields, compatibility with calibration files built other ways.
% 5/16/98 dhb Modifed to allow other gamma fitting methods.
%         dhb Allow 6 or 9 measurements, 9 makes more sense for LCD based devices.
%         dhb Fixed little bugs in saving code, probably these were introduced 4/10/98.
% 7/14/99 dgp At least some parts of VisualGammaDemo assume 8-bit pixelsize, so i forced it to
%             be 8-bits in the Screen('OpenWindow') in GratingNull.m. This doesn't require
%             any alert since it's only temporary.
% 7/24/00 dgp Fixed typo reported by Bosco Tjan, in line 382 changed "cal.dip" to "cal.dpi".
%             Minor update of the help text above.
% 4/06/02 dhb changed reference from cal.who to cal.describe.who
% 4/06/02 awi -Added windows conditional to use Screen 'FrameRate' instead of Screen 'PeekBlanking'
%             to determine the frame rate.  'PeekBlanking' is not yet implemented on Windows.
%             -Turned off warning messages during call to fmins.  Matlab 6 implements implements fmins
%              but reports that it will be dropped from future versions.  
% 4/11/02 awi added comment for dhb on 4/06/02
% 4/13/02 dgp Use FrameRate() function always. Restore old warning state after we suppress warnings.

% Some parameters
theScreen=0;

% Set number of measurements
% Because of the recursive way the
% routine is coded, this can currently
% be either 6 or 9, but nothing else.
%
% It might make more sense to rethink the way
% these values are chosen and thus get more
% flexibility.
nGammaMeasure = 6;

% Specify how gamma function will be fit.
% A value of zero uses classic gamma function.
% Other values are passed on to routine FitGamma.
% Eventually, FitGamma should be modified to 
% incorporate the standard gamma function, but
% I'm in a rush right now.
%
% The problem with a classic gamma function is
% that it is all wrong for an LCD projector.
%
% Use a value of 7 for cubic spline.
gammaFitMethod = 0;

cal=LoadCalFile(theScreen);
if ~isempty(cal)
	fprintf('A calibration file for screen %d already exists:\n',theScreen);
	fprintf('%s measured on %s,\nMonitor %s, using program %s.\n',...
		cal.describe.who,cal.describe.date,cal.describe.monitor,cal.describe.program);
	if (length(cal.describe.program) == length('VisualGammaDemo') & all(cal.describe.program == 'VisualGammaDemo'))
		fprintf('Monitor was displaying %dx%d pixels at %.1f Hz.\n',cal.width,cal.height,cal.hz);
	end
	fprintf('Calibration description: %s\n',cal.describe.comment);
	action=input('Update it (0), Replace it (1), or Cancel (2)? ');
	if (action == 2)
		return;
	end
	if (action == 1)
		fprintf('Are you sure you want to replace?  You will lose\n');
		fprintf('any calibration information that VisualGammaDemo\n');
		fprintf('does not meausure (e.g. phosphor spectra)\n');
		action=input('Update  (0), Replace  (1), or Cancel? ');
		if (action == 2)
			return;
		elseif (action == 1)
			cal=[];
		end
	end
end

if isfield(cal,'monitorGammaGamma')
	fprintf('Old monitorGammaGamma=%.2f\n',cal.monitorGammaGamma);
	fprintf('Old monitorGammaThreshold=%.2f\n',cal.monitorGammaThreshold);
end
if input('Do a nulling task to measure the monitor''s gamma function? yes(1) or no(0): ')
	% Calculate the minimum viewing distance.
	% We want the blend s.f. to be >= 60 c/deg,
	% which cannot pass through a 2 mm pupil.
	theScreen=0;
	dpi=67;
	blendPeriodPix=2;
	distanceM=0.0254*blendPeriodPix*57*dpi/60;
	fprintf('We suggest a viewing distance of at least %.1f m.\n',distanceM);
	fprintf('Measurements at those distances will be pure luminance measurements.\n');
	fprintf('(60 c/deg can''t pass through a 2 mm pupil.) At nearer distances the\n');
	fprintf('spatial frequency would be low enough to reach the retina and\n');
	fprintf('contaminate the optical calibration by visual nonlinearities in\n');
	fprintf('luminance perception. Your mouse cable probably won''t reach that far.\n');
	fprintf('So get a friend to move the mouse for you. Or set up a mirror, at least\n');
	fprintf('%.1f m away, and look at the display''s reflection in the mirror.\n',distanceM/2);
	fprintf('\n');
	distanceM=input('Please enter the actual viewing distance, in meters: ');
	fprintf('\n');
	% fraction is fraction of foreColor in foreColor-backColor mixture
	% weight is scalar multiplier of foreColor for visual match
	% Use each new match color in subsequent mixtures.
	% We are given 2 colors, usually black and white. 
	% Our first iteration bisects the interval between
	% them, by matching to a 50/50 mixture. This yields 3 colors, which we sort. 
	% In each subsequent iteration we bisect every adjoining pair by matching to 50/50 
	% mixture, and then sort again. Each iteration halves the interval size. Doing 2 
	% iterations yields intervals of 1/4, requiring 1+2=3 measurements. We then go on to 
	% make the first measurement of the 3rd, 4th, and 5th iterations, since the lowest 
	% luminances help most in determing the gamma function's threshold. This is a
	% grand total of 6 measurements.
	%
	% 5/16/98 -- DHB.  If nGammaMeasure is set to 9, then the first only break
	% described above happens on the 4th iteration rather than the third, leading
	% to nine measurements in all.  This is a rather kludgy way to adjust the
	% density of the measurements.  Perhaps this code should be rethought to
	% allow more flexibility.
	white=[255 255 255];
	black=[0 0 0];
	foreColor=white;
	backColor=black;
	fraction=[0 1];
	weight=[0 1];
	trials = nGammaMeasure;
	if (nGammaMeasure == 6)
		firstOnlyBreak = 3;
	elseif (nGammaMeasure == 9)
		firstOnlyBreak = 4;
	else
		error('Parameter nGammaMeasure must be either 6 or 9');
	end
	trial=1;
	for j=1:5
		n=length(fraction);
		i=n+1;
		for f=1:n-1
			message=sprintf('%d of %d',trial,trials);
			fraction(i)=(fraction(f)+fraction(f+1))/2;
			color1=weight(f+1)*foreColor+(1-weight(f+1))*backColor;
			color2=weight(f)*foreColor+(1-weight(f))*backColor;
			w=GratingNull(theScreen,message,1/2,2,color1,color2,distanceM,dpi);
			trial=trial+1;
			weight(i)=w*weight(f+1)+(1-w)*weight(f);
			i=i+1;
			if j>=firstOnlyBreak;
				break
			end
		end
		[fraction index]=sort(fraction);
		weight(:)=weight(index);
	end
	
	% Fit and plot gamma function.
	fprintf('Figure No. 1 will now show your monitor''s gamma function.\n');
	if (gammaFitMethod == 0)
		oldWarning=warning;
    warning off
			params=fmins('MonitorGammaError',[.1 2],[],[],weight,fraction);
    warning(oldWarning); 
		threshold=params(1);
		gamma=params(2);
		rmsError=MonitorGammaError(params,weight,fraction);
		fprintf('Monitor gamma=%.1f threshold=%.1f rms error=%.4f\n', ...
			gamma,threshold,rmsError);
		x=0:.01:1;
		lines=plot(x,MonitorGamma(x,threshold,gamma),'-',weight,fraction,'*');
		cal.monitorGammaGamma = gamma;
		cal.monitorGammaThreshold = threshold;
		cal.gammaRmsError = sqrt(MonitorGammaError(params,weight,fraction));
		cal.gammaInput = (0:255)';
		cal.gammaTable = MonitorGamma((0:255)'/255,threshold,gamma)*ones(1,3);
	else
		cal.gammaInput = (0:255)';
		cal.gammaTable = FitGamma(255*weight',fraction',(0:255)',gammaFitMethod)...
			*ones(1,3);
		lines=plot(cal.gammaInput/255,cal.gammaTable(:,1),'-',weight',fraction','*');
	end
	for i=1:length(lines)
		set(lines(i),'LineWidth',2);
		set(lines(i),'MarkerSize',12);
	end
	axis([0 1 0 1])
	axis('square')
	set(gca,'LineWidth',1);
	set(gca,'TickLength',2*get(gca,'TickLength'));
	set(gca,'FontSize',14);
	title(sprintf('Gamma function of screen %d',theScreen),'FontSize',18);
	xlabel('DAC voltage x (normalized)','FontSize',24);
	ylabel('Luminance y (normalized)','FontSize',24);
	if (gammaFitMethod == 0)
		g=sprintf('%.1f',gamma);
		if abs(threshold)>0.005
			equation=sprintf('y=((x%+.2f)/(1%+.2f))^%s^%s^%s ± %.3f',-threshold,-threshold,g(1),g(2),g(3),rmsError);
		else
			equation=sprintf('y=x^%s^%s^%s ± %.3f',g(1),g(2),g(3),rmsError);
		end
		clear text	% we want the built-in function
		text(.05,.9,equation,'FontSize',20);
	end
	drawnow
	figure(gcf); % force figure to the front, so we get to see it.
	
	if (isfield(cal,'rawGammaInput'))
		rmfield(cal,'rawGammaInput');
		rmfield(cal,'rawGammaTable');
	end
end

%calibrate DPI
if isfield(cal,'dpi')
	fprintf('Old dpi=%.1f\n',cal.dpi);
end
if input('Would you like to measure your display''s dots per inch? yes(1) or no(0): ');
	dpi=MeasureDpi(theScreen);
	cal.dpi=dpi;
end

doIt = input('Save the calibration? yes(1) or no(0): ');
if (doIt)
	cal.program = 'VisualGammaDemo';
	cal.date = date;
	cal.gammaMethod = 'visual matching';
	whoString = input('Enter your name, in quotes '''': ');
	cal.who = whoString;
	monitor = input('Enter monitor brand and model, in quotes '''': ');
	cal.monitor = monitor;
	cal.whichScreen = theScreen;
	
	% Get screenRect and measure the frame rate
	screenRect=Screen(theScreen,'Rect');
  hz = FrameRate(theScreen);
	cal.width = RectWidth(screenRect);
	cal.height = RectHeight(screenRect);
	cal.hz = hz;
	fprintf('%s measured gamma on %s by %s.\n',cal.who,cal.date,cal.gammaMethod);
	fprintf('%s, displaying %dx%d pixels at %.1f Hz.',cal.monitor,cal.width,cal.height,cal.hz);
	if isfield(cal,'dpi')
		width=cal.width/cal.dpi;
		height=cal.height/cal.dpi;
		fprintf('%.1f inch display. %.0f dots/inch.\n',sqrt(width^2+height^2),cal.dpi);
	else
		fprintf('\n');
	end
	fprintf('You should further describe the calibration, e.g. ''room lights on''\n');
	description=input('Please describe, in quotes '''': ');
	cal.description = description;
	
	% Save the calibration.	
	% Remove dynamic stuff that should not be in the calibration file
	% but sometimes gets there.
	fprintf('Now saving ''%sscreen%d.mat''\n',CalDataFolder,theScreen);
	if isfield(cal,'gammaMode')
		rmfield(cal,'gammaMode');
	end
	if isfield(cal,'iGammaTable')
		rmfield(cal,'iGammaTable');
	end
	SaveCalFile(cal,theScreen);
end
