1 // This is mul/vil3d/vil3d_trilin_interp.h
2 #ifndef vil3d_trilin_interp_h_
3 #define vil3d_trilin_interp_h_
4 //:
5 // \file
6 // \brief Trilinear interpolation functions for 3D images
7 // \author Tim Cootes
8
9 #include <iostream>
10 #include <cstddef>
11 #include <cassert>
12 #ifdef _MSC_VER
13 # include <vcl_msvc_warnings.h>
14 #endif
15 #include <vil/vil_na.h>
16 #include <vil3d/vil3d_image_view.h>
17
18 //: Compute trilinear interpolation at (x,y,z), no bound checks
19 // Image is nx * ny * nz array of T. x,y,z element is data[z*zstep+ystep*y+x*xstep]
20 // No bound checks are done.
21 template<class T>
vil3d_trilin_interp_raw(double x,double y,double z,const T * data,std::ptrdiff_t xstep,std::ptrdiff_t ystep,std::ptrdiff_t zstep)22 inline double vil3d_trilin_interp_raw(double x, double y, double z,
23 const T* data, std::ptrdiff_t xstep,
24 std::ptrdiff_t ystep, std::ptrdiff_t zstep)
25 {
26 int p1x,p1y,p1z;
27 double normx,normy,normz;
28 p1x=int(x);
29 normx = x-p1x;
30 p1y=int(y);
31 normy = y-p1y;
32 p1z=int(z);
33 normz = z-p1z;
34
35 const T* row11 = data + p1z*zstep+p1y*ystep + p1x*xstep;
36 const T* row21 = row11 + ystep;
37 const T* row12 = row11 + zstep;
38 const T* row22 = row21 + zstep;
39
40 // Bilinear interpolation in plane z=p1z
41 double i11 = (double)row11[0]+(double)(row21[0]-row11[0])*normy;
42 double i21 = (double)row11[xstep]+(double)(row21[xstep]-row11[xstep])*normy;
43 double iz1 = i11+(i21-i11)*normx;
44
45 // Bilinear interpolation in plane z=p1z+1
46 double i12 = (double)row12[0]+(double)(row22[0]-row12[0])*normy;
47 double i22 = (double)row12[xstep]+(double)(row22[xstep]-row12[xstep])*normy;
48 double iz2 = i12+(i22-i12)*normx;
49
50 return iz1+(iz2-iz1)*normz;
51 }
52
53
54 //: Compute trilinear interpolation at (x,y,z), with bound checks
55 // Image is nx * ny * nz array of T. x,y,z element is data[z*zstep+ystep*y+x*xstep]
56 // If (x,y,z) is outside interpolatable image region, returns zero or \a outval
57 template<class T>
58 inline double vil3d_trilin_interp_safe(double x, double y, double z, const T* data,
59 unsigned nx, unsigned ny, unsigned nz,
60 std::ptrdiff_t xstep, std::ptrdiff_t ystep, std::ptrdiff_t zstep,
61 double outval=0)
62 {
63 if (x<0) return static_cast<double>(outval);
64 if (y<0) return static_cast<double>(outval);
65 if (z<0) return static_cast<double>(outval);
66 if (x>=nx-1) return static_cast<double>(outval);
67 if (y>=ny-1) return static_cast<double>(outval);
68 if (z>=nz-1) return static_cast<double>(outval);
69 return vil3d_trilin_interp_raw(x,y,z,data,xstep,ystep,zstep);
70 }
71
72 //: Compute trilinear interpolation at (x,y,z,p), with bound checks
73 // Image is nx * ny * nz array of T. x,y,z element is data[z*zstep+ystep*y+x*xstep]
74 // If (x,y,z) is outside interpolatable image region, returns zero or \a outval
75 template<class T>
76 inline double vil3d_trilin_interp_safe(const vil3d_image_view<T>& image,
77 double x, double y, double z,
78 unsigned p=0,
79 T outval=0)
80 {
81 return vil3d_trilin_interp_safe(x,y,z,&image(0,0,0,p),
82 image.ni(),image.nj(),image.nk(),
83 image.istep(),image.jstep(),image.kstep(),
84 outval);
85 }
86
87 //: Compute trilinear interpolation at (x,y), with bound checks
88 // Image is nx * ny * nz array of Ts. x,y element is data[nx*y+x]
89 // If (x,y) is outside interpolatable image region and NDEBUG is not defined
90 // the code will fail an ASSERT.
91 // The safe interpolatable region is [0,nx)*[0,ny)*[0,nz].
92 template <class T>
93 inline double
vil3d_trilin_interp_assert(double x,double y,double z,const T * data,unsigned nx,unsigned ny,unsigned nz,std::ptrdiff_t xstep,std::ptrdiff_t ystep,std::ptrdiff_t zstep)94 vil3d_trilin_interp_assert(double x, double y, double z, const T *data,
95 unsigned nx, unsigned ny, unsigned nz,
96 std::ptrdiff_t xstep, std::ptrdiff_t ystep,
97 std::ptrdiff_t zstep) {
98 assert(x>=0);
99 assert(y>=0);
100 assert(z>=0);
101 assert(x<nx-1);
102 assert(y<ny-1);
103 assert(z<nz-1);
104 return vil3d_trilin_interp_raw(x,y,z,data,xstep,ystep,zstep);
105 }
106
107 //: Compute trilinear interpolation at (x,y), with bounds checks.
108 // Image is nx * ny array of Ts. x,y element is data[nx*y+x]
109 // If (x,y,z) is outside safe interpolatable image region, nearest pixel value is returned.
110 template<class T>
vil3d_trilin_interp_safe_extend(double x,double y,double z,const T * data,unsigned nx,unsigned ny,unsigned nz,std::ptrdiff_t xstep,std::ptrdiff_t ystep,std::ptrdiff_t zstep)111 inline double vil3d_trilin_interp_safe_extend(double x, double y, double z, const T* data,
112 unsigned nx, unsigned ny, unsigned nz,
113 std::ptrdiff_t xstep, std::ptrdiff_t ystep, std::ptrdiff_t zstep)
114 {
115 if (x<0) x= 0.0;
116 if (y<0) y= 0.0;
117 if (z<0) z= 0.0;
118 if (x>=nx-1) x=(double)nx-1.00000001;
119 if (y>=ny-1) y=(double)ny-1.00000001;
120 if (z>=nz-1) z=(double)nz-1.00000001;
121 return vil3d_trilin_interp_raw(x,y,z,data,xstep,ystep,zstep);
122 }
123
124 //: Compute trilinear interpolation at (x,y), with bounds checks.
125 // Image is nx * ny array of Ts. x,y element is data[nx*y+x]
126 // If (x,y,z) is outside safe interpolatable image region, NA is returned.
127 template<class T>
vil3d_trilin_interp_safe_edgena(double x,double y,double z,const T * data,unsigned nx,unsigned ny,unsigned nz,std::ptrdiff_t xstep,std::ptrdiff_t ystep,std::ptrdiff_t zstep)128 inline double vil3d_trilin_interp_safe_edgena(double x, double y, double z, const T* data,
129 unsigned nx, unsigned ny, unsigned nz,
130 std::ptrdiff_t xstep, std::ptrdiff_t ystep, std::ptrdiff_t zstep)
131 {
132 if (x<0 || y<0 || z<0 ||
133 x>=nx-1 || y>=ny-1 || z>=nz-1) return vil_na(double());
134 return vil3d_trilin_interp_raw(x,y,z,data,xstep,ystep,zstep);
135 }
136 //: Compute trilinear interpolation at (x,y), using the nearest valid value if out of bounds
137 // If (x,y,z) is outside safe interpolatable image region, nearest pixel value is returned.
138 template<class T>
139 inline double vil3d_trilin_interp_safe_extend(const vil3d_image_view<T>& image,
140 double x, double y, double z,
141 unsigned p=0)
142 {
143 return vil3d_trilin_interp_safe_extend(x, y, z, &image(0,0,0,p),
144 image.ni(), image.nj(), image.nk(),
145 image.istep(), image.jstep(), image.kstep());
146 }
147
148
149 #endif // vil3d_trilin_interp_h_
150