如何裁剪凸缺陷？

Question

我正在尝试从轮廓中检测和定位一些物体的位置。我得到的轮廓经常包含一些噪声(也许来自背景，我不知道)。这些物体应该看起来类似于矩形或正方形，如图所示:

使用形状匹配(cv::matchShapes)对包含这些物体的轮廓进行检测时，无论是否有噪声，都能得到非常好的结果，但是在存在噪声的情况下，我在精确定位方面遇到了问题。

噪声看起来像这样:

或 例如。

我的想法是找到凸性缺陷(即非凸出部分)并在它们变得太强时，通过某种方式将导致凹度的部分截掉。检测缺陷没问题，通常每个“不需要的结构”会得到两个缺陷，但我卡在了如何决定从轮廓的哪个位置移除点。

以下是一些轮廓、它们的掩模(这样你就可以轻松提取轮廓)以及包括阈值凸性缺陷的凸壳:

我是否可以顺着轮廓线走，局部地判断轮廓线是否“左转”（如果是逆时针方向），如果是，则删除轮廓点直到再次转弯为止？也许从凸性缺陷开始？

我正在寻找算法或代码，编程语言不应该很重要，算法才更加重要。

- Micka

---

你看过convexityDefects吗？http://docs.opencv.org/2.4/modules/imgproc/doc/structural_analysis_and_shape_descriptors.html#convexitydefects - zeFrenchy

@zeFrenchy 是的，凸包图像中的红点来自于阈值凸缺陷的结果。只是我想不出接下来该如何进行算法。 - Micka

知道了，虽然我从未使用过它，但我只是为了以防万一将其放在那里。 - zeFrenchy

Answer 1

这种方法只适用于点。您不需要为此创建掩码。

主要思想是：

1. 在轮廓上查找缺陷
2. 如果我找到至少两个缺陷，请找到最近的两个缺陷
3. 从轮廓中删除最接近的两个缺陷之间的点
4. 在新轮廓上重新开始从1

我得到了以下结果。如您所见，对于平滑的缺陷（例如第7张图像），它具有一些缺点，但对于清晰可见的缺陷效果很好。我不知道这是否会解决您的问题，但可以作为一个起点。实践中应该很快（您肯定可以优化下面的代码，特别是“removeFromContour”函数）。此外，此方法的唯一参数是凸性缺陷的数量，因此它适用于小型和大型缺陷斑块。

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

---

更新

使用近似轮廓（例如在findContours中使用CHAIN_APPROX_SIMPLE）可能更快，但必须使用arcLength()计算轮廓的长度。

这是在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);
    }
}

- Miki

---

@Micka 可能通过更智能的实现上述代码，针对近似轮廓（类似于CHAIN_APPROX_SIMPLE）进行处理，这可能会非常快速。如果您找到适合您要求的解决方案，请发布一个答案，这将非常有帮助:D - Miki

目前是否交换的决定是由轮廓内的索引距离做出的？这可能就是为什么对于CV_CHAIN_APPROX_SIMPLE有时会裁剪掉错误的部分（错误的方向）的原因。也许arcLength是一个合适的启发式方法？ - Micka

@Micka 检查更新。现在对我来说它像以前一样工作。要计算弧长，您需要点的向量：您可能可以以更聪明的方式安排这些内容:D - Miki

目前已经足够快而且在我的数据上几乎完美地工作，谢谢！ - Micka

需要补充的几件事情：在调用convexityDefects之前，我必须检查轮廓中是否至少存在4个点（否则就会中断）。在removeFromContour中，可能会发生startIdx == endIdx的情况。我不知道为什么会发生这种情况，但我更改了更新minDist条件的条件，以便不更新相同Idx点，并在处理该情况后进行处理（还必须将startIdx和endIdx初始化以测试该情况）。 - Micka

Answer 2

这里提供了一个Python实现，遵循Miki的代码。

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)

- Maxi

Answer 3

我想到了以下检测矩形/正方形边界的方法，它基于以下几个假设：形状为矩形或正方形，居中于图像中，不倾斜。

• 将掩膜（填充）图像沿着x轴分成两半，这样您就得到了两个区域（上半部分和下半部分）
• 对每个区域进行x轴投影
• 获取这些投影的所有非零条目，并取它们的中位数。这些中位数给出了y方向的范围
• 同样地，沿着y轴将图像分成两半，将投影投影到y轴上，然后计算中位数以获得x方向的范围
• 使用这些范围来裁剪该区域

下面显示了样本图像上半部分的中位线和投影。

下面是两个样本的结果边界和裁剪区域：

此代码是Octave/Matlab的，并且我在Octave上测试了它（您需要图像包才能运行此代码）。

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)

- dhanushka

Answer 4

作为起点并假设缺陷相对于您要识别的对象不太大，您可以在使用cv::matchShapes之前尝试简单的腐蚀+膨胀策略，如下所示。

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

- zeFrenchy

---

问题在于对象的大小可能会变化，因此我无法固定 max。您是否对如何根据可提取轮廓的属性（如边界矩形、轮廓面积或类似属性）选择 max 有任何假设？ - Micka

你能否不使用当前测试的 blob 的最大尺寸的百分比呢？这只是一个想法。 - zeFrenchy

除非你所寻找的东西移动得非常快，否则你可以每10帧或任何有效的时间间隔来执行此操作。对图像和时间轴进行子采样可能有助于满足实时约束。 - zeFrenchy

最后需要注意的是：腐蚀/膨胀必须比检测轮廓，然后找到凸缺陷并想办法处理它们更快。 - zeFrenchy

我将这些轮廓作为输入，要使用腐蚀/膨胀，我必须先自己绘制掩模。输入只是轮廓点（很难在SO上分享，所以我发布了绘制的图像）。不确定凸缺陷计算的速度，最后我会进行比较。 - Micka