How to crop away convexity defects?

Question

I'm trying to detect and fine-locate some objects in images from contours. The contours that I get often include some noise (maybe form the background, I don't know). The objects should look similar to rectangles or squares like:

I get very good results with shape matching (cv::matchShapes) to detect contours with those objects in them, with and without noise, but I have problems with the fine-location in case of noise.

Noise looks like:

or for example.

My idea was to find convexity defects and if they become too strong, somehow crop away the part that leads to concavity. Detecting the defects is ok, typically I get two defects per "unwanted structure", but I'm stuck on how to decide what and where I should remove points from the contours.

Here are some contours, their masks (so you can extract the contours easily) and the convex hull including thresholded convexity defects:

Could I just walk through the contour and locally decide whether a "left turn" is performed by the contour (if walking clockwise) and if so, remove contour points until the next left turn is taken? Maybe starting at a convexity defect?

I'm looking for algorithms or code, programming language should not be important, algorithm is more important.

Answer 1

This approach works only on points. You don't need to create masks for this.

The main idea is:

1. Find defects on contour
2. If I find at least two defects, find the two closest defects
3. Remove from the contour the points between the two closest defects
4. Restart from 1 on the new contour

I get the following results. As you can see, it has some drawbacks for smooth defects (e.g. 7th image), but works pretty good for clearly visible defects. I don't know if this will solve your problem, but can be a starting point. In practice should be quite fast (you can surely optimize the code below, specially the removeFromContour function). Also, the only parameter of this approach is the amount of the convexity defect, so it works well with both small and big defecting blobs.

#include <opencv2/opencv.hpp>
using namespace cv;
using namespace std;

int ed2(const Point& lhs, const Point& rhs)
{
    return (lhs.x - rhs.x)*(lhs.x - rhs.x) + (lhs.y - rhs.y)*(lhs.y - rhs.y);
}

vector<Point> removeFromContour(const vector<Point>& contour, const vector<int>& defectsIdx)
{
    int minDist = INT_MAX;
    int startIdx;
    int endIdx;

    // Find nearest defects
    for (int i = 0; i < defectsIdx.size(); ++i)
    {
        for (int j = i + 1; j < defectsIdx.size(); ++j)
        {
            float dist = ed2(contour[defectsIdx[i]], contour[defectsIdx[j]]);
            if (minDist > dist)
            {
                minDist = dist;
                startIdx = defectsIdx[i];
                endIdx = defectsIdx[j];
            }
        }
    }

    // Check if intervals are swapped
    if (startIdx <= endIdx)
    {
        int len1 = endIdx - startIdx;
        int len2 = contour.size() - endIdx + startIdx;
        if (len2 < len1)
        {
            swap(startIdx, endIdx);
        }
    }
    else
    {
        int len1 = startIdx - endIdx;
        int len2 = contour.size() - startIdx + endIdx;
        if (len1 < len2)
        {
            swap(startIdx, endIdx);
        }
    }

    // Remove unwanted points
    vector<Point> out;
    if (startIdx <= endIdx)
    {
        out.insert(out.end(), contour.begin(), contour.begin() + startIdx);
        out.insert(out.end(), contour.begin() + endIdx, contour.end());
    } 
    else
    {
        out.insert(out.end(), contour.begin() + endIdx, contour.begin() + startIdx);
    }

    return out;
}

int main()
{
    Mat1b img = imread("path_to_mask", IMREAD_GRAYSCALE);

    Mat3b out;
    cvtColor(img, out, COLOR_GRAY2BGR);

    vector<vector<Point>> contours;
    findContours(img.clone(), contours, RETR_EXTERNAL, CHAIN_APPROX_NONE);

    vector<Point> pts = contours[0];

    vector<int> hullIdx;
    convexHull(pts, hullIdx, false);

    vector<Vec4i> defects;
    convexityDefects(pts, hullIdx, defects);

    while (true)
    {
        // For debug
        Mat3b dbg;
        cvtColor(img, dbg, COLOR_GRAY2BGR);

        vector<vector<Point>> tmp = {pts};
        drawContours(dbg, tmp, 0, Scalar(255, 127, 0));

        vector<int> defectsIdx;
        for (const Vec4i& v : defects)
        {
            float depth = float(v[3]) / 256.f;
            if (depth > 2) //  filter defects by depth
            {
                // Defect found
                defectsIdx.push_back(v[2]);

                int startidx = v[0]; Point ptStart(pts[startidx]);
                int endidx = v[1]; Point ptEnd(pts[endidx]);
                int faridx = v[2]; Point ptFar(pts[faridx]);

                line(dbg, ptStart, ptEnd, Scalar(255, 0, 0), 1);
                line(dbg, ptStart, ptFar, Scalar(0, 255, 0), 1);
                line(dbg, ptEnd, ptFar, Scalar(0, 0, 255), 1);
                circle(dbg, ptFar, 4, Scalar(127, 127, 255), 2);
            }
        }

        if (defectsIdx.size() < 2)
        {
            break;
        }

        // If I have more than two defects, remove the points between the two nearest defects
        pts = removeFromContour(pts, defectsIdx);
        convexHull(pts, hullIdx, false);
        convexityDefects(pts, hullIdx, defects);
    }


    // Draw result contour
    vector<vector<Point>> tmp = { pts };
    drawContours(out, tmp, 0, Scalar(0, 0, 255), 1);

    imshow("Result", out);
    waitKey();

    return 0;
}

---

UPDATE

Working on an approximated contour (e.g. using CHAIN_APPROX_SIMPLE in findContours) may be faster, but the length of contours must be computed using arcLength().

This is the snippet to replace in the swapping part of removeFromContour:

// Check if intervals are swapped
if (startIdx <= endIdx)
{
    //int len11 = endIdx - startIdx;
    vector<Point> inside(contour.begin() + startIdx, contour.begin() + endIdx);
    int len1 = (inside.empty()) ? 0 : arcLength(inside, false);

    //int len22 = contour.size() - endIdx + startIdx;
    vector<Point> outside1(contour.begin(), contour.begin() + startIdx);
    vector<Point> outside2(contour.begin() + endIdx, contour.end());
    int len2 = (outside1.empty() ? 0 : arcLength(outside1, false)) + (outside2.empty() ? 0 : arcLength(outside2, false));

    if (len2 < len1)
    {
        swap(startIdx, endIdx);
    }
}
else
{
    //int len1 = startIdx - endIdx;
    vector<Point> inside(contour.begin() + endIdx, contour.begin() + startIdx);
    int len1 = (inside.empty()) ? 0 : arcLength(inside, false);


    //int len2 = contour.size() - startIdx + endIdx;
    vector<Point> outside1(contour.begin(), contour.begin() + endIdx);
    vector<Point> outside2(contour.begin() + startIdx, contour.end());
    int len2 = (outside1.empty() ? 0 : arcLength(outside1, false)) + (outside2.empty() ? 0 : arcLength(outside2, false));

    if (len1 < len2)
    {
        swap(startIdx, endIdx);
    }
}

Answer 2

Here is a Python implementation following Miki's code.

import numpy as np
import cv2

def ed2(lhs, rhs):
    return(lhs[0] - rhs[0])*(lhs[0] - rhs[0]) + (lhs[1] - rhs[1])*(lhs[1] - rhs[1])


def remove_from_contour(contour, defectsIdx, tmp):
    minDist = sys.maxsize
    startIdx, endIdx = 0, 0

    for i in range(0,len(defectsIdx)):
        for j in range(i+1, len(defectsIdx)):
            dist = ed2(contour[defectsIdx[i]][0], contour[defectsIdx[j]][0])
            if minDist > dist:
                minDist = dist
                startIdx = defectsIdx[i]
                endIdx = defectsIdx[j]

    if startIdx <= endIdx:
        inside = contour[startIdx:endIdx]
        len1 = 0 if inside.size == 0 else cv2.arcLength(inside, False)
        outside1 = contour[0:startIdx]
        outside2 = contour[endIdx:len(contour)]
        len2 = (0 if outside1.size == 0 else cv2.arcLength(outside1, False)) + (0 if outside2.size == 0 else cv2.arcLength(outside2, False))
        if len2 < len1:
            startIdx,endIdx = endIdx,startIdx     
    else:
        inside = contour[endIdx:startIdx]
        len1 = 0 if inside.size == 0 else cv2.arcLength(inside, False)
        outside1 = contour[0:endIdx]
        outside2 = contour[startIdx:len(contour)]
        len2 = (0 if outside1.size == 0 else cv2.arcLength(outside1, False)) + (0 if outside2.size == 0 else cv2.arcLength(outside2, False))
        if len1 < len2:
            startIdx,endIdx = endIdx,startIdx

    if startIdx <= endIdx:
        out = np.concatenate((contour[0:startIdx], contour[endIdx:len(contour)]), axis=0)
    else:
        out = contour[endIdx:startIdx]
    return out


def remove_defects(mask, debug=False):
    tmp = mask.copy()
    mask = cv2.cvtColor(mask, cv2.COLOR_BGR2GRAY)

    # get contour
    contours, _ = cv2.findContours(
        mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
    assert len(contours) > 0, "No contours found"
    contour = sorted(contours, key=cv2.contourArea)[-1] #largest contour
    if debug:
        init = cv2.drawContours(tmp.copy(), [contour], 0, (255, 0, 255), 1, cv2.LINE_AA)
        figure, ax = plt.subplots(1)
        ax.imshow(init)
        ax.set_title("Initital Contour")

    hull = cv2.convexHull(contour, returnPoints=False)
    defects = cv2.convexityDefects(contour, hull)

    while True:
        defectsIdx = []
        
        for i in range(defects.shape[0]):
            s, e, f, d = defects[i, 0]
            start = tuple(contour[s][0])
            end = tuple(contour[e][0])
            far = tuple(contour[f][0])
            
            depth = d / 256
            if depth > 2:
                defectsIdx.append(f)

        if len(defectsIdx) < 2:
            break

        contour = remove_from_contour(contour, defectsIdx, tmp)
        hull = cv2.convexHull(contour, returnPoints=False)
        defects = cv2.convexityDefects(contour, hull)

    if debug:
      rslt = cv2.drawContours(tmp.copy(), [contour], 0, (0, 255, 255), 1)
      figure, ax = plt.subplots(1)
      ax.imshow(rslt)
      ax.set_title("Corrected Contour")

mask = cv2.imread("a.png")
remove_defects(mask, True)

Answer 3

I came up with the following approach for detecting the bounds of the rectangle/square. It works based on few assumptions: shape is rectangular or square, it is centered in the image, it is not tilted.

• divide the masked(filled) image in half along the x-axis so that you get two regions (a top half and a bottom half)
• take the projection of each region on to the x-axis
• take all the non-zero entries of these projections and take their medians. These medians give you the y bounds
• similarly, divide the image in half along y-axis, take the projections on to y-axis, then calculate the medians to get the x bounds
• use the bounds to crop the region

Median line and the projection for a top half of a sample image is shown below.

Resulting bounds and cropped regions for two samples:

The code is in Octave/Matlab, and I tested this on Octave (you need the image package to run this).

clear all
close all

im = double(imread('kTouF.png'));
[r, c] = size(im);
% top half
p = sum(im(1:int32(end/2), :), 1);
y1 = -median(p(find(p > 0))) + int32(r/2);
% bottom half
p = sum(im(int32(end/2):end, :), 1);
y2 = median(p(find(p > 0))) + int32(r/2);
% left half
p = sum(im(:, 1:int32(end/2)), 2);
x1 = -median(p(find(p > 0))) + int32(c/2);
% right half
p = sum(im(:, int32(end/2):end), 2);
x2 = median(p(find(p > 0))) + int32(c/2);

% crop the image using the bounds
rect = [x1 y1 x2-x1 y2-y1];
cr = imcrop(im, rect);
im2 = zeros(size(im));
im2(y1:y2, x1:x2) = cr;

figure,
axis equal
subplot(1, 2, 1)
imagesc(im)
hold on
plot([x1 x2 x2 x1 x1], [y1 y1 y2 y2 y1], 'g-')
hold off
subplot(1, 2, 2)
imagesc(im2)

Answer 4

As a starting point and assuming the defects are never too big relative to the object you are trying to recognize, you can try a simple erode+dilate strategy before using cv::matchShapes as shown below.

 int max = 40; // depending on expected object and defect size
 cv::Mat img = cv::imread("example.png");
 cv::Mat eroded, dilated;
 cv::Mat element = cv::getStructuringElement(cv::MORPH_ELLIPSE, cv::Size(max*2,max*2), cv::Point(max,max));
 cv::erode(img, eroded, element);
 cv::dilate(eroded, dilated, element);
 cv::imshow("original", img);
 cv::imshow("eroded", eroded);
 cv::imshow("dilated", dilated);