具有正交特征的任意图像旋转角的自动检测

#1 楼

如果我正确理解方法1，并且使用圆形对称区域并围绕该区域的中心进行旋转，则可以消除该区域对旋转角度的依赖性，并通过优点函数获得更公平的比较不同的旋转角度。我将建议一种基本上与之等效的方法，但它使用完整图像并且不需要重复图像旋转，并且将包括低通滤波以消除像素网格各向异性并进行去噪。

渐变各向同性低通滤波图像的图像

首先，让我们计算全尺寸样本图像中绿色通道的每个像素的局部梯度矢量。

我推导了水平和通过将理想的低通滤波器的连续空间脉冲响应与平坦的圆形频率响应区分开来，从而确保垂直微分内核，从而通过确保与水平或垂直对角线不存在不同的细节水平，从而消除了图像轴选择的影响，通过采样结果函数并应用旋转的余弦窗口：

$$ \ begin {gather} h_x [x，y] = \ begin {cases} 0＆\ text {if} x = y = 0，\\-\ displaystyle \ frac {\ omega_c ^ 2 \，x \，J_2 \ left（\ omega_c \ sqrt {x ^ 2 + y ^ 2} \ right）} {2 \ p i \，（x ^ 2 + y ^ 2）}＆\ text {否则，} \ end {cases} \\
h_y [x，y] = \ begin {cases0＆0＆text {if} x = y = 0，\\-\ displaystyle \ frac {\ omega_c ^ 2 \，y \，J_2 \ left（\ omega_c \ sqrt {x ^ 2 + y ^ 2} \ right）} {2 \ pi \，（ x ^ 2 + y ^ 2）}＆\ text {否则，} \ end {cases} \ end {gather} \ tag {1} $$

其中$ J_2 $是二阶Bessel第一种函数，$ \ omega_c $是弧度的截止频率。 Python源（没有等式1的减号）：

import matplotlib.pyplot as plt
import scipy
import scipy.special
import numpy as np

def rotatedCosineWindow(N):  # N = horizontal size of the targeted kernel, also its vertical size, must be odd.
  return np.fromfunction(lambda y, x: np.maximum(np.cos(np.pi/2*np.sqrt(((x - (N - 1)/2)/((N - 1)/2 + 1))**2 + ((y - (N - 1)/2)/((N - 1)/2 + 1))**2)), 0), [N, N])

def circularLowpassKernelX(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size, must be odd.
  kernel = np.fromfunction(lambda y, x: omega_c**2*(x - (N - 1)/2)*scipy.special.jv(2, omega_c*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2))/(2*np.pi*((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2)), [N, N])
  kernel[(N - 1)//2, (N - 1)//2] = 0
  return kernel

def circularLowpassKernelY(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size, must be odd.
  kernel = np.fromfunction(lambda y, x: omega_c**2*(y - (N - 1)/2)*scipy.special.jv(2, omega_c*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2))/(2*np.pi*((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2)), [N, N])
  kernel[(N - 1)//2, (N - 1)//2] = 0
  return kernel

N = 41  # Horizontal size of the kernel, also its vertical size. Must be odd.
window = rotatedCosineWindow(N)

# Optional window function plot
#plt.imshow(window, vmin=-np.max(window), vmax=np.max(window), cmap='bwr')
#plt.colorbar()
#plt.show()

omega_c = np.pi/4  # Cutoff frequency in radians <= pi
kernelX = circularLowpassKernelX(omega_c, N)*window
kernelY = circularLowpassKernelY(omega_c, N)*window

# Optional kernel plot
#plt.imshow(kernelX, vmin=-np.max(kernelX), vmax=np.max(kernelX), cmap='bwr')
#plt.colorbar()
#plt.show()

图1.二维旋转余弦窗口。

图2.窗口水平各向同性低通微分内核，用于不同的截止频率$ \ omega_c $设置。顶部：omega_c = np.pi，中间：omega_c = np.pi/4，底部：omega_c = np.pi/16。等式的负号。 1个被遗漏了。垂直内核看起来相同，但已旋转90度。水平和垂直核的权重和分别为权重$ \ cos（\ phi）$和$ \ sin（\ phi）$，得出了梯度角为$ \ phi $的相同类型的分析核。 br />
脉冲响应的微分不会影响带宽，正如在Python中通过其二维快速傅里叶变换（FFT）可以看出的那样：


图3. $ h_x $的2-d FFT的幅度。在频域中，微分表现为平坦的圆形通带乘以$ \ omega_x $以及90度相移，该幅度在幅度上不可见。

对绿色进行卷积通道并收集二维梯度矢量直方图，以便在Python中进行视觉检查：

# Optional FFT plot
absF = np.abs(np.fft.fftshift(np.fft.fft2(circularLowpassKernelX(np.pi, N)*window)))
plt.imshow(absF, vmin=0, vmax=np.max(absF), cmap='Greys', extent=[-np.pi, np.pi, -np.pi, np.pi])
plt.colorbar()
plt.show()

这也会裁剪数据，从每个边沿丢弃在直方图分析之前被矩形图像边界污染。

$ \ pi $
$ \ frac {\ pi} {2} $
$ \ frac { \ pi} {4} $ $ \ frac {\ pi} {8} $
$ \ frac {\ pi} {16} $
$ \ frac {\ pi} {32} $
$ \ frac {\ pi} {64} $
— $ 0 $图4.梯度矢量的二维直方图，用于不同的低通滤波器截止频率$ \ omega_c $设置。顺序如下：首先使用(N - 1)//2：N=41，omega_c = np.pi，omega_c = np.pi/2（与Python清单相同），omega_c = np.pi/4，omega_c = np.pi/8，然后：omega_c = np.pi/16：N=81，omega_c = np.pi/32：N=161通过低通滤波进行降噪可增强直方图中电路迹线边缘梯度的方向。

矢量长度加权圆均值方向

有Yamartino方法，可以在一遍通过多个风矢量样本中找到“平均”风向。它基于循环量的平均值，该平均值被计算为余弦的位移，该余弦是每个余弦的总和，每个余弦的位移量为周期$ 2 \ pi $。我们可以使用相同方法的向量长度加权版本，但是首先我们需要将所有等于$ \ pi / 2 $取模的方向捆绑在一起。我们可以通过将每个梯度矢量$ [X_k，Y_k] $的角度乘以4并使用复数表示来实现：

$$ Z_k = \ frac {（X_k + Y_k i）^ 4} {\ sqrt {X_k ^ 2 + Y_k ^ 2} ^ 3} = \ frac {X_k ^ 4-6X_k ^ 2Y_k ^ 2 + Y_k ^ 4 +（4X_k ^ 3Y_k-4X_kY_k ^ 3）i} {\ sqrt { X_k ^ 2 + Y_k ^ 2} ^ 3}，\ tag {2} $$

满足$ | Z_k | = \ sqrt {X_k ^ 2 + Y_k ^ 2} $，随后解释为$ Z_k $从$-\ pi $到$ \ pi $的相位表示从$-\ pi / 4 $到$ \ pi /的角度4 $，将计算出的循环平均相位除以4：

$$ \ phi = \ frac {1} {4} \ operatorname {atan2} \ left（\ sum_k \ operatorname {Im}（ Z_k），\ sum_k \ operatorname {Re}（Z_k）\ right）\ tag {3} $$

其中$ \ phi $是估计的图像方向。

可以通过再次遍历数据并计算复数$ Z_k $的相位与估计的圆形平均相位$ 4 \之间的平均加权平方圆距离$ \ text {MSCD} $来评估估计的质量。 phi $，以$ | Z_k | $作为权重：

$$ \ begin {gather} \ text {MSCD} = \ frac {\ sum_k | Z_k | \ bigg（1-\ cos \大（4 \ phi-\ operatorname {atan2} \ big（\ operatorname {Im}（Z_k），\ operatorname {Re}（Z_k）\ big）\ Big）\ bigg）} {\ sum_k | Z_k |} \\\
= \ frac {\ sum_k \ frac {| Z_k |} {2} \ left（\ left（\ cos（4 \ phi）-\ frac {\ operatorname {Re}（Z_k）} {{Z_k | } \ right）^ 2 + \ left（\ sin（4 \ phi）-\ frac {\ o peratorname {Im}（Z_k）} {| Z_k |} \ right）^ 2 \ right）} {\ sum_k | Z_k |} \\
= \ frac {\ sum_k \ big（| Z_k | -\ operatorname {Re}（Z_k）\ cos（4 \ phi）-\ operatorname {Im}（Z_k）\ sin（4 \ phi）\ big）} {\ sum_k | Z_k |}，\ end {gather} \标签{4} $$

根据方程式计算的$ \ phi $将其最小化。 3.在Python中：基于我的omega_c = np.pi/64实验（未显示），我认为即使对于非常大的图像，我们也不会耗尽数值精度。对于不同的滤波器设置（带注释），输出在-45至45度之间报告：

import scipy.ndimage

img = plt.imread('sample.tif').astype(float)
X = scipy.ndimage.convolve(img[:,:,1], kernelX)[(N - 1)//2:-(N - 1)//2, (N - 1)//2:-(N - 1)//2]  # Green channel only
Y = scipy.ndimage.convolve(img[:,:,1], kernelY)[(N - 1)//2:-(N - 1)//2, (N - 1)//2:-(N - 1)//2]  # ...

# Optional 2-d histogram
#hist2d, xEdges, yEdges = np.histogram2d(X.flatten(), Y.flatten(), bins=199)
#plt.imshow(hist2d**(1/2.2), vmin=0, cmap='Greys')
#plt.show()
#plt.imsave('hist2d.png', plt.cm.Greys(plt.Normalize(vmin=0, vmax=hist2d.max()**(1/2.2))(hist2d**(1/2.2))))  # To save the histogram image
#plt.imsave('histkey.png', plt.cm.Greys(np.repeat([(np.arange(200)/199)**(1/2.2)], 16, 0)))

似乎很有用，可以减少以$ \ operatorname {acos}（1-\ text {MSCD}）$计算的均方根圆角距离（RMSCD）等效角。如果没有二维旋转的余弦窗口，则某些结果将偏离大约一个度（未显示），这意味着对分析过滤器进行适当的窗口化非常重要。 RMSCD等效角度不是角度估计中的误差的直接估计值，而误差应该小得多。

另一种平方长权重函数

让我们尝试向量长度作为替代权重函数，方法如下：

$$ Z_k = \ frac {（X_k + Y_k i）^ 4} {\ sqrt {X_k ^ 2 + Y_k ^ 2} ^ 2} = \ frac {X_k ^ 4-6X_k ^ 2Y_k ^ 2 + Y_k ^ 4 +（4X_k ^ 3Y_k-4X_kY_k ^ 3）i} {X_k ^ 2 + Y_k ^ 2}，\ tag {5} $$

在Python中：

absZ = np.sqrt(X**2 + Y**2)
reZ = (X**4 - 6*X**2*Y**2 + Y**4)/absZ**3
imZ = (4*X**3*Y - 4*X*Y**3)/absZ**3
phi = np.arctan2(np.sum(imZ), np.sum(reZ))/4

sumWeighted = np.sum(absZ - reZ*np.cos(4*phi) - imZ*np.sin(4*phi))
sumAbsZ = np.sum(absZ)
mscd = sumWeighted/sumAbsZ

print("rotate", -phi*180/np.pi, "deg, RMSCD =", np.arccos(1 - mscd)/4*180/np.pi, "deg equivalent (weight = length)")

平方长度的重量将RMSCD等效角减小了大约一个角度：

rotate 32.29809399495655 deg, RMSCD = 17.057059965741338 deg equivalent (omega_c = np.pi)
rotate 32.07672617150525 deg, RMSCD = 16.699056648843566 deg equivalent (omega_c = np.pi/2)
rotate 32.13115293914797 deg, RMSCD = 15.217534399922902 deg equivalent (omega_c = np.pi/4, same as in the Python listing)
rotate 32.18444156018288 deg, RMSCD = 14.239347706786056 deg equivalent (omega_c = np.pi/8)
rotate 32.23705383489169 deg, RMSCD = 13.63694582160468 deg equivalent (omega_c = np.pi/16)

这似乎是一个更好的权重功能。我还添加了临界值$ \ omega_c = \ pi / 32 $和$ \ omega_c = \ pi / 64 $。他们使用更大的mpmath导致图像的裁剪不同，而没有严格可比的MSCD值。

1-d直方图

平方长权重函数的好处是更多通过$ Z_k $个相位的一维加权直方图可以明显看出。 Python脚本：

absZ_alt = X**2 + Y**2
reZ_alt = (X**4 - 6*X**2*Y**2 + Y**4)/absZ_alt
imZ_alt = (4*X**3*Y - 4*X*Y**3)/absZ_alt
phi_alt = np.arctan2(np.sum(imZ_alt), np.sum(reZ_alt))/4

sumWeighted_alt = np.sum(absZ_alt - reZ_alt*np.cos(4*phi_alt) - imZ_alt*np.sin(4*phi_alt))
sumAbsZ_alt = np.sum(absZ_alt)
mscd_alt = sumWeighted_alt/sumAbsZ_alt

print("rotate", -phi_alt*180/np.pi, "deg, RMSCD =", np.arccos(1 - mscd_alt)/4*180/np.pi, "deg equivalent (weight = length^2)")

图5.梯度矢量角度的线性插值加权直方图，包裹为$-\ pi / 4 \ ldots \ pi / 4 $并按（从峰顶从下到上的顺序）加权：无权重（黑色），梯度向量长度（红色），梯度向量长度的平方（蓝色）。箱宽度为0.1度。过滤器截止值为N，与Python清单中的相同。底部的图是在峰值处放大的。由等式给出的$ x $和$ y $微分滤波器脉冲响应。可以将图1理解为用于形成可转向微分滤波器的脉冲响应的基本函数，该可微分滤波器从等式的右侧旋转了$ h_x [x，y] $（等式1）。通过转换等式，更容易看出这一点。 1到极坐标：

$$ \ begin {align} h_x（r，\ theta）= h_x [r \ cos（\ theta），r \ sin（\ theta）]＆= \ begin {cases} 0＆\ text {if} r = 0，\\-\ displaystyle \ frac {\ omega_c ^ 2 \，r \ cos（\ theta）\，J_2 \ left（\ omega_c r \ right）} {2 \ pi \，r ^ 2}＆\ text {otherwise} \ end {cases} \\
＆= \ cos（\ theta）f（r），\\
h_y（r，\ theta） = h_y [r \ cos（\ theta），r \ sin（\ theta）]＆= \ begin {cases} 0＆\ text {if} r = 0，\\-\ displaystyle \ frac {\ omega_c ^ 2 \， r \ sin（\ theta）\，J_2 \ left（\ omega_c r \ right）} {2 \ pi \，r ^ 2}＆\ text {other}} end {cases} \\
＆= \ sin（\ theta）f（r），\\
f（r）＆= \ begin {cases} 0＆\ text {if} r = 0，\\-\ displaystyle \ frac {\ omega_c ^ 2 \ ，r \，J_2 \ left（\ omega_c r \ right）} {2 \ pi \，r ^ 2}＆\ text {否则，} \ end {cases} \ end {align} \ tag {6} $$

其中水平和垂直微分滤波器的脉冲响应都具有相同的径向因子函数$ f（r）$。通过转向角$ \ phi $获得的$ h_x（r，\ theta）$的任何旋转版本$ h（r，\ theta，\ phi）$可通过以下方式获得： \ theta，\ phi）= h_x（r，\ theta-\ phi）= \ cos（\ theta-\ phi）f（r）\ tag {7} $$

这个想法是，转向核$ h（r，\ theta，\ phi）$可以构造为$ h_x（r，\ theta）$和$ h_x（r，\ theta）$的加权总和，其中$ \ cos（\ phi）$和$ \ sin（\ phi）$作为权重，的确如此：

$$ \ cos（\ phi）h_x（r，\ theta）+ \ sin（\ phi）h_y（r，\ theta）= \ cos（\ phi）\ cos（\ theta）f（r）+ \ sin（\ phi）\ sin（\ theta）f（r）= \ cos （\ theta-\ phi）f（r）= h（r，\ theta，\ phi）。\ tag {8} $$

如果我们想到各向同性低通滤波后的信号作为输入信号，并针对从坐标$ x $，$ y旋转角度$ \ phi $的第一个旋转坐标$ x_ \ phi $，$ y_ \ phi $构造偏导数$。 （可以将导数视为线性时不变系统。）我们具有：

$$ \ begin {gather} x = \ cos（\ phi）x_ \ phi-\ sin（\ phi） y_ \ phi，\\
y = \ sin（\ phi）x_ \ phi + \ cos（\ phi）y_ \ phi \ end {gather} \ tag {9} $$

使用偏导数的链式规则，相对于$ x_ \ phi $的偏导数算子可以表示为相对于$ x $和$ y $的偏导数的余弦和正弦加权和：

$$ \ begin {gather} \ frac {\ partial} {\ partial x_ \ phi} = \ frac {\ partial x} {\ partial x_ \ phi} \ frac {\ partial} {\ partial x} + \ frac {\ partial y} {\ partial x_ \ phi} \ frac {\ partial} {\ partial y} = \ frac {\ partial \ big（\ cos（\ phi）x_ \ phi-\ sin（\ phi） y_ \ phi \ big）} {\ partial x_ \ phi} \ frac {\ partial} {\ partial x} + \ frac {\ partial \ big（\ sin（\ phi）x_ \ phi + \ cos（\ phi） y_ \ phi \ big）} {\ partial x_ \ phi} \ frac {\ partial} {\ partial y} = \ cos（\ phi）\ frac {\ partial} {\ partial x} + \ sin（\ phi） \ frac {\ partial} {\ partial y} \ end {gather} \ tag {10} $$

还有待探讨的问题是梯度矢量角度的八倍圆均值在某种程度上与“最激活”转向微分滤波器的角度$ \ phi $相关。

可能的改进

为了进一步改善结果，还可以为红色和蓝色通道计算梯度，以将其作为附加数据包括在“平均值”计算中。

我想到了此方法的可能扩展：

1）使用更大的分析过滤器内核集，并检测边缘而不是检测梯度。这需要仔细地设计，以便在所有方向上的边缘都得到同等对待，也就是说，应该通过正交核的加权和获得任何角度的边缘检测器。 （我认为）可以通过应用方程式的微分运算符来获得一组合适的内核。图11，图6（另请参见我的数学堆栈交换文章）介绍了圆形对称低通滤波器的连续空间脉冲响应。
<$ /> $$ \ begin {gather} \ lim_ {h \到0} \ frac {\ sum_ {N = 0} ^ {4N + 1}（-1）^ nf \ bigg（x + h \ cos \ left（\ frac {2 \ pi n} {4N + 2} \\右），y + h \ sin \ left（\ frac {2 \ pi n} {4N + 2} \ right）\ bigg）} {h ^ {2N + 1}}，\\
\ lim_ { h \ to 0} \ frac {\ sum_ {N = 0} ^ {4N + 1}（-1）^ nf \ bigg（x + h \ sin \ left（\ frac {2 \ pi n} {4N + 2 } \ right），y + h \ cos \ left（\ frac {2 \ pi n} {4N + 2} \ right）\ bigg）} {h ^ {2N + 1}} \ end {gather} \ tag { 11} $$

图6.用于构造高阶边缘检测器的微分算子中的狄拉克增量相对位置。

2）计算（加权）平均值循环量的“余量”可以理解为将相同频率的余弦值相乘，并乘以该量的样本（并按权重进行标度），然后求出结果函数的峰值。如果将经过精心选择的相对幅度的，类似的，已移位余弦的已移位和按比例缩放的谐波添加到混合中，形成更清晰的平滑核，则总和中可能会出现多个峰，并且可以报告最大的峰。结合适当的谐波，可以得到一种局部平均值，在很大程度上忽略了分布主峰以外的离群值。

替代方法

也可以通过角度$ \ phi $和角度$ \ phi + \ pi / 2 $旋转的“长边”内核对图像进行卷积，并计算两个卷积像素的均方根图片。将报告最大化均方的角度$ \ phi $。这种方法可能会很好地完善图像方向，因为要大步搜索整个角度$ \ phi $空间是有风险的。

另一种方法是非局部方法，例如cross -将遥远的相似区域相关联，如果您知道水平或垂直走线较长，或者要素在水平或垂直方向上重复多次，则适用。

#2 楼

这里有一个类似的DSP技巧，但我不记得确切的细节。无论方向如何，都必须弄清楚织物图案的匹配情况。因此，您可能需要对此进行研究。

获取一个圆形样本。沿圆的辐条求和以获得圆周轮廓。然后他们对此进行了DFT（毕竟本质上是循环的）。扔掉相位信息（使其独立于方向）并进行比较。

然后他们可以分辨出两种织物的图案是否相同。

您的问题很相似。 />
在我看来，无需先尝试，预DFT轮廓的特征应显示方向。沿辐条而不是总和进行标准偏差处理可能会更好，也许两者都可以。

现在，如果您具有定向的参考图像，则可以使用其技术。

Ced


您对精度的要求非常严格。

我给了这个建议。取每种颜色沿辐条的两个后续点之间的差的绝对值的总和。

此处为圆周图。您的值是用白色标记绘制的。抱歉。


进展报告：一些

我决定了一个三步过程。

1）找到评估点。

2）粗略测量

3）精细测量

当前，第一步是用户投入。它应该是可自动执行的，但是我并不在意。我有第二步的草稿。我想尝试一些调整。最后，对于第三步，我有一些候选人将进行测试，以查看哪种方法最有效。

好消息是照明速度很快。如果您的唯一目的是使图像在网页上看起来水平，则您的公差太严格了，粗略测量应该足够准确。

这是粗略测量。每个像素约为0.6度。 （实际编辑为0.3） />大多数都不是那么好，但是它们很便宜（而且相当本地化），而且对于人类来说，找到易于阅读的地方很容易。蛮力对于程序来说应该可以正常工作。

结果可以大大改善，这是一个简单的基准测试。我还没有准备好做任何解释，也没有发布代码，但是此屏幕快照尚未进行照相处理。


进度报告：代码已发布，我已经完成了

此屏幕快照是处理Marcus 45度镜头的程序。



颜色通道是独立处理的。

选择一个点作为扫描中心。

直径以不连续的角度扫过180度

在每个角度，“波动性”在整个直径上测量。对收集样本的每个通道进行跟踪。样本值是样本点所位于的任何网格平方的四个角值的线性插值。

对于每个通道迹线

样本乘以VonHann窗口函数

对样本进行平滑/差分处​​理

差分的均方根值用作波动率度量

下排图是：

第一个是0到180度的范围，每个像素为0.5度。
第二个是围绕选定角度的范围，每个像素为0.1度。
第三个是围绕选定角度的扫描，每个像素为0.01度。
第四条曲线为Differ曲线

最初的选择是三个通道的最小平均波动率。这将是最佳角度，但通常不是最佳角度。最低处的对称性比最低处更好。在该邻域中最合适的抛物线应该会产生很好的答案。 ：//forum.gambas.one/viewtopic.php？f = 4＆t = 707

这是一个普通的zip文件，因此您无需安装Gambas即可查看源代码。这些文件位于“ .src”子目录中。

删除VonHann窗口可以提高准确性，因为它可以有效地延长跟踪时间，但会增加摆动。如果中心不重要，可能会发现双倍的VonHann会更好，并且会更快地检测到“跷跷板击中地面时”。只要增加图像的走线长度，就可以轻松地提高准确性（是的，这是自动的）。更好的窗口功能sinc？

我在当前设置下所采取的措施证实了3.19值+/-。03 ish。我可以考虑几种策略将其应用于图像。正如他们所说，这是读者的一种练习。或者在这种情况下，就是OP。我稍后再尝试。

算法和程序都有改进的余地，但它们确实非常有用。

线性插值的工作原理如下

'---- Whole Number Portion

        x = Floor(rx)
        y = Floor(ry)

'---- Fractional Portions

        fx = rx - x
        fy = ry - y

        gx = 1.0 - fx
        gy = 1.0 - fy

'---- Weighted Average

        vtl = ArgValues[x, y] * gx * gy         ' Top Left
        vtr = ArgValues[x + 1, y] * fx * gy     ' Top Right
        vbl = ArgValues[x, y + 1] * gx * fy     ' Bottom Left
        vbr = ArgValues[x + 1, y + 1] * fx * fy ' Bottom Rigth

        v = vtl + vtr + vbl + vbr

有人知道它的常规名称吗？

#3 楼

这是我上一个答案的第一个建议扩展。

理想的圆形对称带限滤波器

我们构造了四个滤波器的正交存储区，它们限制在一个圆内频率平面上的半径$ \ omega_c $的半径。这些滤波器的脉冲响应可以线性组合以形成定向边缘检测内核。通过将前两对“沙滩球状”微分算子应用于圆对称理想带限滤波器冲激响应$ h（x，y）的连续空间冲激响应，可以得到任意归一化的正交滤波器冲激响应集。 ）$：

$$ h（x，y）= \ frac {\ omega_c} {2 \ pi \ sqrt {x ^ 2 + y ^ 2}} J_1 \ big（\ omega_c \ sqrt {x ^ 2 + y ^ 2} \ big）\ tag {1} $$

$$ \ begin {align} h_ {0x}（x，y）＆\ propto \ frac {d } {dx} h（x，y），\\
h_ {0y}（x，y）＆\ propto \ frac {d} {dy} h（x，y），\\
h_ {1x}（x，y）＆\ propto \ left（\ left（\ frac {d} {dx} \ right）^ 3-3 \ frac {d} {dx} \ left（\ frac {d} { dy} \ right）^ 2 \ right）h（x，y），\\
h_ {1y}（x，y）＆\ propto \ left（\ left（\ frac {d} {dy} \ right）^ 3-3 \ frac {d} {dy} \ left（\ frac {d} {dx} \ right）^ 2 \ right）h（x，y）\ end {align} \ tag {2} $ $

$$ \ begin {align} h_ {0x}（x，y）＆= \ begin {cases} 0＆\ text {if} x = y = 0，\\-\ displaystyle \ frac {\ omega_c ^ 2 \，x \，J_2 \ left（\ omega_c \ sqrt {x ^ 2 + y ^ 2} \ right）} {2 \ pi \，（x ^ 2 + y ^ 2）}＆\文字{否则，} \ end {cases} \\
h_ {0y}（x，y）＆= h_ {0x} [y，x]，\\
h_ {1x}（x，y）＆= \ begin {cases} 0＆\ text {if} x = y = 0，\\\ frac {\ begin {array} {l} \ Big（ω_cx（3y ^ 2-x ^ 2）\ big（J_0 \ left（ω_c\ sqrt {x ^ 2 + y ^ 2} \ right）ω_c\ sqrt {x ^ 2 + y ^ 2}（ω_c^ 2x ^ 2 +ω_c^ 2y ^ 2-24）\\-8J_1 \ left（ω_c\ sqrt {x ^ 2 + y ^ 2} \ right）（ω_c^ 2x ^ 2 +ω_c^ 2y ^ 2-6）\ big）\ Big）\ end {array}} {2π（x ^ 2 + y ^ 2）^ {7/2}}＆\ text {否则，} \ end {cases} \\
h_ {1y}（x，y）＆= h_ {1x} [y，x]，\ end {align} \ tag {3} $$

其中$ J_ \ alpha $是第一类命令$ \ alpha $的贝塞尔函数，而$ \ propto $表示“与”成比例。我使用Wolfram Alpha查询（（ᵈ/ dx）³;ᵈ/ dx;ᵈ/ dx（ᵈ/ dy）²）进行区分，并简化了结果。

Python中的截断内核：

import matplotlib.pyplot as plt
import scipy
import scipy.special
import numpy as np

def h0x(x, y, omega_c):
  if x == 0 and y == 0:
    return 0
  return -omega_c**2*x*scipy.special.jv(2, omega_c*np.sqrt(x**2 + y**2))/(2*np.pi*(x**2 + y**2))

def h1x(x, y, omega_c):
  if x == 0 and y == 0:
    return 0
  return omega_c*x*(3*y**2 - x**2)*(scipy.special.j0(omega_c*np.sqrt(x**2 + y**2))*omega_c*np.sqrt(x**2 + y**2)*(omega_c**2*x**2 + omega_c**2*y**2 - 24) - 8*scipy.special.j1(omega_c*np.sqrt(x**2 + y**2))*(omega_c**2*x**2 + omega_c**2*y**2 - 6))/(2*np.pi*(x**2 + y**2)**(7/2))

def rotatedCosineWindow(N):  # N = horizontal size of the targeted kernel, also its vertical size, must be odd.
  return np.fromfunction(lambda y, x: np.maximum(np.cos(np.pi/2*np.sqrt(((x - (N - 1)/2)/((N - 1)/2 + 1))**2 + ((y - (N - 1)/2)/((N - 1)/2 + 1))**2)), 0), [N, N])

def circularLowpassKernel(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size, must be odd.
  kernel = np.fromfunction(lambda x, y: omega_c*scipy.special.j1(omega_c*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2))/(2*np.pi*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2)), [N, N])
  kernel[(N - 1)//2, (N - 1)//2] = omega_c**2/(4*np.pi)
  return kernel

def prototype0x(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size, must be odd.
  kernel = np.zeros([N, N])
  for y in range(N):
    for x in range(N):
      kernel[y, x] = h0x(x - (N - 1)/2, y - (N - 1)/2, omega_c)
  return kernel

def prototype0y(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size, must be odd.
  return prototype0x(omega_c, N).transpose()

def prototype1x(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size, must be odd.
  kernel = np.zeros([N, N])
  for y in range(N):
    for x in range(N):
      kernel[y, x] = h1x(x - (N - 1)/2, y - (N - 1)/2, omega_c)
  return kernel

def prototype1y(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size, must be odd.
  return prototype1x(omega_c, N).transpose()

N = 321  # Horizontal size of the kernel, also its vertical size. Must be odd.
window = rotatedCosineWindow(N)

# Optional window function plot
#plt.imshow(window, vmin=-np.max(window), vmax=np.max(window), cmap='bwr')
#plt.colorbar()
#plt.show()

omega_c = np.pi/8  # Cutoff frequency in radians <= pi
lowpass = circularLowpassKernel(omega_c, N)
kernel0x = prototype0x(omega_c, N)
kernel0y = prototype0y(omega_c, N)
kernel1x = prototype1x(omega_c, N)
kernel1y = prototype1y(omega_c, N)

# Optional kernel image save
plt.imsave('lowpass.png', plt.cm.bwr(plt.Normalize(vmin=-lowpass.max(), vmax=lowpass.max())(lowpass)))
plt.imsave('kernel0x.png', plt.cm.bwr(plt.Normalize(vmin=-kernel0x.max(), vmax=kernel0x.max())(kernel0x)))
plt.imsave('kernel0y.png', plt.cm.bwr(plt.Normalize(vmin=-kernel0y.max(), vmax=kernel0y.max())(kernel0y)))
plt.imsave('kernel1x.png', plt.cm.bwr(plt.Normalize(vmin=-kernel1x.max(), vmax=kernel1x.max())(kernel1x)))
plt.imsave('kernel1y.png', plt.cm.bwr(plt.Normalize(vmin=-kernel1y.max(), vmax=kernel1y.max())(kernel1y)))
plt.imsave('kernelkey.png', plt.cm.bwr(np.repeat([(np.arange(321)/320)], 16, 0)))

图1.圆形对称带限滤波器脉冲响应的彩色映射1：1比例图，具有截止频率$ \ omega_c = \ pi / 8 $。颜色键：蓝色：负数，白色：零，红色：最大。
截止频率为$ \ omega_c = \ pi / 8 $，顺序为：$ h_ {0x} $，$ h_ {0y} $，$ h_ {1x} $，$ h_ {0y} $。颜色键：蓝色：最小值，白色：零，红色：最大值。

方向性边缘检测器可以构造为这些值的加权和。在Python中（续）：

composite = kernel0x-4*kernel1x
plt.imsave('composite0.png', plt.cm.bwr(plt.Normalize(vmin=-composite.max(), vmax=composite.max())(composite)))
plt.imshow(composite, vmin=-np.max(composite), vmax=np.max(composite), cmap='bwr')
plt.colorbar()
plt.show()

composite = (kernel0x+kernel0y) + 4*(kernel1x+kernel1y)
plt.imsave('composite45.png', plt.cm.bwr(plt.Normalize(vmin=-composite.max(), vmax=composite.max())(composite)))
plt.imshow(composite, vmin=-np.max(composite), vmax=np.max(composite), cmap='bwr')
plt.colorbar()
plt.show()

图3.定向边缘检测内核构造为图2内核的加权总和。颜色键：蓝色：最小值，白色：零，红色：最大。

与渐变滤镜（图2的前两个滤镜）相比，图3的滤镜应针对连续边缘进行更好的调整。高斯滤波器

由于严格的频带限制，图2的滤波器有很多振荡。像高斯导数滤波器一样，一个更好的凝视点可能是高斯函数。相对而言，它们在数学上要容易得多。让我们尝试一下。我们从高斯“低通”滤波器的脉冲响应定义开始：

$$ h（x，y，\ sigma）= \ frac {e ^ {-\ displaystyle \ frac {x ^ 2 + y ^ 2} {2 \ sigma ^ 2}}} {2 \ pi \ sigma ^ 2}。\ tag {4} $$

我们应用Eq的运算符。 2到$ h（x，y，\ sigma）$并将每个过滤器$ h _ {..} $通过以下方式归一化：

$$ \ int _ {-\ infty} ^ {\ infty} \ int_ {-\ infty} ^ {\ infty} h _ {..}（x，y，\ sigma）^ 2 \，dx \，dy = 1 \\ tag {5} $$

$$ \ begin {align} h_ {0x}（x，y，\ sigma）＆= 2 \ sqrt {2 \ pi}σ^ 2 \ frac {d} {dx} h（x，y，\ sigma）= -\ frac {\ sqrt {2}} {\ sqrt {\ pi}σ^ 2} xe ^ {-\ displaystyle \ frac {x ^ 2 + y ^ 2} {2σ^ 2}}，\\
h_ {0y}（x，y，\ sigma）＆= h_ {0x}（y，x，\ sigma），\\
h_ {1x}（x，y，\ sigma）＆= \ frac {2 \ sqrt {3 \ pi}σ^ 4} {3} \ left（\ left（\ frac {d} {dx} \ right）^ 3-3 \ frac {d} {dx} \ left（\ frac {d} {dy} \ right）^ 2 \ right）h（x，y，\ sigma）=-\ frac {\ sqrt {3}} {3 \ sqrt {\ pi}σ^ 4}（x ^ 3 -3xy ^ 2）e ^ {-\ displaystyle \ frac {x ^ 2 + y ^ 2} {2σ^ 2}}，\\
h_ {1y}（x，y，\ sigma）＆= h_ {1x}（y，x，\ sigma）。\ end {align} \ tag {6} $$

我们想根据它们的加权总和来构造a的脉冲响应最大特异性的垂直边缘检测器滤波器$ S $是相对于可能的边缘旋转角度$ \ beta $和可能的边缘偏移$ s $的平均灵敏度相对于在可能的边缘偏移$ s $的垂直边缘的平均灵敏度：

$$ S = \ frac {2 \ pi \ displaystyle \ int _ {-\ infty} ^ {\ infty} \ Bigg（\ int _ {-\ infty} ^ {\ infty} \ bigg（ \ int _ {-\ infty} ^ { s} h_x（x，y，\ sigma）dx-\ int_ {s} ^ {\ infty} h_x（x，y，\ sigma）dx \ bigg）dy \ Bigg）^ 2ds}
{\ Bigg （\ displaystyle \ int _ {-\ pi} ^ {\ pi} \ int _ {-\ infty} ^ {\ infty} \ bigg（\ int _ {-\ infty} ^ {\ infty} \ Big（\ int _ {-\ infty} ^ {s} h_x \ big（\ cos（\ beta）x- \ sin（\ beta）y，\ sin（\ beta）x + \ cos（\ beta）y \ big）dx \\-\ displaystyle \ int_ {s} ^ {\ infty} h_x \ big（\ cos（\ beta）x-\ sin（\ beta）y，\ sin（\ beta）x + \ cos（\ betay \ big）dx \ Big）dy \ bigg）^ 2ds \，d \ beta \ Bigg）}。\ tag {7} $$

我们只需要加权总和$ h_ {0x} $，且方差$ \带有最佳方差的sigma ^ 2 $和$ h_ {1x} $。事实证明，$ S $通过脉冲响应而最大化：

$$ \ begin {align} h_x（x，y，\ sigma）＆= \ frac {\ sqrt {7625-2440 \ sqrt {5}}} {61} h_ {0x}（x，y，\ sigma）-\ frac {2 \ sqrt {610 \ sqrt {5}-976}} {61} h_ {1x}（x，y ，\ sqrt {5} \ sigma）\\
＆=-\ frac {\ sqrt {（15250-4880 \ sqrt {5}}} {61 \ sqrt {\ pi}σ^ 2} xe ^ {-\ displaystyle \ frac {x ^ 2 + y ^ 2} { 2σ^ 2}} + \ frac {\ sqrt {1830 \ sqrt {5}-2928}} {4575 \ sqrt {\ pi}σ^ 4}（2x ^ 3-6xy ^ 2）e ^ {-\ displaystyle \ frac {x ^ 2 + y ^ 2} {10σ^ 2}} \\
＆= \ frac {2 \ sqrt {\ pi}σ^ 2 \ sqrt {15250-4880 \ sqrt {5}} } {61} \ frac {d} {dx} h（x，y，\ sigma）-\ frac {100 \ sqrt {\ pi}σ^ 4 \ sqrt {1830 \ sqrt {5}-2928}} {183 } \ left（\ left（\ frac {d} {dx} \ right）^ 3-3 \ frac {d} {dx} \ left（\ frac {d} {dy} \ right）^ 2 \ right）h （x，y，\ sqrt {5} \ sigma）\\
＆\约3.8275359956049814 \，\ sigma ^ 2 \ frac {d} {dx} h（x，y，\ sigma）-33.044650082417731 \， \ sigma ^ 4 \ left（\ left（\ left（\ frac {d} {dx} \ right）^ 3-3 \ frac {d} {dx} \ left（\ frac {d} {dy} \ right）^ 2 \ right）h（x，y，\ sqrt {5} \ sigma），\ end {align} \ tag {8} $$

也通过等式5进行了归一化。具有$ S = \ frac {10 \ times5 ^ {1/4}} {9} $的特异性，与之相比，first- $的特异性$ S = 2 $相对于$ x $的高斯导数滤波器。公式8的最后一部分具有规范化兼容w来自Python的scipy.ndimage.gaussian_filter的第2个可分离的2d高斯导数滤波器：

图4.的颜色映射1：1比例图，顺序为：二维高斯函数，相对于$ x $的高斯函数的导数，微分运算符$ \ big（\ frac {d} {dx} \ big）^ 3-3 \ frac {d} {dx} \ big（\ frac {d} {dy} \ big）^ 2 $应用于高斯函数，即最优两分量高斯衍生的垂直边缘检测滤波器$ h_x（x，y，\ sigma）$。 8.每个高斯的标准偏差为$ \ sigma = 8 $，除了最后一个图中的六边形分量的标准偏差为$ \ sqrt {5} \ times8 $。色键：蓝色：最小值，白色：零，红色：最大值。

要继续...

#4 楼

宁可提高性能，但应能获得所需的精度：


边缘检测图像
将图像变换到一个空间，在该空间中您具有足够的像素以达到所需的精度。
因为有足够的正交线；霍夫空间中的图像将包含位于两条线上的最大值。这些很容易检测，并为您提供所需的角度。

#5 楼

我已经进行了基本的调整，将opencv的Hough转换示例调整为您的用例。这个主意不错，但是由于图像的前卫性质已经使您的图像已经有很多边缘，因此边缘检测不会带来太多好处。

因此，我在上面的示例中所做的是


省略边缘检测
将输入图像分解为彩色通道并分别进行处理
计算特定角度下线条的出现（量化角度并取模后） 90°，因为您有足够的直角）
合并色彩通道的计数器
校正这些旋转

可以做些什么来进一步提高估计质量（当您'将在下面看到，最高的猜测是不正确的-第二个是）可能等于将图像转换为代表不同材料之间实际差异的灰度图像的最佳方式-显然，RGB通道并不是最佳的。您是半导体专家，因此请找到一种以最大程度地区别例如金属化和硅。

我的jupyter笔记本在这里。请参阅下面的结果。

要提高角度分辨率，请增加QUANT_STEP变量和hough_transform调用中的角度精度。我没有，因为我希望这段代码能在20分钟内完成编写，因此不想花一分钟时间进行计算。

import cv2
import numpy
from matplotlib import pyplot
import collections

QUANT_STEPS = 360*2

def quantized_angle(line, quant = QUANT_STEPS):
    theta = line[0][1]
    return numpy.round(theta / numpy.pi / 2 * QUANT_STEPS) / QUANT_STEPS * 360 % 90

def detect_rotation(monochromatic_img):
    # edges = cv2.Canny(monochromatic_img, 50, 150, apertureSize = 3) #play with these parameters
    lines = cv2.HoughLines(monochromatic_img, #input
                           1, # rho resolution [px]
                           numpy.pi/180, # angular resolution [radian]
                           200) # accumulator threshold – higher = fewer candidates
    counter = collections.Counter(quantized_angle(line) for line in lines)
    return counter

img = cv2.imread("/tmp/HIKRe.jpg") #Image directly as grabbed from imgur.com
total_count = collections.Counter()
for channel in range(img.shape[-1]):
    total_count.update(detect_rotation(img[:,:,channel]))

most_common = total_count.most_common(5)

for angle,_ in most_common:
    pyplot.figure(figsize=(8,6), dpi=100)
    pyplot.title(f"{angle:.3f}°")
    rotation = cv2.getRotationMatrix2D((img.shape[0]/2, img.shape[1]/2), -angle, 1)
    pyplot.imshow(cv2.warpAffine(img, rotation, img.shape[:2]))