基于多进程的焦点堆栈图像融合方法

Method of Focus Stacks’ Fusion Based on Multi-Processing

摘要: 传统上为了拍摄前景与后景均清晰的图像，相机光圈应当尽可能的小。但是受到衍射极限的影响，过小的光圈会造成图像发生畸变。然而通过将多张相同构图但不同焦距的图像合成，可克服上述困难。本方法包含硬件设备和算法设计，重点在算法设计上，包括高斯滤波、拉普拉斯算子、均值滤波以及探测点的极大值映射，最终实现图像的合成。除此而外，算法中还加入了多进程的思想，使得算法在获取同样锐度图像的条件下，程序运行时间极大缩短，这通过仿真实验得以论证。本论文提出的算法适用于多核CPU的处理，具备普适性。

Abstract: Traditionally, the camera’s aperture should be as small as possible in order to make the foreground and the background clear, respectively. However, affected by the diffraction limit, tiny apertures can cause the image distortion. The above difficulties can be overcome by merging multiple images with the same compositions but with different focus. The thesis includes hardware and algorithm design, emphasizing on algorithm design. The image can be eventually fused including in Gaussian filtering, Laplacian operator, mean filtering and maximum mapping of detection points. Besides, in order to greatly shorten the running time of program, the multi-process is added to the algorithm without deprived of any sharpness, demonstrated by simulation. The algorithm is suitable for the multi-core CPU, universally.

1. 引言

在工业检测领域和摄影领域，当为了获得图像中所有细节部分均清晰的效果，传统方法是选用小光圈的相机来进行拍摄。然而当光圈过小时，光通过光圈时就会发生衍射效应。根据瑞利判据，光学系统存在衍射极限，即艾里斑的半径x [1]：

$x=\frac{1.22\lambda }{d}$ (1)

其中λ为光的波长，d为光圈的直径，f为相机镜头的焦距。当艾里斑尺寸大于像素尺寸时，便会使得图像发生畸变。

为了克服上述困难，考虑采用数字图像处理的技术。通过将同一构图进行多焦距的采样，提取有效信息并将所有图像融合，可最终得到大景深的图像 [2] [3]。正是因为焦点堆栈所具备的优势，因此其应用得到了拓展，斯坦福大学的Marc Levoy等人在显微摄影中成功利用该技术对昆虫的腿部进行成像 [4]，突破了传统显微镜的空间分辨率和角分辨率的衍射极限；Samuel W. Hasinoff等人利用该技术提出了计算3D形状的新方法 [5]，解决了三维重构时对于如毛发等细微物体复杂三维排列的问题 [6] [7] [8]。然而，现今的图像融合算法存在程序运行时间长的问题，这会导致工业检测时间过长，检测效率低的问题；或者是摄影过程中合成时间过长，无法及时地查看图像合成效果，导致图像拍摄过程中存在有部分漏拍的情况时，摄影师无法及时地对拍摄有很好的把控。现对焦点堆栈算法进行算法设计和多进程处理，能够有效的提升图像显示效果和程序运行时间，为工业检测和微距摄影的研究提供了有效的应用拓展。本文的主要思想在于对每一张图像均提取高频信息，通过在每一个像素点位置处选取高频信息的最大值进行图像融合，从而得到合成图像。

2. 项目需求

本项目主要是用来对不同焦距的图像进行快速合成，最终实现在一张图片中所有细节均呈现清晰的效果，即超景深合成图像。运用卷积方法如Sobel算子、Prewitt算子、Roberts算子、Laplacian算子以及Canny算子进行梯度边缘检测，虽然图像检测的精度高，但是程序的运行时间长；而若对图像直接进行简单处理，类似于对图像进行算数平均或者是小波变换，虽然程序运行时间短，但是图像的合成效果不佳。因此本项目考虑在图像检测精度得到保证的情况下，尽可能地压缩程序运行时间。

3. 光路设计

光路设计如图1所示，LD光源发出的光经过待测物，最终被CMOS相机(MV-CA050-20GM)所接收。CMOS相机安置在减震平台的支架上，可做从左往右的平移运动；减震平台的支架可做从上至下的垂直运动以及从前至后的平移运动。CMOS相机先对待测物的最左端聚焦后拍摄(绿色虚线部分)，之后每间隔1 s，通过将减震平台的支架做从下至上的垂直运动，让CMOS相机的焦点往右平移一段距离。此时进行聚焦拍摄，直到CMOS相机运行到待测物最右端，进行最后一次聚焦并拍摄，CMOS相机便停止运行，图像采集完毕。每一次CMOS相机采集到的图像数据均传递给计算机，为焦点堆栈图像的合成做下一步准备。这样在不同焦平面上的物体信息，便通过CMOS记录了下来。

Figure 1. The beam path of schematic setup for focus stacking based on multi-processing

图1. 基于多进程的焦点堆栈系统光路图

4. 软件设计

本系统待处理的图像是通过工业级CMOS相机(MV-CA050-20GM)进行拍摄的，与其配置的镜头(MT05-110-HR-D1.5)为0.5倍，工作距离110 mm。后续图像处理采用python编程语言，因为这能够开源地实现该项目。

4.1. 流程处理

算法是本部分的核心，其流程图如图2所示，包括了图像导入部分、图像处理部分和图像融合部分。在经过图1所示的光路系统对物体信息采集后，计算机中便保存了N张不同焦面，但构图一致的图像。将图像利用多进程进行导入后，便可通过滤波处理以及边缘梯度检测获取到每一张图像的边缘信息，最后通过选取最大梯度信息进行图像融合，便可得到一幅超景深的合成图像。其详细步骤如下。

首先，利用CMOS相机对同一物体从左往右地进行聚焦拍摄，其光路图如图1所示，一共拍摄了N张图片。导入图片时，我们采用opencv的imread函数来对待处理的图像数据进行导入，为后续进行图像处理做准备。

Figure 2. Flow chat of the system for focus stacking based on multi-processing

图2. 基于多进程的焦点堆栈系统流程图

为了避免当光照不均，或高温引起的传感器噪声而带来的高斯噪声影响，我们可通过高斯滤波来减小其对图像处理的影响。若此时我们正准备处理第m张图像 $I_m\left(x,y\right)$，其中 $m\in \left[1,N\right]$，N为图像的总张数，通过卷积操作，我们可得到处理过后的第m张图像 $O_m\left(x,y\right)$，其公式为 [9]

$O_m\left(x,y\right)=I_m\left(x,y\right)\ast K_m\left(x,y,\sigma \right)$ (2)

式中x为输入图像的横向坐标，y为输入图像的纵向坐标， $K_m\left(x,y,\sigma \right)$ 为第m张图像的卷积核， $\ast$ 为卷积运算符。此时的 $K_m\left(x,y,\sigma \right)$ 为第m张图像的高斯卷积核，公式如下：

$K_m\left(x,y,\sigma \right)=\frac{1}{2\pi {\sigma }^2}\exp \left[-\frac{{\left(x-\frac{W-1}{2}\right)}^2+{\left(y-\frac{H-1}{2}\right)}^2}{2{\sigma }^2}\right]$ (3)

其中σ为标准差，W为输入图像 $I_m\left(x,y\right)$ 的宽，H为输入图像 $I_m\left(x,y\right)$ 的高。之后我们可对图像进行梯度边缘检测，通过确定二阶微分为0的点，从而来判断图像的极值点，而这些极值点便为图像中的轮廓，即边缘部分。这里我们采用的是拉普拉斯算子来进行处理，公式如下：

$\begin{array}{c}{\nabla }^2{O}_m=\frac{{\partial }^2{O}_m\left(x,y\right)}{\partial x^2}+\frac{{\partial }^2{O}_m\left(x,y\right)}{\partial y^2}\\ \approx O_m\left(x+1,y\right)+O_m\left(x-1,y\right)+O_m\left(x,y-1\right)+O_m\left(x,y+1\right)-4O_m\left(x,y\right)\end{array}$ (4)

然而由于拉普拉斯算子对噪声并不具备过滤的功能，因此运用该算子来对图像进行处理会造成边缘点的误判，最终导致图像边缘定位不够精准。考虑到后续我们会对所有图像的极大值信息进行比较，而通过拉普拉斯算子检测出来图像的极大值点，这属于高频信息，那么所需处理的噪声便转变为了高频噪声。为了解决上述问题，我们又采用均值滤波来对拉普拉斯算子检测后的图像进行噪点剔除工作，因为均值滤波适合消除图像的高频噪声，具备实现图像模糊、图像平滑的功能。公式如下：

$F_m\left(x,y\right)=O_m\left(x,y\right)\ast K_m\left(x,y\right)$ (5)

其中 $F_m\left(x,y\right)$ 代表了第m张图像经过均值滤波后输出的图像， $O_m\left(x,y\right)$ 代表了第m张图像经过拉普拉斯算子处理过后的图像， $\ast$ 为卷积运算符。 $K_m\left(x,y\right)$ 代表了第m张图像的卷积核，此时的卷积核为均值滤波卷积核，公式如下：

$K_m\left(x,y\right)=\frac{1}{H\times W}{\left[1\right]}_{H\times W}$ (6)

其中W为经过拉普拉斯算子处理过后图像 $O_m\left(x,y\right)$ 的宽，H为经过拉普拉斯算子处理过后图像 $O_m\left(x,y\right)$ 的高。另外均值滤波卷积核是可分离算子，通过卷积核的分离性和卷积的结合律，可减小程序的运算量。可分离的均值滤波卷积核公式为：

$K_m\left(x,y\right)=\frac{1}{1\times W}{\left[1\right]}_{1\times W}\ast \frac{1}{H\times 1}{\left[1\right]}_{H\times 1}=\frac{1}{H\times 1}{\left[1\right]}_{H\times 1}\ast \frac{1}{1\times W}{\left[1\right]}_{1\times W}$ (7)

其中 $\ast$ 为卷积运算符。经过滤波处理和拉普拉斯算子对边缘信息提取后，我们便进入图像融合部分。我们通过把均值滤波处理后的N张图像按照像素点的位置x_i和y_j，进行逐一比较，其中 $\text{1}\le x_i\le W$， $\text{1}\le y_j\le H$， $i,j\in \left[1,N\right]$，从而选取每一个像素点位置所处的最大值点所属的图像标记k， $k\in \left[1,N\right]$。公式如下：

$k=\arg \max {\left(F_m\left(x,y\right)\right)}_{H\times W}$ (8)

之后我们将得到的图像标记k映射到原始图像堆 $I\left(x,y\right)$ 中，按照图像标记k对原始图像堆 $I\left(x,y\right)$ 进行像素点提取，直到提取完整个H × W个像素点。公式如下：

$M{\left(x,y\right)}_{H\times W}=I_k{\left(x,y\right)}_{H\times W}$ (9)

其中 $M\left(x,y\right)$ 为基于多进程的焦点堆栈合成的图像， $I_k\left(x,y\right)$ 为在像素点 $\left(x,y\right)$ 处第k张图像的像素值。这样我们便可实现基于多进程的焦点堆栈图像的合成，合成图像中的每一个像素位置均拥有通过图像处理后的最大像素值，至于多进程的部分我们接下来会谈到。

4.2. 多进程处理

为了加快程序运行速度，我们采用了多进程的方式。如图2所示，多进程方式涵盖了从获取N张原始图像到对图像进行均值滤波处理，基本上涵盖了整个流程的大部分。

线程是操作系统进行运算调度的最小单位，是进程中的实际运行单位。但是由于python有GIL (Global Interpreter Lock)的存在，这会对多线程处理造成很多限制，因此这里不考虑多线程处理。进程指的是可以独立运行的程序单位，当存在多核CPU时，通过多进程处理，就相对于并行处理程序，这可极大提升程序运算效率。

在多进程处理中，进程池的池化处理有着调度的作用。池化处理是指当进程数小于池中规定的最大值时，用户发出的请求可被当作一个新的进程进行并行处理；而当进程数大于或等于池中规定的最大值时，用户发出的请求需被等待，直到某一个进程结束，用户发出的请求方可作为一个新的进程进行并行处理。池化处理的优势在于有效的避免了计算机资源消耗过大，最终导致程序运行效率不高的不利情况。因此，我们采用multiprocessing安装包中的Pool函数来对程序进行调度，从而实现焦点堆栈图像的多进程处理。

5. 实验结果

经过如图1所示的实验装置对待测物拍摄后，最终我们一共获得4张焦距位置不同的图像(2592 pixel × 2048 pixel)。如图3所示我们展示了其中三张图像，可以发现图像中的清晰部分从最左端最终到达最右端，代表了CMOS相机(MV-CA050-20GM)焦点在从左往右的移动。我们采用12核CPU来对图像进行并行处理，CPU型号为Intel酷睿i710750H，GPU型号为GeForce RTX 3060。为了与串行算法相对比，我们采用了另外五种不同算法来作为参照，其效果如图4所示。

(a) 焦点在左方的图像 (b) 焦点在中间的图像 (c) 焦点在右方的图像

Figure 3. Samples

图3. 采集样图

(a) 基于多进程焦点堆栈算法 (b) 拉普拉斯金字塔算法 (c) SIFT算法 + 位运算 (d) SIFT算法 + 最大值映射 (e) 算数平均算法 (f) 小波变换算法

Figure 4. Displays of merging images by different algorithms

图4. 不同算法的合成图像

图4中的图4(a)图是通过本文的算法处理出来的图像；图4(b)图是通过拉普拉斯金字塔算法处理出来的图像，主要思想是先通过对每一张图逐级下采样从而获得高斯金字塔，记录完高斯金字塔每一级下采样后，再上采样与下采样前的差异最终可重构图像；图4(c)图是通过SIFT算法以及位运算来获得的图像，主要思想是在不同的尺度空间上查找关键点，并计算出关键点的方向，在得到每一张图像关键点后，通过把图像关键点的比较用布尔值表示，最终按照布尔值叠加得到重构图像；图4(d)图是通过SIFT算法以及最大值映射来获得的图像，主要思想是获得每一张的关键点后，在每一个像素点的位置确定最大关键点所属的图像编号，最后将图像编号映射到原图像中得到重构图像；图4(e)图是通过算数平均方法来获得图像，主要思想是将所有图像叠加后，除以总张数便可得到重构图像；图4(f)图是通过小波变换来获得图像，主要是通过选择合适小波基来获得图像的高频信息，将处理后的图像叠加后得到重构图像。

为了衡量各算法在程序运行时间、CPU利用率和内存使用率的情况，我们通过表1展示了相应的参数，其中表1展示的算法所处理的图像与图4编号相对应。为了显示不同算法之间程序运行时间的差异，我们通过运行时长耗比C来进行表示，公式如下：

$C=\frac{T_i}{T_a}$ (10)

其中T_i为图4中编号为i所属的算法程序运行时间， $i=\text{b},\text{c},\cdots ,\text{f}$ ；T_a为基于多进程的焦点堆栈算法程序运行时间。

Table 1. Running parameters of different algorithms

表1. 各算法的运行参数

将表1和图4综合来看，我们发现本论文提出的基于多进程的焦点堆栈图像算法在运行时间和图像效果上取得了不错的表现，从程序运行时长耗比来看，只有算数平均算法的程序运行时间比基于多进程的焦点堆栈图像算法的程序运行时间要短，其余算法均超过了基于多进程的焦点堆栈图像算法的程序运行时间，甚至拉普拉斯金字塔算法的程序运行时长耗比是基于多进程的焦点堆栈图像算法的89.4倍。然而从图像效果来看，尽管算数平均算法的程序运行时间比基于多进程的焦点堆栈图像算法的程序运行时间要短，但是最终合成的图像却不够锐利，字迹、图像均有不同程度的模糊。小波变换算法得出的图像也和算数平均算法得出的图像显示效果类似，处理效果不佳。不过，拉普拉斯金字塔算法和以SIFT为主的两种算法均与基于多进程的焦点堆栈图像算法的显示效果差不多，图像足够锐利，易于很好辨认细节。另外，之所以图4中所有合成图像有一定重影或缺失，是因为拍摄图像的过程中CMOS相机有一定振动，导致每一张图像有无法避免的错位。

为了凸显出多进程的优势，我们以本论文的焦点堆栈图像算法来比较多进程与串行的效率，通过利用阿姆达尔定律 [10]，我们可得出并行处理前执行速度与并行处理后执行速度的加速比数值S_p，以此来阐明并行化后的效率提升情况，公式如下 [11]：

$S_p=\frac{T_1}{T_p}$ (11)

其中下标p代表并行处理的CPU核心数，T₁代表并行处理前程序所耗费的时间，T_p代表并行处理后程序所耗费的时间。为了衡量并行处理条件下，多进程的效率数值E_p，我们有如下公式 [12]：

$E_p=\frac{S_p}{p}=\frac{T_1}{T_{p}p}$ (12)

其中p为并行处理的CPU核心数。多进程处理的图像与串行处理的图像如图5所示，串行处理与并行处理的图像对比参数如表2所示。

(a) 基于串行的焦点堆栈算法 (b) 基于多进程的焦点堆栈算法

Figure 5. Merged image by serial algorithm and parallel algorithm, respectively

图5. 串行与并行算法合成的图像

Table 2. Running parameters with serial algorithm and parallel algorithm

表2. 串行算法与并行算法的运行参数

从图5的图像显示效果来看，基于多进程的焦点堆栈算法和基于串行的焦点堆栈算法得到的图像锐度一致。从表2看出当使用串行的焦点堆栈算法时，CPU的利用率会大幅度下降到只有8.6%，而内存使用率会有少许上升至50.5%。从并行加速比来看，串行程序所消耗的时间是并行程序所消耗时间的7.74倍，平均每一个CPU核的运行效率为64.5%。这意味着，通过多进程的改进，尽管CPU利用率会有一定程度的上升，但是程序运行时间会被极大减少，并且最终处理出来的图像与串行程序处理出来的图像无异，因此这是极为具备优势的，尤其是在CPU核心数越多的情况下。

6. 结论

本文提出的焦点堆栈图像算法是为了便于在工业领域对于瑕疵品进行检测或在摄影技术上有所提升，通过将焦距不同的图像进行合成，可最终得到超景深的图像，这便于我们从图像中提取有效且大量的信息。本文根据多核CPU的思路提出了基于多进程的焦点堆栈图像融合方法，通过将不同的算法进行对比，如拉普拉斯金字塔算法、基于SIFT的算法、图像算数平均算法和小波变换算法，若将图像显示效果和程序运行时间综合来看，本文提出的算法具备很大优势。为了探究多进程方法对程序运行时间的提升情况，我们又对多进程算法与串行算法进行对比，发现相比于串行算法，多进程算法的CPU利用率虽然有一定程度的上升，但是在程序运行时间上的确缩短了不少，串行程序运行时间是并行程序运行时间的7.74倍。随着计算机CPU核心数的增加，并行算法与串行算法的运行时间应当还会进一步拉开，这对于处理张数越多的图像更具备优势。然而，本论文提出的算法对合成有一定位移的图像却稍显不足，图像容错性不够，这在未来是应考虑改进的方向。