在二维数组上查找子像素最大值答案

【问题标题】：Find sub-pixel maximum on a 2D array在二维数组上查找子像素最大值
【发布时间】：2020-04-15 09:46:48
【问题描述】：

假设我有一个图像，我想找到一个形状为 3x3 的子数组，其中包含与其他子数组相比的最大总和。

如何在 python 中有效地执行此操作（尽可能快地运行）？如果你能提供一个示例代码那就太好了。

我的具体问题：我想在此热图中提取 blob 中心的位置

我不想只获得最大值，因为这会导致坐标不是很精确。斑点的真实中心实际上可能在 2 个像素之间。因此，最好在多点之间进行加权平均以获得亚像素精度。例如，如果有 2 个点 (x1,y1) 和 (x2,y2) 的值分别为 200 和 100。那么平均坐标将为 x=(200*x1+100*x2)/300 y=(200*y1+100*y2)/300

我的一个解决方案是进行卷积操作。但我认为它不够有效，因为它需要乘法到内核（只包含一个）。我正在寻找一个快速的实现，所以我不能自己循环，因为我不确定它是否会很快。

我想每隔几毫秒对 50 张图像执行此算法。（图像作为一批进来）。具体来说，将这些图像视为输出热图的机器学习模型的输出。为了从这些热图中获取坐标，我需要在高强度坐标之间进行某种加权平均。我的想法是在图像上对 3x3 区域进行加权平均。我也愿意接受其他更快或更优雅的方法。

【问题讨论】：

你已经标记了convolution，所以你显然已经尝试过卷积/过滤。结果如何？您是否能够确定 blob 的范围（可能是边界框？）？这些斑点总是大小相同吗？每张图像只有一个斑点吗？我会尽可能减少搜索区域（或尽可能快地完成）并找到 blob 的加权质心。
OpenCV 或 tensorflow 中的普通卷积实现需要对内核进行乘法运算。它会起作用，但不是最快的。我想要一个不需要做任何乘法的实现，因为只有总和就足够了。另外我认为对于这个问题，也许有比卷积更好的技术，所以我在这里问它。如果您能够在没有内核的情况下实现卷积（并且运行速度很快），那么我也想使用它。 blob 的大小始终相同。理想情况下，应该有一个 blob，但模型可能会输出 2 个 blob（错误）。
加载图像，转换为灰度，Otsu的阈值，执行变形关闭，找到轮廓并选择最大的轮廓区域
试试ndimage.center_of_mass怎么样？
@Alex 看起来可以在 3x3 裁剪的图像周围使用，但不应该在整个图像上使用，因为它会很昂贵（至少理论上，实际上从未尝试过）

标签： image opencv sum 2d convolution

【解决方案1】：

查找“具有最大和的 3x3 形状的子数组”与使用未归一化的 3x3 盒式过滤器过滤后查找图像的最大值相同。因此，它归结为有效地找到图像的最大值，您假设它是基础连续和平滑信号的（可能是“嘈杂”）离散样本 - 因此您希望找到子像素位置。

你真的需要把问题分成两部分：

找到图像最大值的像素位置m=(xm, ym)。这只需要访问图像中的每个像素，并且每个像素进行一次比较，因此它是 O(N)，因此只要您在原始图像分辨率下操作，它就是最佳的。在 OpenCv 中使用 minMaxLoc 函数。
应用您正在使用的图像的任何模型，以在 m 的邻域中找到其（子像素插值）最大值。

澄清第（2）点：你写

我不想只获得最大值，因为这会导致坐标不是很精确。斑点的真实中心实际上可能在 2 个像素之间

虽然直觉上是合理的，但这个断言需要更加精确才能进行计算。也就是说，您需要以数学方式表达您对图像所做的假设，这些假设使您搜索像素采样位置之间的“真实”最大值。

这种假设的一个简单示例是二次平滑度。在这种情况下，您假设在“真实”最大位置的一个小（例如，3x3，5x5）邻域中，图像信号z 很好地近似为二次：

z = A00 dx^2 + A01 dx dy + A11 dy^2 + A02 dx + A12 dy + A22
where:
dx = x - xm; dy = y - ym

由于泰勒级数定理，如果预期基础信号至少是三阶连续且可微的，则此假设是有意义的。在几何上，这意味着您假设（希望？）信号看起来像接近最大值的二次曲面（抛物面或椭圆体）。

然后，您可以为m 邻域中的每个像素计算上述方程，替换z 的实际图像值，从而在未知 Aij 中获得一个线性系统，其中的方程数与那里一样多是相邻像素（因此即使是 3x3 邻域也会产生过度约束的系统）。以最小二乘的方式求解系统可以得到“最优”系数 Aij。该模型预测的理论最大值是一阶偏导数消失的地方：

del z / del dx = 2 A00 dx + A01 dy = 0
del z / del dy = A01 dx + 2 A11 dy = 0

这是两个未知(dx, dy) 中的线性系统，求解它会产生最大值的估计位置，并通过上述z 的方程，得到最大值处的预测图像值。

就计算成本而言，与遍历中等大小的图像相比，所有此类模型估计都非常快。

【讨论】：

我认为您找到最大像素然后围绕该像素取平均值的想法已经足够了。

【解决方案2】：

对不起，我没有完全理解你最后一段的意思，所以我刚刚停在我得到所有坐标都具有最大值的地方。我使用cv2.filter2D 对阈值图像进行卷积，然后使用np.amax 和np.where 找到具有最大值的坐标。

import cv2
import numpy as np
from timeit import default_timer as timer

img = cv2.imread('blob.png', 0)
start = timer()
_, thresh = cv2.threshold(img, 240, 1, cv2.THRESH_BINARY)
mask = np.ones((3, 3), np.uint8)
res = cv2.filter2D(thresh, -1, mask)
result = np.where(res == np.amax(res))
end = timer()
print(end - start)

我不知道它是否像你想要的那样高效，但输出是0.0013461999999435648 s

附：您提供的图像有一个白色边框，我必须为此方法裁剪。

【讨论】：

您的想法与我的想法相似，您尝试对图像进行卷积，然后从新图像中获得最大结果。尽管您用来计算的时间非常长（1 毫秒）。我希望对 50 张图像执行此过程，其中每张图像的形状约为 100x100 大约每秒 50 次。简而言之，每秒将卷积大约 2,500 张图像。
我添加了更多解释来阐明我对坐标加权平均的含义。你可以在问题里再看一遍。

【解决方案3】：

一种方法是二次采样图像并找到所需点的邻域。您可以通过对所有像素进行循环而不是在例如每 5 个像素（row=row+5andcol=col+5 在循环中）。找到附近位置后，考虑该位置周围的特定邻域，并对该特定作物的整个像素进行循环以找到确切位置。

【讨论】：

【解决方案4】：

根据我对图像处理的了解，要获得适用于任何单个 blob 的可靠结果，请执行以下步骤：

将图像设为灰度（如果尚未设置）（像素值 0-255）
对图像进行归一化，使像素强度覆盖 0-255 的整个范围
将图像转换为二进制（像素为 0 或 1） - 这可以通过阈值处理来实现，例如应用强度小于或等于 127 的任何像素的强度为 0 的规则其他任何东西的强度为 1
求所有值为“1”的像素的加权平均值

或

对图像进行腐蚀，直到剩下 2 个像素或 1 个像素。

案例 1

如果你有两个像素，那么你需要找到两个像素的 u 和 v 坐标。斑点的中心将是像素的 u 和 v 坐标之间的中点。

案例 2

如果还剩一个像素，那么该像素的坐标就是中心点。

————

您提到过在 Python 中快速实现这一点：

Python 在设计上是一种解释型语言，因此它是逐行执行的，因此不太适合图像处理等高度迭代的任务。但是，您可以使用诸如 OpenCV (https://docs.opencv.org/2.4/index.html) 之类的库，它是用 C 编写的，除了使您手头的任务变得更容易之外，还可以减轻这种情况。 OpenCV 还为此身份提供了我上面列出的所有步骤的解决方案，因此您应该能够相当快地获得可靠的解决方案，尽管我不能确定它是否会达到您的目标每隔几毫秒就有 50 张图像。其他要考虑的因素是您正在处理的图像的大小。这将成倍增加处理负载。

更新

我刚刚找到了一篇很好的文章，实际上与我的步骤过程相呼应：

https://www.learnopencv.com/find-center-of-blob-centroid-using-opencv-cpp-python/

更重要的是，它还表示在数学上找到质心的公式为：

c = (1/n)sigma(n, i = 1, x_i)

但是这比我在这里写的更好。

【讨论】：