function [varargout] = spm_mni2class2(varargin)
% SPM_MNI2CLASS Anatomical Classification of MNI space
% 
% class_brain = spm_mni2class2(XYZmm) produces classification results of
% input coordinates, with a size of nVOX*5 cell structure.
% Results will be classified as 5 categories: 'Brodmann', 'Gyrus',
% 'Hemisphere', 'Lobe' and 'Type'. XYmm is the input voxel coordinates, and
% it should have the size of nVOX*3 (nVOX is the number of input voxels,
% can be 1 to Inf).
% 
% [class_brain] = spm_mni2class2(XYZmm,'f') automatically look for the
% nearest Gray Matter in the defined search cube, if the specific voxel is
% classified as 'White Matter'. The search cube is 5*5*5 voxels, and if a
% GM is found in this cube, a new classification is applied to replace the
% result of 'White Matter' one and also the distance from the found GM to
% the original voxel will be given in mm. If GM is not found, we will stick
% to the original 'White Matter' voxel.
% 
% I/O Variables:
% XYZmm         - {1 x 3} matrix of coordinates of points (MNI coordinates)
% class_brain   - returned results of the class identification of 5 types
%                 in cell data structure 
%
% Also see spm_mni2class.m by Yuqing Gong, March 2004
%----------------------------------------------------------------------------------
%   07.10.04    Yuqing GONG
%   14.10.04    Last modified
%   For bugs and fix, please contact yqgong@stanford.edu
%----------------------------------------------------------------------------------



% Initialize and set running parameters
search_range = 2; % search range in voxel per side, i.e. (2*search_range+1)^3 search cube
flag_GM = 0; % flag tp indicate whether to look for nearest neighbour.

if (nargin==1),
    XYZmm = varargin{1};    
elseif (nargin==2),
    XYZmm = varargin{1};
    if lower(varargin{2}) == 'f';
        flag_GM = 1;
    else
        error('Incorrect input parameter');
    end
else
    error('Incorrect input parameter');
end

nVOX = size(XYZmm,1); % # input voxels

% LOAD ANATOMIC DATABASE
if exist('CAT_DB.mat'),
    load('CAT_DB.mat');
else
    DB_file  = spm_get(1,'CAT_DB.mat','Select the Database file');close;
    load(DB_file);
end
    
% COORDINATES CONVERSION 
% MNI->CUBE
MNI_coord = XYZmm;
Tal_coord = biac_mni2tal(MNI_coord);
CUBE_coord = inv(aheader.magnet_transform.value)*[MNI_coord,ones(nVOX,1)]';
CUBE_coord = round(CUBE_coord(1:3,:)');% nVox*3

% Initialize for voxel labelling results
CAT_index = zeros( nVOX,nCAT );
Dist_VOX = zeros(1,nVOX);
CAT_result = cell(nVOX,nCAT);
% Process voxel by voxel
for index_VOX = 1:nVOX
    
    if flag_GM == 1, % Finding nearest GM flag
        CAT_index_TP = Limage_mat{5}( CUBE_coord(index_VOX,1),CUBE_coord(index_VOX,2),CUBE_coord(index_VOX,3)  ) ;
        % Results of the lookup tables
        if CAT_index_TP~=0 & strcmp( cell2mat(DB{5}( CAT_index_TP )),'White Matter'),  
            % If the classification result for TD-type is 'White Matter',
            % let's start to look for the nearest Gray Matter.
            
            % Construct Search Cube
            if (CUBE_coord(index_VOX,1)-search_range)>0,       xmin = CUBE_coord(index_VOX,1)-search_range;
            else ,            xmin = 1;        end
            if (CUBE_coord(index_VOX,1)+search_range)>0,            xmax = CUBE_coord(index_VOX,1)+search_range;
            else,            xmax = 91;        end
            if (CUBE_coord(index_VOX,2)-search_range)>0,       ymin = CUBE_coord(index_VOX,2)-search_range;
            else ,            ymin = 1;        end
            if (CUBE_coord(index_VOX,2)+search_range)>0,            ymax = CUBE_coord(index_VOX,2)+search_range;
            else,            ymax = 109;        end
            if (CUBE_coord(index_VOX,3)-search_range)>0,       zmin = CUBE_coord(index_VOX,3)-search_range;
            else ,            zmin = 1;        end
            if (CUBE_coord(index_VOX,3)+search_range)>0,            zmax = CUBE_coord(index_VOX,3)+search_range;
            else,            zmax = 91;        end
            
            Search_Cube = Limage_mat{5}(xmin:xmax,ymin:ymax,zmin:zmax);
            %             Search_Cube = Limage_mat{5}( CUBE_coord(index_VOX,1)-search_range:CUBE_coord(index_VOX,1)+search_range,CUBE_coord(index_VOX,2)-search_range:CUBE_coord(index_VOX,2)+search_range,CUBE_coord(index_VOX,3)-search_range:CUBE_coord(index_VOX,3)+search_range  ); 
            sx = size(Search_Cube,1);
            sy = size(Search_Cube,2);
            sz = size(Search_Cube,3);
            dist =inf; cube_c =[0 0 0];
            % Find the nearest Gray Matter in the search cube and also save
            % the distance between them
            for xi = 1:sx
                for yi = 1:sy
                    for zi = 1:sz
                        if (Search_Cube(xi,yi,zi) ==1) | dist > 1,
                            if sqrt(sum([ (xi-search_range-1)^2 (yi-search_range-1)^2 (zi-search_range-1)^2]))<dist
                                dist = sqrt(sum([ (xi-search_range-1)^2 (yi-search_range-1)^2 (zi-search_range-1)^2]));
                                cube_c = [xi-search_range-1,yi-search_range-1,zi-search_range-1];
                            end
                        end %if (Search_Cube(xi,yi,zi) ==1),
                    end
                end
            end % for xi = 1:sx
            CUBE_coord(index_VOX,:) = CUBE_coord(index_VOX,:) +cube_c;
            Dist_VOX(index_VOX) = dist;
        else
            Dist_VOX(index_VOX) = 0;
        end % if strcmp( cell2mat( CAT_result{index_VOX,5}) ,'White Matter'), 
    end %  if flag_GM == 1, % Finding nearest GM flag
    
    % VOXEL-BASED CLASSIFICATION
    for index_CAT = 1:nCAT
        CAT_index(index_VOX,index_CAT ) = Limage_mat{index_CAT}( CUBE_coord(index_VOX,1),CUBE_coord(index_VOX,2),CUBE_coord(index_VOX,3)  ) ;
        if CAT_index(index_VOX,index_CAT ) ==0,% If unidentified
            CAT_result{index_VOX,index_CAT} = {'N/A'}; 
        else
            CAT_result{index_VOX,index_CAT} = DB{index_CAT}( CAT_index(index_VOX,index_CAT ) );
            
        end
    end
    
    
end % for index_VOX = 1:nVOX

% OUTPUT
if (nargout < 1) ,
    % PRINT OUT CLASSIFICATION RESULTS
    for index_VOX = 1:nVOX
        disp( sprintf( '\n\nVoxel #%d of %d',index_VOX, nVOX) ); 
        disp( sprintf( 'Voxel MNI Coordinate: [%d,%d,%d]',XYZmm(index_VOX,:) ) ); 
        if Dist_VOX(index_VOX),
            disp(sprintf('Found nearest gray matter voxel %1.1f mm from the target voxel. ',2*Dist_VOX(index_VOX) ));
        end
        for index_CAT=1:nCAT        
            disp(   sprintf( '%s: %s',TXT{index_CAT}, cell2mat(CAT_result{index_VOX,index_CAT}) )   );  
            %         disp( sprintf( '%s: %s',TXT{index_CAT},CAT_result{index_VOX,index_CAT} ) );  
            
        end
    end
    disp(sprintf('\nDone!'));
else
    varargout{1} = CAT_result;
end    











function outpoints = biac_mni2tal(inpoints)
% Converts coordinates from MNI brain to best guess
% for equivalent Talairach coordinates
% FORMAT outpoints = mni2tal(inpoints)
% Where inpoints is N by 3 or 3 by N matrix of coordinates
%  (N being the number of points)
% outpoints is the coordinate matrix with Talairach points
% Matthew Brett 10/8/99
% Yuqing Gong, 10/15/04

dimdim = find(size(inpoints) == 3);
if isempty(dimdim)
  error('input must be a N by 3 or 3 by N matrix')
end
if dimdim == 2
  inpoints = inpoints';
end

% Transformation matrices, different zooms above/below AC
upT = biac_spm_matrix([0 0 0 0.05 0 0 0.99 0.97 0.92]);
downT = biac_spm_matrix([0 0 0 0.05 0 0 0.99 0.97 0.84]);

tmp = inpoints(3,:)<0;  % 1 if below AC
inpoints = [inpoints; ones(1, size(inpoints, 2))];
inpoints(:, tmp) = downT * inpoints(:, tmp);
inpoints(:, ~tmp) = upT * inpoints(:, ~tmp);
outpoints = inpoints(1:3, :);
if dimdim == 2
  outpoints = outpoints';
end





function [A] = biac_spm_matrix(P)
% returns an affine transformation matrix
% FORMAT [A] = spm_matrix(P)
% P(1)  - x translation
% P(2)  - y translation
% P(3)  - z translation
% P(4)  - x rotation about - {pitch} (radians)
% P(5)  - y rotation about - {roll}  (radians)
% P(6)  - z rotation about - {yaw}   (radians)
% P(7)  - x scaling
% P(8)  - y scaling
% P(9)  - z scaling
% P(10) - x affine
% P(11) - y affine
% P(12) - z affine
%
% A     - affine transformation matrix
%___________________________________________________________________________
%
% spm_matrix returns a matrix defining an orthogonal linear (translation,
% rotation, scaling or affine) transformation given a vector of
% parameters (P).  The transformations are applied in the following order:
%
% 1) translations
% 2) rotations
% 3) scaling
% 4) affine
%
% SPM uses a PRE-multiplication format i.e. Y = A*X where X and Y are 4 x n
% matrices of n coordinates.
%
%__________________________________________________________________________
% @(#)spm_matrix.m	1.1 95/08/07

% pad P with 'null' parameters
%---------------------------------------------------------------------------
q  = [0 0 0 0 0 0 1 1 1 0 0 0];
P  = [P q((length(P) + 1):12)];
A  = eye(4);

A  = A*[1 	0 	0 	P(1);
        0 	1 	0 	P(2);
        0 	0 	1 	P(3);
        0 	0 	0 	1];

A  = A*[1    0   	0   	   0;
        0    cos(P(4))  sin(P(4))  0;
        0   -sin(P(4))  cos(P(4))  0;
        0    0    	0   	   1];

A  = A*[cos(P(5))  0   	sin(P(5))  0;
        0    	   1    0  	   0;
       -sin(P(5))  0  	cos(P(5))  0;
        0          0    0   	   1];

A  = A*[cos(P(6))   sin(P(6))   0  0;
       -sin(P(6))   cos(P(6))   0  0;
        0           0           1  0;
        0     	    0    	0  1];

A  = A*[P(7) 	0   	0    	0;
        0    	P(8) 	0    	0;
        0    	0    	P(9) 	0;
        0    	0    	0    	1];

A  = A*[1   	P(10)   P(11)   0;
        0   	1 	P(12)   0;
        0   	0   	1	0;
        0    	0    	0    	1];


