Steve on Image Processing

Concepts, algorithms & MATLAB

This is machine translation

Translated by Microsoft
Mouseover text to see original. Click the button below to return to the English version of the page.

Feret Properties – Wrapping Up 17

Posted by Steve Eddins,

Today I want to finish up my long-running discussion of Feret diameters. (Previous posts: 29-Sep-2017, 24-Oct-2017, 20-Feb-2018, and 16-Mar-2018.)

Recall that the Feret diameter measures the extent of an object along a particular direction. The diagram below illustrates the concept. Place the object to be measured inside the jaws of a caliper, with the caliper oriented at a specified angle. Close the jaws tightly on the object while maintaining that angle. The distance between the jaws is the Feret diameter at that angle.

The maximum Feret diameter and minimum Feret diameter measure the maximum and minimum width of an object. In previous posts, I talked about how to identify antipodal vertex pairs from the set the convex hull vertices to speed up the process of finding these maximum and minimum measures. I also considered the various assumptions one can make about the shape of an individual pixel and how those assumptions can affect the results. (In the rest of this post, I'll assume a diamond pixel shape, as suggested by Cris.)

Let's look at these several measurements for a particular object.

bw = imread('shape.png');
imshow(bw)
visualizeFeretProperties(bw)

The first diagram above shows the maximum Feret diameter and its orientation. It also shows the Feret diameter at the angle that is orthogonal to the maximum diameter.

The second diagram is similar, except that it shows the minimum Feret diameter instead of the orthogonal diameter. You can see that they are not the same, and, in general, they won't be.

The third diagram shows the minimum-area bounding box, which can be found using a search procedure similar to the minimum Feret diameter. In this diagram, notice that the bounding box is not exactly aligned with the direction of the maximum Feret dimension. In general, they won't necessarily be aligned, and then the length of the minimum-area bounding box will be less than the maximum Feret dimension.

Below is a function I wrote that adds Feret diameter properties to the table returned by regionprops (using the 'table' option). It makes use of the algorithmic functions that you can find in my previous posts.

T = regionprops('table',bw,'PixelList')
T=1×1 table
        PixelList    
    _________________

    [163875×2 double]

T = feretProperties(T)
T=1×12 table
        PixelList        MaxFeretDiameter    MaxFeretDiameterEndpoints    MaxFeretDiameterOrientation    MinFeretDiameter    MinFeretDiameterTrianglePoints    MinFeretDiameterOrientation    OrthogonalDiameter    OrthogonalDiameterLowerPoints    OrthogonalDiameterUpperPoints    MinAreaBoundingBox    MinAreaBoundingBoxArea
    _________________    ________________    _________________________    ___________________________    ________________    ______________________________    ___________________________    __________________    _____________________________    _____________________________    __________________    ______________________

    [163875×2 double]         781.13               [2×2 double]                     -41.835                   317.08                  [3×2 double]                       -143.8                     331.18                  [1×2 double]                     [1×2 double]                [5×2 double]             2.4321e+05      

I hope that at least a few of you enjoyed this diversion into algorithms associated with object shape measurement.

Appendix: Functions Used

function T = feretProperties(T)
% Copyright 2017-2018 The MathWorks, Inc.

maxfd = zeros(height(T),1);
maxfd_endpoints = cell(height(T),1);
maxfd_orientation = zeros(height(T),1);

minfd = zeros(height(T),1);
minfd_triangle_points = cell(height(T),1);
minfd_orientation = zeros(height(T),1);

minod = zeros(height(T),1);
minod_lower_points = cell(height(T),1);
minod_upper_points = cell(height(T),1);

minbb = cell(height(T),1);
minbb_a = zeros(height(T),1);

for k = 1:height(T)
    pixels = T.PixelList{k};
    V = pixelHull(pixels,'diamond');
    pairs = antipodalPairs(V);
    [maxfd(k),maxfd_endpoints{k}] = maxFeretDiameter(V,pairs);
    points = maxfd_endpoints{k};
    e = points(2,:) - points(1,:);
    maxfd_orientation(k) = atan2d(e(2),e(1));
    
    [minfd(k),minfd_triangle_points{k}] = minFeretDiameter(V,pairs);
    points = minfd_triangle_points{k};
    e = points(2,:) - points(1,:);
    thetad = atan2d(e(2),e(1));
    minfd_orientation(k) = mod(thetad + 180 + 90,360) - 180;
    
    [minod(k),minod_lower_points{k},minod_upper_points{k}] = ...
        feretDiameter(V,maxfd_orientation(k)+90);
    
    [minbb{k},minbb_a(k)] = minAreaBoundingBox(V,pairs);
end

T.MaxFeretDiameter = maxfd;
T.MaxFeretDiameterEndpoints = maxfd_endpoints;
T.MaxFeretDiameterOrientation = maxfd_orientation;
T.MinFeretDiameter = minfd;
T.MinFeretDiameterTrianglePoints = minfd_triangle_points;
T.MinFeretDiameterOrientation = minfd_orientation;
T.OrthogonalDiameter = minod;
T.OrthogonalDiameterLowerPoints = minod_lower_points;
T.OrthogonalDiameterUpperPoints = minod_upper_points;
T.MinAreaBoundingBox = minbb;
T.MinAreaBoundingBoxArea = minbb_a;
end

function [bb,A] = minAreaBoundingBox(V,antipodal_pairs)
% Copyright 2017-2018 The MathWorks, Inc.

if nargin < 2
    antipodal_pairs = antipodalPairs(V);
end

n = size(antipodal_pairs,1);
p = antipodal_pairs(:,1);
q = antipodal_pairs(:,2);

A = Inf;
thetad = [];

for k = 1:n
    if k == n
        k1 = 1;
    else
        k1 = k+1;
    end
    
    pt1 = [];
    pt2 = [];
    pt3 = [];
    
    if (p(k) ~= p(k1)) && (q(k) == q(k1))
        pt1 = V(p(k),:);
        pt2 = V(p(k1),:);
        pt3 = V(q(k),:);
        
    elseif (p(k) == p(k1)) && (q(k) ~= q(k1))
        pt1 = V(q(k),:);
        pt2 = V(q(k1),:);
        pt3 = V(p(k),:);
    end
    
    if ~isempty(pt1)
        % Points pt1, pt2, and pt3 are possibly on the minimum-area
        % bounding box, with points pt1 and pt2 forming an edge coincident with
        % the bounding box. Call the height of the triangle with base
        % pt1-pt2 the height of the bounding box, h.
        
        h = triangleHeight(pt1,pt2,pt3);
        
        pt1pt2_direction = atan2d(pt2(2)-pt1(2),pt2(1)-pt1(1));
        
        w = feretDiameter(V,pt1pt2_direction);
        
        A_k = h*w;
        if (A_k < A)
            A = A_k;
            thetad = pt1pt2_direction;
        end
    end
end

% Rotate all the points so that pt1-pt2 for the minimum bounding box points
% straight up.

r = 90 - thetad;
cr = cosd(r);
sr = sind(r);
R = [cr -sr; sr cr];

Vr = V * R';

xr = Vr(:,1);
yr = Vr(:,2);
xmin = min(xr);
xmax = max(xr);
ymin = min(yr);
ymax = max(yr);

bb = [ ...
    xmin ymin
    xmax ymin
    xmax ymax
    xmin ymax
    xmin ymin ];

% Rotate the bounding box points back.
bb = bb * R;
end

function h = triangleHeight(P1,P2,P3)
% Copyright 2017-2018 The MathWorks, Inc.

h = 2 * abs(signedTriangleArea(P1,P2,P3)) / norm(P1 - P2);
end

function area = signedTriangleArea(A,B,C)
% Copyright 2017-2018 The MathWorks, Inc.

area = ( (B(1) - A(1)) * (C(2) - A(2)) - ...
    (B(2) - A(2)) * (C(1) - A(1)) ) / 2;
end

function [d,V1,V2] = feretDiameter(V,theta)
% Copyright 2017-2018 The MathWorks, Inc.

% Rotate points so that the direction of interest is vertical.

alpha = 90 - theta;

ca = cosd(alpha);
sa = sind(alpha);
R = [ca -sa; sa ca];

% Vr = (R * V')';
Vr = V * R';

y = Vr(:,2);
ymin = min(y,[],1);
ymax = max(y,[],1);

d = ymax - ymin;

if nargout > 1
    V1 = V(y == ymin,:);
    V2 = V(y == ymax,:);
end

function [d,end_points] = maxFeretDiameter(P,antipodal_pairs)
% Copyright 2017-2018 The MathWorks, Inc.

if nargin < 2
    antipodal_pairs = antipodalPairs(P);
end

v = P(antipodal_pairs(:,1),:) - P(antipodal_pairs(:,2),:);
D = hypot(v(:,1),v(:,2));
[d,idx] = max(D,[],1);
P1 = P(antipodal_pairs(idx,1),:);
P2 = P(antipodal_pairs(idx,2),:);

end_points = [P1 ; P2];
end

function [d,triangle_points] = minFeretDiameter(V,antipodal_pairs)
% Copyright 2017-2018 The MathWorks, Inc.

if nargin < 2
    antipodal_pairs = antipodalPairs(V);
end

n = size(antipodal_pairs,1);
p = antipodal_pairs(:,1);
q = antipodal_pairs(:,2);

d = Inf;
triangle_points = [];

for k = 1:n
    if k == n
        k1 = 1;
    else
        k1 = k+1;
    end
    
    pt1 = [];
    pt2 = [];
    pt3 = [];
    
    if (p(k) ~= p(k1)) && (q(k) == q(k1))
        pt1 = V(p(k),:);
        pt2 = V(p(k1),:);
        pt3 = V(q(k),:);
        
    elseif (p(k) == p(k1)) && (q(k) ~= q(k1))
        pt1 = V(q(k),:);
        pt2 = V(q(k1),:);
        pt3 = V(p(k),:);
    end
    
    if ~isempty(pt1)
        % Points pt1, pt2, and pt3 form a possible minimum Feret diameter.
        % Points pt1 and pt2 form an edge parallel to caliper direction.
        % The Feret diameter orthogonal to the pt1-pt2 edge is the height
        % of the triangle with base pt1-pt2.
        
        d_k = triangleHeight(pt1,pt2,pt3);
        if d_k < d
            d = d_k;
            triangle_points = [pt1; pt2; pt3];
        end
    end
end
end      

17 CommentsOldest to Newest

Thomas Close replied on : 1 of 17
Thank you, this was very helpful and insightful! However, there seems to be a very handy function, 'feretDiameter', with arguments of the hull and an angle, that is not provided in the appendix.
adil alazzawi replied on : 3 of 17
hi there is a missing function maxFeretDiameter I could not find it ... thanks good work
FFK replied on : 12 of 17
thanks for your efforts, sir. I'm getting an error in pixelHull() function, in the section below : Q = p + offsets; Error using + "Matrix dimensions must agree." in this line of code, "offsets" is a 1*2*4 matrix and p comes from below lines of code in feretProperties() function. - pixels = T.PixelList{k}; - v = pixelHull(pixels,'diamond'); "p" consists of x and y of true pixels in the binary image. in my example, dimensions of the "p" are 1000*2. I know that your code is correct and I traced code many times but I don't realize my mistake. I need your help to run this code correctly. I'm new in Matlab programming. please accept my sincere apologies for this simple question. thank you.
Steve Eddins replied on : 13 of 17
FFK—I'm guessing you are using a version of MATLAB prior to R2016b. Prior to that version, MATLAB was more strict about matching array sizes for operators such as plus. Try replacing that line of code with Q = bsxfun(@plus,p,offsets);
FFK replied on : 14 of 17
Thank you so much. yes, I was using R2015a version of MATLAB; and with this new line of code "Q = bsxfun(@plus,p,offsets);" the problem is solved.
faride farahnak replied on : 15 of 17
In a segment of my MS final project, I need to rotate letters in text-images (for example 2 degree), and after printing and scanning, measure the amount of rotation. For this aim, I used min area bounding box property in your code, (feretProperties method) and calculate the angle between the bottom side of the bounding box and horizon line. It should be almost 2 degrees .but in some cases like ‘o’, ‘Q’, ‘c’ and … , the degree is the value far from 2 degrees. Please guide me. Is this way (min bbx), a good solution for my aim or not? Do you know another way ( for example, another text-images property ) to help me for measuring rotation angle more accurately?
FFK FFK replied on : 16 of 17
sir, the bounding box in feret properties is not rectangular. It is a parallelogram. Is It True? I need to guide me to modify the bbx code to achieve the rectangular bbx. thank you...
Steve Eddins replied on : 17 of 17
FFK—I expect it to be rectangular. Be sure to call "axis equal" so that the axes' DataAspectRatio is [1 1 1]. Otherwise, a typical axes' aspect ratio will stretch out the x-axis compared to the y-axis, which will make a rectangle look like a parallelogram.

Add A Comment

Your email address will not be published. Required fields are marked *

Preview: hide