1. 背景

在多相CV模型中(https://blog.csdn.net/hit1524468/article/details/79706174)， 我們注意到隨著迭代次數的增加，水平集函式波動的範圍開始逐漸增大，這就是水平集的符號函式重新初始化問題；Li Cunming 提出了對水平集函式進行卷積的方法，避免了水平集的重新初始化問題(參考文獻：Level Set Evolution Without Reinitialization: A New Variational Formulation)

2. 改進多相CV模型

%   This code implements the improved Vese-Chan multiphase level set
%   model in [1] with modification by adding the distance regularizing (DR) term
%   introduced in Li et al's paper [2].
%
%   Reference:
%   [1] A Multiphase Level Set Framework for Image Segmentation Using the Mumford and Shah Model, IJCV 2002.
%   [2] Level Set Evolution Without Reinitialization: A New Variational Formulation", CVPR 2005.
%
%   Although reinitialization in Vese and Chan's original level set methods
%   is not required for all images, the level set function can still be
%   degraded after a certain number of iterations. This not only causes the curve irregular, but also slow down 
%   the curve evolution.By adding the distance regularizing term (the internal energy)in Li's level set method [2], 
%   the need for reinitialization is completely eliminated, and the level set function is smooth during the 
%   evolution, which ensure more accurate computation.
%
%   Note that the distance regularizing term has an effect of maintaining the level set function as an approximate 
%   signed distance function near the zero level set. 
%
%   Note: There may be more sophiscated numerical schemes with better performance than the one in
%   this implementation. We only use a simple difference scheme in this version to demonstrate the desirable 
%   distance regularizing effect.
%   Author: Chunming Li

clear;
delta_t = .1;
lambda_1=1;
lambda_2=1;
h = 1;
epsilon=1;
nu = .001*255*255;
fun_n=2;  % two level set functions for 4 phase level set formulation

Img=imread('fourblock_gray.bmp');  %ԭʼ
% syn = zeros(257);
% [m,n] = size(syn);
% r1 = 30;
% r2 = 60;
% r3 = 100;
% for i = 1:m
%     for j = 1:n
%         if abs(i - 129) <= r1 && abs( j - 129) <= r1
%             syn(i,j) = 40 + 3.*randn(1,1);   
%         elseif abs(i - 129) <= r2 && abs( j - 129) <= r2 && abs(i - 129)> r1 && abs(j - 129) > r1
%             syn(i,j) = 80 + 3.*randn(1,1);
%         elseif abs(i - 129) <= r3 && abs(j - 129) <= r3 &&  abs(i - 129)> r2 && abs(j - 129) > r2
%             syn(i,j) = 120 + 3.*randn(1,1);
%         else 
%             syn(i,j) = 160 + 3.*randn(1,1);
%         end
%     end
% end
% Img = syn;




U=Img(:,:,1);

% get the size
[nrow, ncol] =size(U);

ic=nrow/2;
jc=ncol/2;
r=30;
phi = initial_sdf2circle(nrow,ncol,ic,jc,r,fun_n); 
I=double(U);

figure;
imagesc(uint8(I));colormap(gray)
hold on;
plotLevelSet(phi(:,:,1),0,'r');
plotLevelSet(phi(:,:,2),0,'b');

numIter = 10;

for k=1:70
    mu=0.5;  % coefficient for the distance regularizing term (internal energy) in Li's CVPR05 paper
    phi=EVOLUTION_4PHASE_DR(I, phi, nu, lambda_1, lambda_2, mu, delta_t, epsilon, numIter);   % update level set function
    if mod(k,2)==0
        pause(.1);
        imagesc(uint8(I));colormap(gray)
        hold on;
          phi_1=phi(:,:,1);
          phi_2=phi(:,:,2);
          plotLevelSet(phi_1,0,'r');     
          plotLevelSet(phi_2,0,'b'); 
          hold off;
   end           
end
figure;mesh(phi_1);  
title('\phi_1, improved Vese-Chan model by Li');
figure;mesh(phi_2);
title('\phi_2, improved Vese-Chan model by Li');
 

function phi = EVOLUTION_4PHASE_DR(I, phi0, nu, lambda_1, lambda_2, mu, delta_t, epsilon, numIter);
%  
%   This function updates the level set function in the improved Vese-Chan multiphase level set
%   model in [1] by adding the distance regularizing term introduced in Li et al's paper [2].
%
%   Reference:
%   [1] A Multiphase Level Set Framework for Image Segmentation Using the Mumford and Shah Model, IJCV 2002.
%   [2] Level Set Evolution Without Reinitialization: A New Variational Formulation", CVPR 2005.
% 

phi(:,:,1)=phi0(:,:,1);
phi(:,:,2)=phi0(:,:,2);

for m=1:numIter
    for k=1:2
        tmp_phi=phi(:,:,k);
        tmp_phi=NeumannBoundCond(tmp_phi);
        delta_h=Delta(tmp_phi,epsilon);  
        differenceScheme=0; % use 0, 1, 2 for different schemes to compute curvature
        Curv = CURVATURE_CV(tmp_phi,differenceScheme); 
        [C,mult]=quadrifit(phi,I,epsilon,2);
        b=zeros(size(Curv));

        if k==1
            b=(C(1)-C(2))*(C(1)+C(2)-2*I).*mult(:,:,1,1);
            b=b+(C(3)-C(4))*(C(3)+C(4)-2*I).*mult(:,:,1,2);
        else
            b=(C(1)-C(3))*(C(1)+C(3)-2*I).*mult(:,:,2,1);
            b=b+(C(2)-C(4))*(C(2)+C(4)-2*I).*mult(:,:,2,2);            
        end
        tmp_phi=tmp_phi+delta_t*(delta_h.*(nu*Curv+b) + mu*(4*del2(tmp_phi)-Curv)); % the last term is the distance regularizing term in Li's CVPR05 paper 
        phi(:,:,k)=tmp_phi;        
    end
end

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]);          
 

3 模型效果