function spm_xbrain(P,mode)
% Brain extraction.
% FORMAT spm_xbrain
%
% A rough and ready routine for extracting the brain from segmented
% images.  It begins by taking the white matter, and eroding it a
% couple of times to get rid of any odd voxels.  The algorithm
% continues on to do conditional dilations for several iterations,
% where the condition is based upon gray or white matter being present.
% The function has not yet been extensively validated and tested
% so it may not work every time.
%
% Inputs:
% xxx_seg1.img & xxx_seg2.img - grey and white matter segments extracted
% using spm_segment.
%
% Outputs:
% The extracted brain is written to "brain_xxx.img" in the same
% directory as xxx_seg1.img.
% A "render_xxx.mat" file can also be produced that can be used for
% rendering activations on to.
% In Matlab 5.3 and over, a "surf_xxx.mat" file can also be written.
% This extracted brain surface can be viewed using code something
% like:
%    FV = load(spm_get(1,'surf_*.mat','Select surface data'));
%    fg = spm_figure('GetWin','Graphics');
%    ax = axes('Parent',fg);
%    p  = patch(FV, 'Parent',ax,...
%           'FaceColor', [0.8 0.7 0.7], 'FaceVertexCData', [],...
%           'EdgeColor', 'none',...
%	    'FaceLighting', 'phong',...
%           'SpecularStrength' ,0.7, 'AmbientStrength', 0.1,...
%	    'DiffuseStrength', 0.7, 'SpecularExponent', 10);
%    set(0,'CurrentFigure',fg);
%    set(fg,'CurrentAxes',ax);
%    l  = camlight(-40, 20); 
%    axis image;
%    box on;
%    rotate3d on;
%
%
% Note: that the images should be in a right handed co-ordinate system
% (same as the Talairach system) - otherwise the resulting renderings are
% reflections (mirror images) of the true images.
%
%_______________________________________________________________________
% @(#)spm_xbrain.m	2.14 John Ashburner 99/11/04, general render view support by Tomas Sieger, 02/01/17

if nargin==0,
	global BCH
	SPMid = spm('FnBanner',mfilename,'2.14');
	[Finter,Fgraph,CmdLine] = spm('FnUIsetup','XBrain');
	spm_help('!ContextHelp','spm_xbrain.m');

	if isempty(BCH),
		P    = spm_get(2,'*_seg?.img','Select gray and white matter images');
	else,
		P    = spm_input('batch',{},'images');
	end;

	v    = sscanf(version,'%f');
	v    = v(1);
	if v>= 5.3,
		mode = spm_input('Save','+1','m',...
			['Save Extracted Brain|Save Rendering|Save Extracted Surface|'...
			 'Save Extracted Brain and Rendering|'...
			 'Save Extracted Brain and Surface|'...
			 'Save Rendering and Surface|Save All'],[1 2 3 4 5 6 7],4,...
			 'batch',{},'opt');

	else,
		mode = spm_input('Save','+1','m',...
			['Save Extracted Brain|Save Rendering|'...
			 'Save Extracted Brain and Rendering'],[1 2 4],3,...
			 'batch',{},'opt');
    end;
    
	 if any(mode==[2 4 6 7]),
	    xtract_common_views=spm_input_ui('Extract common orthogonal views?','+1','y/n');
	    n_views = spm_input('number of additional views', 1, 'w', 1);
       if n_views
      	 az=zeros(n_views);
          el=zeros(n_views);
          map_morph=zeros(n_views);
          map_morph_cuts_only=zeros(n_views);
          mni_tal_space=zeros(n_views);
          morphology='';
          for i=1:n_views
             if i>1
	             [az(i),el(i),cnd{i},mni_talspace(i),map_morph(i),morphology,map_morph_cuts_only(i)]=...
                   spm_xbrain_getrotation('draw_figure',[spm('ver') ' (' spm('getuser') '):Specifying view ' num2str(i)],morphology,...
                   cnd{i-1},az(i-1),el(i-1));
             else
	             [az(i),el(i),cnd{i},mni_talspace(i),map_morph(i),morphology,map_morph_cuts_only(i)]=...
                   spm_xbrain_getrotation('draw_figure',[spm('ver') ' (' spm('getuser') '):Specifying view ' num2str(i)],morphology);
             end
          end
          additional_views=struct(...
             'n_views',n_views,...
             'az',az,...
             'el',el,...
             'cnd',{cnd},...
             'mni_tal_space',mni_tal_space,...
             'map_morph',map_morph,...
             'morphology',morphology,...
             'map_morph_cuts_only',map_morph_cuts_only);
   	 elseif xtract_common_views=='n' & mode==2
	       % nothing to do
   	    spm('FigName','Xbrain: done',Finter,CmdLine);
          return;
       else
          additional_views=struct('n_views',0);
       end
    else
       xtract_common_views = 'n';
       additional_views=struct('n_views',0);
    end
    
end;

spm('FigName','Xbrain: working',Finter,CmdLine);
do_it(P,mode,xtract_common_views,additional_views);
spm('FigName','Xbrain: done',Finter,CmdLine);

return;
%_______________________________________________________________________

%_______________________________________________________________________
function do_it(P,mode,xtract_common_views,additional_views)

linfun = inline('fprintf([''%-60s%s''],x,[sprintf(''\b'')*ones(1,60)])','x');

spm('Pointer','Watch')

VG=spm_vol(P(1,:));	% Grey matter
VW=spm_vol(P(2,:));	% White matter

% Brain is initially set to white matter...
%-----------------------------------------------------------------------
br=zeros(VW.dim(1:3));
linfun('Initialising');
for i=1:VW.dim(3),
	br(:,:,i) = spm_slice_vol(VW,spm_matrix([0 0 i]),VW.dim(1:2),1);
end;

% Build a 3x3x3 seperable smoothing kernel
%-----------------------------------------------------------------------
kx=[0.75 1 0.75];
ky=[0.75 1 0.75];
kz=[0.75 1 0.75];
sm=sum(kron(kron(kz,ky),kx))^(1/3);
kx=kx/sm; ky=ky/sm; kz=kz/sm;

% Erosions and conditional dilations
%-----------------------------------------------------------------------

niter = 32;
spm_progress_bar('Init',niter,'Extracting Brain','Iterations completed');
for j=1:niter,
	if j>2, th=0.15; else, th=0.6; end; % Dilate after two its of erosion. 
	linfun(['Iteration ' num2str(j) ' - thresholding and multiplying']);
	for i=1:VW.dim(3),
		w=spm_slice_vol(VW,spm_matrix([0 0 i]),VW.dim(1:2),1);
		g=spm_slice_vol(VG,spm_matrix([0 0 i]),VW.dim(1:2),1);
		br(:,:,i) = (br(:,:,i)>th).*(w+g);
	end;
	linfun(['Iteration ' num2str(j) ' - convolving']);
	spm_conv_vol(br,br,kx,ky,kz,-[1 1 1]);
	spm_progress_bar('Set',j);
end;
spm_progress_bar('Clear');

% Generate output filename
%-----------------------------------------------------------------------
[pth,nam,ext,ver] = fileparts(deblank(P(1,:)));
ind = findstr(nam,'_seg1');
if ~isempty(ind), nam((0:4)+ind(1))=[]; end;
fname = fullfile(pth,['brain_' nam ext ver]);

if any(mode==[2 4 6 7]) & ( xtract_common_views=='y' | additional_views.n_views>0),
	% Produce rendering
	%-----------------------------------------------------------------------
	matname = fullfile(pth,['render_' nam '.mat']);
	tmp = struct('dat',br,'dim',size(br),'mat',VG.mat);
	renviews(tmp,matname,xtract_common_views,additional_views);
end;

if any(mode==[3 5 6 7]),
	% Produce extracted surface
	%-----------------------------------------------------------------------
	linfun(['Extracting surface - please be patient']);
	matname = fullfile(pth,['surf_' nam '.mat']);
	tmp = struct('dat',br,'dim',size(br),'mat',VG.mat);

	% This is done with an eval statement because Matlab5.2
	% complains otherwise...
	eval('[faces,vertices] = isosurface(br,0.5);');

	% Swap around x and y because isosurface does for some
	% wierd and wonderful reason.
	Mat = VG.mat(1:3,:)*[0 1 0 0;1 0 0 0;0 0 1 0; 0 0 0 1];
	vertices = (Mat*[vertices' ; ones(1,size(vertices,1))])';
	save(matname,'faces','vertices');
end;

if any(mode==[1 4 5 7]),
	% Final cleanup
	%-----------------------------------------------------------------------
	linfun(['Final cleanup - thresholding and multiplying']);
	for i=1:VW.dim(3),
		w  = spm_slice_vol(VW,spm_matrix([0 0 i]),VW.dim(1:2),1);
		g  = spm_slice_vol(VG,spm_matrix([0 0 i]),VW.dim(1:2),1);
		th = 0.15;
		br(:,:,i) = ((br(:,:,i)>th).*(w+g))>th;
	end;

	% Write the extracted brain
	%-----------------------------------------------------------------------
	linfun('Writing volume');
	VO=struct('fname',fname,...
		  'dim',[VW.dim(1:3) spm_type('uint8')],...
		  'mat',VW.mat,...
		  'pinfo',[1/255 0 0]',...
		  'descrip','extracted brain');
	spm_write_vol(VO,br);
end;

linfun(' ');
spm('Pointer')

return;
%_______________________________________________________________________

%_______________________________________________________________________
function renviews(V,oname,xtract_common_views,additional_views)
% Produce images for rendering activations to
%
% FORMAT renviews(V,oname,xtract_common_views,AV)
% V     - mapped image to render, or alternatively
%         a structure of:
%         V.dat - 3D array
%         V.dim - size of 3D array
%         V.mat - affine mapping from voxels to millimeters
% oname - the name of the render.mat file.
% xtract_common_views - 0/1 flag specifying whether (not) to render 6 
%			 common views
%    *** following fields require V being a structure, not mapped 
%		   image: ***
% additional_views - structure of:
%         n_views - number of additional views
%   * some of following fields are arrays / cell arrays of 'n_views' 
%	   elements: *
%         az, al - array of azimuths and elevations specifying
%			 		additional views rotations
%         cnd - cell array of string conditions specifying which voxels 
%              to include into a view; condition is evaluated for each voxel,
%					if result of evaluation is 1, the voxel will be included in
%					view, on 0 result, the voxel will be exluced; condition can 
%					refer to x,y,z variables holding MNI/Talairach space voxel 
%					coordinates
%         mni_tal_space - array of flags 0/1 specifying whether the 'cnd' 
%              should be evaluated in MNI (0) or Talairach (1) space
%         map_morph - array of flags 0/1 specifying optional mapping 
%              of morphological scan onto the brain
%         morphology - normalized morphology image filename
%         map_morph_cuts_only - array of flags 0/1 specifying whether morphology
%              should be mapped onto brain cuts only (i.e. not onto brian 
%					surface)
%_______________________________________________________________________
%
% Produces a matrix file "render_xxx.mat" which contains everything that
% "spm_render" is likely to need.
%
% Ideally, the input image should contain values in the range of zero
% and one, and be smoothed slightly.  A threshold of 0.5 is used to
% distinguish brain from non-brain.
%_______________________________________________________________________

linfun = inline('fprintf([''%-30s%s''],x,[sprintf(''\b'')*ones(1,30)])','x');
linfun('Rendering: ');
v      = V.dim(1:3).*sqrt(sum(V.mat(1:3,1:3).^2));
M      = V.mat;
shift0 = inv(spm_matrix(v/2));
shift  = spm_matrix(V.dim(1:3)/2);
zoom   = spm_matrix([0 0 0 0 0 0 sqrt(sum(V.mat(1:3,1:3).^2))]);

if xtract_common_views=='y'
   
   MT0 = zoom;
   
   MT1 = MT0 * shift * spm_matrix([0 0 0 0 pi 0]) * inv(shift);

	MS0 = spm_matrix(v([3 2 1])/2) * spm_matrix([0 0 0 0 pi/2]) ...
		* shift0 * zoom;

	MS1 = MS0 * shift * spm_matrix([0 0 0 0 0 pi]) * inv(shift);

	MC0 = spm_matrix(v([3 1 2])/2) * spm_matrix([0 0 0 0 0 pi/2]) ...
		* spm_matrix([0 0 0 pi/2]) * shift0 * zoom;

	MC1 = MC0 * shift * spm_matrix([0 0 0 0 0 pi]) * inv(shift);
end

% additional views:
if additional_views.n_views>0
	twiddle=[...
   	   [0 0 0]',...
      	[0 0 1]',...
	      [0 1 0]',...
   	   [0 1 1]',...
      	[1 0 0]',...
	      [1 0 1]',...
   	   [1 1 0]',...
	      [1 1 1]'...
   	];
   vertices=twiddle.*repmat(V.dim(1:3)',1,8);
   
   for i=1:additional_views.n_views
 		% rotation: [ around axis normal to sagital, around axis normal to coronal, around axis normal to transverse ]
      M = spm_matrix([[0 0 0] 0 (additional_views.el(i)+90)*(pi/180) (additional_views.az(i))*(pi/180)]) * shift0 * zoom;
		vdim = max(M*[vertices;ones(1,size(vertices,2))],[],2)-min(M*[vertices;ones(1,size(vertices,2))],[],2);
		vdim=vdim(1:3)';
		MO{i} = spm_matrix(vdim/2) * M;
		view_dim{i}=vdim;
   end
end

if max(sqrt(sum(V.mat(1:3,1:3).^2)))<=1,  % 0.5mm resolution
	zm  = diag([2 2 1 1]); % vox vox mm 1
   v   = v*2;
   if xtract_common_views=='y'
	   MT0 = zm*MT0; MT1 = zm*MT1;
		MS0 = zm*MS0; MS1 = zm*MS1;
	   MC0 = zm*MC0; MC1 = zm*MC1;
   end
   for i=1:additional_views.n_views
      view_dim{i}=view_dim{i}*2;
      MO{i} = zm*MO{i};
   end
end;

if det(M)<0,
	flp=spm_matrix([0 (v(2)+1) 0  0 0 0 1 -1 1]); MT0=flp*MT0; MT1=flp*MT1;
	flp=spm_matrix([0 (v(2)+1) 0  0 0 0 1 -1 1]); MS0=flp*MS0; MS1=flp*MS1;
	flp=spm_matrix([0 (v(1)+1) 0  0 0 0 1 -1 1]); MC0=flp*MC0; MC1=flp*MC1;
   for i=1:additional_views.n_views
      flp=spm_matrix([0 (view_dim{i}(2)+1) 0  0 0 0 1 -1 1]); MO{i}=flp*MO{i}; 
   end
end;

if xtract_common_views=='y'
   additional_view=struct(...
      'cnd','1',...
      'm_vox2mni',V.mat,... % not used
      'mni_tal_space',0,...
	   'map_morph',0,...
	   'morphology','',...
      'map_morph_cuts_only',0);
	linfun('Rendering: Transverse 1..');	rend{1} = make_struct(V,MT0,v([1 2]),additional_view);
	linfun('Rendering: Transverse 2..');	rend{2} = make_struct(V,MT1,v([1 2]),additional_view);
	linfun('Rendering: Saggital 1..');	rend{3} = make_struct(V,MS0,v([3 2]),additional_view);
	linfun('Rendering: Saggital 2..');	rend{4} = make_struct(V,MS1,v([3 2]),additional_view);
	linfun('Rendering: Coronal 1..');	rend{5} = make_struct(V,MC0,v([3 1]),additional_view);
   linfun('Rendering: Coronal 2..');	rend{6} = make_struct(V,MC1,v([3 1]),additional_view);
   n_rend=6;
else
   n_rend=0;
end
for i=1:additional_views.n_views
   additional_view=struct(...
      'cnd',additional_views.cnd(i),...
      'm_vox2mni',V.mat,...
      'mni_tal_space',additional_views.mni_tal_space(i),...
	   'map_morph',additional_views.map_morph(i),...
	   'morphology',additional_views.morphology,...
      'map_morph_cuts_only',additional_views.map_morph_cuts_only(i));
   linfun(['Rendering: Additional view ' num2str(i)]);
   rend{n_rend+i} = make_struct(V,MO{i},view_dim{i}([1 2]),additional_view);
end

linfun('Rendering: Save..');
save(oname,'rend');
linfun(' ');
disp_renderings(rend);
spm_print;
return;
%_______________________________________________________________________

%_______________________________________________________________________
function brain=cut_brain(brain,m_vox2mni,cnd,mni_tal_space)
% exclude some voxels from brain
% FORMAT: cut_brain(brain,m_vox2mni,cnd,mni_tal_space)
%
% params:
%		brain - brain volume
%		m_vox2mni - matrix mapping voxels to mm (in MNI space)
%		cnd - condition specifying which voxels to include in the brain
%		mni_tal_space - flag 0/1 specifying MNI/Talairach space
xx=1:size(brain,1);
yy=1:size(brain,2);
zz=1:size(brain,3);
[y0,z0]=ndgrid(yy,zz);

spm_progress_bar('Init',size(brain,1),'Cutting off','slices completed');
for xi=1:size(brain,1)
   spm_progress_bar('Set',xi);
   points=[repmat(xi,1,size(brain,2)*size(brain,3)); y0(:)'; z0(:)'; ones(1,size(brain,2)*size(brain,3))];
   points=m_vox2mni*points;
   points=points(1:3,:);
   if mni_tal_space==1
      points=mni2tal(points);
   end
   x=points(1,:);
   y=points(2,:);
   z=points(3,:);
   msk=eval(['~(' cnd ')']);
   brain(xi,msk)=0;
end
spm_progress_bar('Clear');
%_______________________________________________________________________

%_______________________________________________________________________
function str = make_struct(V,M,D,map_morph)
	[ren,dep]=make_pic(V,M,D,map_morph);
	str=struct('M',M/V.mat,'ren',ren,'dep',dep);
return;
%_______________________________________________________________________

%_______________________________________________________________________
function [ren,zbuf]=make_pic(V,M,D,additional_view)
%
% additional_view: struct containing
%     cnd - condition specifying which voxels will be included in a view 
%				(expression yielding 0/1 value, can refer to x,y,z variables 
%				- voxel coordinates in MNI/Talairach space)
%     m_vox2mni - matrix mapping voxel space to MNI mm space
%     mni_tal_space - flag 0/1 specifying whether x,y,z should be 
%				considered in MNI or Talairach space
% 	   map_morph - flag 0/1 specifying whether to map morphology onto brain
% 	   morphology - (in MNI space) normalized morphology filename
%     map_morph_cuts_only - flag 0/1 specifying whether to map morphology
%				onto cuts only or onto cuts and brain surface
% A bit of a hack to try and make spm_render_vol produce some slightly
% prettier output.  It kind of works...
if isfield(V,'dat'), vv = V.dat; else, vv = V; end;

if strcmp(additional_view.cnd,'1')==0
   vv_whole_brain=vv;
   vv=cut_brain(vv,additional_view.m_vox2mni,additional_view.cnd,additional_view.mni_tal_space);
end

[REN, zbuf, X, Y, Z] = spm_render_vol(vv, M, D, [0.5 1]);
fw        = max(sqrt(sum(M(1:3,1:3).^2)));
msk       = find(zbuf==1024);
brn       = ones(size(X));
brn(msk)  = 0;
brn       = spm_conv(brn,fw);
X(msk)    = 0;
Y(msk)    = 0;
Z(msk)    = 0;
msk       = find(brn<0.5);
tmp       = brn;
tmp(msk)  = 100000;
sX        = spm_conv(X,fw)./tmp;
sY        = spm_conv(Y,fw)./tmp;
sZ        = spm_conv(Z,fw)./tmp;
zbuf      = spm_conv(zbuf,fw)./tmp;
zbuf(msk) = 1024;

vec       = [-1 1 3]; % The direction of the lighting.
vec       = vec/norm(vec);
[t,dx,dy,dz] = spm_sample_vol(vv,sX,sY,sZ,3);
if fw<=1 % if 0.5mm resolution is used (see renviews function)
	IM        = inv(diag([0.5 0.5 1])*M(1:3,1:3))';
else
	IM        = inv(M(1:3,1:3))';
end
ren       = IM(1:3,1:3)*[dx(:)' ; dy(:)' ; dz(:)'];
len       = sqrt(sum(ren.^2,1))+eps;
ren       = [ren(1,:)./len ; ren(2,:)./len ; ren(3,:)./len];
ren       = reshape(vec*ren,[size(dx) 1]);
ren(find(ren<0))  = 0;
ren(msk)  = ren(msk)-0.2;
ren       = ren*0.8+0.2;
mx        = max(ren(:));
ren       = ren/mx;

if additional_view.map_morph==1
	VM=spm_vol(additional_view.morphology);
   vm=spm_sample_vol(VM,sX,sY,sZ,1);
   if additional_view.map_morph_cuts_only==1
      % find voxels on the edge of the cut, not on the brain surface
      
      sX=round(sX);
      sY=round(sY);
      sZ=round(sZ);
      
      if 0
      idx=find(sX~=0);
      sX=sX(idx);sY=sY(idx);sZ=sZ(idx);
      sX_1=max(1,sX-1);
      sX_2=min(size(vv,1),sX+1);
      sY_1=max(1,sY-1);
      sY_2=min(size(vv,2),sY+1);
      sZ_1=max(1,sZ-1);
      sZ_2=min(size(vv,3),sZ+1);
   	end
      
      msk=zeros(size(sX));
		spm_progress_bar('Init',size(msk,1),'Finding cuts','rows completed');
      for i=1:size(sX,1)
         for j=1:size(sX,2)
            if sX(i,j)>0 & ...
               (vv_whole_brain(max(1,sX(i,j)-2),sY(i,j),sZ(i,j))<=.5 | ...
               vv_whole_brain(min(size(vv,1),sX(i,j)+2),sY(i,j),sZ(i,j))<=.5 | ...
	            vv_whole_brain(sX(i,j),max(1,sY(i,j)-2),sZ(i,j))<=.5 | ...
   	         vv_whole_brain(sX(i,j),min(size(vv,2),sY(i,j)+2),sZ(i,j))<=.5 | ...
      	      vv_whole_brain(sX(i,j),sY(i,j),max(1,sZ(i,j)-2))<=.5 | ...
               vv_whole_brain(sX(i,j),sY(i,j),min(size(vv,2),sZ(i,j)+2))<=.5)
      	      msk(i,j)=1;
	         end
   	   end
         spm_progress_bar('Set',i);
      end
      spm_progress_bar('Clear');
   else
      msk=[];
   end
   idx_msk=find(msk);
   vm(idx_msk)=0;
   min_vm=min(vm(:));
   vm=(vm-min_vm)/(max(vm(:))-min_vm); % morphology in range <0,1>
   vm(idx_msk)=1;
	% map morphology by multiplying brain surface lighting in the range 
   % <0,1> by morphology intensity in the range <0,1> 
   ren=ren.*vm;
   if ~isempty(msk)
   %   ren(idx_cuts)=ren(idx_cuts).*vm(idx_cuts);
      ren(find(msk==0))=vm(find(msk==0));
   end
   end

return;
%_______________________________________________________________________
%_______________________________________________________________________
function disp_renderings(rend)
Fgraph = spm_figure('GetWin','Graphics');
spm_results_ui('Clear',Fgraph);
hght = 0.95;
nrow = ceil(length(rend)/2);
ax=axes('Parent',Fgraph,'units','normalized','Position',[0, 0, 1, hght],'Visible','off');
image(0,'Parent',ax);
set(ax,'YTick',[],'XTick',[]);

for i=1:length(rend),
	ren = rend{i}.ren;
	ax=axes('Parent',Fgraph,'units','normalized',...
		'Position',[rem(i-1,2)*0.5, floor((i-1)/2)*hght/nrow, 0.5, hght/nrow],...
		'Visible','off');
	image(ren*64,'Parent',ax);
	set(ax,'DataAspectRatio',[1 1 1], ...
		'PlotBoxAspectRatioMode','auto',...
		'YTick',[],'XTick',[],'XDir','normal','YDir','normal');
end;
drawnow;
return;
%_______________________________________________________________________
