% MakeMovie will create a circle that is color-coded to show a probability 
% distribution.  It will animate possible k-vector transitions on that 
% circle.
%
% Usage:
%	cd <path to MakeMovie.m>
%	mkdir images
%	MakeMovie
%
% Revision History
%	4/29/01 - Initial version, Jason Harris

resolution=256;				% probability density resolution
radialResolution = 256;		% resolution around the rim of the circle
zones = 8;						% number of different probability 'zones'
frames = 30;					% number of transitions to animate
startAngle = 270;				% initial angle of k-vector
lineWidth=5;					% pixel width of the vector

% compute a random number for each probability 'zone'
a=rand(zones,1);
a=[a ; a(1)];

% for each 'zone', compute it's extent along the total probability 
% distribution
b=rand(zones,1);
b=b/sum(b);

% map each zone onto the total probability distribution 'c'.  The space 
% between the zone boundaries is linearly interpolated so that the 
% probability distribution is smooth.  To see the probability distribution, 
% do a 'plot(c)'.
c=[];
for i=1:zones;
	c=[c ; linspace(a(i), a(i+1), round(radialResolution*b(i)))'];
end;

% uncomment this for a very artificial probability distribution that's 
% useful for debugging.
%for i=1:radialResolution;
%   if (i>50 & i<70) | (i >200 & i<250);
%      c=[c ; 1];
%   else
%      c=[c ; 0];
%   end;
%end;

% normalize the probability distribution between zero and one and initialize 
% some convenience variables.
c = c-min(c);
c = c./max(c);
[maxC,arc] = size(c);
[colorMapSize,arc] = size(gray);

% Create a matrix 'circle' that contains an intensity value at every point.  
% This matrix is essentially a grayscale bitmap of the circle.  
circle=zeros(resolution*2+1, resolution*2+1);
for x=-resolution:resolution;
	for y=-resolution:resolution;
		phi = atan2(y,x);
		phi = phi + pi;
		d = sqrt(x^2 + y^2);
		if d <= resolution;
			intensityIndex = round(phi/(2*pi) * (maxC-1)) + 1;
			intensity = c(intensityIndex);
		else
			intensity = 0;
		end;
		intensity = intensity * colorMapSize;
		circle(y+resolution+1, x+resolution+1, 1) = intensity;
	end;
end;
 
% Draw the circle
colormap hot;
image(circle);
axis xy;
axis equal;
axis tight;
axis off;

% Create a vector arrow.  Doing 'patch(x,y,'g')' draws the arrow pointing to
% the left in green.
center = [resolution + 1, resolution + 1];
x=[	center(1)				...
      center(1)				...
	  center(1)-resolution+resolution/15	...
	  center(1)-resolution+resolution/10	...
	  center(1)-resolution ...
	  center(1)-resolution+resolution/10	...
	  center(1)-resolution+resolution/15	...
      center(1)				];
y=[	center(2)-lineWidth	...
      center(2)+lineWidth	...
	  center(2)+lineWidth	...
	  center(2)+lineWidth+resolution/10		...
	  center(2)				...
	  center(2)-lineWidth-resolution/10		...
      center(2)-lineWidth	...
      center(2)-lineWidth	];

% We need an inverse mapping of the probability distribution.  'd' contains the 
% indices of 'c', weighted by their probability.  For example, if c was zero 
% everywhere except between indices 50 and 52 and constant between those 
% indices, 'd' would contain a line from 50 to 52.  Having 'd' lets me 
% 'throw a dart' and hit regions of high probability.
c=c./sum(c);
d=[];
for i=1:maxC;
   for j=1:round(c(i)*radialResolution);
      d = [d ; i];
   end;
end;
[maxD,i]=size(d);

% create the animation.  For each frame, draw an arrow pointing towards 
% 'startAngle'.  Then, compute a probabilistic final orientation and 
% point the arrow in that direction.
%
% MATLAB's movie capabilities are kinda brain-dead on my Windows box for 
% some reason, so I had to run this on Linux.  If I were on a Mac I could 
% use qtexport to save 'm' as a quicktime movie.  But, I'm not, so I'm 
% saving as a bunch of TIFFs and then assembling them into a movie using 
% Adobe Premiere.
%
% Notice that you need to make sure there's a directory called 'images' 
% available and that you've 'cd'ed to it's parent or you'll get an error.
m=moviein(2*frames, gcf);
thisFrame=1;
for i=1:frames;
   index=d(round(rand(1)*(maxD-1))+1);
   phi=index/maxC*360;
		h=patch(x,y,'g');
		rotate(h, center, -startAngle);
   		m(:,thisFrame)=getframe(gcf);
		[im,map]=frame2im(m(:,thisFrame));
		if thisFrame < 10;
			pad='00';
		elseif thisFrame < 100;
			pad='0';
		else pad='';
		end;
		s=sprintf('%s%s%i%s', 'images/image', pad, thisFrame, '.tif');
		thisFrame = thisFrame+1;
		imwrite(im, map, s, 'tif');
		delete(h);
		h=patch(x,y,'g');
		rotate(h, center, -phi);
   		m(:,thisFrame)=getframe(gcf);
		[im,map]=frame2im(m(:,thisFrame));
		if thisFrame < 10;
			pad='00';
		elseif thisFrame < 100;
			pad='0';
		else pad='';
		end;
		s=sprintf('%s%s%i%s', 'images/image', pad, thisFrame, '.tif');
		thisFrame = thisFrame+1;
		imwrite(im, map, s, 'tif');
		delete(h);
end;
%movie(m);
