function [outputNAME] = nic_Image2mask(imageIN,thresholdRange,imageOUT,imageREF)
% Creating a binary mask in which 1 means the value of input image
% (imageIN) is in  the  input range of data value (thresholdRange). The mask
% image will be rearranged and output with the same format as the reference
% image, which is usually an arbitrary image one wish to apply the mask to.
% 
% FORMAT [LaberMask,NUM] = nic_Image2mask(imageP,Pthreshold,imageREF,imageOUT)
% 
% Input:
% 1. imageIN        - full file path and name for the p image, NOT (1-P) image
%                       [path\*.img or .nii]
% 2. thresholdRange - a array of min and max data that we care about in
%                        imageIN. e.g.[-inf -0.05], [-0.05,0.05], or [0.05,inf]
% 3. imageOUT       - the path and filename of image file to output 
%                       [path\*.img or .nii]
% 4. imageREF       - full file path and name for the reference image,
%                 [path\*.img].The reference image is the image which has 
%                 the same orientation, voxel size, dimension, origin, and 
%                 bounding box as the image you wish to apply the mask to. 
%                 If already having the images to apply the mask to, you 
%                 can simply put one of them as the reference image.
% 
% Output:
% 1. outputNAME - full file path and name of the resampled image for output. 
%_______________________________________________________________________
% Copyright (C) 2007 Neuroimage Computing Group, State Key Laboratory of
% Cognitive Neuroscience and Learning
%
% Yujin Zhang, zyjinjin@gmail.com
% @(#)nic_Pmask.m  ver 1.0, 07/11/30

[Data,Head] = nic_ReadNiftiImage(imageIN,1);

V_ref = nic_spm_vol(imageREF);
if any(any(Head.mat - V_ref.mat)) || any((Head.dim-V_ref.dim),2)
    index = strfind(imageIN,'.')-1;
    name = [imageIN(1:index),'_resampled.nii'];
     name = nic_ImageResample(imageIN,name,imageREF);  %trilinear interpolation
     [Data,Head] = nic_ReadNiftiImage(name);
end

Mask = logical(zeros(size(Data)));
Mask (find (Data <= thresholdRange(2) & Data >=thresholdRange(1))) = logical(1);

% [LabelMask,NUM] = nic_spm_bwlabel(Mask,18);

[flag] = nic_WriteMaskImage(Mask,Head,imageOUT);
outputNAME=imageOUT;

return