【瀝血整理】灰度（二值）影象重構演演算法及其應用（morphological reconstruction）。
不記得是怎麼接觸並最終研究這個課題的了，認識我的人都知道我是沒有固定的研究物件的，一切看運氣和當時的興趣。本來研究完了就放在那裡了，一直比較懶的去做總結，但是想一想似乎在網路上就沒有看到關於這個方面的資料，能搜尋到的都是一些關於matlab相關函數的應用，決定還是抽空趁自己對這個演演算法還有點記憶的時候寫點東西吧，畢竟這個演演算法還有一些應用是值得回味和研究的。而且也具有一定的工程價值。
怎麼說呢，其實在很早瀏覽matlab的影象處理工具箱的時候，就無數次的看到過這些函數，但是無奈當時不知道他們有什麼用，就沒怎麼鳥他， 其實M還是很重視他們，這個從他們在工具箱裡佔用的函數列表篇幅裡就能完美的看的出：
這個定義簡單翻譯就是從標記影象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、基本功能
接下來我們來了解下這個演演算法的直接應用，我們先來驗證下前面舉例的那個圖片吧，因為我程式裡認為白色部分為目標，所以影象的結果和上面有點反。
這個演演算法主要針對的是灰度影象，對二值圖沒有什麼意義。
8、其他
這演演算法還有其他的一些應用，比如matlab裡面的imimposemin函數，還比如終極腐蝕點的定位等等，都可以，待將來有時間來再仔細的研究研究。
我已經將這個功能整合到我自己的DEMO中了，速度上划算行，因為這個演演算法不太好用SIMD指令集優化，只能純C實現。
