本次水平集圖像分割并行加速算法設計與實現包含：原理篇、串行實現篇、OpenMP并行實現篇與CUDA GPU并行實現篇四個部分。具體各篇章鏈接如下：

水平集圖像分割并行加速算法設計與實現——原理篇

水平集圖像分割并行加速算法設計與實現——串行實現篇

水平集圖像分割并行加速算法設計與實現——OpenMP并行實現篇

水平集圖像分割并行加速算法設計與實現——CUDA GPU并行實現篇

原理篇主要講解水平集圖像分割的原理與背景。串行實現篇、OpenMP并行實現篇與CUDA GPU并行實現篇主要基于C++與OpenCV實現相應的圖像分割與并行加速任務。本系列屬于圖像處理與并行程序設計結合類文章，希望對你有幫助😊。

串行基本流程設計

水平集圖像分割串行代碼設計流程如上圖所示。演化開始前需要首先對原始圖片進行導入，并對原始圖片相關變量進行初始化，其后根據選定初始化輪廓區域結合符號距離函數進行水平集的初始化。在上述完成后，水平集演化過程開始，根據相應數據依賴關系，首先需要完成對應Heaviside函數、Dirac函數與曲率的計算，其后根據Heaviside函數計算結果對前景與背景均值進行計算。
接著便可將上述計算結果代入式(1)進行演化模擬。具體水平集圖像分割與演化原理，見本人之前發布的原理篇所示。
$$\begin{array}{c}{?}_{i,j}^{n+1}=L\left({?}_{i,j}^n\right)\Delta t+{?}_{i,j}^n\end{array}\text{\hspace{1em}(1)}$$
最后對水平集的收斂性進行判斷，當零水平集所對應輪廓收斂時，則程序結束退出，反之，繼續迭代進行演化。

主要模塊編程實現設計

原始圖片初始化模塊

原始圖片初始化模塊主要用于對原始圖片相關變量進行初始化。其會將原始圖片轉化為單通道的灰度圖像，并根據原始圖像尺寸對水平集矩陣、Heaviside函數矩陣等相關數據進行初始化，具體代碼如下所示：

LevelSet::LevelSet(const Mat& src) {if (src.channels() == 3){cvtColor(src, src_, COLOR_BGR2GRAY);src.copyTo(image_);}else{src.copyTo(src_);cvtColor(src, image_, COLOR_GRAY2BGR);}src_.convertTo(src_, CV_32FC1);phi_ = Mat::zeros(src_.size(), CV_32FC1);dirac_ = Mat::zeros(src_.size(), CV_32FC1);heaviside_ = Mat::zeros(src_.size(), CV_32FC1); }

水平集初始化模塊

根據CV水平集模型，采用倒置的符號距離函數進行水平集的初始化，具體倒置符號距離函數公式如式(2)所示。
$$\begin{array}{c}?(\mathrm{x})=?\begin{array}{cc}?\mathrm{dist}(\mathrm{x},C)& \text{?if?}\mathrm{x}\text{?is?outside?}C\\ 0& \mathrm{x}\in C\\ \mathrm{dist}(\mathrm{x},C)& \text{?if?}\mathrm{x}\text{?is?inside?}C\end{array}\end{array}\text{\hspace{1em}(2)}$$
其中dist(x,C)表示點x到輪廓C上所有點的最短距離，具體轉化示例見下圖所示。

根據上圖可知，其將二維平面的輪廓曲線作為零水平集（水平集曲面與z=0平面的交面），轉化為一個三維水平集曲面，通過三維曲面的演化來模擬二維的輪廓演化。（當迭代收斂后直接獲取對應水平集曲面的零水平集，即為分割完成后區域的輪廓）
具體水平集初始化代碼如下所示：

// 初始化水平集 void LevelSet::initializePhi(Point2f center, float radius) {const float c = 2.0f;float value = 0.0;for (int i = 0; i < src_.rows; i++){for (int j = 0; j < src_.cols; j++){value = -sqrt(pow((j - center.x), 2) + pow((i - center.y), 2)) + radius;if (abs(value) < 1e-3){//在零水平集曲線上phi_.at<float>(i, j) = 0;}else{// 在零水平集內：為正// 在零水平集外：為負phi_.at<float>(i, j) = value;}}} }

曲率模塊

曲率模塊主要用于對曲率進行計算。其中曲面的曲率可以運用梯度除以其模值的散度進行計算，具體計算公式為式(3)所示：
$$\begin{array}{c}{\text{?curvature?}}_{ij}=\mathrm{div}\left(\frac{??}{|??|}\right)\end{array}\text{\hspace{1em}(3)}$$
其具體實現代碼如下所示：

// 計算曲率 void LevelSet::calculateCurvature() {Mat dx, dy;// 計算一階梯度gradient(src_, dx, dy);Mat norm = Mat::zeros(src_.size(), CV_32FC1);// 計算模值for (int i = 0; i < src_.rows; i++){for (int j = 0; j < src_.cols; j++){norm.at<float>(i, j) = pow(dx.at<float>(i, j) * dx.at<float>(i, j) + dy.at<float>(i, j) * dy.at<float>(i, j), 0.5);}}Mat dxx, dxy, dyx, dyy;// 二階梯度與散度計算gradient(dx / norm, dxx, dxy);gradient(dy / norm, dyx, dyy);curv_ = dxx + dyy; }

上述gradient函數用于計算對應矩陣x方向與y方向的梯度，具體可參見本人github倉庫源碼。

前景背景均值模塊

前景背景均值模塊用于對輪廓所劃分的前景與背景均值進行計算，其將調用Heaviside函數計算結果對前景和背景范圍進行劃分，具體代碼如下所示：

// 計算像素點前景與背景均值 void LevelSet::calculatC() {c1_ = 0.0f;c2_ = 0.0f;float sum1 = 0.0f;float h1 = 0.0f;float sum2 = 0.0f;float h2 = 0.0f;float value = 0.0f;float h = 0.0f;for (int i = 0; i < src_.rows; i++){for (int j = 0; j < src_.cols; j++){value = src_.at<float>(i, j);h = heaviside_.at<float>(i, j);h1 += h;sum1 += h * value;h2 += (1 - h);sum2 += (1 - h) * value;}}c1_ = forntpro_ * sum1 / (h1 + 1e-10);c2_ = sum2 / (h2 + 1e-10); }

同時為了更好的對圖像進行分割，為前景均值添加權重項forntpro_，以此可根據不同情況對分割結果進行優化。

演化模擬與收斂性判定模塊

在上述模塊運行完畢后，可對水平集進行迭代演化。當演化后輪廓所包含前景區域與演化前輪廓所包含前景區域相同時(水平集函數大于0的區域沒有變化)，則判定水平集已經收斂，可以停止迭代跳出循環，具體代碼如下所示：

// 具體演化函數 // 運用迭代法解偏微分方程，求解φ對應時間間隔的變化量 int LevelSet::evolving() {showEvolving();showLevelsetEvolving();// 迭代次數int k;bool flag; //是否收斂判定標識for (k = 0; k < iterationnum_; k++){clock_t begin, end; // 計時變量begin = clock();flag = false;heaviside();dirac();calculatC();calculateCurvature();//update phi// 模擬演化過程for (int i = 0; i < src_.rows; i++){for (int j = 0; j < src_.cols; j++){float curv = curv_.at<float>(i, j);float dirac = dirac_.at<float>(i, j);float u0 = src_.at<float>(i, j);// 能量函數中各項計算float lengthTerm = mu_ * dirac * curv;float areamterm = nu_ * dirac;float fittingterm = dirac * (-lambda1_ * pow(u0 - c1_, 2) + lambda2_ * pow(u0 - c2_, 2));float term = lengthTerm + areamterm + fittingterm;float phinew = phi_.at<float>(i, j) + timestep_ * term;float phiold = phi_.at<float>(i, j);phi_.at<float>(i, j) = phinew;if (!flag){if (phinew * phiold < 0){flag = true;}}}}showEvolving();showLevelsetEvolving();// 對是否收斂進行判斷if (!flag){break;}end = clock();}showEvolving();showLevelsetEvolving();return k; }

其中showEvolving()函數用于繪制演化過程中的零水平集輪廓，showLevelsetEvolving()用于繪制水平集函數對應的色溫圖。
上述僅僅展現部分主要代碼，詳細代碼見本人github倉庫所示，具體鏈接如下 。
水平集圖像分割算法串行與并行代碼以及相關測試用例

測試用例說明

本次串行水平集程序的測試用例主要來源于實際生活中存在，或虛擬應用場景中存在的圖形圖像。

實際存在圖片用例

為更加精準測試算法性能，針對實際存在實物方面，選取人臉嘴唇圖片作為測試用例，其尺寸為883像素*594像素，具體如下圖所示。

在上圖中，包含皮膚斑紋、嘴唇反光等噪聲，并含有不規則邊緣，利于對算法的圖像分割效果進行測試。其水平集的初始化區域選擇為一個以點(400,240)為圓心，半徑為200的圓，具體如下圖所示。

虛擬業務場景圖片用例

針對虛擬業務場景方面，本項目選取兩張飛機游戲圖片作為測試用例，其中飛機游戲圖片一尺寸為320像素*570像素，飛機游戲圖片二尺寸為476像素*752像素，具體如下面兩圖所示。

兩幅圖像具有不同背景與圖像特征，利于測試算法針對不同分割對象的通用性。同時其中均含有淺色背景噪聲（星球、隕石等），有利于對算法的分割效果進行更加精準的評估。
其均選擇近乎全圖作為水平集初始化區域。同時為進一步對模型通用進行測試，在飛機游戲圖片一中，初始化區域不包含圖片邊緣四角，而在飛機大戰圖片二中，其包含圖片邊緣中的三個角，具體如下面兩圖所示。

運行環境說明

硬件運行環境

CPU：Intel? Core? i5-9300H 4核CPU處理器

GPU：NVIDIA GeForce GTX 1050顯卡

軟件運行環境

OpenCV: v4.5.0

CUDA: v11.0

運行結果展示

針對人臉嘴唇圖片，以時間步長為1進行測試，最大迭代次數為1000，得到的測試過程與結果展示如下圖所示。

針對飛機游戲圖片一，以時間步長為5進行測試，最大迭代次數為1000，得到的測試過程與結果展示如下所示。

針對飛機游戲圖片二，同樣以時間步長為5進行測試，最大迭代次數為1000，得到的測試過程與結果展示如下所示。

（上圖由于gif大小限制，僅僅展示演化過程開頭部分，完整過程可根據所提供代碼自行運行得到）
其左側圖為對應的水平集演化色溫圖，深藍色表示小于零的背景，其他顏色表示大于零部分的水平集區域，從該圖中可以大致觀察到水平集函數的形態與收斂過程。其右側為對應的演化輪廓曲線（零水平集）。
從上述測試用例分割結果可以看出，水平集圖像分割算法可以實現對所需目標前景的精準分割。同時其也對待分割圖像中的噪聲具有一定的魯棒性。
對上述測試用例進行性能分析，分別取各圖片十次運行時間進行平均，得到測試用例串行時間表，具體如下表所示：

圖片名稱人臉嘴唇飛機游戲一飛機游戲二

串行時間(ms)	342158	34750	220024

從上表中的運行時間可以看出，水平集圖像分割雖然可以得到較好的分割效果，但演化時間較長，針對上述人臉嘴唇圖片需要接近6分鐘才可完成分割任務，極大的限制了該算法的應用，由此對應并行加速算法的設計具有一定的現實意義。

下期將對水平集圖像分割OpenMP并行加速算法的設計與實現進行講解。

總結

以上是生活随笔為你收集整理的水平集图像分割并行加速算法设计与实现(串行、OpenMP、CUDA)——串行实现篇的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。