【沥血整理】灰度（二值）图像重构算法及其应用（morphological reconstruction）。
不记得是怎么接触并最终研究这个课题的了，认识我的人都知道我是没有固定的研究对象的，一切看运气和当时的兴趣。本来研究完了就放在那里了，一直比较懒的去做总结，但是想一想似乎在网络上就没有看到关于这个方面的资料，能搜索到的都是一些关于matlab相关函数的应用，决定还是抽空趁自己对这个算法还有点记忆的时候写点东西吧，毕竟这个算法还有一些应用是值得回味和研究的。而且也具有一定的工程价值。
怎么说呢，其实在很早浏览matlab的图像处理工具箱的时候，就无数次的看到过这些函数，但是无奈当时不知道他们有什么用，就没怎么鸟他， 其实M还是很重视他们，这个从他们在工具箱里占用的函数列表篇幅里就能完美的看的出：
在Intensity and Binary Images功能区域里，除去红线划除掉的2个算法和这个重构没有啥关系，其他的都可以认为是imreconstruct衍生出来的产物。
一、重点函数
我们先重点来看下imrecontruct函数的介绍，重要的文字描述有以下内容：
imreconstruct Morphological reconstruction
Syntax
IM = imreconstruct(marker,mask)
IM = imreconstruct(marker,mask,conn)
Description
IM = imreconstruct(marker,mask) performs morphological reconstruction of the image marker under the image mask. marker and mask can be two intensity images or two binary images with the same size. The returned image IM is an intensity or binary image, respectively. marker must be the same size as mask, and its elements must be less than or equal to the corresponding elements of mask.
By default, imreconstruct uses 8-connected neighborhoods for 2-D images and 26-connected neighborhoods for 3-D images. For higher dimensions, imreconstruct uses conndef(ndims(I),'maximal').
IM = imreconstruct(marker,mask,conn) performs morphological reconstruction with the specified connectivity. conn can have any of the following scalar values.
他的意思是从用户提供的mask图像中重建原图，似乎讲的很模糊啊，有点不知所云。
看看有没有更多的信息呢，在M对应帮助文档的尾部，有这样一段话：
Algorithms
imreconstruct uses the fast hybrid grayscale reconstruction algorithm described in [1].
References
[1] Vincent, L., "Morphological Grayscale Reconstruction in Image Analysis: Applications and Efficient Algorithms," IEEE Transactions on Image Processing, Vol. 2, No. 2, April, 1993, pp. 176-201.
看到没有，有参考资料，我就喜欢matlab和halcon这样比较开明的软件，即使不提供源代码，他也会提供一些非常有用的资料，比如参考资料，比如论文或者计算的数学公式。
拜托几个朋友，FQ终于下载到了这个对应的文章，于是像饥饿的老鼠一样扑在面包上苦心的研究了几天，终于有所成。
要理解这个算法的数学原理呢，我想我是解释不清楚，还是要专心的把那个论文打印出来，然后一个人关在小房间里慢慢的啃里面的数学公式，我呢，看了几天，也就获得这个原理的那么一点点朦脓的感觉，在这里可不敢随便解释。也就拿论文中的一个简单的图来说下事。
以一个简答的二值图的重构来说明下这个算法大概在干什么，以下图为例：
这个定义简单翻译就是从标记图像J中重建图像I的过程为，在I中找到包含至少一个J像素的连续区域。
那么在左侧图中，1、2、3处是我们标记的位置图J，原图就是去掉1、2、3哪些黑色的（对应部分恢复为周边底色），根据这个定义，由1、2、3这个对应的位置去找包含他们的目标，最终找到右侧的结果图，而抛弃掉不含有他们的那些目标。
如果给你一个这样的需求，你如何写代码呢。
这个定义只适合理解的他的意思和需求，但是还是无法从定义中寻找代码的书写方式的。 那么文章里又给出了第二种定义的方式：
其中 
一头雾水是吧，我也是看的一头雾水。
好了，我不装了，我摊牌了，其实就是这么个意思，要从标记的图像中恢复图像，怎么办呢，我们进行迭代，每次迭代中呢，先求Marker图像的3*3领域的最大值（standard dilation of size one），然后再把这个最大值和原始图像求最小值得到一副临时图像，不断的重复这个过程，知道图像没有任何的变换，则结果计算，这个没有任何变化的图像就是我们重构后的图像。针对上面的二值图像，好好的想一想，是不是这个过程确实可以实现刚刚定义1里的需求呢，静下心来想一想哦。
扩展到灰度图像，似乎上述定义1就成了无法理解的行为了，确实是这样的，但是我们如果不管这些，定义的操作从程序的角度来说灰度图也是毫无区别的。那么在论文中也是这样推广到灰度的。
论文里也给出代码的实现的定义：
很多人基础差，写不出代码，我好人做到底，按照上述的意思一个简单的代码如下所示：
// 标准的可并行版本的8领域重建算法,虽可极度优化，但是无赖迭代次数太多
int IM_ReConstruct_Standard_Connected8(unsigned char *Src, unsigned char *Marker, unsigned char *Dest, int Width, int Height, int Stride)
{
int Channel = Stride / Width;
if (Src == NULL) return IM_STATUS_NULLREFRENCE;
if ((Width <= 0) || (Height <= 0)) return IM_STATUS_INVALIDPARAMETER;
if (Channel != 1) return IM_STATUS_NOTSUPPORTED;
int Status = IM_STATUS_OK;
int MaxIteration = 65536, Iteration = 0;
unsigned char *Temp = (unsigned char *)malloc(Height * Stride * sizeof(unsigned char));
if (Temp == NULL) return IM_STATUS_NULLREFRENCE;
// 为了不破坏Marker的数据，对其做个备份
memcpy(Temp, Marker, Height * Stride * sizeof(unsigned char));
while (Iteration < MaxIteration)
{
Iteration++;
for (int Y = 0; Y < Height; Y++)
{
for (int X = 0; X < Width; X++)
{
int X0 = X - 1 >= 0 ? X - 1 : 1;
int X2 = X + 1 < Width ? X + 1 : Width - 2;
int Y0 = Y - 1 >= 0 ? Y - 1 : 1;
int Y2 = Y + 1 < Height ? Y + 1 : Height - 2;
int Index0 = Y0 * Stride;
int Index1 = Y * Stride;
int Index2 = Y2 * Stride;
int V0 = Temp[Index0 + X0];
int V1 = Temp[Index0 + X];
int V2 = Temp[Index0 + X2];
int V3 = Temp[Index1 + X0];
int V4 = Temp[Index1 + X];
int V5 = Temp[Index1 + X2];
int V6 = Temp[Index2 + X0];
int V7 = Temp[Index2 + X];
int V8 = Temp[Index2 + X2];
int Max1 = IM_Max(IM_Max(V0, V1), IM_Max(V2, V3));
int Max2 = IM_Max(IM_Max(V5, V6), IM_Max(V7, V8));
int Max = IM_Max(IM_Max(Max1, Max2), V4);
Dest[Index1 + X] = Max;
}
}
int DiffAmount = 0;
for (int Y = 0; Y < Height; Y++)
{
int Index = Y * Stride;
for (int X = 0; X < Width; X++)
{
int Value = IM_Min(Dest[Index], Src[Index]);
if (Temp[Index] != Value) DiffAmount++;
Temp[Index++] = Value;
}
}
if (DiffAmount == 0) break;
}
memcpy(Dest, Temp, Height * Stride * sizeof(unsigned char));
free(Temp);
return IM_STATUS_OK;
}
是不是很简单。
但是一看这个代码就指导，可能速度不是很好，因为每次迭代基本上也就是处理 Marker图像外围的一个像素左右宽度的图，一般都需要迭代很多次才能收敛。
接着论文有提出了一种改进的方式（序列化方式重建）：
其中额NG+和NG-的意思如下图：
通过上面的代码IM_ReConstruct_Standard_Connected8你能否能写出这个版本的算法呢，我就没有必要提供了吧，不过这个虽然有所提高速度，但还是很慢。
后续论文还给出了2个优化方面的代码，一个叫reconstruction using a queue of pixels，这个的算法基础呢，是什么呢，就是上面的重建工作，其实没有必要针对marker图像J的每一个像素，而只需要针对边缘进行处理，这个边缘要是广义的边缘，对于二值图，就是如果J中一个像素是1，那么主要他3*3领域内有1个像素值不为1，他就是一个边缘，而对于灰度图，这个概念得以扩展，指的是如果一个像素是其3*3领域的最大值，那他距考虑为边缘。 我们找到marker图像中所有的边缘点，并把它们加入到一个叫FIFO（First-In-First-Out ）的数据结构中，好像C++里有这种，似乎是叫dqueue。然后不断的迭代处理，指导收敛，比如灰度的处理方法如下所示：
论文最后提出了一种混合（翻译为杂交总觉得好畜生）的算法，即结合序列化算法和上述FIFO一起实现，第一步先用序列化算法执行一次循环，然后在用FIFO来处理。
matlab里也是使用的这种算法来实现，但是我在用这种算法实现时，发现总是有些图和matlab的结果不一致，但有些图又正常，一直没有找到问题的症结所在，同时我发现第三种写法实际上也没有比最后的混合算法慢多少，而他的结果非常稳定，和matlab基本完全一致，因此，我选择第三种算法的实现。
二、算法的直接和间接应用
理论部分讲完，现在在来谈谈这个算法比较精彩的部分。
1、基本功能
接下来我们来了解下这个算法的直接应用，我们先来验证下前面举例的那个图片吧，因为我程序里认为白色部分为目标，所以图像的结果和上面有点反。
原图I Marker图J 结果图
可以看到结果很完美的体现了论文的需求。
二、清除边界部分的目标
在很多应用中，我们需要清除掉那些和边界连接在一起的目标，要实现这个功能，一个可行的方法是构建一副这样的Marker图像，图像中间部位全部填充为0（就是最小值），而周边区域则为原始图像的值，这样从就从边缘的部位开始向内重构，和边缘连接的目标就找到了，而未和边缘连接的对象则被去除掉，这不是和我们的目标相反了嘛，别急，我们在原图减去这个图不就可以了吗。
如下图所示（未找到有代表性的二值图，下图只是示意功能）。
原图 提取出的Marker图像 提取出的边缘图像 最终得到的结果图
当然这里还有一些问题，有些边界部分可能不需要被去除，那这个时候就需要通过其他的算法对这个结果再次进行补偿了。
这个对应了matlab里的imclearborder。对应的描述如下所示：
IM2 = imclearborder(IM) suppresses structures that are lighter than their surroundings and that are connected to the image border. (In other words, use this function to clear the image border.) I
简单的代码为；
int IM_ClearBorder(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Connectivity)
{
int Channel = Stride / Width;
if (Src == NULL) return IM_STATUS_NULLREFRENCE;
if ((Width <= 0) || (Height <= 0)) return IM_STATUS_INVALIDPARAMETER;
if (Channel != 1) return IM_STATUS_NOTSUPPORTED;
int Status = IM_STATUS_OK;
unsigned char *Temp = (unsigned char *)malloc(Height * Stride);
if (Temp == NULL) return IM_STATUS_NULLREFRENCE;
memset(Temp, 0, Height * Stride); // 先全部赋值为0
memcpy(Temp, Src, Width); // 边缘部位填充为原图的值
memcpy(Temp + (Height - 1) * Stride, Src + (Height - 1) * Stride, Width);
for (int Y = 1; Y < Height - 1; Y++)
{
Temp[Y * Stride] = Src[Y * Stride];
Temp[Y * Stride + Width - 1] = Src[Y * Stride + Width - 1];
}
Status = IM_ReConstruct(Src, Temp, Dest, Width, Height, Stride, Connectivity);
if (Status != IM_STATUS_OK) goto FreeMemory;
Status = IM_SubtractImage(Src, Dest, Dest, Width, Height, Stride);
if (Status != IM_STATUS_OK) goto FreeMemory;
FreeMemory:
free(Temp);
return Status;
}
那么针对灰度图像，这个还真有一些比较牛逼的功能，我们看看下面这个图的处理结果：
原图 提取出的Marker图像 提取出的边缘图像 最终得到的结果图
我们利用这个清除边界的算法成功的提取出了原图中的五个中心十字架的图像。
这个感觉有点像从边缘位置的像素（每个点都作为种子）向内部进行区域生长，这样中间这几个十字架因为基本被黑色的区域包围而没有和周边的圆形环接触，就被独立出来了。
三、填充孔洞
什么是孔洞，针对二值图像，我们的定义为：孔洞指的是那些不和边界连接在一起的最小局部区域（简易的就是黑色区域）。怎么又和边界扯上了呢，呵呵，就是这样。有了这个定义，那和我们的重构算法有什么关系呢。
其实啊，你想啊，如果我把原图反相后（白变黑，黑变白），这个时候我在同样以反相后的图像的边界图像为Marker图像，是不是就那些没有被边界连接起来的最大局部区域（最小的区域已经被反色）就被隔离了呢，这样我把结果再次反色后是不是就得到了想要的结果呢（针对下面的测试图，连通域需要选择为4）。
原图 反色图像 从反色图像提取出的边缘图像 基于反色后的边缘重构 再次反色
matlab对应的函数为：
BW2 = imfill(BW,'holes') fills holes in the binary image BW. A hole is a set of background pixels that cannot be reached by filling in the background from the edge of the image.
我们可以看下这个函数的部分M代码：
marker = mask;
idx = cell(1,ndims(I));
for k = 1:ndims(I)
idx{k} = 2:(size(marker,k) - 1);
end
marker(idx{:}) = Inf;
mask = imcomplement(mask);
marker = imcomplement(marker);
I2 = imreconstruct(marker, mask, conn);
I2 = imcomplement(I2);
I2 = I2(idx{:});
逻辑稍微有点不一样，他是先提取边缘图像后在反色边缘图像，结果是一样的。
当然，这个填充孔洞有个缺点，就是他是填充了所有的孔洞，而不可以运用一些其他的规则连限制孔洞的特性，比如孔洞的大小，圆度等等。这个就需要另外写函数了。
如果对于灰度图像，这个函数也有一些表有意思的结果：
原图 填充孔洞后的结果
可以看到，他把中间那些文字和一些比较黑的地方都去除了，也许这个结果某些场合下比较有用。
四、双阈值图像分割
有些图像比较复杂， 要从复杂的背景图像中分割出目标图像，单个阈值很多情况是难以做到的，如果存在这一种情况：即较小的阈值能分割出目标的主体部分，但是也会带入一些背景，但是背景和主体部分部想连，而较大的阈值侧能分割出目标的部分主体，但是基本没有啥背景图影响，那么这个时候就可以用较大阈值分割后的二值图作为Marker图像，较小阈值得到的那个图作为原图进行重构，这样就能得到较为满意的结果。
5、丢失目标的恢复
我们在进行目标查找的时候，经常进行各种预处理操作，比如开闭操作，顶帽变换等等，这些变换在达到了预处理的目的后，总会多多少少的改变了哪些正常的目标图，这往往不是不是我们想要的。这个时候我们就可以用重构算法来回复他们，比如下图：
原始图 某种处理后去掉了不需要的目标，但改变了正常目标原始形态 使用重建进行回复
6、区域最大值和最小值
这里的区域最大值和最小值不是我们立即的普通意义的最大值和最小值，其严格的定义应该是：
A regional minimum M of an image f at elevation t is a connected component of pixels with the value t whose external boundary pixels have a value strictly greater than t。
也就是这不是指的一个像素，而是一个连续的区域，这些区域具有相同的像素值t,并且其领域的像素都比他或者小。
论文里这样的描述： Image minima and maxima are important morphological features because they often mark relevant image objects: minima for dark objects ami maxima for bright objects. In morphology, the term minimum is used in the sense ofregional minimum, i.e. , a minimum whose extent is not necessarily restricted to a unique pixel.
至于这个东西在形态说如何重要，我还真的不了解。我在halcon的有关函数里也看到这样的算子，比如： local_max local_min
Local_min extracts all points from Image having a gray value smaller than the gray value of all its neighbors and returns them in LocalMinima.
matlab提供了imregionalmax、imregionalmin这样的算子来实现这个功能。
在实现上，其实就是以原图的值-1作为marker图像，对原图进行重构。、
这个算法主要针对的是灰度图像，对二值图没有什么意义。
原图 regional max halcon的local_max 算子（叠加显示）
那么还可以进一步扩展为extendedmax， extendedmin, Hmax, Hmin等算子，不过目前我对他们的应用了解的不多，有点不知道如何让他们产生价值。
7、图像分割封面的辅助功能
这个方面我也不太了解，有兴趣的作者可以去看看相关的论文。仅仅贴两张图予以展示。
8、其他
这算法还有其他的一些应用，比如matlab里面的imimposemin函数，还比如终极腐蚀点的定位等等，都可以，待将来有时间来再仔细的研究研究。
我已经将这个功能集成到我自己的DEMO中了，速度上划算行，因为这个算法不太好用SIMD指令集优化，只能纯C实现。
再次列出参考文献的名称供有需要的朋友了解：
1、Soille, P., Morphological Image Analysis: Principles and Applications, Springer-Verlag, 1999, pp. 172-173.
2、 Vincent, L., "Morphological Grayscale Reconstruction in Image Analysis: Applications and Efficient Algorithms," IEEE Transactions on Image Processing, Vol. 2, No. 2, April, 1993, pp. 176-201.
可执行的DEMO下载地址为：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar?t=1660121429
如果想时刻关注本人的最新文章，也可关注公众号：
