图像上绘制和处理多个roi

【核心代码】

function roi = multiROI(img, nroi)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1) Goal: Draw & process multiple ROIs interactively on an image.                      %
%%                                                                                       %
%% 2) Usage: multiROI(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.     %
%%    You may used img = imread(...) to read an image into your WorkSpace, show it       %
%%    by imshow(img), then use multiROI(img, nroi) to process it;                        %
%%    Alternatively, if there is no image in your WorkSpace, you MUST use square         %
%%    brackets to occupy the argument space for img, for example, multiROI([],nroi)      %
%%    let you open a new image and process correspondingly.                              %
%%                                                                                       %
%% 3) Results: ROI statistics are displayed on screen or output to file (optional).      %
%%                                                                                       %
%% 4) Notes: Since ROI was drawn by Spline interpolation, it is desirable to have        %
%%    more data points at around the sharp corner region; The line/label color was       %
%%    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 process.                                                       %
%%                                                                                       %
%% Shanrong Zhang                                                                        %
%% Department of Radiology                                                               %
%% University of washington                                                              %
%% 02/05/2004                                                                            %
%%                                                                                       %
%% email: zhangs@u.washington.edu                                              %
%%                                                                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin == 2
  
    % get current image handle 
    fighandle = findobj(0, 'Type', 'Figure'); 
    imhandle = findobj(0, 'Type', 'Image');
    
    oldpathname = pwd;

    % if no image opened, open a new one
    if isempty(imhandle)
        if length(img) == 0
            [filename pathname] = uigetfile('*.*','Please select an image file');
            if filename ~= 0
                % cd(pathname);
                img = imread([pathname filename]);
                imhandle = imshow(img);
                cd(oldpathname);
            else
                disp('Cancel by user!')
                return
            end
        else
            imhandle = imshow(img);
        end
    end

    axishandle = gca;
    [nrows, ncols, ncolors] = size(img);
  
    % Save ROIs to a file (optional)
    [outfilename, outpathname] = uiputfile('*', 'Select an output file');
    if outfilename == 0
        disp('ROI results were not saved!');
    else
        fid = fopen([outpathname outfilename], 'w ');
        fprintf(fid, '%20s\t %-50s\n', 'Date\time = ', datestr(now));
    end
    
    % generate a colormap according to number of ROIS
    cmap = jet(nroi);
    rndp = randperm(nroi);

    croi = 1;

    % loop untill total number of ROIs been processed
    while croi <= nroi
        
        x = [];
        y = [];
        
        color = cmap(rndp(croi), :);
        
        if (~isempty(x) & ~isempty(y))
            nx = length(x);
            ny = length(y);
            if (nx<3) | (ny<3)
                disp('  ROI size is too small !')
                return
                x = x(1 : min(nx, ny));
                y = y(1 : min(nx, ny));  
                x( x < 0.5 ) = 0.5;
                x( x > ncols 0.5 ) = ncols   0.5;
                y( y < 0.5 ) = 0.5;
                y( y > nrows 0.5 ) = nrows   0.5;
                hold on; 
                linehandle = plot(x, y, 'Color', color);
            end
        else
            % Get the ROI interactively
            [x , y, linehandle] = getpoints(axishandle, color);
        end;
        
        %Calculate ROI area
        n = length(x);
        diffx = [diff(x) (x(1) - x(n))];
        diffy = [diff(y) (y(1) - y(n))];
        avector = y .* diffx   diffx .* diffy ./2;
        
        % Copy area, vectors and linehandle to roi stucture
        roi.label = croi;
        roi.apix = abs(sum(avector));
        roi.x = x;
        roi.y = y; 
        roi.linehandle = linehandle;
        
        % Change the pointer to something that is familiar to Microsoft users...
        oldpointershape = get(fighandle, 'Pointer');
        set(fighandle, 'Pointer', 'watch');
        
        %Calculate the ROI area in square point units
        XData = get(imhandle, 'XData'); 
        YData = get(imhandle, 'YData');
        pixarea = (diff(XData)   1) * (diff(YData)   1);
        
        % Create the smallest rectangular grid around the ROI
        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);
        mask = zeros(nrows, ncols);
        mask(ygrid, xgrid) = 1;
        
        cdata = get(imhandle, 'CData');
        smallcdata = double(cdata(ygrid, xgrid, :));
        [m, n, ncolors] = size(smallcdata);
        
        % Analyze only the points in the polygon
        k_inside = inpolygon(X, Y, x, y);
        Xin = X(k_inside);
        Yin = Y(k_inside);
        clear X Y
        
        % Determine the center of the polygon and label ROI
        roi.center =  [mean(Xin(:)), mean(Yin(:))];
        text(roi.center(1), roi.center(2), num2str(croi), 'Color', color, 'FontWeight', 'Bold');
        clear Xin Yin
        
        % Calculate the mean, SD, etc... and as fields add to roi structure for each color
        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;
        
        % Add the date and time for future reference (in files)
        roi.timestamp = datestr(now);
        
        % Reset pointer shape
        set(fighandle, 'Pointer', oldpointershape);
        
        % write ROI statistics into file if necessary
        if outfilename ~= 0 
            fprintf(fid, '\n');
            fprintf(fid, '%20s\t %10.0f\n', 'ROI label = ', roi.label);  
            fprintf(fid, '%20s\t %10.2f\n', 'pix area = ', roi.apix);  
            
            fprintf(fid, '%20s\t ', 'roicenter = '); 
            fprintf(fid, '%10.2f\t', roi.center); 
            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');
        end
        
        % dispplay each ROI statistics on screen
        disp(' ');
        disp(fprintf('%20s\t %10.0f', 'ROI label = ', roi.label) );
        disp(fprintf('%20s\t %10.2f', 'pix area = ', roi.apix) );  
        disp(sprintf('%20s\t %10.2f\t %10.2f\t', 'roicenter [x,y] = ', roi.center)); 
        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(' ');
        
        croi = croi   1;
        
    end
    
    if outfilename ~= 0
        disp(['ROI statistics saved into file ', outfilename]);
        disp('But the image with ROIs was not !!!');
        fclose(fid);
    end
    
else
    disp(' ')
    disp('  Number of arguments is incorrect!')
    disp(' ')
    help multiROI
end

%
% LOCAL FUNCTION GETPOINTS
%
function [xs, ys, linehandle] = getpoints(axishandle, color)

% Find parent figure for the argument axishandle
axes(axishandle);
figure(get(axishandle, 'Parent'));

% Change pointer shape
oldpointershape = get(gcf, 'Pointer');

ptrc =  ones(16)   1;
ptrc( 1, :) = 1; 
ptrc(16, :) = 1; 
ptrc(: , 1) = 1; 
ptrc(: ,16) = 1; 
ptrc(1:4,8:9) = 1;
ptrc(8:9,1:4) = 1;
ptrc(13:16, 8:9 ) = 1;
ptrc( 8:9 ,13:16) = 1;
ptrc(5:12, 5:12) = NaN;
set(gcf,'Pointer', 'custom',...
 	    'PointerShapeCData', ptrc,...
 	    'PointerShapeHotSpot', [8 8]);

% Prepare for interactive collection of ROI boundary points
hold on
pointhandles = [];
xpts = [];
ypts = [];
splinehandle = [];
n = 0;
but = 1;
BUTN = 0;
KEYB = 1;
done = 0;

% Loop until right hand mouse button or keayboard is pressed
while ~done;  
  
    % Analyze each buttonpressed event
    keyb_or_butn = waitforbuttonpress;
    if keyb_or_butn == BUTN;
        currpt = get(axishandle, 'CurrentPoint');
        seltype = get(gcf, 'SelectionType');
        switch seltype
            case 'normal',
                but = 1;
            case 'alt',
                but = 2;
            otherwise,
                but = 2;
        end;
    elseif keyb_or_butn == KEYB
        but = 2;
    end;
    
    % Get coordinates of the last buttonpressed event
    xi = currpt(2, 1);
    yi = currpt(2, 2);
    
    % Start a spline throught the points or
    % update the line through the points with a new spline
    
    if but == 1
        if ~isempty(splinehandle)
            delete(splinehandle);
        end;
        pointhandles(n 1) = plot(xi, yi, 'Color', color, 'Marker', 'o');
        n = n   1;
        xpts(n, 1) = xi;
        ypts(n, 1) = yi;
        
        % Draw a spline line through the points
        if n > 1
            t = 1:n;
            ts = 1: 0.1 : n;
            xs = spline(t, xpts, ts);
            ys = spline(t, ypts, ts);
            splinehandle = plot(xs, ys, 'Color', color);
        end;
    elseif but > 1
        % Exit for right hand mouse button or keyboard input
        done = 1;
    end;
end;

% Add first point to the end of the vector for spline 
xpts(n 1, 1) = xpts(1, 1);
ypts(n 1, 1) = ypts(1, 1);

% (re)draw the final spline 
if ~ isempty(splinehandle)
    delete(splinehandle);
end;    

t = 1:n 1;
ts = 1: 0.25 : n 1;
xs = spline(t, xpts, ts);
ys = spline(t, ypts, ts);

linehandle = plot(xs, ys, 'Color', color); 
drawnow;

% Delete the point markers 
if ~isempty(pointhandles)
    delete(pointhandles)
end;

% Reset pointershape 
set(gcf, 'Pointer', oldpointershape);

% END OF LOCAL FUNCTION GETPOINTS 
%