Skip to content

Commit 024bd9d

Browse files
Joey LiJoey Li
Joey Li
authored and
Joey Li
committedMay 6, 2017
添加：原始方法
1 parent b7db5d5 commit 024bd9d

28 files changed

+1568
-150
lines changed
 

‎DrawSegmentedArea2D.m

+48
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
function J=DrawSegmentedArea2D(P,Isize)
2+
% Draw the contour as one closed line in a logical image,
3+
% and make the inside of the contour true with imfill
4+
%
5+
% J=DrawSegmentedArea2D(P,Isize)
6+
%
7+
% inputs,
8+
% P : The list with contour points 2 x N
9+
% Isize : The size of the output image [x y]
10+
%
11+
% outputs,
12+
% J : The binary image with the contour filled
13+
%
14+
% example:
15+
% y=[182 233 251 205 169];
16+
% x=[163 166 207 248 210];
17+
% P=[x(:) y(:)];
18+
%
19+
% J=DrawSegmentedArea(P,[400 400]);
20+
% figure, imshow(J);
21+
% hold on; plot([P(:,2);P(1,2)],[P(:,1);P(1,1)]);
22+
%
23+
% Function is written by D.Kroon University of Twente (July 2010)
24+
25+
26+
J=false(Isize+2);
27+
% Loop through all line coordinates
28+
x=round([P(:,1);P(1,1)]); x=min(max(x,1),Isize(1));
29+
y=round([P(:,2);P(1,2)]); y=min(max(y,1),Isize(2));
30+
for i=1:(length(x)-1)
31+
% Calculate the pixels needed to construct a line of 1 pixel thickness
32+
% between two coordinates.
33+
xp=[x(i) x(i+1)]; yp=[y(i) y(i+1)];
34+
dx=abs(xp(2)-xp(1)); dy=abs(yp(2)-yp(1));
35+
if(dx==dy)
36+
if(xp(2)>xp(1)), xline=xp(1):xp(2); else xline=xp(1):-1:xp(2); end
37+
if(yp(2)>yp(1)), yline=yp(1):yp(2); else yline=yp(1):-1:yp(2); end
38+
elseif(dx>dy)
39+
if(xp(2)>xp(1)), xline=xp(1):xp(2); else xline=xp(1):-1:xp(2); end
40+
yline=linspace(yp(1),yp(2),length(xline));
41+
else
42+
if(yp(2)>yp(1)), yline=yp(1):yp(2); else yline=yp(1):-1:yp(2); end
43+
xline=linspace(xp(1),xp(2),length(yline));
44+
end
45+
% Insert all pixels in the fill image
46+
J(round(xline+1)+(round(yline+1)-1)*size(J,1))=1;
47+
end
48+
J=bwfill(J,1,1); J=~J(2:end-1,2:end-1);

‎ExternalForceImage2D.m

+29
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
function Eextern = ExternalForceImage2D(I,Wline, Wedge, Wterm,Sigma)
2+
% Eextern = ExternalForceImage2D(I,Wline, Wedge, Wterm,Sigma)
3+
%
4+
% inputs,
5+
% I : The image
6+
% Sigma : Sigma used to calculated image derivatives
7+
% Wline : Attraction to lines, if negative to black lines otherwise white
8+
% lines
9+
% Wedge : Attraction to edges
10+
% Wterm : Attraction to terminations of lines (end points) and corners
11+
%
12+
% outputs,
13+
% Eextern : The energy function described by the image
14+
%
15+
% Function is written by D.Kroon University of Twente (July 2010)
16+
17+
Ix=ImageDerivatives2D(I,Sigma,'x');
18+
Iy=ImageDerivatives2D(I,Sigma,'y');
19+
Ixx=ImageDerivatives2D(I,Sigma,'xx');
20+
Ixy=ImageDerivatives2D(I,Sigma,'xy');
21+
Iyy=ImageDerivatives2D(I,Sigma,'yy');
22+
23+
24+
Eline = imgaussian(I,Sigma);
25+
Eterm = (Iyy.*Ix.^2 -2*Ixy.*Ix.*Iy + Ixx.*Iy.^2)./((1+Ix.^2 + Iy.^2).^(3/2));
26+
Eedge = sqrt(Ix.^2 + Iy.^2);
27+
28+
Eextern= (Wline*Eline - Wedge*Eedge -Wterm * Eterm);
29+

‎ExternalForceImage3D.m

+24
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
function Eextern = ExternalForceImage3D(I,Wline, Wedge,Sigma)
2+
% Eextern = ExternalForceImage3D(I,Wline, Wedge,Sigma)
3+
%
4+
% inputs,
5+
% I : The image
6+
% Sigma : Sigma used to calculated image derivatives
7+
% Wline : Attraction to lines, if negative to black lines otherwise white
8+
% lines
9+
% Wterm : Attraction to terminations of lines (end points) and corners
10+
%
11+
% outputs,
12+
% Eextern : The energy function described by the image
13+
%
14+
% Function is written by D.Kroon University of Twente (July 2010)
15+
16+
Ix=ImageDerivatives3D(I,Sigma,'x');
17+
Iy=ImageDerivatives3D(I,Sigma,'y');
18+
Iz=ImageDerivatives3D(I,Sigma,'z');
19+
20+
Eline = imgaussian(I,Sigma);
21+
Eedge = sqrt(Ix.^2 + Iy.^2 + Iz.^2);
22+
23+
Eextern= (Wline*Eline - Wedge*Eedge);
24+

‎GVFOptimizeImageForces2D.m

+44
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
function Fext=GVFOptimizeImageForces2D(Fext, Mu, Iterations, Sigma)
2+
% This function "GVFOptimizeImageForces" does gradient vector flow (GVF)
3+
% on a vector field. GVF gives the edge-forces a larger capature range,
4+
% to make the snake also reach concave regions
5+
%
6+
% Fext = GVFOptimizeImageForces2D(Fext, Mu, Iterations, Sigma)
7+
%
8+
% inputs,
9+
% Fext : The image force vector field N x M x 2
10+
% Mu : Is a trade of scalar between noise and real edge forces
11+
% Iterations : The number of GVF itterations
12+
% Sigma : Used when calculating the Laplacian
13+
%
14+
% outputs,
15+
% Fext : The GVF optimized image force vector field
16+
%
17+
% Function is written by D.Kroon University of Twente (July 2010)
18+
19+
% Squared magnitude of force field
20+
Fx= Fext(:,:,1);
21+
Fy= Fext(:,:,2);
22+
23+
% Calculate magnitude
24+
sMag = Fx.^2+ Fy.^2;
25+
26+
% Set new vector-field to initial field
27+
u=Fx; v=Fy;
28+
29+
% Iteratively perform the Gradient Vector Flow (GVF)
30+
for i=1:Iterations,
31+
% Calculate Laplacian
32+
Uxx=ImageDerivatives2D(u,Sigma,'xx');
33+
Uyy=ImageDerivatives2D(u,Sigma,'yy');
34+
35+
Vxx=ImageDerivatives2D(v,Sigma,'xx');
36+
Vyy=ImageDerivatives2D(v,Sigma,'yy');
37+
38+
% Update the vector field
39+
u = u + Mu*(Uxx+Uyy) - sMag.*(u-Fx);
40+
v = v + Mu*(Vxx+Vyy) - sMag.*(v-Fy);
41+
end
42+
43+
Fext(:,:,1) = u;
44+
Fext(:,:,2) = v;

‎GVFOptimizeImageForces3D.m

+58
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
function Fext=GVFOptimizeImageForces3D(Fext, Mu, Iterations, Sigma)
2+
% This function "GVFOptimizeImageForces" does gradient vector flow (GVF)
3+
% on a vector field. GVF gives the edge-forces a larger capature range,
4+
% to make the snake also reach concave regions
5+
%
6+
% Fext = GVFOptimizeImageForces3D(Fext, Mu, Iterations, Sigma)
7+
%
8+
% inputs,
9+
% Fext : The image force vector field N x M x 3
10+
% Mu : Is a trade of scalar between noise and real edge forces
11+
% Iterations : The number of GVF itterations
12+
% Sigma : Used when calculating the Laplacian
13+
%
14+
% outputs,
15+
% Fext : The GVF optimized image force vector field
16+
%
17+
% Function is written by D.Kroon University of Twente (July 2010)
18+
19+
% Squared magnitude of force field
20+
Fx= Fext(:,:,:,1);
21+
Fy= Fext(:,:,:,2);
22+
Fz= Fext(:,:,:,3);
23+
24+
% Calculate magnitude
25+
sMag = Fx.^2+ Fy.^2 + Fz.^2;
26+
27+
% Set new vector-field to initial field
28+
u=Fx; v=Fy; w=Fz;
29+
30+
% Iteratively perform the Gradient Vector Flow (GVF)
31+
for i=1:Iterations,
32+
% Calculate Laplacian
33+
Uxx=ImageDerivatives3D(u,Sigma,'xx');
34+
Uyy=ImageDerivatives3D(u,Sigma,'yy');
35+
Uzz=ImageDerivatives3D(u,Sigma,'zz');
36+
37+
% Update the vector field
38+
u = u + Mu*(Uxx+Uyy+Uzz) - sMag.*(u-Fx);
39+
clear('Uxx','Uyy','Uzz');
40+
41+
Vxx=ImageDerivatives3D(v,Sigma,'xx');
42+
Vyy=ImageDerivatives3D(v,Sigma,'yy');
43+
Vzz=ImageDerivatives3D(v,Sigma,'zz');
44+
45+
v = v + Mu*(Vxx+Vyy+Vzz) - sMag.*(v-Fy);
46+
clear('Vxx','Vyy','Vzz');
47+
48+
Wxx=ImageDerivatives3D(w,Sigma,'xx');
49+
Wyy=ImageDerivatives3D(w,Sigma,'yy');
50+
Wzz=ImageDerivatives3D(w,Sigma,'zz');
51+
52+
w = w + Mu*(Wxx+Wyy+Wzz) - sMag.*(w-Fz);
53+
clear('Wxx','Wyy','Wzz');
54+
end
55+
56+
Fext(:,:,:,1) = u;
57+
Fext(:,:,:,2) = v;
58+
Fext(:,:,:,3) = w;

‎GetContourNormals2D.m

+33
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
1+
function N=GetContourNormals2D(P)
2+
% This function calculates the normals, of the contour points
3+
% using the neighbouring points of each contour point
4+
%
5+
% N=GetContourNormals2D(P)
6+
%
7+
% inputs,
8+
% P : List with contour coordinates M x 2
9+
%
10+
% outputs,
11+
% N : List with contour normals M x 2
12+
%
13+
% Function is written by D.Kroon University of Twente (July 2010)
14+
15+
% Use the n'th neighbour to calculate the normal (more stable)
16+
a=4;
17+
18+
% From array to separate x,y
19+
xt=P(:,1); yt=P(:,2);
20+
21+
% Derivatives of contour
22+
n=length(xt);
23+
f=(1:n)+a; f(f>n)=f(f>n)-n;
24+
b=(1:n)-a; b(b<1)=b(b<1)+n;
25+
26+
dx=xt(f)-xt(b);
27+
dy=yt(f)-yt(b);
28+
29+
% Normals of contourpoints
30+
l=sqrt(dx.^2+dy.^2);
31+
nx = -dy./l;
32+
ny = dx./l;
33+
N(:,1)=nx; N(:,2)=ny;

‎ImageDerivatives2D.m

+32
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
function J=ImageDerivatives2D(I,sigma,type)
2+
% Gaussian based image derivatives
3+
%
4+
% J=ImageDerivatives2D(I,sigma,type)
5+
%
6+
% inputs,
7+
% I : The 2D image
8+
% sigma : Gaussian Sigma
9+
% type : 'x', 'y', 'xx', 'xy', 'yy'
10+
%
11+
% outputs,
12+
% J : The image derivative
13+
%
14+
% Function is written by D.Kroon University of Twente (July 2010)
15+
16+
% Make derivatives kernels
17+
[x,y]=ndgrid(floor(-3*sigma):ceil(3*sigma),floor(-3*sigma):ceil(3*sigma));
18+
19+
switch(type)
20+
case 'x'
21+
DGauss=-(x./(2*pi*sigma^4)).*exp(-(x.^2+y.^2)/(2*sigma^2));
22+
case 'y'
23+
DGauss=-(y./(2*pi*sigma^4)).*exp(-(x.^2+y.^2)/(2*sigma^2));
24+
case 'xx'
25+
DGauss = 1/(2*pi*sigma^4) * (x.^2/sigma^2 - 1) .* exp(-(x.^2 + y.^2)/(2*sigma^2));
26+
case {'xy','yx'}
27+
DGauss = 1/(2*pi*sigma^6) * (x .* y) .* exp(-(x.^2 + y.^2)/(2*sigma^2));
28+
case 'yy'
29+
DGauss = 1/(2*pi*sigma^4) * (y.^2/sigma^2 - 1) .* exp(-(x.^2 + y.^2)/(2*sigma^2));
30+
end
31+
32+
J = imfilter(I,DGauss,'conv','symmetric');

‎ImageDerivatives3D.m

+43
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,43 @@
1+
function J=ImageDerivatives3D(I,sigma,type)
2+
% Gaussian based image derivatives
3+
%
4+
% J=ImageDerivatives3D(I,sigma,type)
5+
%
6+
% inputs,
7+
% I : The 3D image
8+
% sigma : Gaussian Sigma
9+
% type : 'x', 'y', 'z', 'xx', 'yy', 'zz', 'xy', 'xz', 'yz'
10+
%
11+
% outputs,
12+
% J : The image derivative
13+
%
14+
% Function is written by D.Kroon University of Twente (July 2010)
15+
16+
% Make derivatives kernels
17+
[x,y,z]=ndgrid(floor(-3*sigma):ceil(3*sigma),floor(-3*sigma):ceil(3*sigma),floor(-3*sigma):ceil(3*sigma));
18+
19+
switch(type)
20+
case 'x'
21+
DGauss=-(x./((2*pi)^(3/2)*sigma^5)).*exp(-(x.^2+y.^2+z.^2)/(2*sigma^2));
22+
case 'y'
23+
DGauss=-(y./((2*pi)^(3/2)*sigma^5)).*exp(-(x.^2+y.^2+z.^2)/(2*sigma^2));
24+
case 'z'
25+
DGauss=-(z./((2*pi)^(3/2)*sigma^5)).*exp(-(x.^2+y.^2+z.^2)/(2*sigma^2));
26+
case 'xx'
27+
DGauss = 1/((2*pi)^(3/2)*sigma^5) * (x.^2/sigma^2 - 1) .* exp(-(x.^2 + y.^2 + z.^2)/(2*sigma^2));
28+
case 'yy'
29+
DGauss = 1/((2*pi)^(3/2)*sigma^5) * (y.^2/sigma^2 - 1) .* exp(-(x.^2 + y.^2 + z.^2)/(2*sigma^2));
30+
case 'zz'
31+
DGauss = 1/((2*pi)^(3/2)*sigma^5) * (z.^2/sigma^2 - 1) .* exp(-(x.^2 + y.^2 + z.^2)/(2*sigma^2));
32+
case {'xy','yx'}
33+
DGauss = 1/((2*pi)^(3/2)*sigma^7) * (x .* y) .* exp(-(x.^2 + y.^2 + z.^2)/(2*sigma^2));
34+
case {'xz','zx'}
35+
DGauss = 1/((2*pi)^(3/2)*sigma^7) * (x .* z) .* exp(-(x.^2 + y.^2 + z.^2)/(2*sigma^2));
36+
case {'yz','zy'}
37+
DGauss = 1/((2*pi)^(3/2)*sigma^7) * (y .* z) .* exp(-(x.^2 + y.^2 + z.^2)/(2*sigma^2));
38+
end
39+
40+
K=SeparateKernel(DGauss);
41+
J = imfilter(I,K{1},'conv','symmetric');
42+
J = imfilter(J,K{2},'conv','symmetric');
43+
J = imfilter(J,K{3},'conv','symmetric');

‎InterpolateContourPoints2D.m

+45
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
function K=InterpolateContourPoints2D(P,nPoints)
2+
% This function resamples a few points describing a countour , to a smooth
3+
% contour of uniform sampled points.
4+
%
5+
% K=InterpolateContourPoints(P,nPoints)
6+
%
7+
% input,
8+
% P : Inpute Contour, size N x 2 (with N>=4)
9+
% nPoints : Number of Contour points as output
10+
%
11+
% output,
12+
% K : Uniform sampled Contour points, size nPoints x 2
13+
%
14+
% Function is written by D.Kroon University of Twente (July 2010)
15+
%
16+
% example,
17+
% % Show an image
18+
% figure, imshow(imread('moon.tif'));
19+
% % Select some points with the mouse
20+
% [y,x] = getpts;
21+
% % Make an array with the clicked coordinates
22+
% P=[x(:) y(:)];
23+
% % Interpolate inbetween the points
24+
% Pnew=InterpolateContourPoints2D(P,100)
25+
% % Show the result
26+
% hold on; plot(P(:,2),P(:,1),'b*');
27+
% plot(Pnew(:,2),Pnew(:,1),'r.');
28+
%
29+
% Function is written by D.Kroon University of Twente (July 2010)
30+
31+
% Interpolate points inbetween
32+
O(:,1)=interp([P(end-3:end,1);P(:,1);P(:,1);P(1:4,1)],10);
33+
O(:,2)=interp([P(end-3:end,2);P(:,2);P(:,2);P(1:4,2)],10);
34+
O=O(41:end-39,:);
35+
36+
% Calculate distance between points
37+
dis=[0;cumsum(sqrt(sum((O(2:end,:)-O(1:end-1,:)).^2,2)))];
38+
39+
% Resample to make uniform points
40+
K(:,1) = interp1(dis,O(:,1),linspace(0,dis(end),nPoints*2));
41+
K(:,2) = interp1(dis,O(:,2),linspace(0,dis(end),nPoints*2));
42+
K=K(round(end/4):round(end/4)+nPoints-1,:);
43+
44+
45+

‎MakeContourClockwise2D.m

+16
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
function P=MakeContourClockwise2D(P)
2+
% This function MakeContourClockwise will make a contour clockwise
3+
% contour clockwise. This is done by calculating the area inside the
4+
% contour, if it is positive we change the contour orientation.
5+
%
6+
% P=MakeContourClockwise2D(P);
7+
%
8+
% Function is written by D.Kroon University of Twente (July 2010)
9+
10+
% Area inside contour
11+
O=[P;P(1:2,:)];
12+
area = 0.5*sum((O((1:size(P,1))+1,1) .* (O((1:size(P,1))+2,2) - O((1:size(P,1)),2))));
13+
14+
% If the area inside the contour is positive, change from counter-clockwise to
15+
% clockwise
16+
if(area>0), P=P(end:-1:1,:); end

0 commit comments

Comments
 (0)
Please sign in to comment.