GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-df949ca623
interp.c
Go to the documentation of this file.
1 /*!
2  \file lib/raster/interp.c
3 
4  \brief Raster Library - Interpolation methods
5 
6  (C) 2001-2009,2013 GRASS Development Team
7 
8  This program is free software under the GNU General Public License
9  (>=v2). Read the file COPYING that comes with GRASS for details.
10 
11  \author Original author CERL
12  */
13 
14 #include <math.h>
15 #include <string.h>
16 
17 #include <grass/gis.h>
18 #include <grass/raster.h>
19 #include <grass/glocale.h>
20 
22 {
23  return u * (c1 - c0) + c0;
24 }
25 
26 DCELL Rast_interp_bilinear(double u, double v, DCELL c00, DCELL c01, DCELL c10,
27  DCELL c11)
28 {
29  DCELL c0 = Rast_interp_linear(u, c00, c01);
30  DCELL c1 = Rast_interp_linear(u, c10, c11);
31 
32  return Rast_interp_linear(v, c0, c1);
33 }
34 
35 DCELL Rast_interp_cubic(double u, DCELL c0, DCELL c1, DCELL c2, DCELL c3)
36 {
37  return (u * (u * (u * (c3 - 3 * c2 + 3 * c1 - c0) +
38  (-c3 + 4 * c2 - 5 * c1 + 2 * c0)) +
39  (c2 - c0)) +
40  2 * c1) /
41  2;
42 }
43 
44 DCELL Rast_interp_bicubic(double u, double v, DCELL c00, DCELL c01, DCELL c02,
45  DCELL c03, DCELL c10, DCELL c11, DCELL c12, DCELL c13,
46  DCELL c20, DCELL c21, DCELL c22, DCELL c23, DCELL c30,
47  DCELL c31, DCELL c32, DCELL c33)
48 {
49  DCELL c0 = Rast_interp_cubic(u, c00, c01, c02, c03);
50  DCELL c1 = Rast_interp_cubic(u, c10, c11, c12, c13);
51  DCELL c2 = Rast_interp_cubic(u, c20, c21, c22, c23);
52  DCELL c3 = Rast_interp_cubic(u, c30, c31, c32, c33);
53 
54  return Rast_interp_cubic(v, c0, c1, c2, c3);
55 }
56 
57 DCELL Rast_interp_lanczos(double u, double v, DCELL *c)
58 {
59  double uweight[5], vweight[5], d, d_pi;
60  double usum, vsum;
61  DCELL c0, c1, c2, c3, c4;
62  double sind, sincd1, sincd2;
63 
64  d_pi = u * M_PI;
65  sind = 2 * sin(d_pi);
66  sincd1 = sind * sin(d_pi / 2);
67  uweight[2] = (u == 0 ? 1 : sincd1 / (d_pi * d_pi));
68  usum = uweight[2];
69 
70  d = u + 2;
71  d_pi = d * M_PI;
72  if (d > 2)
73  uweight[0] = 0.;
74  else
75  uweight[0] = (d == 0 ? 1 : -sincd1 / (d_pi * d_pi));
76  usum += uweight[0];
77 
78  d = u + 1.;
79  d_pi = d * M_PI;
80  sincd2 = sind * sin(d_pi / 2);
81  uweight[1] = (d == 0 ? 1 : -sincd2 / (d_pi * d_pi));
82  usum += uweight[1];
83 
84  d = u - 1.;
85  d_pi = d * M_PI;
86  uweight[3] = (d == 0 ? 1 : sincd2 / (d_pi * d_pi));
87  usum += uweight[3];
88 
89  d = u - 2.;
90  d_pi = d * M_PI;
91  if (d < -2)
92  uweight[4] = 0.;
93  else
94  uweight[4] = (d == 0 ? 1 : -sincd1 / (d_pi * d_pi));
95  usum += uweight[4];
96 
97  d_pi = v * M_PI;
98  sind = 2 * sin(d_pi);
99  sincd1 = sind * sin(d_pi / 2);
100  vweight[2] = (v == 0 ? 1 : sincd1 / (d_pi * d_pi));
101  vsum = vweight[2];
102 
103  d = v + 2;
104  d_pi = d * M_PI;
105  if (d > 2)
106  vweight[0] = 0;
107  else
108  vweight[0] = (d == 0 ? 1 : -sincd1 / (d_pi * d_pi));
109  vsum += vweight[0];
110 
111  d = v + 1.;
112  d_pi = d * M_PI;
113  sincd2 = sind * sin(d_pi / 2);
114  vweight[1] = (d == 0 ? 1 : -sincd2 / (d_pi * d_pi));
115  vsum += vweight[1];
116 
117  d = v - 1.;
118  d_pi = d * M_PI;
119  vweight[3] = (d == 0 ? 1 : sincd2 / (d_pi * d_pi));
120  vsum += vweight[3];
121 
122  d = v - 2.;
123  d_pi = d * M_PI;
124  if (d < -2)
125  vweight[4] = 0;
126  else
127  vweight[4] = (d == 0 ? 1 : -sincd1 / (d_pi * d_pi));
128  vsum += vweight[4];
129 
130  c0 = (c[0] * uweight[0] + c[1] * uweight[1] + c[2] * uweight[2] +
131  c[3] * uweight[3] + c[4] * uweight[4]);
132  c1 = (c[5] * uweight[0] + c[6] * uweight[1] + c[7] * uweight[2] +
133  c[8] * uweight[3] + c[9] * uweight[4]);
134  c2 = (c[10] * uweight[0] + c[11] * uweight[1] + c[12] * uweight[2] +
135  c[13] * uweight[3] + c[14] * uweight[4]);
136  c3 = (c[15] * uweight[0] + c[16] * uweight[1] + c[17] * uweight[2] +
137  c[18] * uweight[3] + c[19] * uweight[4]);
138  c4 = (c[20] * uweight[0] + c[21] * uweight[1] + c[22] * uweight[2] +
139  c[23] * uweight[3] + c[24] * uweight[4]);
140 
141  return ((c0 * vweight[0] + c1 * vweight[1] + c2 * vweight[2] +
142  c3 * vweight[3] + c4 * vweight[4]) /
143  (usum * vsum));
144 }
145 
147  DCELL c3)
148 {
149  return (u * (u * (u * (c3 - 3 * c2 + 3 * c1 - c0) +
150  (3 * c2 - 6 * c1 + 3 * c0)) +
151  (3 * c2 - 3 * c0)) +
152  c2 + 4 * c1 + c0) /
153  6;
154 }
155 
156 DCELL Rast_interp_bicubic_bspline(double u, double v, DCELL c00, DCELL c01,
157  DCELL c02, DCELL c03, DCELL c10, DCELL c11,
158  DCELL c12, DCELL c13, DCELL c20, DCELL c21,
159  DCELL c22, DCELL c23, DCELL c30, DCELL c31,
160  DCELL c32, DCELL c33)
161 {
162  DCELL c0 = Rast_interp_cubic_bspline(u, c00, c01, c02, c03);
163  DCELL c1 = Rast_interp_cubic_bspline(u, c10, c11, c12, c13);
164  DCELL c2 = Rast_interp_cubic_bspline(u, c20, c21, c22, c23);
165  DCELL c3 = Rast_interp_cubic_bspline(u, c30, c31, c32, c33);
166 
167  return Rast_interp_cubic_bspline(v, c0, c1, c2, c3);
168 }
169 
170 /*!
171  \brief Get interpolation method from the option.
172 
173  Calls G_fatal_error() on unknown interpolation method.
174 
175  Supported methods:
176  - NEAREST
177  - BILINEAR
178  - CUBIC
179 
180  \code
181  int interp_method
182  struct Option *opt_method;
183 
184  opt_method = G_define_standard_option(G_OPT_R_INTERP_TYPE);
185 
186  if (G_parser(argc, argv))
187  exit(EXIT_FAILURE);
188 
189  interp_method = G_option_to_interp_type(opt_method);
190  \endcode
191 
192  \param option pointer to interpolation option
193 
194  \return interpolation method code
195  */
196 int Rast_option_to_interp_type(const struct Option *option)
197 {
198  int interp_type;
199 
200  interp_type = INTERP_UNKNOWN;
201  if (option->answer) {
202  if (strcmp(option->answer, "nearest") == 0) {
203  interp_type = INTERP_NEAREST;
204  }
205  else if (strcmp(option->answer, "bilinear") == 0) {
206  interp_type = INTERP_BILINEAR;
207  }
208  else if (strcmp(option->answer, "bicubic") == 0) {
209  interp_type = INTERP_BICUBIC;
210  }
211  }
212 
213  if (interp_type == INTERP_UNKNOWN)
214  G_fatal_error(_("Unknown interpolation method: %s"), option->answer);
215 
216  return interp_type;
217 }
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
double DCELL
Definition: gis.h:628
#define M_PI
Definition: gis.h:158
#define _(str)
Definition: glocale.h:10
int Rast_option_to_interp_type(const struct Option *option)
Get interpolation method from the option.
Definition: interp.c:196
DCELL Rast_interp_cubic(double u, DCELL c0, DCELL c1, DCELL c2, DCELL c3)
Definition: interp.c:35
DCELL Rast_interp_bilinear(double u, double v, DCELL c00, DCELL c01, DCELL c10, DCELL c11)
Definition: interp.c:26
DCELL Rast_interp_bicubic_bspline(double u, double v, DCELL c00, DCELL c01, DCELL c02, DCELL c03, DCELL c10, DCELL c11, DCELL c12, DCELL c13, DCELL c20, DCELL c21, DCELL c22, DCELL c23, DCELL c30, DCELL c31, DCELL c32, DCELL c33)
Definition: interp.c:156
DCELL Rast_interp_linear(double u, DCELL c0, DCELL c1)
Definition: interp.c:21
DCELL Rast_interp_cubic_bspline(double u, DCELL c0, DCELL c1, DCELL c2, DCELL c3)
Definition: interp.c:146
DCELL Rast_interp_lanczos(double u, double v, DCELL *c)
Definition: interp.c:57
DCELL Rast_interp_bicubic(double u, double v, DCELL c00, DCELL c01, DCELL c02, DCELL c03, DCELL c10, DCELL c11, DCELL c12, DCELL c13, DCELL c20, DCELL c21, DCELL c22, DCELL c23, DCELL c30, DCELL c31, DCELL c32, DCELL c33)
Definition: interp.c:44
#define INTERP_BILINEAR
Definition: raster.h:21
#define INTERP_NEAREST
Definition: raster.h:20
#define INTERP_UNKNOWN
Interpolation methods.
Definition: raster.h:19
#define INTERP_BICUBIC
Definition: raster.h:22
Structure that stores option information.
Definition: gis.h:556
char * answer
Definition: gis.h:570