水平集图像分割序列——Order LBF模型

1. 参考文献

2. 模型实现

2.1 Order-LBF模型Demo

%demo_Order_LBF.m
%Author: HSW
%Date；2015/4/12
%HARBIN INSTITUTE OF TECHNOLOGY
% Set Matlab
close all;
clear all;
clc;
% demo 编号,需要修改
ii = 1;
% Add path
addpath(genpath('Image\'));
addpath(genpath('Order_LBF solver\'));

% save result path
SaveFilePath = 'Results\';

% Read Image
c0 = 2;
imgID = 6;

Img = imread('Image\vessel.bmp');
Temp = Img;

if ndims(Img) == 3
    Img = rgb2gray(Img);
end
Img = double(Img);
% Initial phi is the level set function
switch imgID
    case 1
        phi= ones(size(Img(:,:,1))).*c0;
        a=43;b=51;c=20;d=28;
        phi(a:b,c:d) = -c0;
        figure;
        imshow(Temp);colormap;
        hold on;
        [c,h] = contour(phi, 0, 'r');
        hold off;
    case 2
        [m,n] = size(Img(:,:,1));
        a=m/2; b=n/2; r=min(m,n)/4;%set the radius
        phi= ones(m,n).*c0;
        phi(a-r:a+r,b-r:b+r) = -c0;
        imshow(Temp);colormap;
        hold on;
        [c,h] = contour(phi, 0, 'r');
        hold off;
    case 3
        figure;
        imshow(Temp);colormap;
        text(6,6,'Left click to get points, right click to get end point','FontSize',12,'Color', 'g');
        BW=roipoly;     %BW is mask
        phi=c0*2*(0.5-BW);
        hold on;
        [c,h] = contour(phi,[0 0],'r');
        hold off;
    case 4
        figure;
        imshow(Temp);colormap;
        rNum = 1;% the cicle number in a row
        cNum = 1;% the cicle number in a colcumn
        [m,n] = size(Img);
        phi= ones(m,n).*c0;
        r = min(m/(2*rNum),n/(2*cNum))/2;
        for i = 1:rNum
            for j = 1:cNum
                px = (2*i-1)*(m/(2*rNum));
                py = (2*j-1)*(n/(2*cNum));%(px,py) is the centre of the initial zero level set cicle
                phi(max(1,px - r):min(size(Img,1),px + r),max(1,py-r):min(size(Img,2),py+r)) = -c0; 
            end%for j
        end%for i
        hold on;
        [c,h] = contour(phi,[0 0],'r');
        hold off;
    case 5
        % 产生随机位置
        figure;
        imshow(Temp);colormap;
        rand('seed',0);
        boardsize = 20; %距离边界的位置
        
        r = 10; %产生圆形时为半径，产生矩形时为(1/2)*边长
        if r > boardsize
            r = boardsize;
        end
        possiblex = (boardsize + 1): (size(Img,1) - boardsize);
        possibley = (boardsize + 1): (size(Img,2) - boardsize);
        labelx = randperm(length(possiblex));
        labely = randperm(length(possibley));
        centrex = possiblex(labelx(1));
        centrey = possibley(labely(1));
        [m,n] = size(Img);
        phi= -ones(m,n).*c0;  
        phi(max(1,centrey-r):min(size(Img,1),centrey+r),max(1,centrex - r):min(size(Img,2),centrex + r)) = c0;
        hold on;
        [c,h] = contour(phi,[0 0],'r');
        hold off;
    case 6
        % 用鼠标获取中心位置
        figure;
        imshow(Temp);colormap;
%         [centrex,centrey] = ginput; % 按右键结束
%         centrex = uint8(centrex(1));
%         centrey = uint8(centrey(1));
        centrex = 52; 
        centrey = 40; 
        boardsize = 20; %距离边界的位置
        iscircle = 1; % 产生圆形,否则产生矩形
        r = 10; %产生圆形时为半径，产生矩形时为(1/2)*边长
        if r > boardsize
            r = boardsize;
        end
        [m,n] = size(Img);
        phi= ones(m,n).*c0;
        phi(max(1,centrey-r):min(size(Img,1),centrey+r),max(1,centrex - r):min(size(Img,2),centrex + r)) = -c0;
        hold on;
        [c,h] = contour(phi,[0 0],'r');
        hold off;
end%switch imgID

iterNum = 300; % the total number of the iteration
lambda1 = 1;   % the weight parameter of the globe term which is inside the level set
lambda2 = 1;   % the weight parameter of the globe term which is ouside the level set
mu = 0.002*255*255; % the weight parameter of the length term
nu = 0; % the weight parameter of the area term
pu = 1.0; %the weight parameter of the penalizing term
timestep = 0.1; % the time step of the level set function evolution
epsilon = 1.0; % the parameter of the Heaviside function and the Dirac delta function

%scale parameter in Gaussian kernel
% this part code is for the LBF model and LIGF model
sigmaG = 3.0;
Klbf = fspecial('gaussian',round(2*sigmaG)*2+1,sigmaG); % Gaussian kernel
KIG = conv2(Img,Klbf,'same'); % KIG is the result of Gaussian kernel Klbf convol Img
KONE = conv2(ones(size(Img)),Klbf,'same'); % KONE is the convolution of Gaussian kernel Klbf convol all one matrix



% all model's initial level set is same
phi_Order_LBF = phi;
phi_star = phi;
%
% figure;
% imshow(Temp); colormap;

%start the level set evolution



% ORDER_LBF
time = cputime;
for iter = 1:iterNum
    numIter=1;
    %level set evolution.
    phi_Order_LBF = EVOL_Order_LBF(phi_Order_LBF,Img,Klbf,KIG,KONE,0.5*mu,timestep,pu,lambda1,lambda2,epsilon,numIter);
    if mod(iter,10) == 0
        pause(0.001);
        imagesc(Img, [0, 255]);colormap(gray);hold on; axis off;
        contour(phi_LBF,[0,0],'r');
        contour(phi_LCV,[0,0],'g');
        contour(phi_LGIF,[0,0],'b');
        contour(phi_LGD,[0,0],'w');
        contour(phi_Order_LBF,[0,0],'k');
        contour(phi_Multi_LBF,[0,0],'c');
    end
end %for
totaltime_Order_LBF = cputime - time;





% Display Results
figure;
imshow(Temp);
hold on;
contour(phi_star,[0,0],'r','linewidth',1);
title('Initial Level set');

figure;
imshow(Temp);
hold on;
contour(phi_Order_LBF,[0,0],'r','linewidth',1);
title('Result of Order LBF model');


% Save Results
OrderLBFFilePath = [SaveFilePath,'Order_LBF\Demo',num2str(ii),'.bmp'];

SaveOrderLBF = phi_Order_LBF >= 0;
imwrite(SaveOrderLBF,OrderLBFFilePath,'bmp');

2.2 Order-LBF模型实现

function u = EVOL_Order_LBF(u0,Img,K,KI,KONE,nu,timestep,mu,lambda1,lambda2,epsilon,numIter)
%  u = EVOL_LBF(u0,Img,Ksigma,KI,KONE,nu,timestep,mu,lambda1,lambda2,epsilon,numIter)

u=u0;
for k1=1:numIter
    u=NeumannBoundCond(u);
    C=curvature_central(u);    % div()  
    HeavU=Heaviside(u,epsilon);
    DiracU=Dirac(u,epsilon);
    
    [f1,f2]=LBF_LocalBinaryFit(K,Img,KI,KONE,HeavU);  
    TempMin = min(f1,f2); 
    f1 = max(f1,f2); 
    f2 = TempMin; 
    LBF=LBF_dataForce(Img,K,KONE,f1,f2,lambda1,lambda2);
  
    areaTerm=-DiracU.*LBF;
    penalizeTerm=mu*(4*del2(u)-C);
    lengthTerm=nu.*DiracU.*C;
    u=u+timestep*(lengthTerm+penalizeTerm+areaTerm);
end

% Make a function satisfy Neumann boundary condition
function g = NeumannBoundCond(f)
[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]);  

function k = curvature_central(u)
% compute curvature for u with central difference scheme
[ux,uy] = gradient(u);
normDu = sqrt(ux.^2+uy.^2+1e-10);
Nx = ux./normDu;
Ny = uy./normDu;
[nxx,junk] = gradient(Nx);
[junk,nyy] = gradient(Ny);
k = nxx+nyy;

function [f1,f2] = LBF_LocalBinaryFit(K,Img,KI,KONE,H)
I=Img.*H;
c1=conv2(H,K,'same');
c2=conv2(I,K,'same');
f1=c2./(c1);
f2=(KI-c2)./(KONE-c1);

function h = Heaviside(x,epsilon)     % function (11)
h=0.5*(1+(2/pi)*atan(x./epsilon));

function f = Dirac(x, epsilon)    % function (12)
f=(epsilon/pi)./(epsilon^2.+x.^2);

function f=LBF_dataForce(Img,K,KONE,f1,f2,lamda1,lamda2)
s1=lamda1.*f1.^2-lamda2.*f2.^2;
s2=lamda1.*f1-lamda2.*f2;
f=(lamda1-lamda2)*KONE.*Img.*Img+conv2(s1,K,'same')-2.*Img.*conv2(s2,K,'same');

3. 模型效果