imfeature.m

有关matlab的电子书籍有一定的帮助希望有用

            r3 = min(tmp);
            r4 = max(tmp);
            c3 = maxC;
            c4 = maxC;

            % Points 5 and 6 are on the bottom row.
            r5 = maxR;
            r6 = maxR;
            % Find the minimum and maximum column coordinates for
            % bottom-row pixels.
            tmp = c(maxRSet);
            c5 = max(tmp);
            c6 = min(tmp);
            
            % Points 7 and 8 are on the left column.
            % Find the minimum and maximum row coordinates for
            % left-column pixels.
            tmp = r(minCSet);
            r7 = max(tmp);
            r8 = min(tmp);
            c7 = minC;
            c8 = minC;
            
            stats(k).Extrema = [c1-0.5 r1-0.5
                c2+0.5 r2-0.5
                c3+0.5 r3-0.5
                c4+0.5 r4+0.5
                c5+0.5 r5+0.5
                c6-0.5 r6+0.5
                c7-0.5 r7+0.5
                c8-0.5 r8-0.5];
        end
    end
    
end
        
%%%
%%% ComputeBoundingBox
%%%
function [stats, computedStats] = ComputeBoundingBox(L, stats, computedStats)
%   [minC minR width height]; minC and minR end in .5.

if ~computedStats.BoundingBox
    computedStats.BoundingBox = 1;
    
    [stats, computedStats] = ComputePixelList(L, stats, computedStats);
    
    for k = 1:length(stats)
        list = stats(k).PixelList;
        if (isempty(list))
            stats(k).BoundingBox = [0.5 0.5 0 0];
        else
            rows = list(:,2);
            cols = list(:,1);
            minR = min(rows) - 0.5;
            maxR = max(rows) + 0.5;
            minC = min(cols) - 0.5;
            maxC = max(cols) + 0.5;
            stats(k).BoundingBox = [minC minR (maxC-minC) (maxR-minR)];
        end
    end
end

%%%
%%% ComputeEllipseParams
%%%
function [stats, computedStats] = ComputeEllipseParams(L, stats, ...
        computedStats)
%   Find the ellipse that has the same 2nd-order moments as the
%   region.  Compute the axes lengths, orientation, and
%   eccentricity of the ellipse.
%   Ref: Haralick and Shapiro, Computer and Robot Vision
%   vol I, Addison-Wesley 1992, Appendix A.


if ~(computedStats.MajorAxisLength & computedStats.MinorAxisLength & ...
            computedStats.Orientation & computedStats.Eccentricity)
    computedStats.MajorAxisLength = 1;
    computedStats.MinorAxisLength = 1;
    computedStats.Eccentricity = 1;
    computedStats.Orientation = 1;
    
    [stats, computedStats] = ComputePixelList(L, stats, computedStats);
    [stats, computedStats] = ComputeCentroid(L, stats, computedStats);
    
    for k = 1:length(stats)
        list = stats(k).PixelList;
        if (isempty(list))
            stats(k).MajorAxisLength = 0;
            stats(k).MinorAxisLength = 0;
            stats(k).Eccentricity = 0;
            stats(k).Orientation = 0;
            
        else
            % Assign X and Y variables so that we're measuring orientation
            % counterclockwise from the horizontal axis.
            
            xbar = stats(k).Centroid(1);
            ybar = stats(k).Centroid(2);
            
            x = list(:,1) - xbar;
            y = -(list(:,2) - ybar);
            
            N = length(x);
            uxx = sum(x.^2)/N + 1/12;
            uyy = sum(y.^2)/N + 1/12;
            uxy = sum(x.*y)/N;
            
            common = sqrt((uxx - uyy)^2 + 4*uxy^2);
            stats(k).MajorAxisLength = 2*sqrt(2)*sqrt(uxx + uyy + common);
            stats(k).MinorAxisLength = 2*sqrt(2)*sqrt(uxx + uyy - common);
            stats(k).Eccentricity = 2*sqrt((stats(k).MajorAxisLength/2)^2 - ...
                    (stats(k).MinorAxisLength/2)^2) / ...
                    stats(k).MajorAxisLength;
            
            if (uyy > uxx)
                stats(k).Orientation = (180/pi)*atan2(uyy - uxx + ...
                        sqrt((uyy - uxx)^2 + 4*uxy^2), 2*uxy);
            else
                stats(k).Orientation = (180/pi)*atan2(2*uxy, ...
                        uxx - uyy + sqrt((uxx - uyy)^2 + 4*uxy^2));
            end
        end
    end
    
end
    
%%%
%%% ComputeSolidity
%%%
function [stats, computedStats] = ComputeSolidity(L, stats, computedStats)
%   Area / ConvexArea

if ~computedStats.Solidity
    computedStats.Solidity = 1;
    
    [stats, computedStats] = ComputeArea(L, stats, computedStats);
    [stats, computedStats] = ComputeConvexArea(L, stats, computedStats);
    
    for k = 1:length(stats)
        if (stats(k).ConvexArea == 0)
            stats(k).Solidity = NaN;
        else
            stats(k).Solidity = stats(k).Area / stats(k).ConvexArea;
        end
    end
end

%%%
%%% ComputeExtent
%%%
function [stats, computedStats] = ComputeExtent(L, stats, computedStats)
%   Area / (BoundingBox(3) * BoundingBox(4))

if ~computedStats.Extent
    computedStats.Extent = 1;
    
    [stats, computedStats] = ComputeArea(L, stats, computedStats);
    [stats, computedStats] = ComputeBoundingBox(L, stats, computedStats);
    
    for k = 1:length(stats)
        if (stats(k).Area == 0)
            stats(k).Extent = NaN;
        else
            stats(k).Extent = stats(k).Area / prod(stats(k).BoundingBox(3:4));
        end
    end
end

%%%
%%% ComputeImage
%%%
function [stats, computedStats] = ComputeImage(L, stats, computedStats)
%   Binary image containing "on" pixels corresponding to pixels
%   belonging to the region.  The size of the image corresponds
%   to the size of the bounding box for each region.

if ~computedStats.Image
    computedStats.Image = 1;

    [stats, computedStats] = ComputePixelList(L, stats, computedStats);
    [stats, computedStats] = ComputeBoundingBox(L, stats, computedStats);

    for k = 1:length(stats)
        M = stats(k).BoundingBox(4);
        N = stats(k).BoundingBox(3);
        firstRow = stats(k).BoundingBox(2) + 0.5;
        firstCol = stats(k).BoundingBox(1) + 0.5;
        
        stats(k).Image = uint8(zeros(M,N));
        
        r = stats(k).PixelList(:,2);
        c = stats(k).PixelList(:,1);
        r = r - firstRow + 1;
        c = c - firstCol + 1;
        
        idx = M * (c - 1) + r;
        stats(k).Image(idx) = 1;
        stats(k).Image = logical(stats(k).Image);
    end
end

%%%
%%% ComputePixelList
%%%
function [stats, computedStats] = ComputePixelList(L, stats, computedStats)
%   A P-by-2 matrix, where P is the number of pixels belonging to
%   the region.  Each row contains the row and column
%   coordinates of a pixel.

if ~computedStats.PixelList
    computedStats.PixelList = 1;
    
    % Form a sparse matrix containing one column per region.  In
    % column P, the location of nonzero values correspond to the
    % linear indices of pixels in L that have value P.  For
    % example, S(100,5) is nonzero if and only L(100) equals 5.
    idx = find(L);
    elementValues = L(idx);
    S = sparse(idx, elementValues, 1);

    % Loop over each column of the sparse matrix.  Finding the
    % row indices of the nonzero entries in S(:,P) is equivalent
    % to finding the linear indices of pixels in L that equal P.
    % Convert the linear indices to row-column indices and store
    % the results in the pixel list.
    [M,N] = size(L);
    for k = 1:length(stats)
        idx = find(S(:,k)) - 1;
        c = floor(idx/M) + 1;
        r = rem(idx,M) + 1;
        stats(k).PixelList = [c(:) r(:)];
    end
end

%%%
%%% ComputePerimeterCornerPixelList
%%%
function [stats, computedStats] = ComputePerimeterCornerPixelList(L, ...
        stats, computedStats)
%   Find the pixels on the perimeter of the region; make a list
%   of the coordinates of their corners; sort and remove
%   duplicates.

if ~computedStats.PerimeterCornerPixelList
    computedStats.PerimeterCornerPixelList = 1;

    [stats, computedStats] = ComputeImage(L, stats, computedStats);
    [stats, computedStats] = ComputeBoundingBox(L, stats, computedStats);

    for k = 1:length(stats)
        perimImage = bwmorph(stats(k).Image, 'perim8');
        firstRow = stats(k).BoundingBox(2) + 0.5;
        firstCol = stats(k).BoundingBox(1) + 0.5;
        [r,c] = find(perimImage);
        % Force rectangular empties.
        r = r(:) + firstRow - 1;
        c = c(:) + firstCol - 1;
        r1 = r - 0.5;
        r2 = r + 0.5;
        c1 = c - 0.5;
        c2 = c + 0.5;
        rr = [r1 ; r1 ; r2 ; r2];
        cc = [c1 ; c2 ; c1 ; c2];
        
        % The r-c order is reversed in the call below so that the
        % result will be suitable for passing to delaunay with the
        % 'sorted' option.
        ccrr = unique([cc rr], 'rows');
        
        stats(k).PerimeterCornerPixelList = ccrr;
    end
    
end

%%%
%%% ComputeConvexHull
%%%
function [stats, computedStats] = ComputeConvexHull(L, stats, computedStats)
%   A P-by-2 array representing the convex hull of the region.
%   The first column contains row coordinates; the second column
%   contains column coordinates.  The resulting polygon goes
%   through pixel corners, not pixel centers.

if ~computedStats.ConvexHull
    computedStats.ConvexHull = 1;

    [stats, computedStats] = ComputePerimeterCornerPixelList(L, stats, ...
            computedStats);
    [stats, computedStats] = ComputeBoundingBox(L, stats, computedStats);

    for k = 1:length(stats)
        list = stats(k).PerimeterCornerPixelList;
        if (isempty(list))
            stats(k).ConvexHull = zeros(0,2);
        else
            rr = list(:,2);
            cc = list(:,1);
            tri = delaunay(rr, cc, 'sorted');
            hullIdx = convhull(rr, cc, tri);
            stats(k).ConvexHull = list(hullIdx,:);
        end
    end
end

%%%
%%% ParseInputs
%%%
function [L,reqStats,n,msg] = ParseInputs(officialStats, varargin)

L = [];
reqStats = [];
n = 8;
msg = '';

if (nargin < 1)
    msg = 'Too few input arguments.';
    return
end

if (nargin >= 2)
    % Check for trailing connectivity argument; strip it off if present.
    if (isequal(varargin{end},4) | isequal(varargin{end},8))
        n = varargin{end};
        varargin(end) = [];
    end
end
   
L = varargin{1};
if (~isnumeric(L))
    msg = 'Invalid input.';
    return;
end

list = varargin(2:end);
if (~isempty(list) & ~iscell(list{1}) & strcmp(lower(list{1}), 'all'))
    reqStats = officialStats;
    reqStatsIdx = 1:length(officialStats);
    
elseif (isempty(list) | (~iscell(list{1}) & strcmp(lower(list{1}),'basic')))
    % Default list
    reqStats = {'Area'
                'Centroid'
                'BoundingBox'};
else
    
    if (iscell(list{1}))
        list = list{1};
    end
    list = list(:);

    officialStatsL = lower(officialStats);

    reqStatsIdx = [];
    for k = 1:length(list)
        if (~isstr(list{k}))
            msg = 'Invalid input argument.';
            return;
        end
        
        idx = strmatch(lower(list{k}), officialStatsL);
        if (isempty(idx))
            msg = sprintf('Unknown measurement: "%s".', list{k});
            return;
        
        elseif (length(idx) > 1)
            msg = sprintf('Ambiguous measurement: "%s".', list{k});
            return;
            
        else
            reqStatsIdx = [reqStatsIdx; idx];
        end
    end
    
    reqStats = officialStats(reqStatsIdx);
end