01

前言

前文我們介紹過demons算法，以及在demons算法的基礎上發展的Active demons算法和Inertial demons算法：

圖像配準算法之demons算法

我們知道，demons算法是一種基于光流理論發展起來的配準方法，因此該算法與其它光流算法一樣，具有小運動的約束條件。如果參考圖像和浮動圖像的差别很大，也即不滿足小運動條件，那麼demons算法的配準效果會很差。

為了解決小運動問題，Vercauteren T等人提出了微分同胚算法[1]，将demons算法計算得到的位移場轉換為微分同胚映射，從而對大運動或大形變都具有較理想的配準效果。本文介紹一種基于Active demons和微分同胚映射的大形變配準算法，其中使用到的微分同胚映射轉換算法參考了Vercauteren T等人的文章。

首先我們來複習一下Active demons算法計算每個像素點位移的公式，其中其中Sx和Sy分别是參考圖像上點(x,y)處x方向和y方向的梯度，Mx和My分别是浮動圖像上點(x,y)處x方向和y方向的梯度，▽f則是點(x,y)處參考圖像與浮動圖像的灰度差：

本文所介紹的微分同胚demons算法，則是在以上Active demons計算公式計算得到的Ux、Uy基礎上，再将Ux、Uy轉換為微分同胚映射。

02

什麼是微分同胚映射？

微分同胚映射，它包含了“同胚”和“微分”兩個條件，也就是說必須同時滿足這兩個條件的映射，才是微分同胚映射。下面我們分别介紹這兩個條件。

• 同胚

介紹同胚映射之前，我們先介紹一下數學中“單射”、“滿射”、“雙射”的概念。假設有函數y=f(x)，函數f也可以稱為從x到y的映射：

1. 如果每個x值都有唯一的y值與之對應，那麼映射f為單射。
2. 如果每個y值都至少有一個x值與之對應，那麼映射f為滿射。
3. 如果映射f同時滿足單射和滿射條件，那麼又稱之為雙射，直觀說就是當所有x值與所有y值一一對應的時候，映射f為雙射。

單射、滿射、雙射的示意圖如下：

其實同胚映射跟雙射是一樣的，都是一一對應的關系，隻是同胚映射通常用于形容拓撲空間的映射關系而已，比如對于兩張圖像，如果它們的所有像素點都滿足一一對應關系，那麼這個對應關系就是同胚映射，如下圖所示：

同胚映射

• 微分

如果映射f滿足以上所說的微分條件，意味着映射f與逆映射f-1都是光滑的映射。什麼又是光滑映射呢？每個x值對應到y值的映射關系，與其鄰域内其它x值對應到y值的映射關系是相似的，即使不完全一樣，也是連續緩慢的變化，而不是突變，那麼稱x到y的映射為光滑映射。如下圖所示，左圖中x到y的映射是連續緩慢變化的，所以是光滑映射，而右圖中x4到y4的映射，大小和方向都突變了，所以右圖的映射不是光滑映射。

此處講的光滑映射可能有點抽象，下文我們再舉例更加直觀地說明。

03

相似度衡量

前文我們使用Active demons算法類配準的時候，每輪叠代都更新了最優位移場。實際上，并不是每輪叠代都是朝着配準效果更好的方向前進的，有的叠代配準效果反而更差，所以每輪叠代都更新最優位移場是有問題的。于是使用相似度來判斷是否要更新最優位移場的方法被提了出來，流程如下圖所示，其中F、M、M1依次為參考圖像、浮動圖像、像素重采樣之後的浮動圖像。

本文我們使用歸一化互相關系數來衡量F、M1的相似度，假設F、M1都是M行N列矩陣，歸一化互相關系數的計算公式如下：

NCC的取值範圍0~1，NCC越接近1表示兩張圖像越相似。

考慮到在叠代過程中，圖像M在像素重采樣之後有可能會出現黑色邊緣（例如下圖），這些黑色邊緣會影響相似度的準确性，所以計算NCC時應該把這些黑色邊緣區域排除在外

我們首先根據計算的位移場獲取黑色邊緣的mask掩碼矩陣，mask矩陣值為0的像素點屬于黑色邊緣，其它則不屬于。假設每輪叠代計算得到的位移場為M行N列的矩陣Ux、Uy，對于每個像素點(x,y)，判斷其像素值是否滿足以下條件（round為四舍五入運算），如果滿足則判定點(x,y)不屬于黑色邊緣，并對mask(x,y)賦值255，反之則屬于黑色邊緣并對mask(x,y)賦值0。

基于以上得到的mask矩陣，NCC計算公式變成下式，也即隻有非黑色邊緣的點才加入計算，黑色邊緣點排除在外：

04

梯度圖的計算

計算位移場的時候，需要使用到參考圖像F、浮動圖像M的梯度。本文我們使用差分梯度來計算位移場。

對于圖像上任意一點(x,y)，其像素值為f(x,y)，那麼該點的x、y方向的梯度gx(x,y)、gy(x,y)可按下式計算：

按理說圖像的邊緣點沒法按照上式計算梯度，所以我們需要對圖像做邊緣填充：在上、下各填充一行，左、右各填充一列。

例如，按照上式計算Lena圖像所有像素點的x、y方向梯度，得到x、y方向的梯度圖如下：

05

将位移場轉換為微分同胚映射

這一步為微分同胚demons算法的核心。參考Vercauteren T等人的文章，可使用複合運算來實現近似微分同胚映射轉換，那麼什麼又是複合運算呢？下面我們詳細道來。

假設x、y方向的位移場分别為Ux、Uy，使用Ux、Uy分别對它們自身進行像素重采樣的操作得到Ux'和Uy'，然後再計算Ux Ux'和Uy Uy'的運算就是複合運算。對于Ux'的任意像素點(x',y')，其對應Ux上的點(x,y)的坐标可按下式計算：

得到(x,y)之後，再把點(x,y)的像素值Ux(x,y)賦值給Ux'上點(x',y')，即完成點(x,y)的像素重采樣操作。

不過由于x、y都是浮點數，無法直接得到Ux(x,y)的值，所以本文使用雙線性插值算法計算Ux(x,y)。關于插值算法與像素重采樣的詳細介紹，可參考前文：

常見圖像插值算法的原理與C 實現

對Ux'的所有像素點都按照上述進行像素插值，則完成了從Ux到Ux'的像素重采樣，從Uy到Uy'的像素重采樣也同理：

最後再計算兩者之和，即完成複合運算：

以上是複合運算的介紹，接下來就是将位移場轉換為近似微分同胚映射的流程：

（1）按下式計算位移場Ux、Uy的平方和矩陣，其中“.*”為點乘運算，也即兩個矩陣中對應位置的像素值相乘。

（2）求矩陣norm2的最大值maxv。

（3）按下式計算n，其中ceil為向上取整運算，比如ceil(1.1)=2.0，max為取較大值運算。

（4）計算矩陣Ux1、Uy1：

（5）對矩陣Ux1、Uy1作n次複合運算，最後結果就是近似的微分同胚映射。

06

微分同胚demons算法的代碼實現

（1）計算梯度代碼：

void get_gradient(Mat src, Mat &Fx, Mat &Fy) { Mat src_board; //邊緣擴充 copyMakeBORDER(src, src_board, 1, 1, 1, 1, BORDER_REPLICATE); Fx = Mat::zeros(src.size(), CV_32FC1); Fy = Mat::zeros(src.size(), CV_32FC1); for (int i = 0; i < src.rows; i ) { float *p_Fx = Fx.ptr<float>(i); float *p_Fy = Fy.ptr<float>(i); for (int j = 0; j < src.cols; j ) { p_Fx[j] = (src_board.ptr<float>(i 1)[j 2] - src_board.ptr<float>(i 1)[j]) / 2.0 ; p_Fy[j] = (src_board.ptr<float>(i 2)[j 1] - src_board.ptr<float>(i)[j 1]) / 2.0; } } }

（2）複合運算實現代碼如下，其中調用Opencv的remap函數可以方便地實現像素重采樣。

//像素重采樣實現 void movepixels_2d2(Mat src, Mat &dst, Mat Tx, Mat Ty, int interpolation) { Mat Tx_map(src.size(), CV_32FC1, 0.0); Mat Ty_map(src.size(), CV_32FC1, 0.0); for (int i = 0; i < src.rows; i ) { for (int j = 0; j < src.cols; j ) { Tx_map.at<float>(i, j) = j Tx.at<float>(i, j); Ty_map.at<float>(i, j) = i Ty.at<float>(i, j); } } remap(src, dst, Tx_map, Ty_map, interpolation); } //複合運算實現 void composite(Mat& vx, Mat& vy) { Mat bxp, byp; movepixels_2d2(vx, bxp, vx, vy, INTER_LINEAR); movepixels_2d2(vy, byp, vx, vy, INTER_LINEAR); vx = vx bxp; vy = vy byp; }

（3）微分同胚映射轉換實現：

void exp_field(Mat &vx, Mat &vy, Mat &vx_out, Mat &vy_out) { Mat normv2 = vx.mul(vx) vy.mul(vy); double minv, maxv; Point pt_min, pt_max; minMaxLoc(normv2, &minv, &maxv, &pt_min, &pt_max); //求最大最小值 float m = sqrt(maxv); float n = ceil(log2(m / 0.5)); n = n > 0.0 ? n : 0.0; float a = pow(2.0, -n); vx_out = vx * a; vy_out = vy * a; //n次複合運算 for (int i = 0; i < (int)n; i ) { composite(vx_out, vy_out); } }

（4）高斯濾波函數實現，調用Opencv的函數GaussianBlur可方便實現高斯濾波：

void imgaussian(Mat src, Mat& dst, float sigma) { int radius = (int)ceil(3 * sigma); int ksize = 2 * radius 1; GaussianBlur(src, dst, Size(ksize, ksize), sigma); }

（5）計算Active demons位移場的代碼實現：

void demons_update(Mat S, Mat M, Mat Sx, Mat Sy, float alpha, Mat& Tx, Mat& Ty) { Mat diff = S - M; Tx = Mat::zeros(S.size(), CV_32FC1); Ty = Mat::zeros(S.size(), CV_32FC1); Mat Mx, My; get_gradient(M, Mx, My); //求M的梯度 float alpha_2 = alpha * alpha; for (int i = 0; i < S.rows; i ) { float* p_sx = Sx.ptr<float>(i); float* p_sy = Sy.ptr<float>(i); float* p_mx = Mx.ptr<float>(i); float* p_my = My.ptr<float>(i); float* p_tx = Tx.ptr<float>(i); float* p_ty = Ty.ptr<float>(i); float* p_diff = diff.ptr<float>(i); for (int j = 0; j < S.cols; j ) { float alpha_diff = alpha_2 * p_diff[j] * p_diff[j]; float a1 = p_sx[j] * p_sx[j] p_sy[j] * p_sy[j] alpha_diff; //分母 float a2 = p_mx[j] * p_mx[j] p_my[j] * p_my[j] alpha_diff; if (a1 > 0.00001 && a2 > 0.00001) { p_tx[j] = p_diff[j] * (p_sx[j] / a1 p_mx[j] / a2); p_ty[j] = p_diff[j] * (p_sy[j] / a1 p_my[j] / a2); } } } }

（6）計算相似度代碼實現：

//獲取黑色邊緣的mask矩陣 void cal_mask(Mat Tx, Mat Ty, Mat &mask) { mask.create(Tx.size(), CV_8UC1); for (int i = 0; i < Tx.rows; i ) { float *pTx = Tx.ptr<float>(i); float *pTy = Ty.ptr<float>(i); uchar *pm = mask.ptr<uchar>(i); for (int j = 0; j < Tx.cols; j ) { int x = (int)(j pTx[j] 0.5); int y = (int)(i pTy[j] 0.5); if (x < 0 || x >= Tx.cols || y < 0 || y >= Tx.rows) { pm[j] = 0; } else { pm[j] = 255; } } } } //計算歸一化互相關系數 double cal_cc_mask(Mat S1, Mat Si, Mat Mask) { double result; double sum1 = 0.0; double sum2 = 0.0; double sum3 = 0.0; for (int i = 0; i < S1.rows; i ) { for (int j = 0; j < S1.cols; j ) { if (Mask.ptr<uchar>(i)[j]) { sum1 = (double)(S1.at<uchar>(i, j)*Si.at<uchar>(i, j)); sum2 = (double)(S1.at<uchar>(i, j)*S1.at<uchar>(i, j)); sum3 = (double)(Si.at<uchar>(i, j)*Si.at<uchar>(i, j)); } } } result = sum1 / sqrt(sum2*sum3); //範圍0~1 return result; }

（7）微分同胚demons算法的最終實現：

void register_diffeomorphic_demons(Mat F0, Mat M0, float alpha, float sigma_fluid, float sigma_diffusion, int niter, Mat &Mp, Mat &sx, Mat &sy) { Mat F, M; F0.convertTo(F, CV_32F); //将參考圖像和浮動圖像都轉換為浮點型矩陣 M0.convertTo(M, CV_32F); //初始化位移場為0位移 Mat vx = Mat::zeros(F.size(), CV_32FC1); Mat vy = Mat::zeros(F.size(), CV_32FC1); float e_min = -9999999999.9; Mat sx_min, sy_min; Mat Sx, Sy; get_gradient(F, Sx, Sy); //計算參考圖像梯度圖 Mat M1 = M.clone(); for (int i = 0; i < niter; i ) { Mat ux, uy; //計算Active demons的位移場 demons_update(F, M1, Sx, Sy, alpha, ux, uy); imgaussian(ux, ux, sigma_fluid); //高斯濾波 imgaussian(uy, uy, sigma_fluid); //高斯濾波 vx = vx ux; //将位移場累加 vy = vy uy; imgaussian(vx, vx, sigma_diffusion); //高斯濾波 imgaussian(vy, vy, sigma_diffusion); //高斯濾波 exp_field(vx, vy, sx, sy); //将累加的位移場轉換為微分同胚映射 Mat mask; cal_mask(sx, sy, mask); //計算黑色邊緣的mask掩碼矩陣 movepixels_2d2(M, M1, sx, sy, INTER_LINEAR); //對M像素重采樣 float e = cal_cc_mask(F, M1, mask); //計算F、M1的相似度 if (e > e_min) //如果相似度提高，則更新最佳位移場 { printf("i=%d, e=%f\n", i, e); e_min = e; sx_min = sx.clone(); sy_min = sy.clone(); } } sx = sx_min.clone(); //得到最優微分同胚映射 sy = sy_min.clone(); movepixels_2d2(M, Mp, sx, sy, INTER_LINEAR); Mp.convertTo(Mp, CV_8U); }

07

測試結果

為了驗證我們實現的微分同胚算法對大形變圖像的配準效果，我們分别使用本文實現的微分同胚demons算法、Active demons算法對以下大形變的心髒圖像進行配準，并比較結果：

由前文可知，Active demons算法中的alpha參數越大，配準步長越小，反之則配準步長越大。由于是大形變圖像，我們統一給alpha取一個較小值0.15，高斯濾波參數sigma_fluid和sigma_diffusion都取0.05，叠代次數設置為3000。配準結果如下：

記錄每輪叠代的配準圖像與參考圖像的相似度，如下圖：

由以上結果可知，微分同胚demons算法的配準結果比Active demons算法的好得多。Active demons算法不适用于大位移，那麼我們嘗試把alpha調大為1.15，使位移減小，得到Active demons算法的以下結果，發現還是不理想。

上文我們講到，微分同胚映射是一種光滑的映射，什麼是光滑映射呢？現在我們分别計算微分同胚demons算法、Active demons算法的最終位移場的位移幅度圖（也即計算每個像素點的位移距離）并作僞彩色處理，同時畫出每一個像素點的位移向量（白色線為位移軌迹，紅點為位移終點），如下圖所示，可以形象地看出來微分同胚映射的位移場是光滑的，相比Active demons算法的位移場就粗糙多了。

本文參考文獻：

[1] Vercauteren T , Pennec X , Perchant A , et al. Diffeomorphic demons: efficient non-parametric image registration.[J]. NeuroImage, 2009, 45( 1):S61-S72.

更多精彩资讯请关注tft每日頭條，我们将持续为您更新最新资讯!