% PhaseDistortDemo
%
% Measures how much spatial phase information the you need
% to discriminate two images.
%
% When you run it, the command window will come forward and
% explain what is going on.
%
% No gamma correction is done on the images before display.

% 6/29/96  dhb,gmb,if  Wrote it.
% 3/4/97   dhb	Bells, whistles, and comments.
% 3/5/97   dhb	More enhancements in response to dgp comments.
% 4/11/97  dhb	Fixed small bug, pThreshold was undefined at time of first use.
% 3/12/98  dgp	Use Ask.
% 3/25/98  dgp	In line 80, delete illegal x argument to Ask.
% 7/9/98   dgp  Cope when GetClut is not available.
% 4/8/02   awi  -Passed window pointer to Screen 'GetClut' instead of screen number. 
%                Windows wants a screen number.
%               -Previosly we opened an onscreen window then positioned the Matlab command
%                window within view and printed instructions.  This doesn't work in Windows
%                because we can't control the command window So instead we close the 
%                onscreen window to display the instructions under Windows.  

% Load in the images.  Convention assumes that images are in
% MAT files with image name equal to file name.  Images are
% also assumed scaled in range 0-255.
LoadTheImagesPD; % and compute magnitude & phase

% Seed random number generator
ClockRandSeed;

% Open up the experimental screen
whichScreen = 0;
[window,screenRect] = Screen(whichScreen,'OpenWindow',0,[]);
offClut = zeros(256,3);
onClut = (0:255)'*[1 1 1];
matlabClut = Screen(window,'GetClut');
if isempty(matlabClut)
	matlabClut=ClutDefault(window);
end
Screen(window,'SetClut',offClut);

% Set up text style for user interface.
Screen(window,'TextFont','Helvetica');
Screen(window,'TextSize',20);
Screen(window,'TextStyle',0);
Screen(window,'TextMode','srcOr');

if strcmp('PCWIN',computer)
    screen(window,'Close')
else
	% Bring up command window and move it to upper left.
	% That way we can print lots of stuff in the command window
	% so the user will know what is going on.
	[cmdDx,cmdDy,bufferWindow] = CmdWinToUpperLeft(window,matlabClut);
end
% Provide introductory remarks.
fprintf('\n**************************************\n\n');
fprintf('When you click the mouse, four images will be\n');
fprintf('displayed.  The top two are images of Reagan and\n');
fprintf('Einstein.  The bottom two are the same images with\n');
fprintf('their phase spectra reversed.  The image under Reagan\n');
fprintf('has the phase spectrum of Reagan with the amplitude\n');
fprintf('spectrum of Einstein and vice-versa\n\n');
fprintf('This is a well-known demonstration of the importance of\n');
fprintf('the importance of spatial phase information.  See for\n');
fprintf('example: Badcock, D. R. (1990) "Phase- or energy-based face\n');
fprintf('discrimination: some problems." JEP:HPP, 16 (1), 217-20.\n\n');
fprintf('Click the mouse now to proceed\n');
GetClicks;

% Bring up screen window to show images.
if strcmp('PCWIN',computer)
    [window,screenRect] = Screen(whichScreen,'OpenWindow',0,[]);
    Screen(window,'SetClut',offClut);
else
	SCREENWinToFront(window,bufferWindow,cmdDx,cmdDy,offClut);
end

HideCursor;

% Show the images
[n,m] = size(image1);
theShift = n/2+20;
imageRect = SetRect(0,0,n,m);
centerRect = CenterRect(imageRect,screenRect);
theRect1 = OffsetRect(centerRect,-theShift,-theShift);
theRect2 = OffsetRect(centerRect,theShift,-theShift);
theRect3 = OffsetRect(centerRect,-theShift,theShift);
theRect4 = OffsetRect(centerRect,theShift,theShift);
Screen(window,'PutImage',255*Scale(image1),theRect1);
Screen(window,'PutImage',255*Scale(image2),theRect2);
Screen(window,'PutImage',...
	255*Scale(real(ifft2(image2Amp.*exp(sqrt(-1)*image1Phase)))),theRect3);
Screen(window,'PutImage',...
	255*Scale(real(ifft2(image1Amp.*exp(sqrt(-1)*image2Phase)))),theRect4);
Screen(window,'SetClut',onClut);
Ask(window,'Click mouse to proceed',255,0);
Screen(window,'FillRect',0);


% Set up QUEST Parameters
grainFactor = 20;
tGuess= 0.35;
beta=3.5;
delta=0.01;
gamma=0.5;
pThreshold=0.82;
tGuessSd=4.0;
trialsDesired=20;
q=QuestCreate(tGuess*grainFactor,tGuessSd*grainFactor,pThreshold,beta,delta,gamma);
wrongRight=str2mat('wrong','right');

ShowCursor;
if strcmp('PCWIN',computer)
    screen(window,'Close')
else
    % Bring command window to the front again.
    CmdWInToUpperLeft(window,matlabClut);
end


% Print out a message for the next stage
% Provide introductory remarks.
fprintf('\n**************************************\n\n');
fprintf('Next you will run an experiment to see how much phase information\n');
fprintf('is required to discriminate between Reagan and Einstein.\n');
fprintf('On each trial of the experiment, you will see two images.\n');
fprintf('Your task is to indicate which is Reagan.  To make the task\n');
fprintf('hard, the images you see will have identical amplitude spectra\n');
fprintf('(the mean of the two original images) and their phase spectra\n');
fprintf('will be titrated (via a parameter p) with a random phase image.\n');
fprintf('When p is 0, the phase of both images becomes completely\n');
fprintf('random.  When p is 1, each image has its proper phase.\n');
fprintf('The Quest algorithm will be used to find your %g%% correct\n',100*pThreshold);
fprintf('threshold on parameter p for discriminating Reagan from Einstein.\n');
fprintf('The program is set to run %g trials.\n\n',trialsDesired);
fprintf('On each trial, use the cursor to click on Reagan or press\n');
fprintf('the key ''q'' to quit early.  A single beep indicates\n');
fprintf('correct, two beeps indicates incorrect.\n\n');
fprintf('Click the mouse now to proceed.\n');
GetClicks;

% Bring up screen window to show images.
if strcmp('PCWIN',computer)
    [window,screenRect] = Screen(whichScreen,'OpenWindow',0,[]);
    Screen(window,'SetClut',offClut);
else
    % Bring up screen window to show images.
    SCREENWinToFront(window,bufferWindow,cmdDx,cmdDy,offClut);
end
HideCursor;

% Set up locations for FC trials.
[n,m] = size(image1);
imageRect = SetRect(0,0,n,m);
centerRect = CenterRect(imageRect,screenRect);
theRect1 = OffsetRect(centerRect,-theShift,0);
theRect2 = OffsetRect(centerRect,theShift,0);
homeX = (centerRect(RectLeft)+centerRect(rectRight))/2;
homeY = (centerRect(RectTop)+centerRect(rectBottom))/2;

% Now run the trials
Screen(window,'DrawText', ...
	['Click on ' targetName ', ''q'' to quit'], ...
	50,100,255);
respString = ' ';
theData = zeros(trialsDesired,4);
for k=1:trialsDesired
	% Get trial level and show the stimulus.  Force it into range
	pTest=QuestQuantile(q)/grainFactor;
	ShowTheStimulusPD;
	
	% Get the response.  There must be a slicker way to code this.
	GetTheResponsePD;
	
	% Update QUEST and keep data
	q=QuestUpdate(q,pTest*grainFactor,correct);
	pEst = QuestMean(q)/grainFactor;
	sdEst =QuestSd(q)/grainFactor;
	theData(k,1) = pTest;
	theData(k,2) = correct;
	theData(k,3) = pEst;
	theData(k,4) = sdEst;
end

% Close up the screen and put command window back where it was.
Screen(window,'SetClut',offClut);
Screen(window,'FillRect',0);
Screen('CloseAll');

% Get the final threshold
pEst = QuestMean(q)/grainFactor;
sdEst = QuestSd(q)/grainFactor;

% Some more printout
fprintf('\n**************************************\n\n');
fprintf('Threshold estimate is %4.2f ± %.2f\n\n',pEst,sdEst);
fprintf('Array theData contains the data.  It has four\n');
fprintf('columns.  The first column provides for each trial\n');
fprintf('the p value tested.  The second column tells\n');
fprintf('whether the trial was right (1) or wrong (0)\n');
fprintf('The third column gives Quest''s running threshold estimate.\n');
fprintf('The fourth column gives Quest''s running estimate standard\n');
fprintf('deviation.\n\n');
fprintf('The plot shows how the threshold estimate\n');
fprintf('converged over trials.\n\n');
fprintf('DHB''s threshold in this task is about 0.25.\n\n');

% Plot the threshold estimates.
plot(theData(:,3),'g');
xlabel('Trial number');
ylabel('Value of p');
title('Threshold estimate vs. Trial Number');

