function STIM = createDotMotionStimulus(win,coherence,direction,numFrames,varargin)
% STIM = createDotMotionStimulus(win,coherence,direction,numFrames,[options])
%
% createDotMotionStimulus: creates a dot-motion stimulus (single trial) on a series of off-screen 
% windows, returning an array of pointers to the windows (STIM). Each entry in STIM is a pointer to
% one frame of the trial.
% inputs:
%
%	win: pointer to onscreen window being used to display stimuli
%
%	coherence: level of coherence for motion in the trial, as a percent. Ranges from 0 
%		(totally random motion) to 100 (all dots moving together).
%
%	direction: if 0, coherent motion will be to the left. Otherwise, motion will be to the right.
%	Other directions are possible using option flags as described below.
%
%	numFrames: number of frames in the trial
%
%	noFixation: Don't add a fixation cross (Default is to add one).
%
% options are: 
%		[comment this in]
%
% created 1/03 by ras (for fMR-A Expt 3)
% 1/31/03 ras: updated to be very flexible, with varargin options incl. setting dot size,
% velocity, stim window size, etc.

%%%%% parameters / defaults %%%%%%%%%%%%%%%%%%%%%%%%
stimSize = 240;
addFixation = 1;
dotSz = 1;
numDots = 500;
velocity = 5; 	% pixels per frame, rate of motion for coherent dots
dotColor = [253 254];	% pointer in screen's lookup table for dot color
dirX = 1;		% defaults to dirX is 1 only -- left/right motion
dirY = 0; dirR = 0; dirT = 0;
radialFlag = 0;
backgndColor = 255;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% a line added for back compatibility with some of my old stuff
if direction==0		dirX = -1;		end

%%%%% parse the option flags
for i = 1:length(varargin)
	switch varargin{i}
	case 'noFixation'
		addFixation = 0;
	case 'left'
		dirX = -1;
	case 'right'
		dirX = +1;
	case 'up'
		dirY = -1; dirX = 0;
	case 'down'
		dirY = +1; dirX = 0;
	case 'cw'
		radialFlag = 1;
		dirY = -1;
	case 'ccw'
		radialFlag = 1;
		dirY = +1;
	case 'in'
		radialFlag = 1;
		dirX = -1;
	case 'out'
		radialFlag = 1;
		dirX = +1;
	case 'radial'
		radialFlag = 1;
	case 'stimSize'
		stimSize = varargin{i+1};
	case 'dotSize'
		dotSz = varargin{i+1};
	case 'velocity'
		velocity = varargin{i+1};
	case 'numDots'
		numDots = varargin{i+1};
	case 'dotColor'
		dotColor = varargin{i+1};
	case 'backgndColor'
		backgndColor = varargin{i+1};
	otherwise,
		if ischar(varargin{i})
			fprintf('\nWARNING: createDotMotionStimulus unrecognized option flag "%s".\n', varargin{i});
		end
	end
end

numCohDots = numDots * (coherence/100);
numNoiseDots = numDots - numCohDots;

radius = floor((stimSize+1)/2);
[X,Y] = meshgrid(1:stimSize);
dX = X-radius; dY = Y-radius;
if radialFlag
	R = sqrt(dX.*dX + dY.*dY);
	Theta = atan(dY./dX);
	% ugly empirical manipulations -- I need to brush up on my trig
	Theta = Theta + pi/2;
	Theta(1:stimSize,1:stimSize/2) = Theta(1:stimSize,1:stimSize/2) + pi;
end

rect = screen(win,'Rect');
szX = rect(3); szY = rect(4);
minX = round(szX/2 - stimSize);
maxX = round(szX/2 + stimSize);
minY = round(szY/2 - stimSize);
maxY = round(szY/2 + stimSize);
r = radius;
imgRect = [minX minY maxX maxY];

%%%%% Main loop 
% fprintf('Creating dot stimuli.');
rand('state',sum(100*clock));
numColors = length(dotColor);
numCohDots = numCohDots / numColors;
numNoiseDots = numNoiseDots / numColors;
for frame  = 1:numFrames
	STIM(frame) = screen(win,'OpenOffScreenWindow',backgndColor);
	currimg = zeros(stimSize,stimSize);
	for color = 1:numColors	

		%%%%% draw coherent dots on current image, taking into account
		%%%%% their location on the previous frame, direction and velocity
		if frame==1
			loc = round((stimSize^2)*rand(numCohDots,1));
		elseif radialFlag==0	% dots moving l/r/up/down
			loc = lastloc{color} + dirX*velocity*stimSize + dirY*velocity;
		else					% dots moving in/out/ccw/cw
			locR = R(lastloc{color}) + dirX*velocity;
			locTheta = Theta(lastloc{color}) + dirY*velocity;
			ind = find(locR > r | locR <= 0);
			while ~isempty(ind)
				locR(ind) = round((r-1)*rand(size(ind))) + 1;
				ind = find(locR > r | locR <= 0);
			end
			locX = locR .* cos(locTheta);
			locY = locR .* sin(locTheta);
			locX = round(locX + r);
			locY = round(locY + r);
			ind = find(locX > stimSize | locY > stimSize | locX <= 0 | locY <= 0);
			while ~isempty(ind)
				locX(ind) = round((stimSize-1)*rand(size(ind))) + 1;
				locY(ind) = round((stimSize-1)*rand(size(ind))) + 1;
				ind = find(locX > stimSize | locY > stimSize);
			end
		
			% let's also randomly remove a subset of the moving dots
			ind = randi(length(locX),[1 round(0.1*length(locX))]);
			while ~isempty(ind)
				locX(ind) = round((stimSize-1)*rand(size(ind))) + 1;
				locY(ind) = round((stimSize-1)*rand(size(ind))) + 1;
				ind = find(locX > stimSize | locY > stimSize);
			end
			
			loc = sub2ind([stimSize, stimSize],locX,locY);	
			ind2 = find(isnan(loc));
			loc(ind2) = round(stimSize*rand(size(ind2)));
		end

		% constrain to a circle of radius r
		ind = find( (mod(loc,stimSize)-r).^2 + (floor(loc/stimSize)-r).^2 > r^2);
		while ~isempty(ind)
			loc(ind) = round((stimSize^2) * rand(size(ind)));
			ind = find( (mod(loc,stimSize)-r).^2 + (floor(loc/stimSize)-r).^2 > r^2);
		end
	
		% record positions for next frame
		lastloc{color} = loc;

		% expand to dotSz 
		coh = [];
		for row = 1:dotSz
			for col = 1:dotSz
				coh = [coh; loc + row - stimSize*col];
			end
		end	
		coh = coh(find(coh > 0 & coh < stimSize^2));
		
% 		% convert to radial units if necssary
% 		if radialFlag
% 			locX = R(loc).*cos(Theta(loc));
% 			locY = R(loc).*sin(Theta(loc));
% 			locX = round(locX + r);
% 			locY = round(locY + r);
% 			loc = sub2ind([stimSize, stimSize],locX,locY);
% 		end		
	
		currimg(coh) = dotColor(color);
		
		%%%%% draw noise dots on current image
		locN = round((stimSize^2)*rand(numNoiseDots,1));
		
		% constrain to a circle of radius r
		ind = find( (mod(locN,stimSize)-r).^2 + (floor(locN/stimSize)-r).^2 > r^2);
		while ~isempty(ind)
			locN(ind) = round((stimSize^2) * rand(size(ind)));
			ind = find( (mod(locN,stimSize)-r).^2 + (floor(locN/stimSize)-r).^2 > r^2);
		end
	
		% expand to dotSz
		noise = [];
		for row = 1:dotSz
			for col = 1:dotSz
				noise = [noise; locN + row - stimSize*col];
			end
		end	
		noise = noise(find(noise > 0 & noise < stimSize^2));
		
		currimg(noise) = dotColor(color);
		
		%%%%% put up next frame of movie
		screen(STIM(frame),'PutImage',currimg,imgRect,'srcCopyQuickly');
		
		if addFixation
			PutFixation(STIM(frame),szX,szY);
		end
		
	% 	fprintf('.');
	end		
end
% fprintf('\n');

return



function PutFixation(win,X,Y);
%%%%% apply fixation point
SCREEN(win, 'FillRect', 255, centerrect([0 0 2 6], [0 0 X Y]));
SCREEN(win, 'FillRect', 255, centerrect([0 0 6 2], [0 0 X Y]));
return
