function ROI(img, nroi)

% 1) Goal: Draw & process multiple ROIs interactively within an image.
%
% 2) Usage: ROI(img, nroi), where 'img' is your image, and 'nroi' is a
% total number to ROIs to be processed.  The opened image will be processed
% BY DEFAULT.  A prefered procedure is as follows: a) img = imread(...); b)
% imh = imshow(img) and c) ROI(img, nroi). 
% Alternatively, if there is no image in your WorkSpace, you MUST use
% square brackets to occupy the argument space for img, such as, ROI([],5),
% will let you open a new image and process with 5 ROIs.
%
% 3) Since getline('closed') is used to get the polygon interactively,
% please click left mouse button to select, and right button to finish up a
% ROI (Backspace to delete the lastest click). You may repeat this process
% till all ROIs were processed.  For more infos, please may see help
% getline. 
%
% 4) Results: ROI statistics are displayed on screen or output to a text
% file (optional).
%
% 5) The colors of line/text label are generated by 'jet' colormap,
% therefore, certain color may be too close to tell, especially when you
% select too many ROIs. In that case, you may need to edit the color after
% ROI processing. For more infos, please see help jet. 
%
% 6) For your convenience, you may use the following m-scripts to clear
% your enviroments: 'dt', 'dl' and 'df' to delete all text characters,
% lines and figures, respectively. 
%
% 7) This package includes 5 m-files: 'ROI.m', 'MultiROI', 'df.m', 'dl.m' and 'dt.m'.
%    MultiROI is similar to ROI. The only difference is that it uses spline
%    interprate to smooth the ROI edge. 
%
% Shanrong Zhang
% Department of Radiology
% University of Washington
% 02/09/2004
%
% email: zhangs@u.washington.edu

if nargin == 2
  
    imh = findobj(0, 'Type', 'Image');
    
    % if no image opened, open a new one
    if isempty(imh)
        if length(img) == 0
            [infn inpn] = uigetfile('*.*','Please select an image file');
            if infn ~= 0
                img = imread([inpn infn]);
                imh = imshow(img);
                axis image; axis off;
            else
                disp('Cancel by user!')
                return
            end
        else
            imh = imshow(img);
        end
    end

    [nrows, ncols, ncolors] = size(img);
  
    % Save ROIs to a file (optional)
    SaveIt = questdlg('Do you want to save ROI outputs ?');
    switch SaveIt
        case 'Yes'
            [outfn, outpn] = uiputfile('*', 'Select an output file');
            if outfn == 0
                disp('Cancel by user !');
            else
                fid = fopen([outpn outfn], 'w+');
                fprintf(fid, '%20s\t %-50s\n', 'Date\time = ', datestr(now));
            end
        otherwise
            outfn = 0;
    end
    
    % generate a jet colormap according to nroi
    cmap = jet(nroi);
    rndp = randperm(nroi);
    
    hold on;
    
    croi = 1;
    
    while croi <= nroi
        [x,y] =getline('closed');
        
        XData = get(imh, 'XData');
        YData = get(imh, 'YData');
        xmingrid = max( XData(1), floor(min(x)) );
        xmaxgrid = min( XData(2),  ceil(max(x)) );
        ymingrid = max( YData(1), floor(min(y)) );
        ymaxgrid = min( YData(2),  ceil(max(y)) );
        xgrid = xmingrid : xmaxgrid;
        ygrid = ymingrid : ymaxgrid;
        [X, Y] = meshgrid(xgrid, ygrid);
        k_inside = inpolygon(X, Y, x, y);
        Xin = X(k_inside);
        Yin = Y(k_inside);
        
        cdata = get(imh, 'CData');
        smallcdata = double(cdata(ygrid, xgrid, :));
        
        roi.index = croi;
        roi.area = polyarea(x,y);
        roi.center =  [mean(Xin(:)), mean(Yin(:))];
        
        for i=1:ncolors
            roicidata     = smallcdata(:, :, i);
            roi.mean(i)   =   mean(roicidata(k_inside));
            roi.std(i)    =    std(roicidata(k_inside));
            roi.min(i)    =    min(roicidata(k_inside));
            roi.max(i)    =    max(roicidata(k_inside));
            roi.median(i) = median(roicidata(k_inside));
        end;
        
        plot(x,y,'Color',cmap(rndp(croi), :));
        text(roi.center(1), roi.center(2), num2str(croi), 'Color', cmap(rndp(croi), :), 'FontWeight','Bold');
        
        % write ROI statistics into file if necessary
        if outfn ~= 0 
            fprintf(fid, '\n');
            fprintf(fid, '%20s\t %10.0f\n', 'ROI Index = ', roi.index);  
            
            fprintf(fid, '%20s\t ', 'area = '); 
            fprintf(fid, '%10.2f\t', roi.area); 
            fprintf(fid, '\n');
            
            fprintf(fid, '%20s\t ', 'mean = ');  
            fprintf(fid, '%10.2f\t', roi.mean);  
            fprintf(fid, '\n');
            
            fprintf(fid, '%20s\t ', 'std = ');  
            fprintf(fid, '%10.2f\t', roi.std);  
            fprintf(fid, '\n');  
            
            fprintf(fid, '%20s\t ', 'min = ');  
            fprintf(fid, '%10.2f\t', roi.min);  
            fprintf(fid, '\n');  
            
            fprintf(fid, '%20s\t ', 'max = ');  
            fprintf(fid, '%10.2f\t', roi.max);  
            fprintf(fid, '\n');  
            
            fprintf(fid, '%20s\t ', 'median = ');  
            fprintf(fid, '%10.2f\t', roi.median);  
            fprintf(fid, '\n');
                                
            fprintf(fid, '%20s\t ', 'roicenter = '); 
            fprintf(fid, '%10.2f\t', roi.center); 
            fprintf(fid, '\n'); 
        end

        disp(' ');
        disp(sprintf('%20s\t %10.0f', 'ROI index = ', roi.index));
        disp(sprintf('%20s\t %10.2f', 'area = ', roi.area));
        disp(sprintf('%20s\t %10.2f\t %10.2f\t %10.2f', 'mean = ', roi.mean) );  
        disp(sprintf('%20s\t %10.2f\t %10.2f\t %10.2f', 'std = ',  roi.std) );  
        disp(sprintf('%20s\t %10.2f\t %10.2f\t %10.2f', 'min = ', roi.min) );  
        disp(sprintf('%20s\t %10.2f\t %10.2f\t %10.2f', 'max = ', roi.max) );   
        disp(sprintf('%20s\t %10.2f\t %10.2f\t %10.2f', 'median = ', roi.median) );
        disp(sprintf('%20s\t %10.2f\t %10.2f\t', 'roicenter [x,y] = ', roi.center));
        disp(' ');
        
        croi = croi + 1;
        
    end
    
    if outfn == 0
        disp('Done, but ROI statistics were not saved !!!');
        disp(' ' );
    else
        disp(['Done, ROI statistics been output to ', outfn]);
        disp('But the image with ROI lines/labels was not saved yet !!!');
        fclose(fid);
    end
    
else
    disp(' ')
    disp('  Number of arguments is incorrect !!!')
    disp(' ')
    help ROI
end

% end of code
