vil_bicub_interp.hxx
Go to the documentation of this file.
1 // core/vil/vil_bicub_interp.hxx
2 #ifndef vil_bicub_interp_hxx_
3 #define vil_bicub_interp_hxx_
4 //:
5 // \file
6 // \brief Bicubic interpolation functions for 2D images
7 //
8 // The vil bicub source files were derived from the corresponding
9 // vil bilin files, thus the vil bilin/bicub source files are very
10 // similar. If you modify something in this file, there is a
11 // corresponding bilin file that would likely also benefit from
12 // the same change.
13 //
14 // In this particular case, there is no corresponding
15 // vil_bilin_interp.hxx file, see vil_bilin_interp.h instead.
16 
17 #include "vil_bicub_interp.h"
18 #ifdef _MSC_VER
19 # include <vcl_msvc_warnings.h>
20 #endif
21 
22 // vil_bilin_interp.h defines only inline functions, but some of the
23 // corresponding vil_bicub_interp functions are a little big to be
24 // inline. Plus, on one platform, msvc 6.0 with /O2 optimization
25 // compiled the vil_bicub_interp functions without a peep but gave
26 // incorrect numerical results when these functions were inline and
27 // defined in vil_bicub_interp.h.
28 
29 template<class T>
30 double vil_bicub_interp_unsafe(double x, double y, const T* data,
31  std::ptrdiff_t xstep, std::ptrdiff_t ystep)
32 {
33  int p1x=int(x);
34  double normx = x-p1x;
35  int p1y=int(y);
36  double normy = y-p1y;
37 
38  const T* pix1 = data + p1y*ystep + p1x*xstep;
39 
40  // like bilinear interpolation, use separability.
41  // the s's are for the x-direction and the t's for the y-direction.
42  double s0 = ((2-normx)*normx-1)*normx; // -1
43  double s1 = (3*normx-5)*normx*normx+2; // 0
44  double s2 = ((4-3*normx)*normx+1)*normx; // +1
45  double s3 = (normx-1)*normx*normx; // +2
46 
47  double t0 = ((2-normy)*normy-1)*normy;
48  double t1 = (3*normy-5)*normy*normy+2;
49  double t2 = ((4-3*normy)*normy+1)*normy;
50  double t3 = (normy-1)*normy*normy;
51 
52 #define vil_I(dx,dy) (pix1[(dx)*xstep+(dy)*ystep])
53 
54  double xi0 = s0*vil_I(-1,-1) + s1*vil_I(+0,-1) + s2*vil_I(+1,-1) + s3*vil_I(+2,-1);
55  double xi1 = s0*vil_I(-1,+0) + s1*vil_I(+0,+0) + s2*vil_I(+1,+0) + s3*vil_I(+2,+0);
56  double xi2 = s0*vil_I(-1,+1) + s1*vil_I(+0,+1) + s2*vil_I(+1,+1) + s3*vil_I(+2,+1);
57  double xi3 = s0*vil_I(-1,+2) + s1*vil_I(+0,+2) + s2*vil_I(+1,+2) + s3*vil_I(+2,+2);
58 
59 #undef vil_I
60 
61  double val = 0.25 * ( xi0*t0 + xi1*t1 + xi2*t2 + xi3*t3 );
62 
63  return val;
64 }
65 
66 template<class T>
67 double vil_bicub_interp_raw(double x, double y, const T* data,
68  std::ptrdiff_t xstep, std::ptrdiff_t ystep)
69 {
70  int p1x=int(x);
71  double normx = x-p1x;
72  int p1y=int(y);
73  double normy = y-p1y;
74 
75  const T* pix1 = data + p1y*ystep + p1x*xstep;
76 
77  // special boundary cases can be handled more quickly first; also
78  // avoids accessing an invalid pix1[t] which is going to have
79  // weight 0.
80 
81  if (normx == 0.0 && normy == 0.0) return pix1[0];
82 
83  // coefficients for interpolation
84  double s0=-1.0, s1=-1.0, s2=-1.0, s3=-1.0; // in the x-direction
85  double t0=-1.0, t1=-1.0, t2=-1.0, t3=-1.0; // in the y-direction
86 
87  if (normx != 0.0) {
88  s0 = ((2-normx)*normx-1)*normx; // -1
89  s1 = (3*normx-5)*normx*normx+2; // 0
90  s2 = ((4-3*normx)*normx+1)*normx; // +1
91  s3 = (normx-1)*normx*normx; // +2
92  }
93 
94  if (normy != 0.0) {
95  t0 = ((2-normy)*normy-1)*normy; // -1
96  t1 = (3*normy-5)*normy*normy+2; // 0
97  t2 = ((4-3*normy)*normy+1)*normy; // +1
98  t3 = (normy-1)*normy*normy; // +2
99  }
100 
101 #define vil_I(dx,dy) (pix1[(dx)*xstep+(dy)*ystep])
102 
103  if (normy == 0.0) {
104  double val = 0.0;
105  val += s0*vil_I(-1,+0);
106  val += s1*vil_I(+0,+0);
107  val += s2*vil_I(+1,+0);
108 
109  val += s3*vil_I(+2,+0);
110  val *= 0.5;
111  return val;
112  }
113 
114  if (normx == 0.0) {
115  // The computation of 'val' in this section seems to compile
116  // fine, even when the very similar computation just above has
117  // trouble.
118  double val = t0*vil_I(+0,-1) + t1*vil_I(+0,+0) + t2*vil_I(+0,+1) + t3*vil_I(+0,+2);
119  val *= 0.5;
120  return val;
121  }
122 
123  double xi0 = s0*vil_I(-1,-1) + s1*vil_I(+0,-1) + s2*vil_I(+1,-1) + s3*vil_I(+2,-1);
124  double xi1 = s0*vil_I(-1,+0) + s1*vil_I(+0,+0) + s2*vil_I(+1,+0) + s3*vil_I(+2,+0);
125  double xi2 = s0*vil_I(-1,+1) + s1*vil_I(+0,+1) + s2*vil_I(+1,+1) + s3*vil_I(+2,+1);
126  double xi3 = s0*vil_I(-1,+2) + s1*vil_I(+0,+2) + s2*vil_I(+1,+2) + s3*vil_I(+2,+2);
127 
128 #undef vil_I
129 
130  double val = 0.25 * ( xi0*t0 + xi1*t1 + xi2*t2 + xi3*t3 );
131 
132  return val;
133 }
134 
135 #define VIL_BICUB_INTERP_INSTANTIATE(T) \
136 template double \
137 vil_bicub_interp_unsafe (double x, double y, const T* data, \
138  std::ptrdiff_t xstep, std::ptrdiff_t ystep); \
139 template double \
140 vil_bicub_interp_raw (double x, double y, const T* data, \
141  std::ptrdiff_t xstep, std::ptrdiff_t ystep)
142 
143 #endif // vil_bicub_interp_hxx_
vil_bicub_interp_raw
double vil_bicub_interp_raw(double x, double y, const T *data, std::ptrdiff_t xstep, std::ptrdiff_t ystep)
Compute bicubic interpolation at (x,y), no bound checks.
Definition: vil_bicub_interp.hxx:67
vil_bicub_interp.h
Bicubic interpolation functions for 2D images.
vil_bicub_interp_unsafe
double vil_bicub_interp_unsafe(double x, double y, const T *data, std::ptrdiff_t xstep, std::ptrdiff_t ystep)
Compute bicubic interpolation at (x,y), no bound checks. Requires 1<x<ni-3, 1<y<nj-3.
Definition: vil_bicub_interp.hxx:30
vil_I
#define vil_I(dx, dy)
Generated by   doxygen 1.8.15