function [flag]=nic_Coordinate2mask(imageREF,center,radius,imageOUT)
% Create spheric ROIs (indexed by 1,2,3...) with specified Talairach 
% coordinates and radius (center,radius), and output to a mask image.
% The mask image is rearranged and output (imageOUT) with the  same 
% format as the reference image (imageREF). 
%
% FORMAT [flag]=nic_Coordinate2mask(imageREF,center,radius,imageOUT)
% 
% Input:
% 1. imageREF - full file path and name for the reference image
%                i.e., the post-processed image [path\*.img]
% 2. center - N*3 matrix of Talairach coordinates, which are the spheric
%             ROI centers in row, e.g. [-10,34,20;56,10,33]   
% 3. radius - the radius for all spheric ROIs (millimeter), e.g. [6]
% 4. imageOUT - full file path and name for the output mask [path\*.img]
% 
% Output:
% 1. flag - a flag for all done, 1: successful, 0: fail
%
% Note: If you want to use MNI coordinate to create mask, do convert MNI to
%       Talairach coordinate with tal2mni.m before using nic_Coordinate2mask
%_______________________________________________________________________
% Copyright (C) 2007 Neuroimage Computering Group, State Key Laboratory of
% Cognitive Neuroscience and Learning
%
% Chun Meng, charlieofmeng@gmail.com
% @(#)nic_Coordinate2mask.m  ver 1.0, 07/11/27

% -------------------------check input---------------------------
if nargin ~=4
    error('Providing sufficient parameters again');
end

dim_center=size(center);
if dim_center(2)~=3
    error('Coordinate of center is not in the form of (x,y,z)');
end
% ---------------------------------------------------------------

num_ROI=dim_center(1);
clear('dim_center');

if num_ROI>1
    % check distance between centers >2*radius, in case of overlapped ROI
    distance=[];
    c=nchoosek(1:num_ROI,2);
    for i=1:size(c,1)
        distance(i)=norm(center(c(i,1),:)-center(c(i,2),:));
        if distance(i)<=2*radius
            error('ROIs overlap');
        end
    end

    clear('distance','c');
end

% read imageREF, nifti img, to get volume size and affine matrix
[Data,v] = nic_ReadNiftiImage(imageREF,1);
clear('Data');

dimsize=v.dim;% e.g. [61 73 61]
mat=v.mat;
voxelsize=sum(abs(mat(1:3,1:3)));

% % convert Talairach to MNI
TTcenter=center;
center=tal2mni(TTcenter);

% -------------------------------------------------------
% 创建ROI  [MNI 1]=[col_vector_x,y,z,constant]*[Voxel 1]
% -------------------------------------------------------
t1 = clock;% 记录计算费时
Data=zeros(dimsize);
num_voxel=[];%记录ROI中有多少个voxel
for n=1:num_ROI
    fprintf('--------create ROI: %d\nTalairach:(%s),radius:%d mm\n'...
        ,n,num2str(TTcenter(n,:)),radius);
    [mask,num_voxel(n)]=seed2ROI(center(n,:),radius,mat,dimsize,voxelsize,n);
    Data=Data+mask;
    fprintf('MNI:(%s)\tROI size:%d voxels \n',...
        num2str(center(n,:)),num_voxel(n));
end

v.fname=imageOUT;% 转换为输出mask文件名

[flag] = nic_WriteMaskImage(Data,v,imageOUT);
if flag == 1
    t2 = etime(clock, t1);
    fprintf('--------Finish mask creation and output to %s\n',imageOUT);
    fprintf('--------Elapsed time: %s s\n',num2str(t2));
else
    error('Please carefully check the whole processing step, do again')
end

% --------------------------sub function-----------------------------------
% -------------------------------------------------------------------------
% -----------给定一个MNI坐标作为球心，以及半径，返回指定标号的mask
    function [mask,voxelnumber]=seed2ROI(singleCENTER,singleR,mat,dimsize,voxelsize,index)
        % 检查给定球心坐标是否超出volume
        vx_center=inv(mat)*([singleCENTER 1]');
        vx_center=round(vx_center(1:3));% 四舍五入，保持整数；列向量代表i j k
        if any(vx_center>dimsize')
            fprintf('第%d个坐标 , out of box!\n',index);
            error('The specified center of ROI is out of the box');
        end
        % 遍历以该点为中心，2倍半径以内范围的cubic中的voxel
        delta_i=ceil(singleR/voxelsize(1))+2;% 注意 voxelsize为正
        delta_j=ceil(singleR/voxelsize(2))+2;
        delta_k=ceil(singleR/voxelsize(3))+2;% ＋2是增大cubic防止中心点四舍五入可能造成的误差
        i1=max([0,vx_center(1)-delta_i]);
        i2=min([vx_center(1)+delta_i,dimsize(1)]);
        j1=max([0,vx_center(2)-delta_j]);
        j2=min([vx_center(2)+delta_j,dimsize(2)]);
        k1=max([0,vx_center(3)-delta_k]);
        k2=min([vx_center(3)+delta_k,dimsize(3)]);
        mask=zeros(dimsize);
        voxelnumber=0;
        for i=i1:i2
            for j=j1:j2
                for k=k1:k2
                    tmp=mat*[i j k 1]';                     
                    if norm(tmp(1:3)'-singleCENTER)<=singleR
                        mask(i,j,k)=index;
                        voxelnumber=voxelnumber+1;
                        % 判断该点i j k 在脑外的代码暂无，因为脑内模板未定，目前不判断ROI是否在脑外
                    end
                end
            end
        end

% ------------------------------------从宋晓伟REST拷贝来-------Talairach2MNI
    function outpoints = tal2mni(inpoints)
        % Converts coordinates to MNI brain best guess
        % from Talairach coordinates
        % FORMAT outpoints = tal2mni(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 MNI points
        % Matthew Brett 2/2/01

        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
        rotn = rest_spm_matrix([0 0 0 0.05]);
        upz = rest_spm_matrix([0 0 0 0 0 0 0.99 0.97 0.92]);
        downz = rest_spm_matrix([0 0 0 0 0 0 0.99 0.97 0.84]);

        inpoints = [inpoints; ones(1, size(inpoints, 2))];
        % Apply inverse translation
        inpoints = inv(rotn)*inpoints;

        tmp = inpoints(3,:)<0;  % 1 if below AC
        inpoints(:, tmp) = inv(downz) * inpoints(:, tmp);
        inpoints(:, ~tmp) = inv(upz) * inpoints(:, ~tmp);
        outpoints = inpoints(1:3, :);
        if dimdim == 2
            outpoints = outpoints';
        end

        function [A] = rest_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
            % (i.e., the opposite to which they are specified):
            %
            % 1) shear
            % 2) scale
            % 3) rotation - yaw, roll & pitch
            % 4) translation
            %
            % SPM uses a PRE-multiplication format i.e. Y = A*X where X and Y are 4 x n
            % matrices of n coordinates.
            %
            %__________________________________________________________________________
            % Copyright (C) 2005 Wellcome Department of Imaging Neuroscience

            % Karl Friston
            % $Id: spm_matrix.m 112 2005-05-04 18:20:52Z john $


            % 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)];

            T  =   [1 	0 	0 	P(1);
                0 	1 	0 	P(2);
                0 	0 	1 	P(3);
                0 	0 	0 	1];

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

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

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

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

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

            A = T*R1*R2*R3*Z*S;
