1 // This is mul/vil3d/vil3d_resample_trilinear.hxx
2 #ifndef vil3d_resample_trilinear_hxx_
3 #define vil3d_resample_trilinear_hxx_
4 //:
5 // \file
6 // \brief Resample a 3D image by a different factor in each dimension
7 // \author Kevin de Souza, Ian Scott
8
9 #include "vil3d_resample_trilinear.h"
10
11 #include <vil/vil_convert.h>
12 #include <vil3d/vil3d_trilin_interp.h>
13 #include <vil3d/vil3d_plane.h>
14 #include <cassert>
15 #ifdef _MSC_VER
16 # include <vcl_msvc_warnings.h>
17 #endif
18
19
vil3dresample_trilin_corner_in_image(double x0,double y0,double z0,const vil3d_image_view_base & image)20 inline bool vil3dresample_trilin_corner_in_image(double x0, double y0, double z0,
21 const vil3d_image_view_base& image)
22 {
23 if (x0<0.0) return false;
24 if (x0>image.ni()-1.0) return false;
25 if (y0<0.0) return false;
26 if (y0>image.nj()-1.0) return false;
27 if (z0<0.0) return false;
28 if (z0>image.nk()-1.0) return false;
29 return true;
30 }
31
32
33 // Sample grid of points in one image and place in another, using trilinear interpolation.
34 // dest_image(i,j,k,p) is sampled from the src_image at
35 // (x0+i.dx1+j.dx2+k.dx3, y0+i.dy1+j.dy2+k.dy3, z0+i.dz1+j.dz2+k.dz3),
36 // where i=[0..n1-1], j=[0..n2-1], k=[0..n3-1]
37 // dest_image resized to (n1,n2,n3,src_image.nplanes())
38 // Points outside image return value of nearest pixel.
39 template <class S, class T>
vil3d_resample_trilinear_edge_extend(const vil3d_image_view<S> & src_image,vil3d_image_view<T> & dest_image,double x0,double y0,double z0,double dx1,double dy1,double dz1,double dx2,double dy2,double dz2,double dx3,double dy3,double dz3,int n1,int n2,int n3)40 void vil3d_resample_trilinear_edge_extend(const vil3d_image_view<S>& src_image,
41 vil3d_image_view<T>& dest_image,
42 double x0, double y0, double z0,
43 double dx1, double dy1, double dz1,
44 double dx2, double dy2, double dz2,
45 double dx3, double dy3, double dz3,
46 int n1, int n2, int n3)
47 {
48 bool all_in_image =
49 vil3dresample_trilin_corner_in_image(x0,
50 y0,
51 z0,
52 src_image) &&
53 vil3dresample_trilin_corner_in_image(x0 + (n1-1)*dx1,
54 y0 + (n1-1)*dy1,
55 z0 + (n1-1)*dz1,
56 src_image) &&
57 vil3dresample_trilin_corner_in_image(x0 + (n2-1)*dx2,
58 y0 + (n2-1)*dy2,
59 z0 + (n2-1)*dz2,
60 src_image) &&
61 vil3dresample_trilin_corner_in_image(x0 + (n1-1)*dx1 + (n2-1)*dx2,
62 y0 + (n1-1)*dy1 + (n2-1)*dy2,
63 z0 + (n1-1)*dz1 + (n2-1)*dz2,
64 src_image) &&
65 vil3dresample_trilin_corner_in_image(x0 + (n3-1)*dx3,
66 y0 + (n3-1)*dy3,
67 z0 + (n3-1)*dz3,
68 src_image) &&
69 vil3dresample_trilin_corner_in_image(x0 + (n1-1)*dx1 + (n3-1)*dx3,
70 y0 + (n1-1)*dy1 + (n3-1)*dy3,
71 z0 + (n1-1)*dz1 + (n3-1)*dz3,
72 src_image) &&
73 vil3dresample_trilin_corner_in_image(x0 + (n2-1)*dx2 + (n3-1)*dx3,
74 y0 + (n2-1)*dy2 + (n3-1)*dy3,
75 z0 + (n2-1)*dz2 + (n3-1)*dz3,
76 src_image) &&
77 vil3dresample_trilin_corner_in_image(x0 + (n1-1)*dx1 + (n2-1)*dx2 + (n3-1)*dx3,
78 y0 + (n1-1)*dy1 + (n2-1)*dy2 + (n3-1)*dy3,
79 z0 + (n1-1)*dz1 + (n2-1)*dz2 + (n3-1)*dz3,
80 src_image);
81
82 vil_convert_round_pixel<double, T> cast_and_possibly_round;
83 // Rounds only if T is int, unsigned or kind of. Doesn't round if T
84 // is float, double or kind of.
85
86 const unsigned ni = src_image.ni();
87 const unsigned nj = src_image.nj();
88 const unsigned nk = src_image.nk();
89 const unsigned np = src_image.nplanes();
90 const std::ptrdiff_t istep = src_image.istep();
91 const std::ptrdiff_t jstep = src_image.jstep();
92 const std::ptrdiff_t kstep = src_image.kstep();
93 const std::ptrdiff_t pstep = src_image.planestep();
94 const S* plane0 = src_image.origin_ptr();
95
96 dest_image.set_size(n1,n2,n3,np);
97 const std::ptrdiff_t d_istep = dest_image.istep();
98 const std::ptrdiff_t d_jstep = dest_image.jstep();
99 const std::ptrdiff_t d_kstep = dest_image.kstep();
100 const std::ptrdiff_t d_pstep = dest_image.planestep();
101 T* d_plane0 = dest_image.origin_ptr();
102
103 if (all_in_image)
104 {
105 if (np==1)
106 {
107 double xk=x0, yk=y0, zk=z0;
108 T *slice = d_plane0;
109 for (int k=0; k<n3; ++k, xk+=dx3, yk+=dy3, zk+=dz3, slice+=d_kstep)
110 {
111 double xj=xk, yj=yk, zj=zk; // Start of k-th slice
112 T *row = slice;
113 for (int j=0; j<n2; ++j, xj+=dx2, yj+=dy2, zj+=dz2, row+=d_jstep)
114 {
115 double x=xj, y=yj, z=zj; // Start of j-th row
116 T *dpt = row;
117 for (int i=0; i<n1; ++i, x+=dx1, y+=dy1, z+=dz1, dpt+=d_istep)
118 cast_and_possibly_round( vil3d_trilin_interp_raw( x, y, z,
119 plane0,
120 istep, jstep, kstep ),
121 *dpt);
122 }
123 }
124 }
125 else
126 {
127 double xk=x0, yk=y0, zk=z0;
128 T *slice = d_plane0;
129 for (int k=0; k<n3; ++k, xk+=dx3, yk+=dy3, zk+=dz3, slice+=d_kstep)
130 {
131 double xj=xk, yj=yk, zj=zk; // Start of k-th slice
132 T *row = slice;
133 for (int j=0; j<n2; ++j, xj+=dx2, yj+=dy2, zj+=dz2, row+=d_jstep)
134 {
135 double x=xj, y=yj, z=zj; // Start of j-th row
136 T *dpt = row;
137 for (int i=0; i<n1; ++i, x+=dx1, y+=dy1, z+=dz1, dpt+=d_istep)
138 {
139 for (unsigned int p=0; p<np; ++p)
140 cast_and_possibly_round( vil3d_trilin_interp_raw( x, y, z,
141 plane0+p*pstep,
142 istep, jstep, kstep),
143 dpt[p*d_pstep]);
144 }
145 }
146 }
147 }
148 }
149 else
150 {
151 // Use safe interpolation with edge-extension
152 if (np==1)
153 {
154 double xk=x0, yk=y0, zk=z0;
155 T *slice = d_plane0;
156 for (int k=0; k<n3; ++k, xk+=dx3, yk+=dy3, zk+=dz3, slice+=d_kstep)
157 {
158 double xj=xk, yj=yk, zj=zk; // Start of k-th slice
159 T *row = slice;
160 for (int j=0; j<n2; ++j, xj+=dx2, yj+=dy2, zj+=dz2, row+=d_jstep)
161 {
162 double x=xj, y=yj, z=zj; // Start of j-th row
163 T *dpt = row;
164 for (int i=0; i<n1; ++i, x+=dx1, y+=dy1, z+=dz1, dpt+=d_istep)
165 cast_and_possibly_round( vil3d_trilin_interp_safe_extend( x, y, z,
166 plane0,
167 ni, nj, nk,
168 istep, jstep, kstep),
169 *dpt );
170 }
171 }
172 }
173 else
174 {
175 double xk=x0, yk=y0, zk=z0;
176 T *slice = d_plane0;
177 for (int k=0; k<n3; ++k, xk+=dx3, yk+=dy3, zk+=dz3, slice+=d_kstep)
178 {
179 double xj=xk, yj=yk, zj=zk; // Start of k-th slice
180 T *row = slice;
181 for (int j=0; j<n2; ++j, xj+=dx2, yj+=dy2, zj+=dz2, row+=d_jstep)
182 {
183 double x=xj, y=yj, z=zj; // Start of j-th row
184 T *dpt = row;
185 for (int i=0; i<n1; ++i, x+=dx1, y+=dy1, z+=dz1, dpt+=d_istep)
186 {
187 for (unsigned int p=0; p<np; ++p)
188 cast_and_possibly_round( vil3d_trilin_interp_safe_extend( x, y, z,
189 plane0+p*pstep,
190 ni, nj, nk,
191 istep, jstep, kstep),
192 dpt[p*d_pstep] );
193 }
194 }
195 }
196 }
197 }
198 }
199
200
201 // Sample grid of points in one image and place in another, using trilinear interpolation.
202 // dest_image(i,j,k,p) is sampled from the src_image at
203 // (x0+i.dx1+j.dx2+k.dx3, y0+i.dy1+j.dy2+k.dy3, z0+i.dz1+j.dz2+k.dz3),
204 // where i=[0..n1-1], j=[0..n2-1], k=[0..n3-1]
205 // dest_image resized to (n1,n2,n3,src_image.nplanes())
206 // Points outside image return zero or \a outval
207 template <class S, class T>
vil3d_resample_trilinear(const vil3d_image_view<S> & src_image,vil3d_image_view<T> & dest_image,double x0,double y0,double z0,double dx1,double dy1,double dz1,double dx2,double dy2,double dz2,double dx3,double dy3,double dz3,int n1,int n2,int n3,T outval,double edge_tol)208 void vil3d_resample_trilinear(const vil3d_image_view<S>& src_image,
209 vil3d_image_view<T>& dest_image,
210 double x0, double y0, double z0,
211 double dx1, double dy1, double dz1,
212 double dx2, double dy2, double dz2,
213 double dx3, double dy3, double dz3,
214 int n1, int n2, int n3,
215 T outval/*=0*/, double edge_tol/*=0*/)
216 {
217 bool all_in_image =
218 vil3dresample_trilin_corner_in_image(x0,
219 y0,
220 z0,
221 src_image) &&
222 vil3dresample_trilin_corner_in_image(x0 + (n1-1+edge_tol)*dx1,
223 y0 + (n1-1+edge_tol)*dy1,
224 z0 + (n1-1+edge_tol)*dz1,
225 src_image) &&
226 vil3dresample_trilin_corner_in_image(x0 + (n2-1+edge_tol)*dx2,
227 y0 + (n2-1+edge_tol)*dy2,
228 z0 + (n2-1+edge_tol)*dz2,
229 src_image) &&
230 vil3dresample_trilin_corner_in_image(x0 + (n1-1+edge_tol)*dx1 + (n2-1+edge_tol)*dx2,
231 y0 + (n1-1+edge_tol)*dy1 + (n2-1+edge_tol)*dy2,
232 z0 + (n1-1+edge_tol)*dz1 + (n2-1+edge_tol)*dz2,
233 src_image) &&
234 vil3dresample_trilin_corner_in_image(x0 + (n3-1+edge_tol)*dx3,
235 y0 + (n3-1+edge_tol)*dy3,
236 z0 + (n3-1+edge_tol)*dz3,
237 src_image) &&
238 vil3dresample_trilin_corner_in_image(x0 + (n1-1+edge_tol)*dx1 + (n3-1+edge_tol)*dx3,
239 y0 + (n1-1+edge_tol)*dy1 + (n3-1+edge_tol)*dy3,
240 z0 + (n1-1+edge_tol)*dz1 + (n3-1+edge_tol)*dz3,
241 src_image) &&
242 vil3dresample_trilin_corner_in_image(x0 + (n2-1+edge_tol)*dx2 + (n3-1+edge_tol)*dx3,
243 y0 + (n2-1+edge_tol)*dy2 + (n3-1+edge_tol)*dy3,
244 z0 + (n2-1+edge_tol)*dz2 + (n3-1+edge_tol)*dz3,
245 src_image) &&
246 vil3dresample_trilin_corner_in_image(x0 + (n1-1+edge_tol)*dx1 + (n2-1+edge_tol)*dx2 + (n3-1+edge_tol)*dx3,
247 y0 + (n1-1+edge_tol)*dy1 + (n2-1+edge_tol)*dy2 + (n3-1+edge_tol)*dy3,
248 z0 + (n1-1+edge_tol)*dz1 + (n2-1+edge_tol)*dz2 + (n3-1+edge_tol)*dz3,
249 src_image);
250
251 vil_convert_round_pixel<double, T> cast_and_possibly_round;
252
253 const unsigned ni = src_image.ni();
254 const unsigned nj = src_image.nj();
255 const unsigned nk = src_image.nk();
256 const unsigned np = src_image.nplanes();
257 const std::ptrdiff_t istep = src_image.istep();
258 const std::ptrdiff_t jstep = src_image.jstep();
259 const std::ptrdiff_t kstep = src_image.kstep();
260 const std::ptrdiff_t pstep = src_image.planestep();
261 const S* plane0 = src_image.origin_ptr();
262
263 dest_image.set_size(n1,n2,n3,np);
264 const std::ptrdiff_t d_istep = dest_image.istep();
265 const std::ptrdiff_t d_jstep = dest_image.jstep();
266 const std::ptrdiff_t d_kstep = dest_image.kstep();
267 const std::ptrdiff_t d_pstep = dest_image.planestep();
268 T* d_plane0 = dest_image.origin_ptr();
269
270 if (all_in_image)
271 {
272 if (np==1)
273 {
274 double xk=x0, yk=y0, zk=z0;
275 T *slice = d_plane0;
276 for (int k=0; k<n3; ++k, xk+=dx3, yk+=dy3, zk+=dz3, slice+=d_kstep)
277 {
278 double xj=xk, yj=yk, zj=zk; // Start of k-th slice
279 T *row = slice;
280 for (int j=0; j<n2; ++j, xj+=dx2, yj+=dy2, zj+=dz2, row+=d_jstep)
281 {
282 double x=xj, y=yj, z=zj; // Start of j-th row
283 T *dpt = row;
284 for (int i=0; i<n1; ++i, x+=dx1, y+=dy1, z+=dz1, dpt+=d_istep)
285 cast_and_possibly_round( vil3d_trilin_interp_raw( x, y, z,
286 plane0,
287 istep, jstep, kstep),
288 *dpt);
289 }
290 }
291 }
292 else
293 {
294 double xk=x0, yk=y0, zk=z0;
295 T *slice = d_plane0;
296 for (int k=0; k<n3; ++k, xk+=dx3, yk+=dy3, zk+=dz3, slice+=d_kstep)
297 {
298 double xj=xk, yj=yk, zj=zk; // Start of k-th slice
299 T *row = slice;
300 for (int j=0; j<n2; ++j, xj+=dx2, yj+=dy2, zj+=dz2, row+=d_jstep)
301 {
302 double x=xj, y=yj, z=zj; // Start of j-th row
303 T *dpt = row;
304 for (int i=0; i<n1; ++i, x+=dx1, y+=dy1, z+=dz1, dpt+=d_istep)
305 {
306 for (unsigned int p=0; p<np; ++p)
307 cast_and_possibly_round( vil3d_trilin_interp_raw( x, y, z,
308 plane0+p*pstep,
309 istep, jstep, kstep),
310 dpt[p*d_pstep]);
311 }
312 }
313 }
314 }
315 }
316 else
317 {
318 // Use safe interpolation
319 if (np==1)
320 {
321 double xk=x0, yk=y0, zk=z0;
322 T *slice = d_plane0;
323 for (int k=0; k<n3; ++k, xk+=dx3, yk+=dy3, zk+=dz3, slice+=d_kstep)
324 {
325 double xj=xk, yj=yk, zj=zk; // Start of k-th slice
326 T *row = slice;
327 for (int j=0; j<n2; ++j, xj+=dx2, yj+=dy2, zj+=dz2, row+=d_jstep)
328 {
329 double x=xj, y=yj, z=zj; // Start of j-th row
330 T *dpt = row;
331 for (int i=0; i<n1; ++i, x+=dx1, y+=dy1, z+=dz1, dpt+=d_istep)
332 cast_and_possibly_round( vil3d_trilin_interp_safe( x, y, z,
333 plane0,
334 ni, nj, nk,
335 istep, jstep, kstep, outval),
336 *dpt);
337 }
338 }
339 }
340 else
341 {
342 double xk=x0, yk=y0, zk=z0;
343 T *slice = d_plane0;
344 for (int k=0; k<n3; ++k, xk+=dx3, yk+=dy3, zk+=dz3, slice+=d_kstep)
345 {
346 double xj=xk, yj=yk, zj=zk; // Start of k-th slice
347 T *row = slice;
348 for (int j=0; j<n2; ++j, xj+=dx2, yj+=dy2, zj+=dz2, row+=d_jstep)
349 {
350 double x=xj, y=yj, z=zj; // Start of j-th row
351 T *dpt = row;
352 for (int i=0; i<n1; ++i, x+=dx1, y+=dy1, z+=dz1, dpt+=d_istep)
353 {
354 for (unsigned int p=0; p<np; ++p)
355 cast_and_possibly_round( vil3d_trilin_interp_safe( x, y, z,
356 plane0+p*pstep,
357 ni, nj, nk,
358 istep, jstep, kstep, outval),
359 dpt[p*d_pstep]);
360 }
361 }
362 }
363 }
364 }
365 }
366
367
368 // Resample image to a specified dimensions (n1 * n2 * n3)
369 template <class S, class T>
vil3d_resample_trilinear(const vil3d_image_view<S> & src_image,vil3d_image_view<T> & dest_image,int n1,int n2,int n3)370 void vil3d_resample_trilinear(const vil3d_image_view<S>& src_image,
371 vil3d_image_view<T>& dest_image,
372 int n1, int n2, int n3)
373 {
374 double f= 0.9999999; // so sampler doesn't go off edge of image
375 double x0=0;
376 double y0=0;
377 double z0=0;
378
379 double dx1=f*(src_image.ni()-1)*1.0/(n1-1);
380 double dy1=0;
381 double dz1=0;
382
383 double dx2=0;
384 double dy2=f*(src_image.nj()-1)*1.0/(n2-1);
385 double dz2=0;
386
387 double dx3=0;
388 double dy3=0;
389 double dz3=f*(src_image.nk()-1)*1.0/(n3-1);
390
391 vil3d_resample_trilinear(src_image, dest_image, x0, y0, z0, dx1, dy1, dz1,
392 dx2, dy2, dz2, dx3, dy3, dz3, n1, n2, n3);
393 }
394
395
396 // Resample a 3D image by a different factor in each dimension.
397 // Note: The upper image boundaries are extended.
398 template <class T>
vil3d_resample_trilinear(const vil3d_image_view<T> & src_image,vil3d_image_view<T> & dst_image,const double dx,const double dy,const double dz)399 void vil3d_resample_trilinear(const vil3d_image_view<T>& src_image,
400 vil3d_image_view<T>& dst_image,
401 const double dx,
402 const double dy,
403 const double dz)
404 {
405 assert (dx > 0.0 && dy > 0.0 && dz > 0.0);
406
407 // Assume planes are the same for both images
408 const unsigned np = src_image.nplanes();
409
410 const unsigned sni = src_image.ni();
411 const unsigned snj = src_image.nj();
412 const unsigned snk = src_image.nk();
413 const std::ptrdiff_t s_istep = src_image.istep();
414 const std::ptrdiff_t s_jstep = src_image.jstep();
415 const std::ptrdiff_t s_kstep = src_image.kstep();
416 const std::ptrdiff_t s_pstep = src_image.planestep();
417 const T* s_plane = src_image.origin_ptr();
418
419 const unsigned dni = static_cast<unsigned>(sni*dx);
420 const unsigned dnj = static_cast<unsigned>(snj*dy);
421 const unsigned dnk = static_cast<unsigned>(snk*dz);
422 dst_image.set_size(dni, dnj, dnk, np);
423 const std::ptrdiff_t d_istep = dst_image.istep();
424 const std::ptrdiff_t d_jstep = dst_image.jstep();
425 const std::ptrdiff_t d_kstep = dst_image.kstep();
426 const std::ptrdiff_t d_pstep = dst_image.planestep();
427 T* d_plane = dst_image.origin_ptr();
428
429 // Use this to convert from double to T with appropriate rounding
430 vil_convert_round_pixel<double, T> cast_and_possibly_round;
431
432 // Loop over all voxels in the destination image
433 // and sample from the corresponding point in the source image
434 // (except near upper boundaries).
435 for (unsigned p=0; p<np; ++p, s_plane+=s_pstep, d_plane+=d_pstep)
436 {
437 T* d_slice = d_plane;
438 for (unsigned k=0; k<static_cast<unsigned>(dnk-dz); ++k, d_slice+=d_kstep)
439 {
440 T* d_row = d_slice;
441 for (unsigned j=0; j<static_cast<unsigned>(dnj-dy); ++j, d_row+=d_jstep)
442 {
443 T* d_pix = d_row;
444 for (unsigned i=0; i<static_cast<unsigned>(dni-dx); ++i, d_pix+=d_istep)
445 {
446 cast_and_possibly_round( vil3d_trilin_interp_raw( i/dx, j/dy, k/dz,
447 s_plane,
448 s_istep, s_jstep, s_kstep),
449 *d_pix);
450 }
451 // Process the pixels near the upper i boundary - safe_extend interpolation
452 for (unsigned i=static_cast<unsigned>(dni-dx); i<dni; ++i, d_pix+=d_istep)
453 {
454 cast_and_possibly_round( vil3d_trilin_interp_safe_extend(i/dx, j/dy, k/dz,
455 s_plane,
456 sni, snj, snk,
457 s_istep, s_jstep, s_kstep),
458 *d_pix);
459 }
460 }
461
462 // Process the pixels near the upper j boundary - safe_extend interpolation
463 for (unsigned j=static_cast<unsigned>(dnj-dy); j<dnj; ++j, d_row+=d_jstep)
464 {
465 T* d_pix = d_row;
466 for (unsigned i=0; i<dni; ++i, d_pix+=d_istep)
467 {
468 cast_and_possibly_round( vil3d_trilin_interp_safe_extend( i/dx, j/dy, k/dz,
469 s_plane,
470 sni, snj, snk,
471 s_istep, s_jstep, s_kstep),
472 *d_pix);
473 }
474 }
475 }
476
477 // Process the pixels near the upper k boundary - safe_extend interpolation
478 for (unsigned k=static_cast<unsigned>(dnk-dz); k<dnk; ++k, d_slice+=d_kstep)
479 {
480 T* d_row = d_slice;
481 for (unsigned j=0; j<dnj; ++j, d_row+=d_jstep)
482 {
483 T* d_pix = d_row;
484 for (unsigned i=0; i<dni; ++i, d_pix+=d_istep)
485 {
486 cast_and_possibly_round( vil3d_trilin_interp_safe_extend( i/dx, j/dy, k/dz,
487 s_plane,
488 sni, snj, snk,
489 s_istep, s_jstep, s_kstep),
490 *d_pix);
491 }
492 }
493 }
494 }
495 }
496
497
498 template <class T>
vil3d_resample_trilinear_scale_2(const vil3d_image_view<T> & src_im,vil3d_image_view<T> & dest_im)499 void vil3d_resample_trilinear_scale_2(
500 const vil3d_image_view<T>& src_im,
501 vil3d_image_view<T>& dest_im)
502 {
503 // Assume planes are the same for both images
504 const unsigned np = src_im.nplanes();
505
506 const unsigned ni = src_im.ni();
507 const unsigned nj = src_im.nj();
508 const unsigned nk = src_im.nk();
509 dest_im.set_size(ni*2-1, nj*2-1, nk*2-1, np);
510
511 for (unsigned p = 0; p<np; ++p)
512 {
513 const vil3d_image_view<T> src = vil3d_plane(src_im, p);
514 vil3d_image_view<T> dest = vil3d_plane(dest_im, p);
515
516 for (unsigned k=0; k<nk-1; ++k)
517 {
518 for (unsigned j=0; j<nj-1; ++j)
519 {
520 // Do all except last slice.
521 for (unsigned i=0; i<ni-1; ++i)
522 {
523 // s0-s7 are the values at 8 neighbouring positions (on a cube) in src.
524 // d0-d7 are the values at 8 neighbouring positions (on a cube of 1/8 the above cube's size) in dest.
525 // They are aligned so that s0 = d0.
526 // d6t2 is the value of d6 time 2, etc.
527 T s0 = src(i,j,k);
528 T s1 = src(i+1, j , k );
529 T s2 = src(i, j+1, k );
530 T s3 = src(i+1, j+1, k );
531 T s4 = src(i , j , k+1);
532 T s5 = src(i+1, j , k+1);
533 T s6 = src(i , j+1, k+1);
534 T s7 = src(i+1, j+1, k+1);
535 dest(2*i , 2*j , 2*k ) = s0;
536 T d1t2 = s0+s1; dest(2*i+1, 2*j , 2*k ) = d1t2/2;
537 T d2t2 = s0+s2; dest(2*i , 2*j+1, 2*k ) = d2t2/2;
538 T d4t2 = s0+s4; dest(2*i , 2*j , 2*k+1) = d4t2/2;
539 T d3t4 = d2t2 + s1+s3; dest(2*i+1, 2*j+1, 2*k ) = d3t4/4;
540 T d5t4 = d4t2 + s1+s5; dest(2*i+1, 2*j , 2*k+1) = d5t4/4;
541 T d6t4 = d4t2 + s2+s6; dest(2*i , 2*j+1, 2*k+1) = d6t4/4;
542 dest(2*i+1, 2*j+1, 2*k+1) = (d6t4 + s1+s3+s5+s7)/8;
543 }
544 T s0 = src(ni-1, j , k );
545 T s2 = src(ni-1, j+1, k );
546 T s4 = src(ni-1, j , k+1);
547 T s6 = src(ni-1, j+1, k+1);
548 dest(ni*2-2, j*2 , k*2 ) = s0;
549 dest(ni*2-2, j*2+1, k*2 ) = (s0 + s2)/2;
550 dest(ni*2-2, j*2 , k*2+1) = (s0 + s4)/2;
551 dest(ni*2-2, j*2+1, k*2+1) = (s0 + s2 + s4 + s6)/4;
552 }
553 // Do last plane
554 for (unsigned i=0; i<ni-1; ++i)
555 {
556 T s0 = src(i , nj-1, k );
557 T s1 = src(i+1, nj-1, k );
558 T s4 = src(i , nj-1, k+1);
559 T s5 = src(i+1, nj-1, k+1);
560 dest(i*2 , nj*2-2, k*2 ) = s0;
561 dest(i*2+1, nj*2-2, k*2 ) = (s0 + s1)/2;
562 dest(i*2 , nj*2-2, k*2+1) = (s0 + s4)/2;
563 dest(i*2+1, nj*2-2, k*2+1) = (s0 + s1 + s4 + s5)/4;
564 }
565 T s0 = src(ni-1, nj-1, k);
566 dest(ni*2-2, nj*2-2, k*2 ) = s0;
567 dest(ni*2-2, nj*2-2, k*2+1) = (s0 + src(ni-1, nj-1, k+1))/2;
568 }
569 // Do last plane
570 for (unsigned j=0; j<nj-1; ++j)
571 {
572 for (unsigned i=0; i<ni-1; ++i)
573 {
574 T s0 = src(i , j , nk-1);
575 T s1 = src(i+1, j , nk-1);
576 T s2 = src(i , j+1, nk-1);
577 T s3 = src(i+1, j+1, nk-1);
578 dest(i*2 , j*2 , nk*2-2) = s0;
579 dest(i*2+1, j*2 , nk*2-2) = (s0 + s1)/2;
580 dest(i*2 , j*2+1, nk*2-2) = (s0 + s2)/2;
581 dest(i*2+1, j*2+1, nk*2-2) = (s0 + s1 + s2 + s3)/4;
582 }
583 T s0 = src(ni-1, j, nk-1);
584 dest(ni*2-2, j*2 , nk*2-2) = s0;
585 dest(ni*2-2, j*2+1, nk*2-2) = (s0 + src(ni-1, j+1, nk-1))/2;
586 }
587 for (unsigned i=0; i<ni-1; ++i)
588 {
589 T s0 = src(i, nj-1, nk-1);
590 dest(i*2 , nj*2-2, nk*2-2) = s0;
591 dest(i*2+1, nj*2-2, nk*2-2) = (s0 + src(i+1, nj-1, nk-1))/2;
592 }
593 dest(ni*2-2, nj*2-2, nk*2-2) = src(ni-1, nj-1, nk-1);
594 }
595 }
596
597
598 #define VIL3D_RESAMPLE_TRILINEAR_INSTANTIATE( T, S ) \
599 template void vil3d_resample_trilinear(const vil3d_image_view< S >& src_image, \
600 vil3d_image_view< T >& dest_image, \
601 double x0, double y0, double z0, \
602 double dx1, double dy1, double dz1, \
603 double dx2, double dy2, double dz2, \
604 double dx3, double dy3, double dz3, \
605 int n1, int n2, int n3, \
606 T outval, double edge_tol); \
607 template void vil3d_resample_trilinear_edge_extend(const vil3d_image_view< S >& src_image, \
608 vil3d_image_view< T >& dest_image, \
609 double x0, double y0, double z0, \
610 double dx1, double dy1, double dz1, \
611 double dx2, double dy2, double dz2, \
612 double dx3, double dy3, double dz3, \
613 int n1, int n2, int n3); \
614 template void vil3d_resample_trilinear(const vil3d_image_view< S >& src_image, \
615 vil3d_image_view< T >& dest_image, \
616 int n1, int n2, int n3); \
617 template void vil3d_resample_trilinear(const vil3d_image_view< T >& src_image, \
618 vil3d_image_view< T >& dst_image, \
619 const double dx, \
620 const double dy, \
621 const double dz); \
622 template void vil3d_resample_trilinear_scale_2(const vil3d_image_view< T >& src_image, \
623 vil3d_image_view< T >& dst_image)
624
625 #endif // vil3d_resample_trilinear_hxx_
626