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.

Underclassman

Have you looked at convexityDefects? docs.opencv.org/2.4/modules/imgproc/doc/… – Conover

@Conover yes, the red dots in the convex-hull images are from thresholded convexityDefects' result. I just can't think of an algorithm on how to go on from there. – Underclassman

Got you, never used it but I just dropped that in there just in case :) – Conover

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);
    }
}

Singer

@Underclassman Probably with a smarter implementation of the above code, working on an approximated contour (similar to CHAIN_APPROX_SIMPLE), this may actually be really fast. Please post an answer if you find something that works for your requirements, it can be very helpful :D – Singer

currently the decision whether to swap is made by the index-distance within the contour? That's probably the reason, why for CV_CHAIN_APPROX_SIMPLE sometimes the wrong parts are cropped away (wrong direction)? Could arcLength be an appropriate heuristic instead? – Underclassman

@Underclassman check update. Now it works as before for me. To compute arc length you need the vector of points: you can probably arrange stuff in a smarter way :D – Singer

currently already fast enough without optimizations and works nearly perfectly for my data, thank you! – Underclassman

a few things to add: Before calling convexityDefects I have to check whether at least 4 points are present in the contour (I just break otherwise). And within removeFromContour it can happen, that startIdx == endIdx. I don't know WHY that happens, but I changed the the condition for updating minDist to not update for same Idx points and handle that case afterwards (had additionally to initialize startIdx and endIdx to test for that case). – Underclassman

@miki I am trying to translate this code in to python, I don't have much background in C++ specially vector data structure. Can you please explain what is happening in code snippet 1, in function removeFromContour lines starts with out.insert(). I tried to look up for the details of vector but couldn't find reasonable detail. It would be great if you can explain. Thanks. – Frechette

@MuazUsmani: were you succeeded in creating the python verison of it? – Graybeard

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)

Marcelmarcela

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)

Herwin

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);

Conover

the problem is that the object size can vary, so I can't fix max. Do you have any assumptions about how to choose max depending on some extractable contour's properties like boundind rectangle, contour area or similar? – Underclassman

Can you not use a percentage of the maximum dimension of the blob you are currently testing? Just a thought. – Conover

You could also try increasing amount of erosion/dilation until you find what you're looking for or nothing is left from erosion. – Conover

Unfortunately I have to do this with about 100 fps, maybe I can use subsampling though. I'll try to find a property to estimate the erosion size, thx. First results look ok, the "corners" are cropped etc but that won't be a real problem. However I guess/hope there could be a more elegant (and maybe faster) way, working with the contours and the assumption that the object is basically convex. – Underclassman

Unless the stuff you looking for is moving really fast you might get away with doing this every 10 frame or whatever works. Sub-sampling the images and the timeline may help meeting your real-time constraints. – Conover

On a (final) note: erosion/dilation must be faster than detecting contour, then finding convexity defects, then finding a way to deal with those. – Conover

I get those contours as input, To use erision/dilation I have to even draw the mask myself first. Input are just the contour points (which are hard to share on SO, so I posted the drawn images. Not sure about the speed of convexity defect computation, I'll compare in the end. – Underclassman