1 #ifndef vimt_resample_bilin_h_
2 #define vimt_resample_bilin_h_
3 //:
4 // \file
5 // \brief Sample grid of points in one image and place in another
6 // \author Tim Cootes
7
8 #include <cassert>
9 #ifdef _MSC_VER
10 # include <vcl_msvc_warnings.h>
11 #endif
12 #include <vil/vil_resample_bilin.h>
13 #include <vil/algo/vil_gauss_reduce.h>
14 #include <vgl/vgl_point_2d.h>
15 #include <vgl/vgl_vector_2d.h>
16 #include <vimt/vimt_image_2d_of.h>
17
18 //: Sample grid of points in one image and place in another, using bilinear interpolation.
19 // dest_image(i,j,p) is sampled from the src_image at
20 // p+i.u+j.v, where i=[0..n1-1], j=[0..n2-1] in world co-ordinates.
21 //
22 // dest_image resized to (n1,n2,src_image.nplanes())
23 //
24 // dest_image.world2im() set up so that the world co-ordinates in src and dest match
25 //
26 // Points outside image return zero.
27 // \relatesalso vimt_image_view
28 template <class sType, class dType>
vimt_resample_bilin(const vimt_image_2d_of<sType> & src_image,vimt_image_2d_of<dType> & dest_image,const vgl_point_2d<double> & p,const vgl_vector_2d<double> & u,const vgl_vector_2d<double> & v,int n1,int n2)29 inline void vimt_resample_bilin(
30 const vimt_image_2d_of<sType>& src_image,
31 vimt_image_2d_of<dType>& dest_image,
32 const vgl_point_2d<double>& p,
33 const vgl_vector_2d<double>& u,
34 const vgl_vector_2d<double>& v,
35 int n1, int n2)
36 {
37 // Not implemented for projective yet
38 assert(src_image.world2im().form()!=vimt_transform_2d::Projective);
39
40 const vimt_transform_2d& s_w2i = src_image.world2im();
41 vgl_point_2d<double> im_p = s_w2i(p);
42 vgl_vector_2d<double> im_u = s_w2i.delta(p, u);
43 vgl_vector_2d<double> im_v = s_w2i.delta(p, v);
44
45 vil_resample_bilin(src_image.image(),dest_image.image(),
46 im_p.x(),im_p.y(), im_u.x(),im_u.y(),
47 im_v.x(),im_v.y(), n1,n2);
48
49 // Point (i,j) in dest corresponds to p+i.u+j.v,
50 // an affine transformation for image to world
51 vimt_transform_2d d_i2w;
52 d_i2w.set_affine(p,u,v);
53 dest_image.set_world2im(d_i2w.inverse());
54 }
55
56
57 //: Sample grid of points in one image and place in another, using bilinear interpolation.
58 // dest_image(i,j,p) is sampled from the src_image at
59 // p+i.u+j.v, where i=[0..n1-1], j=[0..n2-1] in world co-ordinates.
60 //
61 // dest_image resized to (n1,n2,src_image.nplanes())
62 //
63 // dest_image.world2im() set up so that the world co-ordinates in src and dest match
64 //
65 // Points outside image return the value of the nearest valid pixel.
66 // \relatesalso vimt_image_view
67 template <class sType, class dType>
vimt_resample_bilin_edge_extend(const vimt_image_2d_of<sType> & src_image,vimt_image_2d_of<dType> & dest_image,const vgl_point_2d<double> & p,const vgl_vector_2d<double> & u,const vgl_vector_2d<double> & v,int n1,int n2)68 inline void vimt_resample_bilin_edge_extend(
69 const vimt_image_2d_of<sType>& src_image,
70 vimt_image_2d_of<dType>& dest_image,
71 const vgl_point_2d<double>& p,
72 const vgl_vector_2d<double>& u,
73 const vgl_vector_2d<double>& v,
74 int n1, int n2)
75 {
76 // Not implemented for projective yet
77 assert(src_image.world2im().form()!=vimt_transform_2d::Projective);
78
79 const vimt_transform_2d& s_w2i = src_image.world2im();
80 vgl_point_2d<double> im_p = s_w2i(p);
81 vgl_vector_2d<double> im_u = s_w2i.delta(p, u);
82 vgl_vector_2d<double> im_v = s_w2i.delta(p, v);
83
84 vil_resample_bilin_edge_extend(src_image.image(),dest_image.image(),
85 im_p.x(),im_p.y(), im_u.x(),im_u.y(),
86 im_v.x(),im_v.y(), n1,n2);
87
88 // Point (i,j) in dest corresponds to p+i.u+j.v,
89 // an affine transformation for image to world
90 vimt_transform_2d d_i2w;
91 d_i2w.set_affine(p,u,v);
92 dest_image.set_world2im(d_i2w.inverse());
93 }
94
95
96 //: Resample an image using appropriate smoothing if the resolution changes significantly.
97 // dest_image(i,j,p) is sampled from the src_image at
98 // p+i.u+j.v, where i=[0..n1-1], j=[0..n2-1] in world co-ordinates.
99 //
100 // dest_image resized to (n1,n2,src_image.nplanes())
101 //
102 // dest_image.world2im() set up so that the world co-ordinates in src and dest match
103 //
104 // Points outside image return the value of the nearest valid pixel.
105 // \relatesalso vimt_image_view
106 template <class sType, class dType>
vimt_resample_bilin_smoothing_edge_extend(const vimt_image_2d_of<sType> & src_image,vimt_image_2d_of<dType> & dest_image,const vgl_point_2d<double> & p,const vgl_vector_2d<double> & u,const vgl_vector_2d<double> & v,int n1,int n2)107 inline void vimt_resample_bilin_smoothing_edge_extend(
108 const vimt_image_2d_of<sType>& src_image,
109 vimt_image_2d_of<dType>& dest_image,
110 const vgl_point_2d<double>& p,
111 const vgl_vector_2d<double>& u,
112 const vgl_vector_2d<double>& v,
113 int n1, int n2)
114 {
115 // Not implemented for projective yet
116 assert(src_image.world2im().form()!=vimt_transform_2d::Projective);
117
118 vimt_transform_2d scaling;
119 scaling.set_zoom_only(0.5,0,0);
120
121 vimt_image_2d_of<sType> im = src_image;
122 vgl_vector_2d<double> im_d = im.world2im().delta(p, u) + im.world2im().delta(p, v);
123
124 // If step length (in pixels) >> 1, smooth and reduce image x2.
125 // Don't use strict Nyqusit limit.
126 // Since we are just using factor 2 smoothing and reduction,
127 // we have to tradeoff aliasing with information loss.
128 while (im_d.length() > 1.33*1.414 && im.image().ni() > 5 && im.image().nj() > 5)
129 {
130 vimt_image_2d_of<sType> dest;
131 vil_image_view<sType> work;
132 vil_gauss_reduce(im.image(), dest.image(), work);
133
134 dest.set_world2im(scaling * im.world2im());
135
136 // re-establish loop invariant
137 im = dest;
138 im_d = im.world2im().delta(p, u) + im.world2im().delta(p, v);
139 }
140
141 vimt_resample_bilin_edge_extend(im, dest_image, p, u, v, n1, n2);
142 }
143
144
145 //: Resample an image using appropriate smoothing if the resolution changes significantly.
146 // dest_image(i,j,p) is sampled from the src_image at
147 // p+i.u+j.v, where i=[0..n1-1], j=[0..n2-1] in world co-ordinates.
148 //
149 // dest_image resized to (n1,n2,src_image.nplanes())
150 //
151 // dest_image.world2im() set up so that the world co-ordinates in src and dest match
152 //
153 // Points outside image return 0.
154 // \relatesalso vimt_image_view
155 template <class sType, class dType>
vimt_resample_bilin_smoothing(const vimt_image_2d_of<sType> & src_image,vimt_image_2d_of<dType> & dest_image,const vgl_point_2d<double> & p,const vgl_vector_2d<double> & u,const vgl_vector_2d<double> & v,int n1,int n2)156 inline void vimt_resample_bilin_smoothing(
157 const vimt_image_2d_of<sType>& src_image,
158 vimt_image_2d_of<dType>& dest_image,
159 const vgl_point_2d<double>& p,
160 const vgl_vector_2d<double>& u,
161 const vgl_vector_2d<double>& v,
162 int n1, int n2)
163 {
164 // Not implemented for projective yet
165 assert(src_image.world2im().form()!=vimt_transform_2d::Projective);
166
167 vimt_transform_2d scaling;
168 scaling.set_zoom_only(0.5,0,0);
169
170 vimt_image_2d_of<sType> im = src_image;
171 vgl_vector_2d<double> im_d = im.world2im().delta(p, u) + im.world2im().delta(p, v);
172
173 // If step length (in pixels) >> 1, smooth and reduce image x2.
174 // Don't use strict Nyqusit limit.
175 // Since we are just using factor 2 smoothing and reduction,
176 // we have to tradeoff aliasing with information loss.
177 while (im_d.length() > 1.33*1.414 && im.image().ni() > 5 && im.image().nj() > 5)
178 {
179 vimt_image_2d_of<sType> dest;
180 vil_image_view<sType> work;
181 vil_gauss_reduce(im.image(), dest.image(), work);
182
183 dest.set_world2im(scaling * im.world2im());
184
185 // re-establish loop invariant
186 im = dest;
187 im_d = im.world2im().delta(p, u) + im.world2im().delta(p, v);
188 }
189
190 vimt_resample_bilin(im, dest_image, p, u, v, n1, n2);
191 }
192
193 #endif // vimt_resample_bilin_h_
194