GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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,
27  DCELL c00, DCELL c01, DCELL c10, 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)) + (c2 - c0)) +
39  2 * c1) / 2;
40 }
41 
42 DCELL Rast_interp_bicubic(double u, double v,
43  DCELL c00, DCELL c01, DCELL c02, DCELL c03,
44  DCELL c10, DCELL c11, DCELL c12, DCELL c13,
45  DCELL c20, DCELL c21, DCELL c22, DCELL c23,
46  DCELL c30, DCELL c31, DCELL c32, DCELL c33)
47 {
48  DCELL c0 = Rast_interp_cubic(u, c00, c01, c02, c03);
49  DCELL c1 = Rast_interp_cubic(u, c10, c11, c12, c13);
50  DCELL c2 = Rast_interp_cubic(u, c20, c21, c22, c23);
51  DCELL c3 = Rast_interp_cubic(u, c30, c31, c32, c33);
52 
53  return Rast_interp_cubic(v, c0, c1, c2, c3);
54 }
55 
56 DCELL Rast_interp_lanczos(double u, double v, DCELL *c)
57 {
58  double uweight[5], vweight[5], d, d_pi;
59  double usum, vsum;
60  DCELL c0, c1, c2, c3, c4;
61  double sind, sincd1, sincd2;
62 
63  d_pi = u * M_PI;
64  sind = 2 * sin(d_pi);
65  sincd1 = sind * sin(d_pi / 2);
66  uweight[2] = (u == 0 ? 1 : sincd1 / (d_pi * d_pi));
67  usum = uweight[2];
68 
69  d = u + 2;
70  d_pi = d * M_PI;
71  if (d > 2)
72  uweight[0] = 0.;
73  else
74  uweight[0] = (d == 0 ? 1 : -sincd1 / (d_pi * d_pi));
75  usum += uweight[0];
76 
77  d = u + 1.;
78  d_pi = d * M_PI;
79  sincd2 = sind * sin(d_pi / 2);
80  uweight[1] = (d == 0 ? 1 : -sincd2 / (d_pi * d_pi));
81  usum += uweight[1];
82 
83  d = u - 1.;
84  d_pi = d * M_PI;
85  uweight[3] = (d == 0 ? 1 : sincd2 / (d_pi * d_pi));
86  usum += uweight[3];
87 
88  d = u - 2.;
89  d_pi = d * M_PI;
90  if (d < -2)
91  uweight[4] = 0.;
92  else
93  uweight[4] = (d == 0 ? 1 : -sincd1 / (d_pi * d_pi));
94  usum += uweight[4];
95 
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]) / (usum * vsum));
143 }
144 
146 {
147  return (u * (u * (u * (c3 - 3 * c2 + 3 * c1 - c0) +
148  (3 * c2 - 6 * c1 + 3 * c0)) +
149  (3 * c2 - 3 * c0)) +
150  c2 + 4 * c1 + c0) / 6;
151 }
152 
154  DCELL c00, DCELL c01, DCELL c02, DCELL c03,
155  DCELL c10, DCELL c11, DCELL c12, DCELL c13,
156  DCELL c20, DCELL c21, DCELL c22, DCELL c23,
157  DCELL c30, DCELL c31, DCELL c32, DCELL c33)
158 {
159  DCELL c0 = Rast_interp_cubic_bspline(u, c00, c01, c02, c03);
160  DCELL c1 = Rast_interp_cubic_bspline(u, c10, c11, c12, c13);
161  DCELL c2 = Rast_interp_cubic_bspline(u, c20, c21, c22, c23);
162  DCELL c3 = Rast_interp_cubic_bspline(u, c30, c31, c32, c33);
163 
164  return Rast_interp_cubic_bspline(v, c0, c1, c2, c3);
165 }
166 
167 /*!
168  \brief Get interpolation method from the option.
169 
170  Calls G_fatal_error() on unknown interpolation method.
171 
172  Supported methods:
173  - NEAREST
174  - BILINEAR
175  - CUBIC
176 
177  \code
178  int interp_method
179  struct Option *opt_method;
180 
181  opt_method = G_define_standard_option(G_OPT_R_INTERP_TYPE);
182 
183  if (G_parser(argc, argv))
184  exit(EXIT_FAILURE);
185 
186  interp_method = G_option_to_interp_type(opt_method);
187  \endcode
188 
189  \param option pointer to interpolation option
190 
191  \return interpolation method code
192 */
193 int Rast_option_to_interp_type(const struct Option *option)
194 {
195  int interp_type;
196 
197  interp_type = INTERP_UNKNOWN;
198  if (option->answer) {
199  if (strcmp(option->answer, "nearest") == 0) {
200  interp_type = INTERP_NEAREST;
201  }
202  else if (strcmp(option->answer, "bilinear") == 0) {
203  interp_type = INTERP_BILINEAR;
204  }
205  else if (strcmp(option->answer, "bicubic") == 0) {
206  interp_type = INTERP_BICUBIC;
207  }
208  }
209 
210  if (interp_type == INTERP_UNKNOWN)
211  G_fatal_error(_("Unknown interpolation method: %s"), option->answer);
212 
213  return interp_type;
214 }
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
DCELL Rast_interp_lanczos(double u, double v, DCELL *c)
Definition: interp.c:56
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:153
DCELL Rast_interp_linear(double u, DCELL c0, DCELL c1)
Definition: interp.c:21
double DCELL
Definition: gis.h:614
#define M_PI
Definition: gis.h:144
int Rast_option_to_interp_type(const struct Option *option)
Get interpolation method from the option.
Definition: interp.c:193
DCELL Rast_interp_cubic_bspline(double u, DCELL c0, DCELL c1, DCELL c2, DCELL c3)
Definition: interp.c:145
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:42
DCELL Rast_interp_bilinear(double u, double v, DCELL c00, DCELL c01, DCELL c10, DCELL c11)
Definition: interp.c:26
char * answer
Definition: gis.h:555
#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:542
#define _(str)
Definition: glocale.h:10
#define INTERP_BILINEAR
Definition: raster.h:21
DCELL Rast_interp_cubic(double u, DCELL c0, DCELL c1, DCELL c2, DCELL c3)
Definition: interp.c:35