function [varargout] = spm_findGM(XYZmm)
% Format 1: [class_brain, GMcoord] = spm_findGM(XYZmm)
% Format 2: [class_brain, class_brain_s, GMcoord] = spm_findGM(XYZmm)
tic;
% SearchRange = [4 4 4]; % In mm
SearchRange = [3 3 3]; % In mm

SearchStep = 2; % In mm

% Get the XYZ MNI coordinates for the target voxel
x_tar = XYZmm(1);
y_tar = XYZmm(2);
z_tar = XYZmm(3);
% Initialize
% GMlist = {};
GMcounter = 0;
GMdist = 100;
GMcoord = zeros(1,3);
% temp = 0;

for index_x = sort_abs( -1*SearchRange(1):SearchStep:SearchRange(1) ),
    for index_y = sort_abs( -1*SearchRange(2):SearchStep:SearchRange(2) ),
        for index_z = sort_abs( -1*SearchRange(3):SearchStep:SearchRange(3) ),
            if (GMdist > SearchStep),
                curr_coord = XYZmm + [index_x,index_y,index_z];
                if sum( abs([index_x,index_y,index_z]) ) ~= 0 ...
                        &  sqrt(   sum( ([index_x,index_y,index_z]).^2 )  ) < GMdist ,              
                    [class_brain, class_brain_s] = spm_mni2class( curr_coord  );
%                     temp = temp + 1
%                     index_x
                    if strcmp(   cell2mat( class_brain_s{5}),'gray'  ) ==1,
                        % GMlist = {GMlist,curr_coord};
                        % GMdist = [ GMdist, sqrt(   sum( ([index_x,index_y,index_z]).^2 )  ) ]
                        GMcoord = curr_coord;
                        GMdist = sqrt(   sum( ([index_x,index_y,index_z]).^2 )  );
                        GMcounter = GMcounter + 1;
                    end % if strcmp(   cell2mat( class_brain_s{5}),'gray'  ) ==1,
                end % if sum( abs([index_x,index_y,index_z]) ) ~= 0 
            end % if (GMdist > 1),
        end
    end
end

if GMcounter == 0,
    class_brain=cell(5,1);
    class_brain_s=cell(5,1);
%     for index=1:5
%         class_brain{index}={'None'};
%         class_brain_s{index}={'None'};
%     end %for index=1:5
    GMcoord = [inf,inf,inf];
    if (nargout ==0),
        disp(sprintf('No GM is found even in %dx%dx%dmm cube!',(SearchRange*2+2)));
%         spm_popup_class(class_name,class_brain);
    elseif (nargout ==2),
        varargout{1} = class_brain;
        varargout{2} = GMcoord;
    elseif (nargout ==3),
        varargout{1}= class_brain;
        varargout{2}= class_brain_s;
        varargout{3} = GMcoord;
    end
    
else % if GMcounter == 0,
    [class_brain, class_brain_s] = spm_mni2class(GMcoord);
    
    
    
    if (nargout ==0),
        disp(sprintf('Found nearest gray matter voxel at [%2.1f,%2.1f,%2.1f] (MNI),',GMcoord));
        disp(sprintf('and it is %1.1f mm from the target voxel.\n',GMdist));    
        class_brain = spm_mni2class(GMcoord);
    elseif (nargout ==2),
        varargout{1} = class_brain;
        varargout{2} = GMcoord;
    elseif (nargout ==3),
        varargout{1}= class_brain;
        varargout{2}= class_brain_s;
        varargout{3} = GMcoord;
    end    
    
end

toc
return

% ======================================
function B = sort_abs(A)

A_abs = abs(A);
[D,IX] = sort(A_abs,2);
B = A(IX);
return