core/vil/algo/vil_cartesian_differential_invariants.txx Source File

Go to the documentation of this file.
00001 #ifndef vil_cartesian_differential_invariants_txx_
00002 #define vil_cartesian_differential_invariants_txx_
00003 //:
00004 // \file
00005 // \brief Find Cartesian differential invariants
00006 // \author Ian Scott
00007 // Based on some Matlab code by K. Walker 1999.
00008 
00009 #include "vil_cartesian_differential_invariants.h"
00010 #include <vil/vil_image_view.h>
00011 #include <vil/algo/vil_gauss_filter.h>
00012 #include <vil/algo/vil_convolve_1d.h>
00013 #include <vil/vil_transpose.h>
00014 #include <vil/vil_plane.h>
00015 #include <vcl_cassert.h>
00016 
00017 //: Compute 1st, 2nd, and 3rd order C.d.i.s of an image.
00018 // The input must be 1 plane, the output will be 8 planes.
00019 template <class S, class T>
00020 inline void vil_cartesian_differential_invariants_3_1plane(
00021   const vil_image_view<S>& src, vil_image_view<T>& dest, double scale,
00022   unsigned max_kernel_width /*=0*/)
00023 {
00024   assert(src.nplanes()==1);
00025 
00026   unsigned filt_n=int(3*scale + 0.5)*2+1;
00027 
00028   if (max_kernel_width!=0 && filt_n > max_kernel_width)
00029     filt_n = (max_kernel_width | 1); // make sure it is an odd number
00030 
00031   const vcl_ptrdiff_t filt_c=filt_n/2;
00032 
00033   vcl_vector<double> filt_0(filt_n), filt_2(filt_n), filt_1(filt_n), filt_3(filt_n);
00034   vil_gauss_filter_gen_ntap(scale, 0, filt_0);
00035   vil_gauss_filter_gen_ntap(scale, 1, filt_1);
00036   vil_gauss_filter_gen_ntap(scale, 2, filt_2);
00037   vil_gauss_filter_gen_ntap(scale, 3, filt_3);
00038 
00039   const vil_convolve_boundary_option bo = vil_convolve_constant_extend;
00040 
00041   // Filter all the x directions
00042   vil_image_view<T> LX, LXx, LXxx, LXxxx;
00043   vil_convolve_1d(src, LX, &filt_0[filt_c], -filt_c, filt_c, double(), bo, bo);
00044   vil_convolve_1d(src, LXx, &filt_1[filt_c], -filt_c, filt_c, double(), bo, bo);
00045   vil_convolve_1d(src, LXxx, &filt_2[filt_c], -filt_c, filt_c, double(), bo, bo);
00046   vil_convolve_1d(src, LXxxx, &filt_3[filt_c], -filt_c, filt_c, double(), bo, bo);
00047 
00048   // Now calculate the full values.
00049   vil_image_view<T> Lx, Ly, Lxx, Lxy, Lyy, Lxxx, Lxxy, Lxyy, Lyyy;
00050 
00051   // construct first order partial derivatives
00052   vil_convolve_1d(vil_transpose(LXx), Lx,
00053                   &filt_0[filt_c], -filt_c, filt_c, double(), bo, bo);
00054   vil_convolve_1d(vil_transpose(LX), Ly,
00055                   &filt_1[filt_c], -filt_c, filt_c, double(), bo, bo);
00056 
00057   // construct second order partial derivatives
00058   vil_convolve_1d(vil_transpose(LXxx), Lxx,
00059                   &filt_0[filt_c], -filt_c, filt_c, double(), bo, bo);
00060   vil_convolve_1d(vil_transpose(LXx), Lxy,
00061                   &filt_1[filt_c], -filt_c, filt_c, double(), bo, bo);
00062   vil_convolve_1d(vil_transpose(LX), Lyy,
00063                   &filt_2[filt_c], -filt_c, filt_c, double(), bo, bo);
00064 
00065   // construct third order partial derivatives
00066   vil_convolve_1d(vil_transpose(LXxxx), Lxxx,
00067                   &filt_0[filt_c], -filt_c, filt_c, double(), bo, bo);
00068   vil_convolve_1d(vil_transpose(LXxx), Lxxy,
00069                   &filt_1[filt_c], -filt_c, filt_c, double(), bo, bo);
00070   vil_convolve_1d(vil_transpose(LXx), Lxyy,
00071                   &filt_2[filt_c], -filt_c, filt_c, double(), bo, bo);
00072   vil_convolve_1d(vil_transpose(LX), Lyyy,
00073                   &filt_3[filt_c], -filt_c, filt_c, double(), bo, bo);
00074 
00075   dest.set_size(src.ni(), src.nj(), 8);
00076   for (unsigned j=0;j<dest.ni();++j)
00077   {
00078     for (unsigned i=0;i<dest.nj();++i)
00079     {
00080       // I1 = LiLi=LxLx+LyLy
00081       // cdi(1,:,:)=(Lx.*Lx)+(Ly.*Ly);
00082       dest(j,i,0)= Lx(i,j)*Lx(i,j)
00083                   +Ly(i,j)*Ly(i,j);
00084 
00085       // LiLijLj=LxxLxLx+2LxLxyLy+LyyLyLy
00086       // cdi(2,:,:)=(Lxx.*Lx.*Lx)+(2.*Lx.*Lxy.*Ly)+(Lyy.*Ly.*Ly);
00087       dest(j,i,1)= Lxx(i,j)*Lx(i,j)*Lx(i,j)
00088                   +T(2)*Lx(i,j)*Lxy(i,j)*Ly(i,j)
00089                   +Lyy(i,j)*Ly(i,j)*Ly(i,j);
00090 
00091       // %LiiLjLj-LijLiLj=LxxLyLy-2LxyLxLy+LyyLxLx
00092       // cdi(3,:,:)=(Lxx.*Ly.*Ly)-(2.*Lxy.*Lx.*Ly)+(Lyy.*Lx.*Lx);
00093       dest(j,i,2)= Lxx(i,j)*Ly(i,j)*Ly(i,j)
00094                   -T(2)*Lxy(i,j)*Lx(i,j)*Ly(i,j)
00095                   +Lyy(i,j)*Lx(i,j)*Lx(i,j);
00096 
00097 
00098       // %-eijLjkLiLk=LxxLxLy+LxyLyLy-LyyLxLy-LxyLxLx
00099       // cdi(4,:,:)=(Lxx.*Lx.*Ly)+(Lxy.*Ly.*Ly)-(Lyy.*Lx.*Ly)-(Lxy.*Lx.*Lx);
00100       dest(j,i,3)= Lxx(i,j)*Lx(i,j)*Ly(i,j)
00101                   +Lxy(i,j)*Ly(i,j)*Ly(i,j)
00102                   -Lyy(i,j)*Lx(i,j)*Ly(i,j)
00103                   -Lxy(i,j)*Lx(i,j)*Lx(i,j);
00104 
00105       // %eij(LjklLiLkLl-LjkkLiLlLl)=(2LxyyLxLxLy)-(2LxxyLxLyLy)-(LxxyLxLyLy)-(LyyyLxLxLx)+(LxyyLyLxLx)+(LxxxLyLyLy)
00106       // cdi(5,:,:) = (2.*Lxyy.*Lx.*Lx.*Ly)-(2.*Lxxy.*Lx.*Ly.*Ly)-(Lxxy.*Lx.*Ly.*Ly)
00107       //             -(Lyyy.*Lx.*Lx.*Lx)+(Lxyy.*Ly.*Lx.*Lx)+(Lxxx.*Ly.*Ly.*Ly);
00108       dest(j,i,4)= T(2)*Lxyy(i,j)*Lx(i,j)*Lx(i,j)*Ly(i,j)
00109                   -T(2)*Lxxy(i,j)*Lx(i,j)*Ly(i,j)*Ly(i,j)
00110                   -Lxxy(i,j)*Lx(i,j)*Ly(i,j)*Ly(i,j)
00111                   -Lyyy(i,j)*Lx(i,j)*Lx(i,j)*Lx(i,j)
00112                   +Lxyy(i,j)*Ly(i,j)*Lx(i,j)*Lx(i,j)
00113                   +Lxxx(i,j)*Ly(i,j)*Ly(i,j)*Ly(i,j);
00114 
00115 
00116       // %LiijLjLkLk-LijkLiLjLk=(LxxxLxLyLy)-(2LxxyLxLxLy)+(LxxyLyLyLy)+(LxyyLxLxLx)-(2LxyyLxLyLy)+(LyyylxLxLy)
00117       // cdi(6,:,:) = (Lxxx.*Lx.*Ly.*Ly)-(2.*Lxxy.*Lx.*Lx.*Ly)+(Lxxy.*Ly.*Ly.*Ly)
00118       //             +(Lxyy.*Lx.*Lx.*Lx)-(2.*Lxyy.*Lx.*Ly.*Ly)+(Lyyy.*Lx.*Lx.*Ly);
00119       dest(j,i,5)= Lxxx(i,j)*Lx(i,j)*Ly(i,j)*Ly(i,j)
00120                   -T(2)*Lxxy(i,j)*Lx(i,j)*Lx(i,j)*Ly(i,j)
00121                   +Lxxy(i,j)*Ly(i,j)*Ly(i,j)*Ly(i,j)
00122                   +Lxyy(i,j)*Lx(i,j)*Lx(i,j)*Lx(i,j)
00123                   -T(2)*Lxyy(i,j)*Lx(i,j)*Ly(i,j)*Ly(i,j)
00124                   +Lyyy(i,j)*Lx(i,j)*Lx(i,j)*Ly(i,j);
00125 
00126 
00127       // %-eijLjklLiLkLl=(LxxxLxLxLy)+(2LxxyLxLyLy)+(LxyyLyLyLy)-(LyyyLxLyLy)-(2LxyyLxLxLy)-(LxxyLxLxLx)
00128       // cdi(7,:,:) = (Lxxx.*Lx.*Lx.*Ly)+(2.*Lxxy.*Lx.*Ly.*Ly)+(Lxyy.*Ly.*Ly.*Ly)
00129       //             -(Lyyy.*Lx.*Ly.*Ly)-(2.*Lxyy.*Lx.*Lx.*Ly)-(Lxxy.*Lx.*Lx.*Lx);
00130       dest(j,i,6)= Lxxx(i,j)*Lx(i,j)*Lx(i,j)*Ly(i,j)
00131                   +T(2)*Lxxy(i,j)*Lx(i,j)*Ly(i,j)*Ly(i,j)
00132                   +Lxyy(i,j)*Ly(i,j)*Ly(i,j)*Ly(i,j)
00133                   -Lyyy(i,j)*Lx(i,j)*Ly(i,j)*Ly(i,j)
00134                   -T(2)*Lxyy(i,j)*Lx(i,j)*Lx(i,j)*Ly(i,j)
00135                   -Lxxy(i,j)*Lx(i,j)*Lx(i,j)*Lx(i,j);
00136 
00137 
00138       // %LijkLiLjLk=(LxxxLxLxLx)+(LyyyLyLyLy)+(3LxxyLxLxLy)+(3LxyyLxLyLy)
00139       // cdi(8,:,:)=(Lxxx.*Lx.*Lx.*Lx)+(Lyyy.*Ly.*Ly.*Ly)+(3.*Lxxy.*Lx.*Lx.*Ly)+(3.*Lxyy.*Lx.*Ly.*Ly);
00140       dest(j,i,7)= Lxxx(i,j)*Lx(i,j)*Lx(i,j)*Lx(i,j)
00141                   +Lyyy(i,j)*Ly(i,j)*Ly(i,j)*Ly(i,j)
00142                   +T(3)*Lxxy(i,j)*Lx(i,j)*Lx(i,j)*Ly(i,j)
00143                   +T(3)*Lxyy(i,j)*Lx(i,j)*Ly(i,j)*Ly(i,j);
00144     }
00145   }
00146 }
00147 
00148 
00149 //: Compute 1st, 2nd, and 3rd order C.d.i.s of an image.
00150 template <class S, class T>
00151 void vil_cartesian_differential_invariants_3(
00152   const vil_image_view<S>& src, vil_image_view<T>& dest, double scale, unsigned max_kernel_width /*=0*/)
00153 {
00154   dest.set_size(src.ni(), src.nj(), src.nplanes()*8);
00155   assert(dest.planestep() >= int(dest.ni() * dest.nj()));
00156 
00157   for (unsigned p=0; p < src.nplanes(); ++p)
00158   {
00159     vil_image_view<T> destplanes = vil_planes(dest, p*8, 1, 8);
00160 #ifndef NDEBUG
00161     vil_image_view<T> check = destplanes;
00162 #endif
00163     vil_cartesian_differential_invariants_3_1plane(vil_plane(src, p), destplanes, scale, max_kernel_width);
00164     assert(destplanes == check);
00165   }
00166 }
00167 
00168 #undef VIL_CARTESIAN_DIFFERENTIAL_INVARIANTS_INSTANTIATE
00169 #define VIL_CARTESIAN_DIFFERENTIAL_INVARIANTS_INSTANTIATE(S, T) \
00170 template void \
00171 vil_cartesian_differential_invariants_3( const vil_image_view< S >& src, \
00172                                          vil_image_view< T >& dest, double, unsigned )
00173 
00174 #endif // vil_cartesian_differential_invariants_txx_