function img = lyngby_talaimg(volume,mask,sx,sy,sz,x1,x2,y1,y2,z1,z2)

% lyngby_talaimg		- Returns Talairach brain plot in the style of SPM
%
%	function img = lyngby_talaimg(volume,mask,numVoxel,x1,x2,y1,y2,z1,z2)
%
%	Input:	volume: 	vector of length sx*sy*sz ordered in sz number of slices sx*sy
%         			that is, the first sx*sy elements in volume is slice1 etc.
%		mask:		
% 		numVoxel	Size x of volume
%		x1,x2,y1,y2,z1,z2:	Volume to be picked out
%
%	Output:	img: 		Image returned
%
%	Plot the 3 projections of a 3D brain volume in a Talairach grid
%	The brain volume contained in 'volume' must have been normalized
%	to Talaraich.

%	Date		Author			Version
%	1996		Blund, Bent Lundsager	Initial version
%			Benny Kristensen
%	1996-02-19	fnielsen		
%	1996-09-18	fnielsen

% cvs :  $Revision: 1.1.1.1 $

	sx = numVoxel(1);
	sy = numVoxel(2);
	sz = numVoxel(3);


	% Rescale volume
		%disp('Rescaling voxel values...')
		volume = volume - min(volume');

	% First do the 3 projections
		disp('Doing projections..')

		if 1==1
			proj1=reshape(sum(reshape(volume,sx,sy*sz)),sy,sz)';
			proj1=(proj1-min(min(proj1)))/(max(max(proj1))-min(min(proj1)))*158/160+1/160;
			proj2=reshape(sum(reshape(reshape(volume,sx,sy*sz)',sy,sz*sx)),sz,sx);
			proj2=(proj2-min(min(proj2)))/(max(max(proj2))-min(min(proj2)))*158/160+1/160;
			proj3=reshape(sum(reshape(volume,sx*sy,sz)'),sx,sy);
			proj3=(proj3-min(min(proj3)))/(max(max(proj3))-min(min(proj3)))*158/160+1/160;
		else
			%themin1=3.7958e+04;
			%themin2=3.7474e+04;
			%themin3=1.2376e+04;
			%themax1=2.6153e+05;
			%themax2=2.5677e+05;
			%themax3=1.0743e+05;
	
			%proj1=reshape(sum(reshape(volume,sx,sy*sz)),sy,sz)';
			%proj1=(proj1-themin1)/(themax1-themin1)*158/160+1/160;
			%proj2=reshape(sum(reshape(reshape(volume,sx,sy*sz)',sy,sz*sx)),sz,sx);
			%proj2=(proj2-themin2)/(themax2-themin2)*158/160+1/160;
			%proj3=reshape(sum(reshape(volume,sx*sy,sz)'),sx,sy);
			%proj3=(proj3-themin3)/(themax3-themin3)*158/160+1/160;
		end

	% Add mask
		disp('Masking...')
		proj=reshape(sum(reshape(mask,sx,sy*sz)),sy,sz)';
		proj1=proj1.*(proj>0);
		proj=reshape(sum(reshape(reshape(mask,sx,sy*sz)',sy,sz*sx)),sz,sx);
		proj2=proj2.*(proj>0);
		proj=reshape(sum(reshape(mask,sx*sy,sz)'),sx,sy);
		proj3=proj3.*(proj>0);
		clear mask

	% Pick out box of right size
		disp('Picking out volume with voxels containing brain...')
		sx=x2-x1+1;
		sy=y2-y1+1;
		sz=z2-z1+1;
		[x1 x2 y1 y2 z1 z2]
		proj1=proj1(z1:z2,y1:y2);
		proj2=proj2(z1:z2,x1:x2);
		proj3=proj3(x1:x2,y1:y2);

	% Flip projections as needed
		proj1=flipud(proj1);
		proj2=flipud(proj2);
		proj1=fliplr(proj1);
		proj3=fliplr(proj3);

	% Get the Talairach grid
		load MIP.mat
		mip=flipud(mip') & 1;

	% Important corners in MIP grid:
		x1=41;
		y1=22;
		x2=144;
		y2=192;
		xx1=41;
		yy1=213;
		xx2=144;
		yy2=340;
		xxx1=185;
		yyy1=22;
		xxx2=312;
		yyy2=192;

	% Stretch the projections to match grid
		disp('Stretching projections...')
		projimage=zeros(size(mip));
		for x=x1:x2
			for y=y1:y2
				projimage(x,y)=proj1(round((x-x1)/(x2-x1)*(sz-1)+1),round((y-y1)/(y2-y1)*(sy-1)+1));
  			end
		end
		clear proj1
		for x=xx1:xx2
			for y=yy1:yy2
				projimage(x,y)=proj2(round((x-xx1)/(xx2-xx1)*(sz-1)+1),round((y-yy1)/(yy2-yy1)*(sx-1)+1));
  			end
		end
		clear proj2
		for x=xxx1:xxx2
			for y=yyy1:yyy2
				projimage(x,y)=proj3(round((x-xxx1)/(xxx2-xxx1)*(sx-1)+1),round((y-yyy1)/(yyy2-yyy1)*(sy-1)+1));
			end
		end
		clear proj3

	% Add grid and projections
		disp('Getting ready to show image...')
		img=(~mip & projimage).*projimage+mip;
		clear mip	
		clear projimage

	% Show image
		mymap=jet(200);
		mymap=mymap(41:200,:);
		mymap(1,:)=[1 1 1];
		mymap(160,:)=[0 0 0];
		colormap(mymap)
		imagesc(img);
		axis('off')





