MUSICA（多尺度映像對比增強）演算法的簡要原理及VC實現

MUSICA的專利文檔：MUSICA_patent - Original document.pdf

    程式源碼下載：IPLab_MUSICA.rar

 

    演算法原理：
    映像增強一般方法是對比展開和長條圖均衡，這兩種方法對於灰階級過於集中， 還有大量的灰階級沒有充分利用的情況下比較適用，並且這兩種方法都是基於長條圖的變換，和像素的位置資訊無關。假設有原始訊號如(a)，可看出細節訊號（高頻，例如指紋、衣服褶縐這樣的細節）集中在較窄的灰階範圍內，人眼很難分辨。並且細節訊號集中的灰階範圍在整個灰階級空間內，使用對比增強和長條圖均衡無法對不同的地區內像素進行不同程度的放大或者縮小。因此就有了MUSICA（Multi-Scale Image Contrast Amplification），可以翻譯成多尺度映像對比增強。
    基本的原理是先提取細節資訊，然後對細節訊號進行增強放大，然後再重構到原圖中。局部細節資訊的提取現在研究的熱點一直都是小波分解。
    暫且不管理論，先從實際運用的角度來看。（b），假設對訊號進行平滑濾波，可以得到新的藍色訊號。可以認為藍色訊號是對原訊號的低解析度近似，保留了整體的特徵而丟棄了高頻的細節訊號，例如要分析映像的整體特徵就要使用近似訊號以排除高頻（可能是雜訊）的幹擾。
    那麼（c），原訊號減去近似訊號，得到的就是高頻訊號，也就是細節訊號（當然也可能是雜訊）。
    然後對（c）進行增強，最簡單例如放大2倍，或者進行log變換得到（d）。
    然後把（d）加回到近似訊號上，重構出增強後的原訊號，如（e）所示。可以看到，映像的整體特性沒有改變，兩個細節集團的細節訊號被放大了，適合人眼辨認。

圖1

    這個算是從最直觀的角度來考慮演算法原理，只是一種近似的抽象，實際過程需要更多更具體的考慮，同時這也是個人的理解，會有偏頗的地方，留待完善。

    然後，總得給出點具體的原理和演算法的實際步驟。詳細的說明可以從MUSICA的專利文檔裡獲得，這裡給出個人角度的理解和說明。
    首先要從小波說起，直觀地看，小波分解將（w*h）的映像分解成4個（w/2 * h/2） 的子圖，其中1個為近似映像，3個為小波係數矩陣。然後對近似映像進行下一級的分解。如，A為原圖的低解析度近似映像，D為小波係數，D1、D2、D3的物理意義可以看作是水平、垂直、以及對角方向的細節矩陣。小波的重構是這個分解過程的逆過程，使用4個低解析度的子圖重組成一個高解析度的圖。這部分知識最好看看小波的教科書，可以比較全面瞭解小波的特點。

圖2

    直觀地看小波分解，近似映像就是對原圖的一個低解析度的近似，保留映像的整體特徵，而小波係數矩陣，就是低解析度映像和原圖的之間的差異，就是高頻的訊號，也就是細節。通過這樣的分解，可以得到各個解析度上的細節，目的是使某種解析度下難以區分的特徵在另一種解析度下將很容易被區分和檢測。人的眼睛可以自動聚焦來選擇最佳的解析度對眼前的景物進行分辨，而小波也就是因為類似人眼而受到很多重視。還有一個重要的概念是小波分解的小波，可以理解為用來分解映像而使用的小波基。直觀地看，不同特性的小波基分解同一個映像，得到近似映像應該差不多，而提取得到的細節矩陣則反映的映像細節則有所差異。關於小波基的選擇以及各種小波基特性的物理含義一直沒有尋找到比較好的解釋。有一種解釋是，某種小波基對某類特徵的映像得到的小波係數矩陣比較有代表性，例如xxx1小波對圖形組合簡單的映像效果比較好，而xxx2小波對噪雜的自然映像好，xxx3對條紋狀映像好...這僅是一個同學間討論的假設，有待查證。
    但是MUSICA裡使用的不是這種分解，而是使用高斯影像金字塔，在Rafael C. Gonzalez.的 Digital Image Processing.（電子工業出版社. 2004）一書裡有介紹，這裡引用其中的兩個圖來說明，前面提到的下波分解作為數學的基礎，教材比較多，容易深入。
    高斯影像金字塔也是小波分解的一種，它產生一個w/2*h/2的低解析度近似，和一個與原圖等大的殘差距陣（也就是高頻細節）。圖中的抽樣（2↓）指對w*h的映像進行隔行採樣，得到w/2*h/2的子圖。插入器（2↑）則將w/2*h/2的映像行列間插入一個值擴充成w*h的映像，插入的值可以是0，或者是左右行的均值等方式。近似濾波和插入濾波就是所選擇的小波基轉化的濾波器組形式，這裡選擇的是高斯濾波器，濾波器的參數在專利文檔裡提供了一組經驗參數，但是到具體裝置上還需要根據實際成像的雜訊等情況重新計算過才能達到較優的結果。

圖3

    專利文檔提供得整體演算法框圖如下：

圖4

    2為輸入映像，30為j級影像金字塔的分解過程，得到31為所有的殘差矩陣，31'為映像的第j級近似。61緩衝映像（對具體裝置有意義吧）和62並對31進行s型函數的查表變換。34為重構過程，將變換後的殘差距陣33和j級近似31'重新組合成增強後的映像4，輸出。
    根據這個流程，就可以開始考慮演算法的程式實現了。
    在開始在pc機上編程之前，需要說明幾點：
    首先是演算法原理中影響效果的幾個部分。一是62裡s型函數的選擇。增強演算法都有會帶來雜訊也一同增強副作用，專利文檔裡討論了這個問題，選擇y=m*(x/m)p的指數變換，p的取值會對增強頻段產生影響，合理選擇p、進行分段選擇不同的p，可以抑制雜訊的增強。二是30分解過程和34重構過程，這個其實屬於小波分解的部分，分解的級數、分解步驟和小波基的選擇（這裡就是濾波器參數的選擇），會對分解得到殘差距陣產生影響，得到的殘差距陣是否具有有效細節的代表性，能否有效過濾雜訊，是這裡的關鍵。編程中可以將一些參數設定成變數，使用配製檔案輸入，則可以將演算法的研究和實現分開。專利文檔中對這兩個方面也進行了討論，可以參考深入。
    然後還需要說明的一點的是，在pc上實現這個演算法可能並不適合，一般都是在醫療設備上直接實現這個演算法（嵌入式系統或者是某種特殊韌體），使用關鍵字“AGFA MUSICA”可以在google上找到些相關資料。但是在這裡並不考慮這點，僅把MUSICA作為一道有趣的演算法題目，結合一些影像處理的知識努力去得到最優解。有相關背景的朋友也希望可以介紹一下相關專業的知識，那樣就更有趣了。

    程式源碼下載：IPLab_MUSICA.rar

    整體的流程還是比較簡單的，問題是影像處理中的一些細節，例如濾波時邊界的處理，映像抽樣時w或者h為奇數時w/2（h/2）的取值，也就是離散整數的取整誤差問題。還有些編程時記憶體配置的問題，臨時緩衝的分配等。但是不考慮得太複雜都使用最常用的方式，以上僅將大概的問題列出，碰到時候需要注意。
    入口函數虛擬碼：
void GDIPlusImage::IPFuncMUSICA()
{
    O = 原始映像;

    for（分解級數）
      {計算每級的長寬儲存到矩陣S;}

    ImageGrey24To8(O);//灰階化，僅對灰階圖象進行處理

    X = MUSICA_Decomposition(O,S);//分解，得到X為所有的映像按S裡的順序排列，30

    X = MUSICA_Mapping(X,S);//對分解映像進行變換，32

    O = MUSICA_Reconstruction(X,S);//重構，得到增強後的映像，34

    ImageGrey8To24(O);//轉化成24為圖，顯示、後繼處理起來方便，僅此而已
}

    MUSICA_Mapping(X,S)函數相對簡單，只需要對S中的殘差距陣進行灰階級的查表替換，替換的表根據s型函數事先構造好，其過程類似進行反色處理，可以參考前面的文章。可以考慮將替換表的構造獨立出來，通過參數傳遞，方便改進。

    MUSICA_Decomposition分解和MUSICA_Reconstruction重構過程相對複雜，而且專利文檔中提供了幾種改進後的分解重構過程。對於編程實現來說其實本質都相通，因此所附程式中實現了專利文檔中的4.a和4.b，也就是最簡單的一種。僅對4.a分解的實現進行說明，4.c重構為逆過程，參考原始碼。

圖5

    一次分解流程：
1) 對X進行一級分解，X拷貝到buffer
2) 對buffer進行平滑濾波
3) buffer隔行抽樣得到下一級的近似圖X+1
4) 近似圖X+1隔行插值0到buffer
5) 進行濾波參數乘4（對置0行列的補償）的平滑濾波
6) X減去平滑後的buffer得到殘差映像儲存於X

    這樣就完成一次分解，然後再對X+1進行下一級的分解。

圖6

    所附程式在xp+vs.net2003下編譯運行通過，程式運行後開啟工程目錄下的圖片，點擊快捷按鈕1或者功能表項目“影像處理-〉MUSICA”，然後要稍等一會，處理速度比較慢。工程目錄下有兩張樣本圖片，“心臟-交響樂處理前.jpg”為原始成像圖片，“心臟-交響樂處理後.jpg”是一位老師給我的在專業裝置上處理過後的結果圖，作為和程式運行結果的比較。當然，程式啟動並執行結果比它差很多，但是基本方向是對的，專利文檔裡還有很多改進的步驟和參數的擇優都未實現，有興趣的朋友自己動手。：）