FCMLSM模糊聚类水平集

FCMLSM 模糊聚类 变分法 水平集

function [imglsoptssls]=FCLSM（imgimgfcmbeta）

img=double（img）;

se=5;       %template radius for spatial filtering
sigma=2;    %spatial filter weight
d0=.5;      %fuzzy thresholding
epsilon=1.5;    %Dirac regulator

%adaptive definition of penalizing item mu
u=d0<=imgfcm;
bwa=bwarea（u）;  %area of initial contour
bw2=bwperim（u）;
bwp=sum（sum（bw2））;  %peripherium of initial contour
mu=bwp/bwa;     %Coefficient of the internal （penalizing） energy term P（\phi）;
opts.mu=mu;
timestep=0.2/mu; %The product timestep*mu must be less than 0.25 for stability
opts.timestep=timestep;
%end

% fs=fspecial（‘disk‘5）;
% imgs=imfilter（imgfs‘replicate‘）;
fs=fspecial（‘gaussian‘sesigma）;
img_smooth=conv2（double（img）double（fs）‘same‘）;
[IxIy]=gradient（img_smooth）;
f=Ix.^2+Iy.^2;
g=1./（1+f）;  % edge indicator function.

% define initial level set function as -c0 c0 
%   at points outside and inside of a region R respectively.
u=u-0.5;
u=4*epsilon*u;
sls（1::）=double（u）;

lambda=1/mu;
opts.lambda=lambda;
% nu=2*imgfcm-1;     % Coefficient of the weighted area term Ag（\phi）;
nu=-2*（2*beta*imgfcm-（1-beta））;
opts.alf=nu;
%Note: Choose a positive（negative） alf if the initial contour is
% outside（inside） the object.

% start level set evolution
bGo=1;
nTi=0;
while bGo
    u=EVOLUTION（u g lambda mu nu epsilon timestep 100）;
    nTi=nTi+1;
    sls（nTi+1::）=u;
    
    pause（0.001）;
    imshow（img[]）;hold on;
    [ch] = contour（u[0 0]‘m‘）;
    title（sprintf（‘Time Step: %d‘nTi*100））;
    hold off
    if ~strcmp（questdlg（‘Continue or not?‘）‘Yes‘）bGo=0;end
    pause（0.1）;
end

imgls=u;

imshow（img[]）;
hold on
imgt（::）=sls（1::）;
[ch] = contour（imgt[0 0]‘m‘）;
[ch] = contour（u[0 0]‘g‘‘linewidth‘2）; 
totalIterNum=[num2str（nTi*100） ‘ iterations‘];  
title（[‘Magenta: Initial; Green: Final after ‘ totalIterNum]）;
hold off

%% core functions of level set methods
function u = EVOLUTION（u0 g lambda mu nu epsilon delt numIter）
u=u0;
[vxvy]=gradient（g）;
 
for k=1:numIter
    u=NeumannBoundCond（u）;
    [uxuy]=gradient（u）; 
    normDu=sqrt（ux.^2 + uy.^2 + 1e-10）;
    Nx=ux./normDu;
    Ny=uy./normDu;
    diracU=Dirac（uepsilon）;
    K=curvature_central（NxNy）;
    weightedLengthTerm=lambda*diracU.*（vx.*Nx + vy.*Ny + g.*K）;
    penalizingTerm=mu*（4*del2（u）-K）;
    weightedAreaTerm=nu.*diracU.*g;
    u=u+delt*（weightedLengthTerm + weightedAreaTerm + penalizingTerm）; 
    % update the level set function
end

% the following functions are called by the main function EVOLUTION
function f = Dirac（x sigma）
f=（1/2/sigma）*（1+cos（pi*x/sigma））;
b = （x<=sigma） & （x>=-sigma）;
f = f.*b;

function K = curvature_central（nxny）
[nxxjunk]=gradient（nx）;  
[junknyy]=gradient（ny）;
K=nxx+nyy;

function g = NeumannBoundCond（f）
% Make a function satisfy Neumann boundary condition
[nrowncol] = size（f）;
g = f;
g（[1 nrow][1 ncol]） = g（[3 nrow-2][3 ncol-2]）;  
g（

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           0  2011-04-16 23:04  FCMLSM\
     文件     1150012  2011-04-16 01:35  FCMLSM\CBM2011.pdf
     文件       51080  2009-10-30 05:41  FCMLSM\ctlivertumor.bmp
     文件        3102  2009-12-04 15:16  FCMLSM\FCLSM.m
     文件        3564  2011-04-16 23:03  FCMLSM\fuzzyLSM.m
     文件       51078  2009-06-22 01:39  FCMLSM\mrihead.bmp
     文件         980  2011-04-16 23:04  FCMLSM\runtest.m
     文件        2455  2009-10-30 20:23  FCMLSM\SFCM2D.m
     文件        1332  2011-07-16 01:02  license.txt