-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathSnake3D.m
151 lines (139 loc) · 5.48 KB
/
Snake3D.m
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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
function FV=Snake3D(I,FV,Options)
% This function SNAKE implements the basic snake segmentation. A snake is an
% active (moving) contour, in which the points are attracted by edges and
% other boundaries. To keep the contour smooth, an membrame and thin plate
% energy is used as regularization.
%
% OV=Snake3D(I,FV,Options)
%
% inputs,
% I : An Image of type double preferable ranged [0..1]
% FV : Structure with triangulated mesh, with list of faces FV.faces N x 3
% and list of vertices M x 3
% Options : A struct with all snake options
%
% outputs,
% OV : Structure with triangulated mesh of the final surface
%
% options (general),
% Option.Verbose : If true show important images, default false
% Options.Gamma : Time step, default 1
% Options.Iterations : Number of iterations, default 100
%
% options (Image Edge Energy / Image force))
% Options.Sigma1 : Sigma used to calculate image derivatives, default 2
% Options.Wline : Attraction to lines, if negative to black lines otherwise white
% lines , default 0.04
% Options.Wedge : Attraction to edges, default 2.0
% Options.Sigma2 : Sigma used to calculate the gradient of the edge energy
% image (which gives the image force), default 2
%
% options (Gradient Vector Flow)
% Options.Mu : Trade of between real edge vectors, and noise vectors,
% default 0.2. (Warning setting this to high >0.5 gives
% an instable Vector Flow)
% Options.GIterations : Number of GVF iterations, default 0
% Options.Sigma3 : Sigma used to calculate the laplacian in GVF, default 1.0
%
% options (Snake)
% Options.Alpha : Membrame energy (first order), default 0.2
% Options.Beta : Thin plate energy (second order), default 0.2
% Options.Delta : Baloon force, default 0.1
% Options.Kappa : Weight of external image force, default 2
% Options.Lambda : Weight which changes the direction of the image
% potential force to the direction of the surface
% normal, value default 0.8 (range [0..1])
% (Keeps the surface from self intersecting)
%
% Literature:
% - Michael Kass, Andrew Witkin and Demetri TerzoPoulos "Snakes : Active
% Contour Models", 1987
% - Christoph Lurig, Leif Kobbelt, Thomas Ertl, "Hierachical solutions
% for the Deformable Surface Problem in Visualization"
%
% Example,
%
% load testvolume
% load SphereMesh
% Options=struct;
% Options.Verbose=1;
% Options.Wedge=0;
% Options.Wline=-1;
% Options.Alpha=0.2;
% Options.Beta=0.2;
% Options.Kappa=0.5;
% Options.Delta=0.1000;
% Options.Gamma=0.1000;
% Options.Iterations=20;
% Options.Sigma1=2;
% Options.Sigma2=2;
% Options.Lambda=0.8;
% FV.vertices(:,1)=FV.vertices(:,1)+35;
% FV.vertices(:,2)=FV.vertices(:,2)+25;
% FV.vertices(:,3)=FV.vertices(:,3)+20;
% OV=Snake3D(I,FV,Options)
%
% Function is written by D.Kroon University of Twente (July 2010)
% Process inputs
defaultoptions=struct('Verbose',false,'Wline',0.04,'Wedge',2,'Sigma1',2,'Sigma2',2,'Alpha',0.2,'Beta',0.2,'Delta',0.1,'Gamma',1,'Kappa',2,'Iterations',100,'GIterations',0,'Mu',0.2,'Sigma3',1,'Lambda',0.8);
if(~exist('Options','var')),
Options=defaultoptions;
else
tags = fieldnames(defaultoptions);
for i=1:length(tags)
if(~isfield(Options,tags{i})), Options.(tags{i})=defaultoptions.(tags{i}); end
end
if(length(tags)~=length(fieldnames(Options))),
warning('snake:unknownoption','unknown options found');
end
end
% Convert input to single if xintxx
if(~strcmpi(class(I),'single')&&~strcmpi(class(I),'double'))
I=single(I);
end
% The surface faces must always be clockwise (because of the balloon force)
FV=MakeContourClockwise3D(FV);
% Transform the Image into an External Energy Image
Eext = ExternalForceImage3D(I,Options.Wline, Options.Wedge,Options.Sigma1);
% Make the external force (flow) field.
Fx=ImageDerivatives3D(Eext,Options.Sigma2,'x');
Fy=ImageDerivatives3D(Eext,Options.Sigma2,'y');
Fz=ImageDerivatives3D(Eext,Options.Sigma2,'z');
Fext(:,:,:,1)=-Fx*2*Options.Sigma2^2;
Fext(:,:,:,2)=-Fy*2*Options.Sigma2^2;
Fext(:,:,:,3)=-Fz*2*Options.Sigma2^2;
% Do Gradient vector flow, optimalization
Fext=GVFOptimizeImageForces3D(Fext, Options.Mu, Options.GIterations, Options.Sigma3);
% Show the image, contour and force field
if(Options.Verbose)
drawnow; pause(0.1);
h=figure; set(h,'render','opengl')
subplot(2,3,1),imshow(squeeze(Eext(:,:,round(end/2))),[]);
subplot(2,3,2),imshow(squeeze(Eext(:,round(end/2),:)),[]);
subplot(2,3,3),imshow(squeeze(Eext(round(end/2),:,:)),[]);
subplot(2,3,4),imshow(squeeze(Fext(:,:,round(end/2),:))+0.5);
subplot(2,3,5),imshow(squeeze(Fext(:,round(end/2),:,:))+0.5);
subplot(2,3,6),imshow(squeeze(Fext(round(end/2),:,:,:))+0.5);
h=figure; set(h,'render','opengl'); hold on;
%patch(i,'facecolor',[1 0 0],'facealpha',0.1);
ind=find(I(:)>0);
[ix,iy,iz]=ind2sub(size(Eext),ind);
plot3(ix,iy,iz,'b.');
hold on;
h=patch(FV,'facecolor',[1 0 0],'facealpha',0.1);
drawnow; pause(0.1);
end
% Make the interal force matrix, which constrains the moving points to a
% smooth contour
S=SnakeInternalForceMatrix3D(FV,Options.Alpha,Options.Beta,Options.Gamma);
for i=1:Options.Iterations
FV=SnakeMoveIteration3D(S,FV,Fext,Options.Gamma,Options.Kappa,Options.Delta,Options.Lambda);
% Show current contour
if(Options.Verbose)
if(ishandle(h));
delete(h);
h=patch(FV,'facecolor',[1 0 0],'facealpha',0.1);
drawnow; %pause(0.1);
end
end
end