如何检测图像中的光点？

Question

为了从某个地理区域检测inselbergs，我下载了该区域的地形图像（地形和倾斜度）。倾斜度图像似乎最适合这项任务。

应用高斯模糊（或更快的普通模糊，使用 ImageMagick）后，图像似乎已准备好进行自动检测。

现在我想知道在黑色背景上检测这些白色污点的最佳/最快方法。我的第一个想法是使用一个简单的函数（没有外部库），它的工作原理类似于绘图程序的“桶画”功能，计算比预定义阈值更轻的对象的面积。问题是：由于图像非常大，5400x3600 像素，通常的递归函数可能会导致堆栈溢出，特别是在像下面这样的大山上。

那么，有什么建议吗？我认为理想的语言可能是 C++（也许是 JavaScript）。我不习惯 Python。也许使用像 OpenCV 这样的库会变得更容易，但也许问题太小而无法请求外部库（包括隐含的学习曲线）。

---

TIFF 图像来自here（我可以将它们转换为其他格式进行处理）。示例图像位于 17S42 四边形（“Declividade”选项），靠近坐标 18°S 41°W。在中间图像（上图）中，每个“白球”都是一个 inselberg。精度取决于所选的灰度阈值。

Answer 1

一种方法是使用 ImageMagick 连接组件并在适当的阈值后从每个白色区域中提取质心。

这是 ImageMagick 6 中的命令行。但如果需要，您可以使用 C+ API。

毛刺图像：

convert blur.png -auto-level -threshold 65% -type bilevel \
-define connected-components:verbose=true \
-define connected-components:mean-color=true \
-connected-components 4 binary.png |\
grep "gray(255)" | awk '{print $3}'

质心（像素坐标）：

915.9,367.4
873.4,76.8
936.4,269.5
664.8,212.0
598.2,554.3
954.0,59.3
732.6,446.0
353.5,696.4
69.7,848.0
750.4,357.6
537.5,609.9
832.0,230.9
517.9,539.1

阈值图像：

Answer 2

简而言之，使用你的方法你会得到一个非常糟糕的大小估计。模糊区域的边缘，无论您选择什么阈值，都与您要测量的 inselbergs 的边缘不对应。 相反，我建议你按照下面的食谱。为此，我在 Python 中使用了DIPlib（免责声明：我是作者）。 Python 绑定是 C++ 库上的一个薄层，将下面的 Python 代码转换为 C++ 相当简单（我更容易在 Python 中以交互方式开发它）。

我从您提供的链接（而不是倾斜度）下载了原始高度数据。 DIPlib 可以直接读取浮点值的 TIFF 文件，因此不需要任何特殊的转换。我裁剪了一个类似于 OP 用于此演示的区域，但没有理由不将该方法应用于整个图块。

import PyDIP as dip

height = dip.ImageRead('17S42_ZN.tif')
height.SetPixelSize(0.000278, 'rad')  # not really radian, but we don't have degrees
height = height[3049:3684, 2895:3513];

代码还根据 TIFF 文件中的数据设置像素大小（使用弧度单位，因为 DIPlib 不做度数）。

接下来，我应用具有特定直径（25 像素）的顶帽滤镜。这将隔离所有直径为 25 像素或更小的峰。根据您认为的 inselberg 的最大宽度调整此尺寸。

local_height = dip.Tophat(height, 25)

实际上，结果是一个局部高度，高于某个基线的高度由过滤器的大小决定。

接下来，我应用滞后阈值（双阈值）。这会产生一个二值图像，在基线以上 100m 处进行阈值处理，其中地形高于该基线 200m。也就是说，我决定一个 inselberg 应该至少高于基线 200m，但在 100m 处将它们中的每一个都切断。在这个高度，我们将测量尺寸（面积）。再次，根据您的需要调整阈值。

inselbergs = dip.HysteresisThreshold(local_height, 100, 200)

现在剩下的就是测量我们找到的区域：

labels = dip.Label(inselbergs)
result = dip.MeasurementTool.Measure(labels, features=['Size', 'Center'])
print(result)

这个输出：

   |       Size |                  Center | 
-- | ---------- | ----------------------- | 
   |            |       dim0 |       dim1 | 
   |     (rad²) |      (rad) |      (rad) | 
-- | ---------- | ---------- | ---------- | 
 1 |  1.863e-05 |     0.1514 |    0.01798 | 
 2 |  4.220e-05 |     0.1376 |    0.02080 | 
 3 |  6.214e-05 |    0.09849 |    0.04429 | 
 4 |  6.492e-06 |     0.1282 |    0.04710 | 
 5 |  3.022e-05 |     0.1354 |    0.04925 | 
 6 |  4.274e-05 |     0.1510 |    0.05420 | 
 7 |  2.218e-05 |     0.1228 |    0.05802 | 
 8 |  1.932e-05 |     0.1420 |    0.05689 | 
 9 |  7.690e-05 |     0.1493 |    0.06960 | 
10 |  3.285e-05 |     0.1120 |    0.07089 | 
11 |  5.248e-05 |     0.1389 |    0.07851 | 
12 |  4.637e-05 |     0.1096 |    0.09016 | 
13 |  3.787e-05 |    0.07146 |     0.1012 | 
14 |  2.133e-05 |    0.09046 |    0.09908 | 
15 |  3.895e-05 |    0.08553 |     0.1064 | 
16 |  3.308e-05 |    0.09972 |     0.1143 | 
17 |  3.277e-05 |    0.05312 |     0.1174 | 
18 |  2.581e-05 |    0.07298 |     0.1167 | 
19 |  1.955e-05 |    0.04038 |     0.1304 | 
20 |  4.846e-05 |    0.03657 |     0.1448 | 

（请记住，“rad”实际上是度数。）以平方度为单位的面积有点奇怪，但您可以将其转换为平方米，因为您知道地球上的位置。实际上，在计算之前将像素大小转换为米可能更容易。

这里给出的“中心”值是相对于左上角的像素，如果我们没有裁剪瓷砖开始，我们可以添加瓷砖的坐标（可以从相应的TIFF 文件中的标记）：（-42.0，-17.0）。

---

在 C++ 中，代码应如下所示：

#include <diplib/simple_file_io.h>
#include <diplib/morphology.h>
#include <diplib/segmentation.h>
#include <diplib/regions.h>
#include <diplib/measurement.h>

//...

dip::Image height = dip::ImageRead("17S42_ZN.tif");
height.SetPixelSize(0.000278 * dip::Units::Radian());
height = height.At(dip::Range(3049, 3684), dip::Range(2895, 3513));

dip::Image local_height = dip::Tophat(height, 25);

dip::Image inselbergs = dip::HysteresisThreshold(local_height, 100, 200);

dip::Image labels = dip::Label(inselbergs);
dip::MeasurementTool measurementTool;
dip::Measurement result = measurementTool.Measure(labels, {}, {"Size", "Center"});
std::cout << result;