reduce side-count of a polygon

Question

Suppose we have image of a simple graphic, and we know it is a polygon, slightly distorted. Is there an image processing way to approximate the original parameters of the graphic object?

The matrix below was created by code and then reduced in size to show the fifth area of interest:

EnumeratedMask=bwlabel(Zdata<-.06);

0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
0   0   0   0   5   5   5   0   0   0   0   0   0   0   0   0   0   0   0
0   0   0   0   5   5   5   5   0   0   0   0   0   0   0   0   0   0   0
0   0   0   0   5   5   5   5   5   5   5   0   0   0   0   0   0   0   0
0   0   0   5   5   5   5   5   5   5   5   5   0   0   0   0   0   0   0
0   0   0   5   5   5   5   5   5   5   5   5   5   5   0   0   0   0   0
0   0   0   5   5   5   5   5   5   5   5   5   5   5   5   5   0   0   0
0   0   0   5   5   5   5   5   5   5   5   5   5   5   5   5   5   0   0
0   0   0   5   5   5   5   5   5   5   5   5   5   5   5   0   0   0   0
0   0   0   5   5   5   5   5   5   5   5   5   0   0   0   0   0   0   0
0   0   0   5   5   5   5   5   5   5   5   0   0   0   0   0   0   0   0
0   0   0   0   5   5   5   5   5   5   5   0   0   0   0   0   0   0   0
0   0   0   0   5   5   5   5   5   5   5   0   0   0   0   0   0   0   0
0   0   0   5   5   5   5   5   5   5   0   0   0   0   0   0   0   0   0
0   0   0   5   5   5   5   0   0   0   0   0   0   0   0   0   0   0   0
0   0   0   5   5   5   5   0   0   0   0   0   0   0   0   0   0   0   0
0   0   0   5   5   5   5   5   0   0   0   0   0   0   0   0   0   0   0
0   0   0   5   5   0   0   0   0   0   0   0   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0

For next step I need the ABC/ABCD coordinates to get z-profile across lines further defined by those points.

[Maybe this?](https://en.m.wikipedia.org/wiki/Ramer–Douglas–Peucker_algorithm) — Cris Luengo, Jun 18 '19 at 02:33
I think some of what's done [here](https://stackoverflow.com/q/41988945/52738) and [here](https://stackoverflow.com/q/45073798/52738) might be related/useful, at least in regard to the initial steps needed. — gnovice, Jun 19 '19 at 15:04
@gnovice the links you have posted address just the first part of the simplification. Cris' comment gave a much more useful tool. — Crowley, Jun 19 '19 at 15:08

Patrick Happel · Accepted Answer · 2019-06-19T09:40:46.513

Here is an implemetation of the Ramer–Douglas–Peucker algorithm as suggested in the comment by Cris Luengo above.

This is a full edit of the first version of the answer, which used edge to find the boundaries of the object. As Cris Luengo pointed out in a comment, bwboundaries is a better choice for binary images. bwboundaries returns the points sorted, which simplifies the code a lot.

The following code does the following:

1) find the edges of your object using bwboundaries. They are already sorted.

2) Use the Ramer–Douglas–Pecker algorithm to simplify the point list.

Since I needed some visual cues for debugging, the code opens a figure that shows what is going on.

Please note that the code is far from being properly tested.

img = [...
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 5 5 5 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0
    0 0 0 5 5 5 5 5 5 5 5 5 0 0 0 0 0 0 0
    0 0 0 5 5 5 5 5 5 5 5 5 5 5 0 0 0 0 0
    0 0 0 5 5 5 5 5 5 5 5 5 5 5 5 5 0 0 0
    0 0 0 5 5 5 5 5 5 5 5 5 5 5 5 5 5 0 0
    0 0 0 5 5 5 5 5 5 5 5 5 5 5 5 0 0 0 0
    0 0 0 5 5 5 5 5 5 5 5 5 0 0 0 0 0 0 0
    0 0 0 5 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0
    0 0 0 0 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0
    0 0 0 0 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0
    0 0 0 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0 0
    0 0 0 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 5 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 5 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];

watch = true;

if watch
    f = figure;
    subplot(1,2,1);
    imagesc(img);
end


% first, find the edges

b = bwboundaries(img);

% note that the first and the last point are the same.
% However they are already sorted.
x = b{1}(:,1);
y = b{1}(:,2);

edges = zeros(size(img));
edges(sub2ind(size(e), x,y)) = 1;

if watch
    ax = subplot(1,2,2);

    img_h = imagesc(edges);
    hold on;

end 

title('Performing Douglas-Peucker algorithm');

% Omit the last point for the algorithm.

points = l_DouglasPeucker( [x(1:end-1), y(1:end-1)], 1, watch);

title('Final result');
plot([points(:,2); points(1,2)], [points(:,1); points(1,1)]);



function res = l_DouglasPeucker( points, ep, watch )

    if nargin < 3
        watch = false;
    end

    if watch
        subplot(1,2,2);
        hold on;
        hp = plot(points(:,2), points(:,1), 'ko-');
        hp2 = plot([points(1,2) points(end,2)], [points(1,1) points(end,1)], 'r-');
    end
    distances = zeros(size(points,1),1);
    for i = 1:size(points,1)
        distances(i) = l_distance_to_line(points(1,:), points(end,:), points(i,:));
    end

    idx = find(distances == max(distances),1);


    if watch
        hp4 = plot(points(idx,2), points(idx,1), 'mo', 'MarkerFaceColor', [1,1,1]);
        pause(.5);
        delete(hp);
        delete(hp2);
        delete(hp4);
    end

    if max(distances) > ep
        res = [l_DouglasPeucker(points(1:idx,:), ep, watch); l_DouglasPeucker(points(idx:end,:), ep, watch)];
    else
        res = [points(1,:); points(end,:)];
    end


end

function d = l_distance_to_line(p1,p2,p)
% computes the distance of p to the line through by p1 and p2
% There might be much better implementations of this...

% we need 3-dimensional data for this

p1 = [p1(1), p1(2), 0];
p2 = [p2(1), p2(2), 0];
p = [p(1,1) p(1,2) 0];

a = p2 - p1;
b = p - p1;

d = norm(cross(a,b)) ./ norm(a);
end

`edge` is to detect edges in a grey-value image. We’re dealing with a binary image here, no detection is necessary. Just use [`bwboundaries`](https://www.mathworks.com/help/images/ref/bwboundaries.html) to obtain a polygon of the outline. — Cris Luengo, Jun 19 '19 at 01:50

Crowley · Answer 2 · 2019-06-20T16:15:18.593

I have started to work before Patrick posted his educational answer and I have faced one "issue" with the Ramer-Douglas-Peucker algorithm: It, by definition, keeps the first and last points. Both convhull and boundary functions starts somewhere and not always in the corner, which would trigger false positive. Third and fourth steps address this issue - The other point is more likely to be the true corner.

Algorithm:

detect convex hull/outer points.
Apply the Ramer-Douglas-Peucker algorithm (RDP) with tighter fit
Rewind the loop
Apply RDP algorithm to remove false positive corner with looser fit

Code:

eps1=2;
eps2=4;
img=[...
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 5 5 5 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0
0 0 0 5 5 5 5 5 5 5 5 5 0 0 0 0 0 0 0
0 0 0 5 5 5 5 5 5 5 5 5 5 5 0 0 0 0 0
0 0 0 5 5 5 5 5 5 5 5 5 5 5 5 5 0 0 0
0 0 0 5 5 5 5 5 5 5 5 5 5 5 5 5 5 0 0
0 0 0 5 5 5 5 5 5 5 5 5 5 5 5 0 0 0 0
0 0 0 5 5 5 5 5 5 5 5 5 0 0 0 0 0 0 0
0 0 0 5 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0
0 0 0 0 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0
0 0 0 0 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0
0 0 0 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0 0
0 0 0 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 5 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0
0 0 0 5 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];

[X,Y]=find(matr==5); % // find image coordinates of points of interest
is=boundary(X,Y,1);  % // find points on the boundary

xy=[X(is),Y(is)];    % // boundary points (step 1)

% // apply RDP algorithm (step 2)
xy=RDP(xy,eps1);

% // rewind the loop (step 3)
xy=xy([2:end-1 1:2],:);

% // apply the DRP algorithm (step 4)
xy=RDP(xy,eps2);

function[hull]=RDP(hull,eps)
    sp=hull(1,:);
    ep=hull(end,:);
    ip=hull(2:end-1,:);
    % // calculate distances of inner points from the first-last line
    dst=PerpDist(sp,ep,ip);
    % // find the point furthest from the f-l line
    [mx,mi]=max(dst);
    if mx>eps % // furthest point does not fit in - split and apply RDP recursively
        lp=[sp;ip(1:mi,:)];
        if size(lp,1)>2 % // there are points left to assess
            lp=RDP(lp,eps);
        end
        rp=[ip(mi:end,:);ep];
        if size(rp,1)>2 % // there are points left to assess
            rp=RDP(rp,eps);
        end
        hull=[lp;rp(2:end,:)]; % // concatenate the branches
    else % // the furthest poit fits in the limit, drop all inner points
        hull=[sp;ep];
    end
end

function[D]=PerpDist(A,B,C)
    D=nan(size(C,1),1);
    if A==B % // edge is defined by one point, use euclidean distance between points
        for PDi=1:size(C,1)
            D(PDi)=norm(C(PDi,:)-A);
        end
    else % // edge is a line, use eucleidean distance from a line
        for PDi=1:size(C,1)
            D(PDi)=abs(A(1)*(B(2)-C(PDi,2))+B(1)*(C(PDi,2)-A(2))+C(PDi,1)*(A(2)-B(2)))/norm(B-A);
        end
    end
end

Dots: points within the shape.
Red square: First and last point returned from boundary function.
green line: Result from first RDP simlification.
Magenta line: final triangle (original shape was triangle).

Edits:

Leaving the idea of finding two furthest points to split the loop because it triggered false positives even for triangles. Starting from arbitrary point and using RDP twice do the trick.

Credits:

Steve Eddins for max_feret diameter
Cris Luengo for Ramer-Douglas-Peucker algorithm comment

If you know the original shape (triangle in your case) and just want to find its vertices, another way would be to fit a triangle to the data. — Patrick Happel, Jun 19 '19 at 16:29
@PatrickHappel In this case it was a triangle, more generally it can be triangle or tetragon. In general I aimed for "simpler" polygons. — Crowley, Jun 20 '19 at 12:20

reduce side-count of a polygon

2 Answers2