1 // This is mul/vil3d/algo/vil3d_world_gradients.hxx
2 #ifndef vil3d_world_gradients_hxx_
3 #define vil3d_world_gradients_hxx_
4 //:
5 // \file
6 // \brief Given image gradients compute world gradients and gradient magnitude
7 // \author Tim Cootes
8
9 #include <iostream>
10 #include <cmath>
11 #include "vil3d_world_gradients.h"
12 #include <vil3d/algo/vil3d_fill_border.h>
13 #include <vil3d/vil3d_transform.h>
14 #include <vil3d/vil3d_plane.h>
15 #ifdef _MSC_VER
16 # include <vcl_msvc_warnings.h>
17 #endif
18 #include <cassert>
19
20 //: Functor class to scale by s
21 class vil3d_math_scale_functor
22 {
23 private:
24 double s_;
25 public:
vil3d_math_scale_functor(double s)26 vil3d_math_scale_functor(double s) : s_(s) {}
operator ()(vxl_byte x) const27 float operator()(vxl_byte x) const { return float(s_*x); }
operator ()(unsigned x) const28 float operator()(unsigned x) const { return float(s_*x); }
operator ()(short x) const29 float operator()(short x) const { return float(s_*x); }
operator ()(int x) const30 float operator()(int x) const { return float(s_*x); }
operator ()(float x) const31 float operator()(float x) const { return float(s_*x); }
operator ()(double x) const32 float operator()(double x) const { return float(s_*x); }
33 };
34
35 //: Given image gradients compute world gradients and gradient magnitude
36 // Input gradient images are assumed to be un-normalised pixel gradients
37 // (ie no scaling has been done to take account of world pixel widths).
38 // Divides each by corresponding pixel dimension to give gradient in world units
39 // (ie intensity change per unit world length) in world_grad (3 plane image)
40 // The gradient magnitude output is in units of intensity change per world length
41 // (ie it does take account of voxel sizes).
42 //
43 // Note: Currently assumes single plane only.
44 // 1 pixel border around output set to zero.
45 //
46 // \relatesalso vil3d_image_view
47 template<class srcT, class destT>
vil3d_world_gradients(const vil3d_image_view<srcT> & grad_i,const vil3d_image_view<srcT> & grad_j,const vil3d_image_view<srcT> & grad_k,double voxel_width_i,double voxel_width_j,double voxel_width_k,vil3d_image_view<destT> & world_grad,vil3d_image_view<destT> & grad_mag)48 void vil3d_world_gradients(const vil3d_image_view<srcT>& grad_i,
49 const vil3d_image_view<srcT>& grad_j,
50 const vil3d_image_view<srcT>& grad_k,
51 double voxel_width_i,
52 double voxel_width_j,
53 double voxel_width_k,
54 vil3d_image_view<destT>& world_grad,
55 vil3d_image_view<destT>& grad_mag)
56 {
57 assert(grad_i.nplanes()==grad_j.nplanes());
58 assert(grad_i.nplanes()==grad_k.nplanes());
59 assert(grad_i.nplanes()==1);
60 unsigned ni = grad_i.ni(), nj = grad_i.nj(), nk = grad_i.nk();
61 assert(ni>2 && nj>2 && nk>2);
62 assert(grad_j.ni()==ni && grad_j.nj()==nj && grad_j.nk()==nk);
63 assert(grad_k.ni()==ni && grad_k.nj()==nj && grad_k.nk()==nk);
64 world_grad.set_size(ni,nj,nk,3);
65 grad_mag.set_size(ni,nj,nk,1);
66
67 vil3d_image_view<destT> w_grad_i=vil3d_plane(world_grad,0);
68 vil3d_image_view<destT> w_grad_j=vil3d_plane(world_grad,1);
69 vil3d_image_view<destT> w_grad_k=vil3d_plane(world_grad,2);
70
71 vil3d_transform(grad_i,w_grad_i,vil3d_math_scale_functor(1.0/voxel_width_i));
72 vil3d_transform(grad_j,w_grad_j,vil3d_math_scale_functor(1.0/voxel_width_j));
73 vil3d_transform(grad_k,w_grad_k,vil3d_math_scale_functor(1.0/voxel_width_k));
74
75 // Fill 1 voxel border with zero
76 vil3d_fill_border(grad_mag,1,1,1,destT(0));
77
78 const std::ptrdiff_t gi_istep = w_grad_i.istep(), gi_jstep = w_grad_i.jstep(),
79 gi_kstep = w_grad_i.kstep();
80 const std::ptrdiff_t gj_istep = w_grad_j.istep(), gj_jstep = w_grad_j.jstep(),
81 gj_kstep = w_grad_j.kstep();
82 const std::ptrdiff_t gk_istep = w_grad_k.istep(), gk_jstep = w_grad_k.jstep(),
83 gk_kstep = w_grad_k.kstep();
84 const std::ptrdiff_t gm_istep = grad_mag.istep(), gm_jstep = grad_mag.jstep(),
85 gm_kstep = grad_mag.kstep();
86
87 unsigned ihi=ni-2;
88 unsigned jhi=nj-2;
89 unsigned khi=nk-2;
90
91 // Scaling to allow for relative size of voxels
92 double c2i=1.0/(voxel_width_i*voxel_width_i);
93 double c2j=1.0/(voxel_width_j*voxel_width_j);
94 double c2k=1.0/(voxel_width_k*voxel_width_k);
95
96 // Compute gradient magnitude at every point
97 const srcT * gi_slice = &w_grad_i(1,1,1);
98 const srcT * gj_slice = &w_grad_j(1,1,1);
99 const srcT * gk_slice = &w_grad_k(1,1,1);
100 destT * gm_slice = &grad_mag(1,1,1);
101
102 for (unsigned k=1; k<=khi; ++k, gi_slice+=gi_kstep, gj_slice+=gj_kstep,
103 gk_slice+=gk_kstep, gm_slice+=gm_kstep)
104 {
105 const srcT* gi_row=gi_slice;
106 const srcT* gj_row=gj_slice;
107 const srcT* gk_row=gk_slice;
108 destT *gm_row = gm_slice;
109
110 for (unsigned j=2; j<=jhi; ++j, gi_row+=gi_jstep, gj_row+=gj_jstep,
111 gk_row+=gk_jstep, gm_row+=gm_jstep)
112 {
113 const srcT* pgi = gi_row;
114 const srcT* pgj = gj_row;
115 const srcT* pgk = gk_row;
116 destT *pgm = gm_row;
117 for (unsigned i=2; i<=ihi; ++i, pgi+=gi_istep, pgj+=gj_istep,
118 pgk+=gk_istep, pgm+=gm_istep)
119 {
120 *pgm=destT(std::sqrt(double(c2i*pgi[0]*pgi[0] + c2j*pgj[0]*pgj[0] + c2k*pgk[0]*pgk[0])));
121 }
122 }
123 }
124 }
125
126 #undef VIL3D_WORLD_GRADIENTS_INSTANTIATE
127 #define VIL3D_WORLD_GRADIENTS_INSTANTIATE(srcT, destT) \
128 template void vil3d_world_gradients(const vil3d_image_view<srcT >& grad_i,\
129 const vil3d_image_view<srcT >& grad_j,\
130 const vil3d_image_view<srcT >& grad_k,\
131 double voxel_width_i,\
132 double voxel_width_j,\
133 double voxel_width_k,\
134 vil3d_image_view<destT >& world_grad,\
135 vil3d_image_view<destT >& grad_mag)
136
137 #endif // vil3d_world_gradients_hxx_
138