1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
| %-----------------------------------------------------------------------------------------------
% Loccally One-dimensional (LOD) implicit scheme for closed geometric
% active contour model+minimal variation+GAC+robust alignement
% I=Input image matrix, Phi= Initial contour matrix (implicit form)
% Balloon= Weighted area factor (scalar)
% Align = Alignement force factro (scalar)
% Max_Lloyd= The Max-Lloyd/ Chan-Vese threshold factor (scalar)
% k= the time step (tau scalar)
% iter= Maximal number of iterations (scalar)
function Phi=LOD_Active_Contour(I,Phi,Balloon,Align,Max_Lloyd,k,iter)
D2I=Dxx(I)+Dyy(I); P=Dx(I); Q=Dy(I);
g=1./sqrt(1+P.^2+Q.^2); % example for computing g(x,y)
delta=2; count=1;
while and(delta>0.0001,count<iter) % check if converged
Phi_old=Phi;
threshold=LloydMax(Phi,I); % Max-Lloyd/ Chan-Vese term
alignement=-sign(P.*Dx(Phi)+Q.*Dy(Phi)).*D2I; % Laplacien term
Phi=Phi+k*(Balloon*g+Align*alignement+Max_Lloyd*thershold);
for i=1:2, % i=1 => (I-tau*Ax) i=2 => (I-tau*Ay)
Phi=Implicit(g(:),k,Phi); % (1/(I-tau*Ai))Phi
Phi=Phi';g=g'; % Transpose for Ay
end % for i
Phi=redistance(Phi); % Use fast marching for re-distancing
delta=sum(sum((Phi_old-Phi).^2)) % Compute L2 norm
count=count+1;
imshow(I,[]); hold on; contour(Phi,[00],'r'); hold off; drawnow;
end % while and funtion
%-----------------------------------------------------------------------------------------
% Compute (1/(I-k*Al)) Phi using Thomas algorithm
% k=time step, g=g(x,y) in column stack form
function u=Implicit(g,k,Phi)
gm= -k.*(g+g([end 1:end-1]))/2; % lower diag
gc= 1-k.*(-2*g-g([end 1:end-1])-g([2:end 1]))/2; % main diag
gp= -k.*(g+g([2:end 1]))/2; % upper diag
u= Thomas(gc,gp(1:end-1),gm(2:end),Phi(:));
u= reshape(u,size(Phi));
%-------------------------------------------------------------------------------------------
% Compute the Lloyd-Max/ Chan-Vese thersholding
function force=LloydMax(Phi,I)
mask_in=(Phi<0); inside the contour
mask_out= 1-mask_in; % rest of the domain
I_in=sum(sum(mask_in.*I))/sum(mask_in(:)); % mean value
I_out=sum(sum(mask_out.*I))/sum(mask_out(:));
force=(I_out-I_in).*(I-(I_in+I_out)/2);
%---------------------------------------------------------------------------------------------
% 'Roughly' correct Phi to be a distance map of its zero set
function u=redistance(Phi);
u=(sign(Phi)+1)*999999; % set to infinity all positive
for i=1:2,
l2=2;
if i>1 u=(1-sign(Phi))*999999; end % set to infinity all negative
while l2>1,
v=Update(u,1);
l2=sum(sum((u-v).^2));
u=v;;
end % while
if i>1 u=u-up; else up=u; end % if
end % for
%-----------------------------------------------------------------------------------------------
%Solve |grad u|=F(x,y) parallel version of the FMM
function res=Update(u,F)
mx=min(u([2:end end],:), u([1 1:end-1],:));
my= min(u(:,[2:end end]), u(:,[1 1:end-1]));
delm=(mx-my);
mask=(abs(delm)<F);
res=min(mask.*(mx+my+sqrt(2.*F.^2-delm.^2))./2+...
(1-mask).*(min(mx,my)+F),u);
%-----------------------------------------------------------------------------------------------
function f=Dmx(P)
f=P-P([1 1:end-1],:);
%-----------------------------------------------------------------------------------------------
function f=Dpx(P)
f=P([2:end end],:)-P;
%-----------------------------------------------------------------------------------------------
function f=Dx(P)
f=(Dpx(P)+Dmx(P))/2;
%-----------------------------------------------------------------------------------------------
function f=Dy(P)
f=(Dx(P'))';
%-----------------------------------------------------------------------------------------------
function f=Dxx(P)
f=Dpx(P)-Dmx(P);
%-----------------------------------------------------------------------------------------------
function f=Dyy(P)
f=(Dxx(P'))';
%-----------------------------------------------------------------------------------------------
% Thomas Algorithm for trilinear diagonally dominant system:
% B u=d(solve for u); B is given by its 3 diagonals:
% alpha(1:N)=main diagonal, beta(1:N-1)=upper diagonal,
% Gamma(1:N-1)=lower diagonal, (Compile first!)
function u=Thomas(alpha,beta,gamma,d)
N=length(alpha); r=beta;
l=zeros(N-1,1); u=zeros(N,1);
m=u; %zero
m(1)=alpha(1);
for i=1:N-1,
l(i)=gamma(i)/m(i);
m(i+1)=alpha(i+1)-l(i)*beta(i);
end % for
y=u; %zero
y(1)=d(1);
for i=2:N,
y(i)=d(i)-l(i-1)*y(i-1);
end % for
u(N)=y(N)/m(N);
for i=N-1:-1:1,
u(i)=(y(i)-beta(i)*u(i+1))/m(i);
end % for |
