从紧耦合线和噪声曲线中寻找直线

Question

我有这张树线作物的图像。我需要找到作物对齐的大致方向。我正在尝试获取图像的霍夫线，然后找到角度的分布模式。

我一直在关注this tutorial的裁剪线，但是在那一个中​​，裁剪线很稀疏。在这里，它们是密集的，经过灰度化、模糊化和使用精明的边缘检测，这就是我得到的

import cv2
import numpy as np
import matplotlib.pyplot as plt

img = cv2.imread('drive/MyDrive/tree/sample.jpg')
gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
gauss = cv2.GaussianBlur(gray, (3,3), 3)

plt.figure(figsize=(15,15))
plt.subplot(1,2,1)
plt.imshow(gauss)

gscale = cv2.Canny(gauss, 80, 140)
plt.subplot(1,2,2)
plt.imshow(gscale)
plt.show()

（左侧没有canny的模糊图像，留下一张经过canny预处理的图像）

之后，我按照教程，将预处理后的图像“骨架化”

size = np.size(gscale)

skel = np.zeros(gscale.shape, np.uint8)

ret, gscale = cv2.threshold(gscale, 128, 255,0)
element = cv2.getStructuringElement(cv2.MORPH_CROSS, (3,3))
done = False

while not done:
  eroded = cv2.erode(gscale, element)
  temp = cv2.dilate(eroded, element)
  temp = cv2.subtract(gscale, temp)
  skel = cv2.bitwise_or(skel, temp)
  gscale = eroded.copy()
  
  zeros = size - cv2.countNonZero(gscale)
  if zeros==size:
    done = True

给我

如您所见，仍然有一堆曲线。在上面使用HoughLines算法时，有11k条线散布在各处

lines = cv2.HoughLinesP(skel,1,np.pi/180,130)
a,b,c = lines.shape
for i in range(a):
    rho = lines[i][0][0]
    theta = lines[i][0][1]    
    a = np.cos(theta)
    b = np.sin(theta)
    x0 = a*rho
    y0 = b*rho
    x1 = int(x0 + 1000*(-b))
    y1 = int(y0 + 1000*(a))
    x2 = int(x0 - 1000*(-b))
    y2 = int(y0 - 1000*(a))
    cv2.line(img,(x1,y1),(x2,y2),(0,0,255),2, cv2.LINE_AA)#showing the results:

plt.figure(figsize=(15,15))
plt.subplot(121)#OpenCV reads images as BGR, this corrects so it is displayed as RGB
plt.plot()
plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB)) 
plt.title('Row Detection')
plt.xticks([])
plt.yticks([])
plt.subplot(122)
plt.plot()
plt.imshow(skel,cmap='gray')
plt.title('Skeletal Image')
plt.xticks([])
plt.yticks([])
plt.show()

我是 cv2 的新手，所以我不知道该怎么做。搜索并尝试了一堆东西，但没有一个有效。如何去除稍微大一点的点，去除弯曲的线条？

Answer 1

您可以使用 2D FFT 来找到作物对齐的大致方向（正如 mozway 在 cmets 中提出的那样）。这个想法是，当输入包含许多相同方向的线时，可以很容易地从出现在幅度谱中的中心光束中提取大致方向。您可以在this previous post 中找到有关其工作原理的更多信息。它直接处理输入图像，但最好应用 Gaussian + Canny 过滤器。

这是过滤后的灰度图像幅度谱中有趣的部分：

可以很容易地看到主光束。您可以通过以增加的角度迭代多条线来提取其角度，并将每条线上的幅度值求和，如下图所示：

这里是每条线相对于线的角度（弧度）绘制的幅度和：

基于此，您只需要找到使计算和最大化的角度。

这是生成的代码：

def computeAngle(arr):
    # Naive inefficient algorithm
    n, m = arr.shape
    yCenter, xCenter = (n-1, m//2-1)
    lineLen = m//2-2
    sMax = 0.0
    bestAngle = np.nan
    for angle in np.arange(0, math.pi, math.pi/300):
        i = np.arange(lineLen)
        y, x = (np.sin(angle) * i + 0.5).astype(np.int_), (np.cos(angle) * i + 0.5).astype(np.int_)
        s = np.sum(arr[yCenter-y, xCenter+x])
        if s > sMax:
            bestAngle = angle
            sMax = s
    return bestAngle

# Load the image in gray
img = cv2.imread('lines.jpg')
gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)

# Apply some filters
gauss = cv2.GaussianBlur(gray, (3,3), 3)
gscale = cv2.Canny(gauss, 80, 140)

# Compute the 2D FFT of real values
freqs = np.fft.rfft2(gscale)

# Shift the frequencies (centering) and select the low frequencies
upperPart = freqs[:freqs.shape[0]//4,:freqs.shape[1]//2]
lowerPart = freqs[-freqs.shape[0]//4:,:freqs.shape[1]//2]
filteredFreqs = np.vstack((lowerPart, upperPart))

# Compute the magnitude spectrum
magnitude = np.log(np.abs(filteredFreqs))

# Correct the angle
magnitude = np.rot90(magnitude).copy()

# Find the major angle
bestAngle = computeAngle(magnitude)

Answer 2

为了完整起见，我也想发布 Sobel Gradient Angle 方式。

总体思路：

1. 为每个像素计算 X 和 Y 梯度（例如使用 Sobel）
2. 根据 X 和 Y 梯度信息及其分布计算角度
3. 找到模式，例如带有直方图并选择最高的一个

用 C++ 编写，但可能很容易转换为 python：

int main()
{
    try
    {
        cv::Mat img = cv::imread("C:/data/StackOverflow/cropLines/lines.jpg", cv::IMREAD_GRAYSCALE);

        // tests with artificial lines:
        //img = cv::Mat::zeros(img.size(), CV_8UC1);
        //img = cv::Mat(img.size(), CV_8UC1, cv::Scalar::all(255));
        //cv::line(img, cv::Point(0, img.rows), cv::Point(img.cols, 0), cv::Scalar::all(255), 10); // sample to check angles
        //cv::line(img, cv::Point(img.cols, img.rows), cv::Point(0, 0), cv::Scalar::all(255), 10); // sample to check angles
        //cv::line(img, cv::Point(img.cols, img.rows/2), cv::Point(0, img.rows/2), cv::Scalar::all(255), 10); // sample to check angles
        //cv::line(img, cv::Point(img.cols/2, img.rows), cv::Point(img.cols/2, 0), cv::Scalar::all(255), 10); // sample to check angles
        //cv::line(img, cv::Point(img.cols / 2, img.rows), cv::Point(img.cols / 2, 0), cv::Scalar::all(255), 10); // sample to check angles
        //cv::line(img, cv::Point(img.cols / 2, img.rows), cv::Point(img.cols / 2, 0), cv::Scalar::all(0), 10); // sample to check angles
        cv::imshow("img", img);

        cv::Mat gradX, gradY, mag, angle;
        cv::Sobel(img, gradX, CV_32F, 1, 0, 3);
        cv::Sobel(img, gradY, CV_32F, 0, 1, 3);

        cv::cartToPolar(gradX, gradY, mag, angle, true);

        cv::Mat magMask = mag > 0; // dont use pixels where angle is 0 just because there is no gradient.

        float scaleX = 3;
        float scaleY = 0.15;
        float maxValueY = 3000;
        cv::Mat chart = cv::Mat(maxValueY * scaleY, 360 * scaleX, CV_8UC3);

        cv::Point prev(-1, -1);
        double window = 5.0; // window size 1 is much more noisy but still works
        // this loop can probably be optimized with an optimized histogram compuation from any library
        for (int i = 0; i <= 360; i = i + 1)
        {
            double amount = cv::countNonZero((angle >= i-window/2 & angle < i + window/2) & (magMask));
            std::cout << i << "°: " << amount << std::endl;

            cv::Point current(i*scaleX, chart.rows - amount*scaleY/window);
            if (prev.x >= 0) cv::line(chart, prev, current, cv::Scalar(0, 0, 255), 1);
            prev = current;
        }

        cv::line(chart, cv::Point(45 * scaleX, 0), cv::Point(45 * scaleX, 100 * scaleY), cv::Scalar(255, 0, 0), 1);
        cv::line(chart, cv::Point(90 * scaleX, 0), cv::Point(90 * scaleX, 100 * scaleY), cv::Scalar(255, 0, 0), 1);
        cv::line(chart, cv::Point(135 * scaleX, 0), cv::Point(135 * scaleX, 100 * scaleY), cv::Scalar(255, 0, 0), 1);
        cv::line(chart, cv::Point(180 * scaleX, 0), cv::Point(180 * scaleX, 100 * scaleY), cv::Scalar(255, 0, 0), 1);
        cv::line(chart, cv::Point(225 * scaleX, 0), cv::Point(225 * scaleX, 100 * scaleY), cv::Scalar(255, 0, 0), 1);
        cv::line(chart, cv::Point(270 * scaleX, 0), cv::Point(270 * scaleX, 100 * scaleY), cv::Scalar(255, 0, 0), 1);
        cv::line(chart, cv::Point(315 * scaleX, 0), cv::Point(315 * scaleX, 100 * scaleY), cv::Scalar(255, 0, 0), 1);
        cv::line(chart, cv::Point(360 * scaleX, 0), cv::Point(360 * scaleX, 100 * scaleY), cv::Scalar(255, 0, 0), 1);

        cv::imshow("chart", chart);
        cv::imwrite("C:/data/StackOverflow/cropLines/chart.png", chart);

        cv::imwrite("C:/data/StackOverflow/cropLines/input.png", img);

        cv::waitKey(0);
    }
    catch (std::exception& e)
    {
        std::cout << e.what() << std::endl;
    }
}

给出这个结果，图像顶部的每条蓝线都是 45°。最大值为 52°（及其 180° 的倍数），其中梯度相对于线方向旋转了 90°。因此，该图像中的线方向是从 x 轴顺时针方向 142° 或逆时针方向 38°。