function [flag]=nic_MNI2mask(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_MNI2mask(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 MNI coordinates. Each row vector is a MNI coordinate
%             , which is the center of a spheric ROI, 
%             e.g. [-10,34,20;34,56,1;-46,29,-33;-1,67,20]   
% 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: 
% 1. If you want to use Talairach coordinate to create mask, do convert 
% Talairach to MNI coordinate with tal2mni.m before using nic_MNI2mask
%_______________________________________________________________________
% Copyright (C) 2007 Neuroimage Computing Group
%
% Chun Meng, charlieofmeng@gmail.com
% @(#)nic_MNI2mask.m  ver 1.0, 07/11/29, checked on 07/12/19

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

% -------------------------------------------------------
% Create ROIs  [MNI 1]=[col_vector_x,y,z,constant]*[Voxel 1]
% -------------------------------------------------------
t1 = clock;% start time 
Data=zeros(dimsize);
num_voxel=[];% how many voxels in one ROI
for n=1:num_ROI
    fprintf('--------create ROI: %d\nMNI:(%s),radius:%d mm\t'...
        ,n,num2str(center(n,:)),radius);
    [mask,num_voxel(n)]=seed2ROI(center(n,:),radius,mat,dimsize,voxelsize,n);
    Data=Data+mask;
    fprintf('ROI size:%d voxels \n',num_voxel(n));
end

v.fname=imageOUT;% output name

[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-----------------------------------
% Create one ROI with specified MNI coordinate, radius and index-----------
% -------------------------------------------------------------------------
    function [mask,voxelnumber]=seed2ROI(singleCENTER,singleR,mat,dimsize,voxelsize,index)
        % check whether center is out of box
        vx_center=inv(mat)*([singleCENTER 1]');
        vx_center=round(vx_center(1:3));
        if any(vx_center>dimsize')
            fprintf('Coordinate %d , out of box!\n',index);
            error('The specified center of ROI is out of the box');
        end

        delta_i=ceil(singleR/voxelsize(1))+2;
        delta_j=ceil(singleR/voxelsize(2))+2;
        delta_k=ceil(singleR/voxelsize(3))+2;
        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;
                        
                    end
                end
            end
        end

