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