% ZonePlateDemo
%
% Little program to put up a zone plate image
% and then play with the monitor gamma.  When
% the display is not linearized, you can see
% distortion products.

% 5/24/97  dhb, eha  Wrote it.
% 6/4/97   dhb	More comments.
% 1/12/98	 dgp	Restrict array size, if necessary, for compatibility with Student Matlab.
% 4/13/02  awi  Changed "DeviceToSettings" to new name "PrimaryToSettings".

[c,maxElements]=computer;	% check for limitation of Student Matlab
maxSize=round(sqrt(maxElements));

% Print out some information.
fprintf('A zone plate is basically a radially symmetric\n');
fprintf('sweep frequency grating.  At each location, the\n');
fprintf('local spatial frequency and orientation matches\n');
fprintf('what you''d find in the same location in the fourier\n');
fprintf('plane.  So what you should see is concentric rings.\n');
fprintf('Any small circles are distortions of one sort or\n');
fprintf('another.  What Ted Adelson told me is that a zone plate\n');
fprintf('is a good test image because it is very sensitive to\n'); 
fprintf('linearization failures.  He also said he''s never seen\n');
fprintf('one on a monitor without distortion.  The appearance\n');
fprintf('is surprisingly insensitive to changes in gamma\n');
fprintf('between about 1.8 and 2.3, suggesting perhaps\n');
fprintf('that other monitor distortions dominate the gamma\n');
fprintf('correction after a certain degree of linearity has been\n');
fprintf('obtained.\n\n');

% Print out message about Eero's tools
fprintf('\nThis demo depends on routines in Eero Simoncelli''s\n');
fprintf('image pyramid tools.  The necessary routines are\n');
fprintf('included here with permission.  To get the whole toolbox\n');
fprintf('see Eero''s home page: http://www.cns.nyu.edu/~eero/\n');

% Set up some useful parameters
imageSize = 256;
if imageSize>maxSize
	imageSize=maxSize;
end
whichScreen = 0;

% Create a zone plate image in a matrix
fprintf(1,'\nComputing zone plate image ...');
theImage =round(255*((mkZonePlate(imageSize)+1)/2));
fprintf(1,' done\n');

% Load calibration information
if (exist('calZP') ~= 1)
	fprintf(1,'Loading calibration data ...');
	calZP = LoadCalFile(whichScreen);
	calZP = SetGamma(calZP,1);
	fprintf(' done\n');
end
% Use Psychophysics Toolbox to show the image uncalibrated
[theWindow,theRect] = Screen(whichScreen,'OpenWindow',0);
HideCursor;
Screen(theWindow,'SetClut',(0:255)'*ones(1,3));
Ask(theWindow,'Click mouse for uncalibrated zone plate',255,0);
Screen(theWindow,'PutImage',theImage);

% Now we'll load a calibration file and try to linearize
% the display
Ask(theWindow,'Click mouse for linearized zone plate',255,0);

linearValues = linspace(0,1,256)';
linearClut = PrimaryToSettings(calZP,[linearValues linearValues linearValues]')';
Screen(theWindow,'SetClut',linearClut);

% Now we'll linearize based on an exponent and cycle through a series
% of values.
Ask(theWindow,'Click mouse to see effect of varying gamma',255,0);
for gamma = 1:0.5:4
	offset = 0;
	oneClut = 255*((linearValues.^(1/gamma))*(1-offset) + offset);
	Screen(theWindow,'SetClut',[oneClut oneClut oneClut]);
	Ask(theWindow,sprintf('Corrected for gamma = %g, click to proceed',gamma),255,0);	
end

% Close up the window
ShowCursor;
Screen(theWindow,'Close');
