不記得是怎么接觸并最終研究這個課題的了，認識我的人都知道我是沒有固定的研究物件的，一切看運氣和當時的興趣，本來研究完了就放在那里了，一直比較懶的去做總結，但是想一想似乎在網路上就沒有看到關于這個方面的資料，能搜索到的都是一些關于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 =https://www.cnblogs.com/Imageshop/archive/2022/08/11/= 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 =https://www.cnblogs.com/Imageshop/archive/2022/08/11/ 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 =https://www.cnblogs.com/Imageshop/archive/2022/08/11/= 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
如果想時刻關注本人的最新文章，也可關注公眾號：
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/501627.html
標籤：其他
上一篇：SOAP 和 REST 的區別