1 // This is core/vil/vil_bilin_interp.h
2 #ifndef vil_bilin_interp_h_
3 #define vil_bilin_interp_h_
4 //:
5 // \file
6 // \brief Bilinear interpolation functions for 2D images
7 // \author Tim Cootes
8 //
9 // The vil bicub source files were derived from the corresponding
10 // vil bilin files, thus the vil bilin/bicub source files are very
11 // similar. If you modify something in this file, there is a
12 // corresponding bicub file that would likely also benefit from
13 // the same change.
14
15 #include <cstddef>
16 #include <cassert>
17 #ifdef _MSC_VER
18 # include <vcl_msvc_warnings.h>
19 #endif
20 #include "vil_image_view.h"
21 #include "vil_na.h"
22
23 //: Compute bilinear interpolation at (x,y), no bound checks. Requires 0<x<ni-2, 0<y<nj-2
24 // Image is nx * ny array of Ts. x,y element is data[xstep*x+ystep*y]
25 // No bound checks are done.
26 template<class T>
vil_bilin_interp_unsafe(double x,double y,const T * data,std::ptrdiff_t xstep,std::ptrdiff_t ystep)27 inline double vil_bilin_interp_unsafe(double x, double y, const T* data,
28 std::ptrdiff_t xstep, std::ptrdiff_t ystep)
29 {
30 int p1x=int(x);
31 double normx = x-p1x;
32 int p1y=int(y);
33 double normy = y-p1y;
34
35 const T* pix1 = data + p1y*ystep + p1x*xstep;
36
37 double i1 = pix1[0 ]+(pix1[ ystep]-pix1[0 ])*normy;
38 double i2 = pix1[xstep]+(pix1[xstep+ystep]-pix1[xstep])*normy;
39
40 return i1+(i2-i1)*normx;
41 }
42
43
44 //: Compute bilinear interpolation at (x,y), no bound checks. Requires 0<x<ni-2, 0<y<nj-2
45 // Image is nx * ny array of Ts. x,y element is data[xstep*x+ystep*y]
46 // No bound checks are done.
47 // This is a version of vil_bilin_interp_unsafe with the same function
48 // signature as vil_bilin_interp_safe.
49 template<class T>
vil_bilin_interp_unsafe(double x,double y,const T * data,int,int,std::ptrdiff_t xstep,std::ptrdiff_t ystep)50 inline double vil_bilin_interp_unsafe(double x, double y, const T* data,
51 int /*nx*/, int /*ny*/,
52 std::ptrdiff_t xstep, std::ptrdiff_t ystep)
53 {
54 return vil_bilin_interp_unsafe(x, y, data, xstep, ystep);
55 }
56
57 //: Compute bilinear interpolation at (x,y), no bound checks
58 // Image is nx * ny array of Ts. x,y element is data[xstep*x+ystep*y]
59 // No bound checks are done.
60 template<class T>
vil_bilin_interp_raw(double x,double y,const T * data,std::ptrdiff_t xstep,std::ptrdiff_t ystep)61 inline double vil_bilin_interp_raw(double x, double y, const T* data,
62 std::ptrdiff_t xstep, std::ptrdiff_t ystep)
63 {
64 int p1x=int(x);
65 double normx = x-p1x;
66 int p1y=int(y);
67 double normy = y-p1y;
68
69 const T* pix1 = data + p1y*ystep + p1x*xstep;
70
71 // special boundary cases can be handled more quickly first;
72 // also avoids accessing an invalid pix1[t] which is going to have weight 0.
73 if (normx == 0 && normy == 0) return pix1[0];
74 if (normx == 0) return pix1[0]+(pix1[ystep]-pix1[0])*normy;
75 if (normy == 0) return pix1[0]+(pix1[xstep]-pix1[0])*normx;
76
77 double i1 = pix1[0 ]+(pix1[ ystep]-pix1[0 ])*normy;
78 double i2 = pix1[xstep]+(pix1[xstep+ystep]-pix1[xstep])*normy;
79
80 return i1+(i2-i1)*normx;
81 }
82
83 //: Compute bilinear interpolation at (x,y), no bound checks.
84 // Image is nx * ny array of Ts. x,y element is data[xstep*x+ystep*y]
85 // No bound checks are done.
86 // This is a version of vil_bilin_interp_raw with the same function
87 // signature as vil_bilin_interp_safe.
88 template<class T>
vil_bilin_interp_raw(double x,double y,const T * data,int,int,std::ptrdiff_t xstep,std::ptrdiff_t ystep)89 inline double vil_bilin_interp_raw(double x, double y, const T* data,
90 int /*nx*/, int /*ny*/,
91 std::ptrdiff_t xstep, std::ptrdiff_t ystep)
92 {
93 return vil_bilin_interp_raw(x, y, data, xstep, ystep);
94 }
95
96 //: Compute bilinear interpolation at (x,y), with bound checks
97 // Image is nx * ny array of Ts. x,y element is data[xstep*x+ystep*y]
98 // If (x,y) is outside interpolatable image region, zero is returned.
99 // The safe interpolatable region is [0,nx-1]*[0,ny-1].
100 template<class T>
vil_bilin_interp_safe(double x,double y,const T * data,int nx,int ny,std::ptrdiff_t xstep,std::ptrdiff_t ystep)101 inline double vil_bilin_interp_safe(double x, double y, const T* data,
102 int nx, int ny,
103 std::ptrdiff_t xstep, std::ptrdiff_t ystep)
104 {
105 if (x<0) return 0.0;
106 if (y<0) return 0.0;
107 if (x>nx-1) return 0.0;
108 if (y>ny-1) return 0.0;
109 return vil_bilin_interp_raw(x,y,data,xstep,ystep);
110 }
111
112
113 //: Compute bilinear interpolation at (x,y), with bound checks
114 // Image is nx * ny array of Ts. x,y element is data[xstep*x+ystep*y]
115 // If (x,y) is outside interpolatable image region, NA is returned.
116 // The safe interpolatable region is [0,nx-1]*[0,ny-1].
117 template<class T>
vil_bilin_interp_safe_edgena(double x,double y,const T * data,int nx,int ny,std::ptrdiff_t xstep,std::ptrdiff_t ystep)118 inline double vil_bilin_interp_safe_edgena(double x, double y, const T* data,
119 int nx, int ny,
120 std::ptrdiff_t xstep, std::ptrdiff_t ystep)
121 {
122 if (x<0 || y<0 || x>nx-1 || y>ny-1) return vil_na(double());
123 return vil_bilin_interp_raw(x,y,data,xstep,ystep);
124 }
125
126 //: Compute bilinear interpolation at (x,y), with bound checks
127 // If (x,y) is outside interpolatable image region, zero is returned.
128 // The safe interpolatable region is [0,view.ni()-1]*[0,view.nj()-1].
129 // \relatesalso vil_image_view
130 template<class T>
131 inline double vil_bilin_interp_safe(const vil_image_view<T> &view,
132 double x, double y, unsigned p=0)
133 {
134 return vil_bilin_interp_safe(x, y, &view(0,0,p),
135 view.ni(), view.nj(),
136 view.istep(), view.jstep());
137 }
138
139
140 //: Compute bilinear interpolation at (x,y), with bound checks
141 // If (x,y) is outside interpolatable image region, NA is returned.
142 // The safe interpolatable region is [0,view.ni()-1]*[0,view.nj()-1].
143 // \relatesalso vil_image_view
144 template<class T>
145 inline double vil_bilin_interp_safe_edgena(const vil_image_view<T> &view,
146 double x, double y, unsigned p=0)
147 {
148 return vil_bilin_interp_safe_edgena(x, y, &view(0,0,p),
149 view.ni(), view.nj(),
150 view.istep(), view.jstep());
151 }
152
153 //: Compute bilinear interpolation at (x,y), with minimal bound checks
154 // Image is nx * ny array of Ts. x,y element is data[ystep*y+xstep*x]
155 // If (x,y) is outside interpolatable image region and NDEBUG is not defined
156 // the code will fail an ASSERT.
157 // The safe interpolatable region is [0,nx-1]*[0,ny-1].
158 template <class T>
vil_bilin_interp(double x,double y,const T * data,int nx,int ny,std::ptrdiff_t xstep,std::ptrdiff_t ystep)159 inline double vil_bilin_interp(double x, double y, const T *data, int nx,
160 int ny, std::ptrdiff_t xstep,
161 std::ptrdiff_t ystep) {
162 assert (x>=0);
163 assert (y>=0);
164 assert (x<=nx-1);
165 assert (y<=ny-1);
166 return vil_bilin_interp_raw(x,y,data,xstep,ystep);
167 }
168
169 //: Compute bilinear interpolation at (x,y), with minimal bound checks
170 // If (x,y) is outside interpolatable image region and NDEBUG is not defined
171 // the code will fail an ASSERT.
172 // The safe interpolatable region is [0,view.ni()-1]*[0,view.nj()-1].
173 // \relatesalso vil_image_view
174 template<class T>
175 inline double vil_bilin_interp(const vil_image_view<T> &view,
176 double x, double y, unsigned p=0)
177 {
178 return vil_bilin_interp(x, y, &view(0,0,p),
179 view.ni(), view.nj(),
180 view.istep(), view.jstep());
181 }
182
183
184 //: Compute bilinear interpolation at (x,y), with bound checks
185 // Image is nx * ny array of Ts. x,y element is data[nx*y+x]
186 // If (x,y) is outside safe interpolatable image region, nearest pixel value is returned.
187 // The safe interpolatable region is [0,nx-1]*[0,ny-1].
188 template<class T>
vil_bilin_interp_safe_extend(double x,double y,const T * data,int nx,int ny,std::ptrdiff_t xstep,std::ptrdiff_t ystep)189 inline double vil_bilin_interp_safe_extend(double x, double y, const T* data,
190 int nx, int ny,
191 std::ptrdiff_t xstep, std::ptrdiff_t ystep)
192 {
193 if (x<0) x= 0.0;
194 if (y<0) y= 0.0;
195 if (x>nx-1) x=nx-1.0;
196 if (y>ny-1) y=ny-1.0;
197 return vil_bilin_interp_raw(x,y,data,xstep,ystep);
198 }
199
200 //: Compute bilinear interpolation at (x,y), with bound checks
201 // If (x,y) is outside safe interpolatable image region, nearest pixel value is returned.
202 // The safe interpolatable region is [0,view.ni()-1]*[0,view.nj()-1].
203 // \relatesalso vil_image_view
204 template<class T>
205 inline double vil_bilin_interp_safe_extend(const vil_image_view<T> &view,
206 double x, double y, unsigned p=0)
207 {
208 return vil_bilin_interp_safe_extend(x, y, &view(0,0,p),
209 view.ni(), view.nj(),
210 view.istep(), view.jstep());
211 }
212
213 #endif // vil_bilin_interp_h_
214