1 // This is mul/vimt3d/vimt3d_resample_trilinear.h
2 #ifndef vimt3d_resample_trilinear_h_
3 #define vimt3d_resample_trilinear_h_
4 //:
5 // \file
6 // \brief Resample a 3D image by a different factor in each dimension
7 // \author Kevin de Souza
8
9 #include <vgl/vgl_point_3d.h>
10 #include <vgl/vgl_vector_3d.h>
11 #include <vil3d/vil3d_image_view.h>
12 #include <vil3d/algo/vil3d_gauss_reduce.h>
13 #include <vil3d/vil3d_resample_trilinear.h>
14 #include <vimt3d/vimt3d_image_3d_of.h>
15
16 //: Resample a 3D image by a factor of 2 in each dimension.
17 // \p dst_image has size 2*src.image().n?()-1 in each direction.
18 // Transform is modified by an appropriate scaling of 0.5
19 // Interpolated values are truncated when the type T is smaller than double.
20 // \sa vil3d_resample_trilinear_scale_2()
21 template <class T>
vimt3d_resample_trilinear_scale_2(const vimt3d_image_3d_of<T> & src,vimt3d_image_3d_of<T> & dst)22 void vimt3d_resample_trilinear_scale_2(
23 const vimt3d_image_3d_of<T>& src,
24 vimt3d_image_3d_of<T>& dst)
25 {
26 vil3d_resample_trilinear_scale_2(src.image(), dst.image());
27
28 vimt3d_transform_3d scaling;
29 scaling.set_zoom_only(2.0, 2.0, 2.0, 0.0, 0.0, 0.0);
30 dst.set_world2im(scaling * src.world2im());
31 }
32
33
34 //: Sample grid of points in one image and place in another, using trilinear interpolation.
35 // dest_image(i,j,k,p) is sampled from the src_image at p+i.u+j.v+k.w,
36 // where i=[0..n1-1], j=[0..n2-1], k=[0..n3-1] in world co-ordinates.
37 //
38 // dest_image resized to (n1,n2,n3,src_image.nplanes())
39 //
40 // dest_image.world2im() set up so that the world co-ordinates in src and dest match
41 //
42 // Points outside image return zero or \a outval
43 template <class sType, class dType>
44 inline void vimt3d_resample_trilinear(
45 const vimt3d_image_3d_of<sType>& src_image,
46 vimt3d_image_3d_of<dType>& dest_image,
47 const vgl_point_3d<double>& p,
48 const vgl_vector_3d<double>& u,
49 const vgl_vector_3d<double>& v,
50 const vgl_vector_3d<double>& w,
51 int n1, int n2, int n3,
52 dType outval=0, double edge_tol=0)
53 {
54 const vimt3d_transform_3d& s_w2i = src_image.world2im();
55 vgl_point_3d<double> im_p = s_w2i(p);
56 vgl_vector_3d<double> im_u = s_w2i.delta(p, u);
57 vgl_vector_3d<double> im_v = s_w2i.delta(p, v);
58 vgl_vector_3d<double> im_w = s_w2i.delta(p, w);
59
60 vil3d_resample_trilinear(src_image.image(),dest_image.image(),
61 im_p.x(), im_p.y(), im_p.z(),
62 im_u.x(), im_u.y(), im_u.z(),
63 im_v.x(), im_v.y(), im_v.z(),
64 im_w.x(), im_w.y(), im_w.z(),
65 n1, n2, n3,
66 outval, edge_tol);
67
68 // Point (i,j,k) in dest corresponds to p+i.u+j.v+k.w,
69 // an affine transformation for image to world
70 vimt3d_transform_3d d_i2w;
71 d_i2w.set_affine(p,u,v,w);
72 d_i2w.simplify();
73 dest_image.set_world2im(d_i2w.inverse());
74 }
75
76
77 //: Sample grid of points in one image and place in another, using trilinear interpolation.
78 // dest_image(i,j,k,p) is sampled from the src_image at
79 // p+i.u+j.v+k.w, where i=[0..nk-1], j=[0..nj-1], k=[0..nk-1] in world co-ordinates.
80 //
81 // dest_image resized to (ni,nj,nk,src_image.nplanes())
82 //
83 // dest_image.world2im() set up so that the world co-ordinates in src and dest match
84 //
85 // Points outside image return the value of the nearest valid pixel.
86 // \relatesalso vimt3d_image_3d_of
87 template <class sType, class dType>
vimt3d_resample_trilin_edge_extend(const vimt3d_image_3d_of<sType> & src_image,vimt3d_image_3d_of<dType> & dest_image,const vgl_point_3d<double> & p,const vgl_vector_3d<double> & u,const vgl_vector_3d<double> & v,const vgl_vector_3d<double> & w,int ni,int nj,int nk)88 inline void vimt3d_resample_trilin_edge_extend(
89 const vimt3d_image_3d_of<sType>& src_image,
90 vimt3d_image_3d_of<dType>& dest_image,
91 const vgl_point_3d<double>& p,
92 const vgl_vector_3d<double>& u,
93 const vgl_vector_3d<double>& v,
94 const vgl_vector_3d<double>& w,
95 int ni, int nj, int nk)
96 {
97 const vimt3d_transform_3d& s_w2i = src_image.world2im();
98 vgl_point_3d<double> im_p = s_w2i(p);
99 vgl_vector_3d<double> im_u = s_w2i.delta(p, u);
100 vgl_vector_3d<double> im_v = s_w2i.delta(p, v);
101 vgl_vector_3d<double> im_w = s_w2i.delta(p, w);
102
103 vil3d_resample_trilinear_edge_extend(src_image.image(),dest_image.image(),
104 im_p.x(),im_p.y(),im_p.z(), im_u.x(),im_u.y(),im_u.z(),
105 im_v.x(),im_v.y(),im_v.z(), im_w.x(),im_w.y(),im_w.z(), ni,nj,nk);
106
107 // Point (i,j,k) in dest corresponds to p+i.u+j.v+k.w,
108 // an affine transformation for image to world
109 vimt3d_transform_3d d_i2w;
110 d_i2w.set_affine(p,u,v,w);
111 d_i2w.simplify();
112 dest_image.set_world2im(d_i2w.inverse());
113 }
114
115
116 //: Resample an image using appropriate smoothing if the resolution changes significantly.
117 // dest_image(i,j,k,p) is sampled from the src_image at
118 // p+i.u+j.v+k.w, where i=[0..ni-1], j=[0..nj-1], k=[0..nk-1] in world co-ordinates.
119 //
120 // dest_image resized to (ni,nj,nk,src_image.nplanes())
121 //
122 // dest_image.world2im() set up so that the world co-ordinates in src and dest match
123 //
124 // Points outside image return the value of the nearest valid pixel.
125 // \relatesalso vimt3d_image_3d_of
126 template <class sType, class dType>
vimt3d_resample_trilin_smoothing_edge_extend(const vimt3d_image_3d_of<sType> & src_image,vimt3d_image_3d_of<dType> & dest_image,const vgl_point_3d<double> & p,const vgl_vector_3d<double> & u,const vgl_vector_3d<double> & v,const vgl_vector_3d<double> & w,int ni,int nj,int nk)127 inline void vimt3d_resample_trilin_smoothing_edge_extend(
128 const vimt3d_image_3d_of<sType>& src_image,
129 vimt3d_image_3d_of<dType>& dest_image,
130 const vgl_point_3d<double>& p,
131 const vgl_vector_3d<double>& u,
132 const vgl_vector_3d<double>& v,
133 const vgl_vector_3d<double>& w,
134 int ni, int nj, int nk)
135 {
136 vimt3d_transform_3d scaling;
137 scaling.set_zoom_only(0.5,0,0,0);
138
139 vimt3d_image_3d_of<sType> im = src_image;
140 vgl_vector_3d<double> im_d = im.world2im().delta(p, u) + im.world2im().delta(p, v)
141 + im.world2im().delta(p, w);
142
143 // If step length (in pixels) >> 1, smooth and reduce image x2.
144 // Don't use strict Nyqusit limit.
145 // Since we are just using factor 2 smoothing and reduction,
146 // we have to tradeoff aliasing with information loss.
147 while (im_d.length() > 1.33*1.732 && im.image().ni() > 5 && im.image().nj() > 5 && im.image().nk() > 5)
148 {
149 vimt3d_image_3d_of<sType> dest;
150 vil3d_image_view<sType> work1, work2;
151 vil3d_gauss_reduce(im.image(), dest.image(), work1, work2);
152
153 dest.set_world2im(scaling * im.world2im());
154
155 // re-establish loop invariant
156 im = dest;
157 im_d = im.world2im().delta(p, u) + im.world2im().delta(p, v) + im.world2im().delta(p, w);
158 }
159
160 vimt3d_resample_trilin_edge_extend(im, dest_image, p, u, v, w, ni, nj, nk);
161 }
162
163 //: Resample src, using the grid defined by dest.
164 // Smooths appropriately if the resolution changes significantly.
165 // dest(i,j,k,p) is sampled from src at the wc grid defined by
166 // the world co-ords of the pixel centres in dest.
167 //
168 // dest is not resized, nor has its world2im transform modified.
169 //
170 // Points outside src return the value of the nearest valid pixel.
171 // \relatesalso vimt3d_image_3d_of
172 template <class sType, class dType>
vimt3d_resample_trilin_smoothing_edge_extend(const vimt3d_image_3d_of<sType> & src,vimt3d_image_3d_of<dType> & dest)173 inline void vimt3d_resample_trilin_smoothing_edge_extend(
174 const vimt3d_image_3d_of<sType>& src,
175 vimt3d_image_3d_of<dType>& dest)
176 {
177 vimt3d_transform_3d orig_w2i = dest.world2im();
178 vimt3d_transform_3d i2w = orig_w2i.inverse();
179
180 vimt3d_resample_trilin_smoothing_edge_extend(
181 src, dest,
182 i2w.origin(),
183 i2w.delta(vgl_point_3d<double>(0,0,0), vgl_vector_3d<double>(1,0,0)),
184 i2w.delta(vgl_point_3d<double>(0,0,0), vgl_vector_3d<double>(0,1,0)),
185 i2w.delta(vgl_point_3d<double>(0,0,0), vgl_vector_3d<double>(0,0,1)),
186 dest.image().ni(), dest.image().nj(), dest.image().nk() );
187
188 dest.world2im() = orig_w2i;
189 }
190
191
192 #endif // vimt3d_resample_trilinear_h_
193