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
Rast_interp_linear(double u,DCELL c0,DCELL c1)21 DCELL Rast_interp_linear(double u, DCELL c0, DCELL c1)
22 {
23 return u * (c1 - c0) + c0;
24 }
25
Rast_interp_bilinear(double u,double v,DCELL c00,DCELL c01,DCELL c10,DCELL c11)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
Rast_interp_cubic(double u,DCELL c0,DCELL c1,DCELL c2,DCELL c3)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
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)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
Rast_interp_lanczos(double u,double v,DCELL * c)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
Rast_interp_cubic_bspline(double u,DCELL c0,DCELL c1,DCELL c2,DCELL c3)145 DCELL Rast_interp_cubic_bspline(double u, DCELL c0, DCELL c1, DCELL c2, DCELL c3)
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
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)153 DCELL Rast_interp_bicubic_bspline(double u, double v,
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 */
Rast_option_to_interp_type(const struct Option * option)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 }
215