function [lower,upper]=threshold_image(I, cutoff)
% [lower_image,upper_image]=threshold_image(I, cutoff)
%
% I      : input image. It's either 2D or 3D images
% cutoff : valueof threshold
%
% if will return thresholded images seperated by cutoff
%
% Moo K. Chung (c)
% mkchung@wisc.edu
% University of Wisconsin-Madison
%
% Update History:
% Dec 25, 2003 
% Nov 17, 2007 
% Feb 14, 2010 
% Oct 10, 2010


if size(I,3)== 1
    d=size(I);


    dim1=d(1);dim2=d(2); dim3=1;


    temp=reshape(I,dim1*dim2*dim3,1);
    lower_image_ind=find(temp < cutoff);
    upper_image_ind=find(temp > cutoff);

    lower_image=zeros(dim1*dim2*dim3,1)+ cutoff;
    upper_image=zeros(dim1*dim2*dim3,1)+ cutoff; % constant image

    lower_image(lower_image_ind)=temp(lower_image_ind);
    upper_image(upper_image_ind)=temp(upper_image_ind);

    lower=reshape(lower_image,dim1,dim2);
    upper=reshape(upper_image,dim1,dim2);

else
    d=size(I);


    dim1=d(1);dim2=d(2); dim3=d(3);


    temp=reshape(I,dim1*dim2*dim3,1);
    lower_image_ind=find(temp < cutoff);
    upper_image_ind=find(temp > cutoff);

    lower_image=zeros(dim1*dim2*dim3,1)+ cutoff;
    upper_image=zeros(dim1*dim2*dim3,1)+ cutoff;

    lower_image(lower_image_ind)=temp(lower_image_ind);
    upper_image(upper_image_ind)=temp(upper_image_ind);

    lower=reshape(lower_image,dim1,dim2,dim3);
    upper=reshape(upper_image,dim1,dim2,dim3);


end;
    


