1 /*
2 Gri - A language for scientific graphics programming
3 Copyright (C) 2008 Daniel Kelley
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; version 3 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License along
16 with this program; if not, write to the Free Software Foundation, Inc.,
17 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 */
19
20 #include <string>
21 #include <math.h>
22 #include <stdio.h>
23 #include "private.hh"
24 #include "gr.hh"
25 extern double _grMissingValue;/* defined in gr.c */
26 int number_good_xyz(const std::vector<double>& x,
27 const std::vector<double>& y,
28 const std::vector<double>& f,
29 int n);
30 #define root_of_2 1.414213562
31
32 int
number_good_xyz(const std::vector<double> & x,const std::vector<double> & y,const std::vector<double> & f,int n)33 number_good_xyz(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& f, int n)
34 {
35 int good = 0;
36 for (int i = 0; i < n; i++)
37 if (!gr_missingx(x[i]) && !gr_missingy(y[i]) && !gr_missing(f[i]))
38 good++;
39 return good;
40 }
41
42 /*
43 * gr_grid1: objective analysis of 2d field
44 *
45 * DESCRIPTION: Given f[i] defined at points (x[i],y[i]), where 0<=i<n, this
46 * routine uses objective analysis to get a predicted value fOut. The points
47 * f[i] used to calculate fOut come from a specified search rectangle xr wide
48 * by yr high.
49 *
50 * Also does boxcar gridding if 'method' is 0.
51 *
52 * Third method exists to .... better clean this doc up!
53 *
54 * If there are `neighbors' or more points in a rectangular neighborhood 2*xr
55 * wide and 2*yr high, then only the data inside that region are used, with x
56 * and y distances weighted by xr and yr. If there are fewer than `neighbors'
57 * points in the neighborhood, then the search rectangle is enlarged in area
58 * [xr and yr are multiplied by sqrt(2)] and the process starts over. This
59 * enlargement may be repeated up to 'enlargements' times, until either
60 * `neighbors' data are located in the search region, or all the data are
61 * used. If `enlargements' is set to a negative value, the enlargement
62 * process will be continued until all the data are used, if necessary.
63 *
64 * When all relevant points are located, the output is calculated using a
65 * weighting formula for the data within the rectangle:
66 *
67 * fOut = sum(f[i]*w[i])/sum(w[i]
68 *
69 * Various weightings 'w' are provided, each as a function of
70 * the scaled distance 'r = sqrt([(x1-x2)/xr]^2 + [(y1-y2)/yr]^2)'.
71 *
72 * method = 0: w = 1 (boxcar smoothing)
73 *
74 * method = 1: w = (2 - r^2) / (2 + r^2) (Levy + Brown method "1")
75 *
76 * method = 2: w = [(2 - r^2) / (2 + r^2)]^2 (Levy + Brown method "2")
77 *
78 * method = 3: w = [(2 - r^2) / (2 + r^2)]^4 (Levy + Brown method "3")
79 *
80 * REFERENCE: Levy & Brown [1986 Journal of
81 * Geophysical Research 91C4: 5153-5158]. Note that the present method
82 * extends for unequal weightings in the (x,y) direction, so that nonisotropic
83 * data can be handled.
84 *
85 * RETURN VALUE: The number of points finally used in the sum, or 0 if there was
86 * an error (if n<=0).
87 */
88 int
gr_grid1(const std::vector<double> & x,const std::vector<double> & y,const std::vector<double> & f,double x0,double y0,double xr,double yr,int method,unsigned int neighbors,int enlargements,double * fOut)89 gr_grid1(const std::vector<double> &x,
90 const std::vector<double> &y,
91 const std::vector<double> &f,
92 double x0, double y0,
93 double xr, double yr,
94 int method,
95 unsigned int neighbors,
96 int enlargements,
97 double *fOut)
98 {
99 unsigned int n = x.size();
100 double dx, dy;
101 int enlarge = 0;
102 unsigned int i;
103 unsigned int nn, num_in_rect;
104 double d2; /* squared distance */
105 double sumw; /* sum of weights */
106 double sumfw; /* sum of weighted values */
107 double w; /* weight of f[i] */
108 /* Check for obvious errors */
109 if (n <= 0 || neighbors == 0) {
110 *fOut = _grMissingValue;
111 return 0;
112 }
113 nn = number_good_xyz(x, y, f, n);
114 if (neighbors > nn)
115 neighbors = nn;
116 /*
117 * Search the rectangle, increasing its size if necessary.
118 */
119 do {
120 num_in_rect = 0;
121 sumw = sumfw = 0.0;
122 switch (method) {
123 case 0:
124 for (i = 0; i < n; i++) {
125 dx = x0 - x[i];
126 if (-xr <= dx && dx <= xr) {
127 dy = y0 - y[i];
128 if (-yr <= dy && dy <= yr) {
129 sumw += 1.0;
130 sumfw += f[i];
131 num_in_rect++;
132 }
133 }
134 }
135 break;
136 case 1:
137 for (i = 0; i < n; i++) {
138 dx = GRI_ABS(x0 - x[i]);
139 if (dx <= xr) {
140 dy = GRI_ABS(y0 - y[i]);
141 if (dy <= yr) {
142 dx /= xr;
143 dy /= yr;
144 d2 = dx * dx + dy * dy; /* note 0<=d2<=2 */
145 w = (2.0 - d2) / (2.0 + d2);
146 sumw += w;
147 sumfw += f[i] * w;
148 num_in_rect++;
149 }
150 }
151 }
152 break;
153 case 2:
154 for (i = 0; i < n; i++) {
155 dx = GRI_ABS(x0 - x[i]);
156 if (dx <= xr) {
157 dy = GRI_ABS(y0 - y[i]);
158 if (dy <= yr) {
159 dx /= xr;
160 dy /= yr;
161 d2 = dx * dx + dy * dy; /* note 0<=d2<=2 */
162 w = (2.0 - d2) / (2.0 + d2);
163 w *= w;
164 sumw += w;
165 sumfw += f[i] * w;
166 num_in_rect++;
167 }
168 }
169 }
170 break;
171 case 3:
172 for (i = 0; i < n; i++) {
173 dx = GRI_ABS(x0 - x[i]);
174 if (dx <= xr) {
175 dy = GRI_ABS(y0 - y[i]);
176 if (dy <= yr) {
177 dx /= xr;
178 dy /= yr;
179 d2 = dx * dx + dy * dy; /* note 0<=d2<=2 */
180 w = (2.0 - d2) / (2.0 + d2);
181 w *= w;
182 w *= w;
183 sumw += w;
184 sumfw += f[i] * w;
185 num_in_rect++;
186 }
187 }
188 }
189 break;
190 default:
191 return 0; /* unknown method! */
192 }
193 xr *= root_of_2;
194 yr *= root_of_2;
195 } while ((++enlarge <= enlargements || enlargements < 0)
196 && num_in_rect < neighbors);
197 if (num_in_rect > 0) {
198 *fOut = sumfw / sumw;
199 return num_in_rect;
200 } else {
201 *fOut = _grMissingValue;
202 return 0;
203 }
204 }
205