【图像分割】基于计算机视觉实现道路提取分割附matlab代码_start level set evolution代码

1 内容介绍

MATLAB作为功能强大的仿真分析软件,被广泛用于科学研究和解决各种具体问题,而它在图像处理方面也发挥了不可小觑的优势.基于此介绍了图像分割的基本理论和常用方法,包括基于阈值,区域特性图像分割,并介绍了一种基于计算机视觉的图像分割方法.然后借助MATLAB平台,分别对这些方法进行了图像仿真,并分析了仿真效率与效果.最后提出了多种分割方法整合的观点.

2 仿真代码

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:iter
    phi=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 zero
    Ny=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.
    else
        disp('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 force
    edgeTerm=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]);  

This code is used to compare zhenzhou threshold selection method
%%%% with others. This code should be run many times since the synthesized
%%%% images vary a little bit every time
close all;
clear all;
clc;

%%%%Synthesizing the image
S=30;
C1=40+S*randn(128,128);
Center1=90+S*randn(128,8);
Center0=140+S*randn(128,8);
Center2=180+S*randn(8,128);
Center3=230+S*randn(8,128);
C1(:,41:48)=Center1;
C1(:,81:88)=Center0;
C1(41:48,:)=Center2;
C1(81:88,:)=Center3;
C=double(C1);
figure,imagesc(C);colormap gray;

IS=size(C1);
Seg1=double(C1);
S=size(Seg1);
    for l=1:10
        for i=2:IS(1)-1
            for j=2:IS(2)-1
                meanc=(Seg1(i+1,j+1)+Seg1(i+1,j-1)+Seg1(i-1,j+1)+Seg1(i-1,j-1)+Seg1(i+1,j)+Seg1(i,j+1)+Seg1(i-1,j)+Seg1(i,j-1)+Seg1(i,j))/9;
                Seg1(i,j)=round(meanc);
            end;
        end;
    end;
    Seg1(1:2,:)=Seg1(3:4,:);
    Seg1((IS(1)-1):IS(1),:)=Seg1((IS(1)-3):(IS(1)-2),:);
    Seg1(:,1:2)=Seg1(:,3:4);
    Seg1(:,(IS(2)-1):IS(2))=Seg1(:,(IS(2)-3):(IS(2)-2));

Seg=Seg1;
figure,imagesc(Seg);colormap gray;
%%%%%%%%%%%End of Synthesizing the image

%%%Zhenzhou threshold selection method
T=zhenzhou_threshold_selection_updated(double(Seg),2,15,1);
bSeg=Seg>T;
figure,imagesc(bSeg);colormap gray;

      
 
%%%%EM Results
[Segs,mu,v,p]=EMSeg(Seg,2);
figure,imagesc(Segs);colormap gray;

%%%%%%%%%%%%%%kmeans results 

[mu,Segs]=kmeans(uint8(Seg),2);
figure,imagesc(Segs);colormap gray;

%  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.
% clear all;
% close all;

    tic;
    Img=double(Seg);
    %% parameter setting
    timestep=5;  % time step
    mu=0.2/timestep;  % coefficient of the distance regularization term R(phi)
%     iter_inner=5;
%     iter_outer=40;
    iter_inner=20;
    iter_outer=80;
    lambda=5; % coefficient of the weighted length term L(phi)
    alfa=1.5;  % coefficient of the weighted area term A(phi)
    epsilon=1.5; % papramater that specifies the width of the DiracDelta function
    
    sigma=1.5;     % scale parameter in Gaussian kernel
    G=fspecial('gaussian',15,sigma);
    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 a rectangle
    NS=size(Img);
    % initialLSF(10:55, 10:75)=-c0;
    initialLSF(10:NS(1)-10,10:NS(2)-10)=-c0;
%     initialLSF(round(NS(1)/2-10):round(NS(1)/2+10),round(NS(2)/2-10):round(NS(2)/2+10))=-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 ==1
        potentialFunction = 'single-well';  % use single well potential p1(s)=0.5*(s-1)^2, which is good for region-based model
    elseif potential == 2
        potentialFunction = 'double-well';  % use double-well potential in Eq. (16), which is good for both edge and region based models
    else
        potentialFunction = 'double-well';  % default choice of potential function
    end
    
    
    % start level set evolution
    for n=1:iter_outer
        phi = drlse_edge(phi, g, lambda, mu, alfa, epsilon, timestep, iter_inner, potentialFunction);
        if mod(n,2)==0
            figure(2);
            imagesc(Img,[0, 255]); axis off; axis equal; colormap(gray); hold on;  contour(phi, [0,0], 'r','LineWidth',3);
        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','LineWidth',3);
    str=['Final zero level contour, ', num2str(iter_outer*iter_inner+iter_refine), ' iterations'];
    title(str);
    
    pause(1);
    figure;
    mesh(-finalLSF); % for a better view, the LSF is displayed upside down
    hold on;  contour(phi, [0,0], 'r-','LineWidth',3);
%     legend('Dis-LS');
%     str=['Final level set function, ', num2str(iter_outer*iter_inner+iter_refine), ' iterations'];
%     title(str);
    axis on;
    
    toc;

% Description: This code implements the paper: "Active Contours Without
% Edges" By Chan Vese.

I=uint8(Seg);
tic;
%     figure,imagesc(I);colormap gray;
m = zeros(size(I,1),size(I,2));          %-- create initial mask
NS=size(I);
 m(30:NS(1)-10,2:NS(2)-10) = 1;
seg = region_seg(I, m, 1200); %-- Run segmentation
figure,imagesc(seg);colormap gray;
toc;

% Shawn Lankton's IEEE TIP 2008 paper 'Localizing Region-Based Active
% Contours'.

% clc;clear all;close all;
imgID = 1; % 1,2,3  % choose one of the five test images
    
%     path='Images\test\100007.jpg';
%     path='Images\test\8068.jpg';%%C(:,:,1)
%     
%     path='Images\test\388067.jpg';%%C(:,:,1)
%     path='Images\val\3096.jpg';%%C(:,:,1)
%     path='Images\val\86016.jpg';%%C(:,:,3)
%     Seg=im2double(imread(path));
%     C=Seg(:,:,3)*255;
Img=(Seg);
NS=size(Img);
 
imgID = 3;
tic;
epsilon = 1;
switch imgID

    case 1
        num_it =1000;
        rad = 8;
        alpha = 0.3;% coefficient of the length term
        mask_init  = zeros(size(Img(:,:,1)));
        mask_init(10:NS(1)-10,10:NS(2)-10) = 1;
        seg = local_AC_MS(Img,mask_init,rad,alpha,num_it,epsilon);
    case 2
        num_it =1800;
        rad = 9;
        alpha = 0.003;% coefficient of the length term
        mask_init = zeros(size(Img(:,:,1)));
        mask_init(20:NS(1)-20,20:NS(2)-20) = 1;
        seg = local_AC_UM(Img,mask_init,rad,alpha,num_it,epsilon);
    case 3
        num_it = 2000;
        rad = 5;
        alpha = 0.001;% coefficient of the length term
        mask_init  = zeros(size(Img(:,:,1)));
        mask_init(10:NS(1)-10,10:NS(2)-10) = 1;
        seg = local_AC_UM(Img,mask_init,rad,alpha,num_it,epsilon);
end
toc;

%%%%%%%%%%%Otsu
X=(Seg);

tic;
for n = 2:2
    IDX = otsu(X,n);   
end
toc;
figure,imagesc(IDX);colormap gray;

 

%%%%%%%%%%%%ISO data
% Reference :T.W. Ridler, S. Calvard, Picture thresholding using an iterative selection method, 
%    IEEE Trans. System, Man and Cybernetics, SMC-8 (1978) 630-632.
    
 C=Seg;   
level = isodata(uint8(C));
 BW = im2bw(uint8(C),level);
 figure,imagesc(BW);colormap gray;
      
 

 
 %%%%maximum entropy    

 C=Seg;   
    
[threshold_max I1]=maxentropie(uint8(C));
figure,imagesc(I1);colormap gray;

 

 %%%cross entropy

C=Seg;

[threshold_cross I1]=minCE(uint8(C));
figure,imagesc(I1);colormap gray;

 

 %%%fuzzy entropy
C=Seg;   
    
NS=size(C);
T_fuzzy= entropy_fuzzy_threshold(C)
BC=zeros(NS(1),NS(2));
for i=1:NS(1)
    for j=1:NS(2)
        if C(i,j)>T_fuzzy
            BC(i,j)=1;
        end;
    end;
end;
figure,imagesc(BC);colormap gray;

%%%%%%%%%%%2015 thresholding
C=Seg;   
MG=fth(C,2,1,1.5);  
figure,imagesc(MG);colormap gray;
% MG1=MG<2;
% figure,imagesc(MG1);colormap gray;
% I1=MG-1;

%%%fuzzy Cmeans
C=Seg;    
fim=mat2gray(C);
level=graythresh(fim);
bwfim=im2bw(fim,level);

[bwfim0,level0]=fcmthresh(fim,0);
[bwfim1,level1]=fcmthresh(fim,1);
figure,imagesc(fim);colormap gray;
figure,imagesc(bwfim);colormap gray;

 

3 运行结果

4 参考文献

[1]马永慧. "基于MATLAB的图像分割的技术研究." 山西电子技术 4(2012):2.​

博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。

部分理论引用网络文献，若有侵权联系博主删除。