當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

水平集分割

發布時間：2025/3/17 编程问答 30 豆豆

生活随笔收集整理的這篇文章主要介紹了水平集分割小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

基于距離正則的水平集分割MATLAB代碼，無需初始化

% This Matlab code demonstrates an edge-based active contour model as an application of % the Distance Regularized Level Set Evolution (DRLSE) formulation in the following paper: % % C. Li, C. Xu, C. Gui, M. D. Fox, "Distance Regularized Level Set Evolution and Its Application to Image Segmentation", % IEEE Trans. Image Processing, vol. 19 (12), pp. 3243-3254, 2010. % % Author: Chunming Li, all rights reserved % E-mail: lchunming@gmail.com % li_chunming@hotmail.com % URL: http://www.imagecomputing.org/~cmli// clear all; close all; Img=imread('gourd.bmp'); Img=double(Img(:,:,1)); %% parameter setting timestep=1; % time step mu=0.2/timestep; % coefficient of the distance regularization term R(phi) iter_inner=5; iter_outer=20; lambda=5; % coefficient of the weighted length term L(phi) alfa=-3; % coefficient of the weighted area term A(phi) epsilon=1.5; % papramater that specifies the width of the DiracDelta function sigma=.8; % scale parameter in Gaussian kernel G=fspecial('gaussian',15,sigma); % Caussian kernel Img_smooth=conv2(Img,G,'same'); % smooth image by Gaussiin convolution [Ix,Iy]=gradient(Img_smooth); f=Ix.^2+Iy.^2; g=1./(1+f); % edge indicator function. % initialize LSF as binary step function c0=2; initialLSF = c0*ones(size(Img)); % generate the initial region R0 as two rectangles initialLSF(25:35,20:25)=-c0; initialLSF(25:35,40:50)=-c0; phi=initialLSF; figure(1); mesh(-phi); % for a better view, the LSF is displayed upside down hold on; contour(phi, [0,0], 'r','LineWidth',2); title('Initial level set function'); view([-80 35]); figure(2); imagesc(Img,[0, 255]); axis off; axis equal; colormap(gray); hold on; contour(phi, [0,0], 'r'); title('Initial zero level contour'); pause(0.5); potential=2; if potential ==1potentialFunction = 'single-well'; % use single well potential p1(s)=0.5*(s-1)^2, which is good for region-based model elseif potential == 2potentialFunction = 'double-well'; % use double-well potential in Eq. (16), which is good for both edge and region based models elsepotentialFunction = 'double-well'; % default choice of potential function end % start level set evolution for n=1:iter_outerphi = drlse_edge(phi, g, lambda, mu, alfa, epsilon, timestep, iter_inner, potentialFunction); if mod(n,2)==0figure(2);imagesc(Img,[0, 255]); axis off; axis equal; colormap(gray); hold on; contour(phi, [0,0], 'r');end end % refine the zero level contour by further level set evolution with alfa=0 alfa=0; iter_refine = 10; phi = drlse_edge(phi, g, lambda, mu, alfa, epsilon, timestep, iter_inner, potentialFunction); finalLSF=phi; figure(2); imagesc(Img,[0, 255]); axis off; axis equal; colormap(gray); hold on; contour(phi, [0,0], 'r'); hold on; contour(phi, [0,0], 'r'); str=['Final zero level contour, ', num2str(iter_outer*iter_inner+iter_refine), ' iterations']; title(str); figure; mesh(-finalLSF); % for a better view, the LSF is displayed upside down hold on; contour(phi, [0,0], 'r','LineWidth',2); view([-80 35]); str=['Final level set function, ', num2str(iter_outer*iter_inner+iter_refine), ' iterations']; title(str); axis on; [nrow, ncol]=size(Img); axis([1 ncol 1 nrow -5 5]); set(gca,'ZTick',[-3:1:3]); set(gca,'FontSize',14) % This Matlab code demonstrates an edge-based active contour model as an application of % the Distance Regularized Level Set Evolution (DRLSE) formulation in the following paper: % % C. Li, C. Xu, C. Gui, M. D. Fox, "Distance Regularized Level Set Evolution and Its Application to Image Segmentation", % IEEE Trans. Image Processing, vol. 19 (12), pp. 3243-3254, 2010. % % Author: Chunming Li, all rights reserved % E-mail: lchunming@gmail.com % li_chunming@hotmail.com % URL: http://www.imagecomputing.org/~cmli// clear all; close all; Img=imread('gourd.bmp'); Img=double(Img(:,:,1)); %% parameter setting timestep=1; % time step mu=0.2/timestep; % coefficient of the distance regularization term R(phi) iter_inner=5; iter_outer=20; lambda=5; % coefficient of the weighted length term L(phi) alfa=-3; % coefficient of the weighted area term A(phi) epsilon=1.5; % papramater that specifies the width of the DiracDelta function sigma=.8; % scale parameter in Gaussian kernel G=fspecial('gaussian',15,sigma); % Caussian kernel Img_smooth=conv2(Img,G,'same'); % smooth image by Gaussiin convolution [Ix,Iy]=gradient(Img_smooth); f=Ix.^2+Iy.^2; g=1./(1+f); % edge indicator function. % initialize LSF as binary step function c0=2; initialLSF = c0*ones(size(Img)); % generate the initial region R0 as two rectangles initialLSF(25:35,20:25)=-c0; initialLSF(25:35,40:50)=-c0; phi=initialLSF; figure(1); mesh(-phi); % for a better view, the LSF is displayed upside down hold on; contour(phi, [0,0], 'r','LineWidth',2); title('Initial level set function'); view([-80 35]); figure(2); imagesc(Img,[0, 255]); axis off; axis equal; colormap(gray); hold on; contour(phi, [0,0], 'r'); title('Initial zero level contour'); pause(0.5); potential=2; if potential ==1potentialFunction = 'single-well'; % use single well potential p1(s)=0.5*(s-1)^2, which is good for region-based model elseif potential == 2potentialFunction = 'double-well'; % use double-well potential in Eq. (16), which is good for both edge and region based models elsepotentialFunction = 'double-well'; % default choice of potential function end % start level set evolution for n=1:iter_outerphi = drlse_edge(phi, g, lambda, mu, alfa, epsilon, timestep, iter_inner, potentialFunction); if mod(n,2)==0figure(2);imagesc(Img,[0, 255]); axis off; axis equal; colormap(gray); hold on; contour(phi, [0,0], 'r');end end % refine the zero level contour by further level set evolution with alfa=0 alfa=0; iter_refine = 10; phi = drlse_edge(phi, g, lambda, mu, alfa, epsilon, timestep, iter_inner, potentialFunction); finalLSF=phi; figure(2); imagesc(Img,[0, 255]); axis off; axis equal; colormap(gray); hold on; contour(phi, [0,0], 'r'); hold on; contour(phi, [0,0], 'r'); str=['Final zero level contour, ', num2str(iter_outer*iter_inner+iter_refine), ' iterations']; title(str); figure; mesh(-finalLSF); % for a better view, the LSF is displayed upside down hold on; contour(phi, [0,0], 'r','LineWidth',2); view([-80 35]); str=['Final level set function, ', num2str(iter_outer*iter_inner+iter_refine), ' iterations']; title(str); axis on; [nrow, ncol]=size(Img); axis([1 ncol 1 nrow -5 5]); set(gca,'ZTick',[-3:1:3]); set(gca,'FontSize',14) function phi = drlse_edge(phi_0, g, lambda,mu, alfa, epsilon, timestep, iter, potentialFunction) % This Matlab code implements an edge-based active contour model as an % application of the Distance Regularized Level Set Evolution (DRLSE) formulation in Li et al's paper: % % C. Li, C. Xu, C. Gui, M. D. Fox, "Distance Regularized Level Set Evolution and Its Application to Image Segmentation", % IEEE Trans. Image Processing, vol. 19 (12), pp.3243-3254, 2010. % % Input: % phi_0: level set function to be updated by level set evolution % g: edge indicator function % mu: weight of distance regularization term % timestep: time step % lambda: weight of the weighted length term % alfa: weight of the weighted area term % epsilon: width of Dirac Delta function % iter: number of iterations % potentialFunction: choice of potential function in distance regularization term. % As mentioned in the above paper, two choices are provided: potentialFunction='single-well' or % potentialFunction='double-well', which correspond to the potential functions p1 (single-well) % and p2 (double-well), respectively.% % Output: % phi: updated level set function after level set evolution % % Author: Chunming Li, all rights reserved % E-mail: lchunming@gmail.com % li_chunming@hotmail.com % URL: http://www.imagecomputing.org/~cmli/ phi=phi_0; [vx, vy]=gradient(g); for k=1:iterphi=NeumannBoundCond(phi);[phi_x,phi_y]=gradient(phi);s=sqrt(phi_x.^2 + phi_y.^2);smallNumber=1e-10; Nx=phi_x./(s+smallNumber); % add a small positive number to avoid division by zeroNy=phi_y./(s+smallNumber);curvature=div(Nx,Ny);if strcmp(potentialFunction,'single-well')distRegTerm = 4*del2(phi)-curvature; % compute distance regularization term in equation (13) with the single-well potential p1.elseif strcmp(potentialFunction,'double-well');distRegTerm=distReg_p2(phi); % compute the distance regularization term in eqaution (13) with the double-well potential p2.elsedisp('Error: Wrong choice of potential function. Please input the string "single-well" or "double-well" in the drlse_edge function.');end diracPhi=Dirac(phi,epsilon);areaTerm=diracPhi.*g; % balloon/pressure forceedgeTerm=diracPhi.*(vx.*Nx+vy.*Ny) + diracPhi.*g.*curvature;phi=phi + timestep*(mu*distRegTerm + lambda*edgeTerm + alfa*areaTerm); end function f = distReg_p2(phi) % compute the distance regularization term with the double-well potential p2 in eqaution (16) [phi_x,phi_y]=gradient(phi); s=sqrt(phi_x.^2 + phi_y.^2); a=(s>=0) & (s<=1); b=(s>1); ps=a.*sin(2*pi*s)/(2*pi)+b.*(s-1); % compute first order derivative of the double-well potential p2 in eqaution (16) dps=((ps~=0).*ps+(ps==0))./((s~=0).*s+(s==0)); % compute d_p(s)=p'(s)/s in equation (10). As s-->0, we have d_p(s)-->1 according to equation (18) f = div(dps.*phi_x - phi_x, dps.*phi_y - phi_y) + 4*del2(phi); function f = div(nx,ny) [nxx,junk]=gradient(nx); [junk,nyy]=gradient(ny); f=nxx+nyy; function f = Dirac(x, sigma) f=(1/2/sigma)*(1+cos(pi*x/sigma)); b = (x<=sigma) & (x>=-sigma); f = f.*b; function g = NeumannBoundCond(f) % Make a function satisfy Neumann boundary condition [nrow,ncol] = size(f); g = f; g([1 nrow],[1 ncol]) = g([3 nrow-2],[3 ncol-2]); g([1 nrow],2:end-1) = g([3 nrow-2],2:end-1); g(2:end-1,[1 ncol]) = g(2:end-1,[3 ncol-2]);

LeveSet 水平集方法主要的思想是利用三維（高維）曲面的演化來表示二維曲線的演化過程。在計算機視覺領域，利用水平集方法可以實現很好的圖像分割效果。

1.數學原理

根據維基百科的定義，在數學上一個包含n個變量的實值函數其水平集可以表示為下面的公式：
$Lc(f)=(x1,x2,...,xn)∣f(x1,x2,...,xn)=cL_c(f) = {(x_1,x_2,...,x_n)|f(x_1,x_2,...,x_n) = c}Lc?(f)=(x1?,x2?,...,xn?)∣f(x1?,x2?,...,xn?)=c$
可以看出，水平集指的是這個函數的取值為一個給定的常數c.那么當變量個數為2時，這個函數的水平集就變味了一條曲線，也可以成為等高線。這時函數f就可以描述一個曲面。（同樣，三個變量時就能得到一個等值面，大于3的就是水平超面了）
下面來看一個簡單的曲面的水平集（等值線的例子）：

上圖中，左邊是Himmelblau方程所描述的曲面即 $z = f (x, y) z = f (x, y) z = f (x, y)$ ，右邊就是這個曲面的等高線集合，其中每一條就對應著一個常數c的水平集。

2.水平集方法

在計算機視覺中，利用水平集方法的優勢在于，可以利用這種方法在固定的坐標系中，計算曲線、曲面的演化，而無需知道曲線曲面的參數，所以這種方法又稱為幾何法。
具體來講，在圖像分割任務中如果要得到每個物體準確的包絡曲線，就需要對去描述這個曲線在xy坐標系下的演化。最初曲線在二維圖像平面上是這樣的：
$y = g (x) y = g (x) y = g (x)$
圖割過程就是要找出一條比較好的曲面來包圍物體。
但是這個函數求解和計算比較困難，那么就可以用水平集的思想和方法來解決：

1.首先將這個曲線看成是某個三維曲面下的某一條等高線：
$z = f (x, y) z = f (x, y) z = f (x, y)$
$f (x, y) f (x, y) f (x, y)$ 就可以看做是一個xy的隱函數方程。

2.特別的，在二維圖像領域一般將上面講到的函數 $z = f (x, y) z = f (x, y) z = f (x, y)$ 設置為0，剛剛所說的曲面的零水平集就是圖像上的邊緣包絡。
$z = f (x, y) = 0 z = f (x, y) = 0 z = f (x, y) = 0$
或者更為標準的形式如下：
$Γ=(x,y)∣z(x,y)=0Γ=(x,y)∣z(x,y)=0Γ=(x,y)∣z(x,y)=0Γ=(x,y)∣z(x,y)=0\Gamma = {(x,y)| z(x,y)=0}Γ=(x,y)∣z(x,y)=0$
這時，z(x,y)=f(x,y)就是輔助函數，用三維的曲面來輔助表示二維的曲線。

這里還有一個假設， $ΓΓΓΓ\GammaΓ$ 零水平集內部的（包絡內的）為正外面為負。

這里有個重要的性質，基于高維曲面水平集方法得到的這個零水平集曲線是封閉的、連續的、處處可導的，這就為后面的圖像分割提供了基礎。

3.零水平集演化

為了獲得與圖像邊緣相適應的結果，z需要進行演化才能使得其0水平集所描述的曲線很好的包裹住需要分割的物體，那么此時就涉及到方程的演化：
$Γ(t)=(x,y)∣z(x,y)=0Γ(t)=(x,y)∣z(x,y)=0Γ(t)=(x,y)∣z(x,y)=0Γ(t)=(x,y)∣z(x,y)=0\Gamma(t) = {(x,y)| z(x,y)=0}Γ(t)=(x,y)∣z(x,y)=0$
$?z?t=v∣Δz∣?z?t=v∣Δz∣?t?z?=v∣Δz∣?z?t=v∣Δz∣\dfrac{\partial z}{\partial t} = v |\Delta z|?t?z?=v∣Δz∣$
這里是一個曲面演化后，其零水平集演化的過程：

要促使邊緣演化成包絡object的形式，需要一個所謂的圖像力來促使其變化，這時候就需要利用圖像中的梯度來作為這個圖像力的組成部分驅動曲線去靠近邊緣。下圖是一個例子，作圖是原圖，右圖是梯度圖。

<font color=dodgerblue<那么這個圖像力就需要在原理邊緣的時候大，靠近邊緣的時候縮小，直到貼合邊緣。

在梯度的驅動下，其演化過程如上圖所示。

4.如何求解最小值

上面的問題在數學中可以歸結為能量泛函問題，求解邊緣的過程就是這個水平集能量泛函最小化的過程。
著名的方法來自于下面這篇文章：
Level Set Evolution Without Re-initialization: A New Variational Formulation
作者(李明純)通過引入了內部能量項來懲罰由信號距離函數造成的水平集函數的偏差,以及外部能量項來驅動零水平集向期望的圖像特征運動，實現了很好的分割效果。

5.Demo

下圖是matlab 中，稍作修改運行作者給出的工具函數demo：
原圖和分割效果圖：

pics contenx ref from:
https://profs.etsmtl.ca/hlombaert/algorithms.php
https://math.berkeley.edu/~sethian/2006/Explanations/level_set_explain.html
https://www.zhihu.com/question/22608763
https://blog.csdn.net/github_35768306/article/details/64129197
https://blog.csdn.net/songzitea/article/details/46385271
李明純：http://www.engr.uconn.edu/~cmli/
paper:http://www.imagecomputing.org/~cmli/paper/levelset_cvpr05.pdf
matlab：https://www.mathworks.com/matlabcentral/profile/authors/870631-chunming-li
https://blog.csdn.net/yutianxin123/article/details/69802364
cells: https://www.histology.leeds.ac.uk/blood/blood_wbc.php

總結

以上是生活随笔為你收集整理的水平集分割的全部內容，希望文章能夠幫你解決所遇到的問題。

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

水平

上一篇：在Browser Application
下一篇： TensorFlow生成.mat文件