]> AND Private Git Repository - these_gilles.git/blob - THESE/codes/snake/basic_code/ImageDerivatives3D.m
Logo AND Algorithmique Numérique Distribuée

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