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 <math.h>
21 #include <stdio.h>
22 #include "gr.hh"
23 #include "GMatrix.hh"
24
25
26 // gr_smootharray() -- smooth z(x,y) data
27 //
28 // DESCRIPTION Smooths data z[i][j] defined on a rectangular grid x[i] and y[j]
29 // (where 0<=i<nx and 0<=j<ny). The resultant smoothed field is put into
30 // zS[][]. The array legit[i][j] is 1 if z[i][j] is a good data point; 0
31 // otherwise. When computing the smoothed value at a point, neighboring
32 // points having legit[][]=0 are ignored. So are points outside the region
33 // 0<=i<nx and 0<=j<ny. The distance between grid points in the x direction
34 // is dx; in the y direction it's dy. dt indicates the strength of smoothing,
35 // as described below. Holes in the grid (i.e., missing values) will be
36 // smoothed over if fill_in is false.
37 //
38 // Returns true if all was OK; false if unknown method supplied.
39 //
40 // METHOD: 'method' specifies the smoothing algorithm as follows:
41 //
42 // method = 0
43 //
44 // Smoothing is by equal weights over points in a + shaped window about the
45 // datum: zS[i][j] = z[i][j]/2 + (z[i-1][j]+z[i+1][j])/8 +
46 // (z[i][j-1]+z[i][j+1])/8 The parameters dx, dy, and dt are ignored.
47 //
48 // method = 1
49 //
50 // Smoothing is across x: zS[i][j] = z[i][j]/2 + (z[i-1][j] + z[i+1][j]) / 4;
51 //
52 // method = 2
53 //
54 // Smoothing is across y: zS[i][j] = z[i][j]/2 + (z[i][j-1] + z[i][j+1]) / 4;
55 //
56 // BUGS: None known.
57 bool
gr_smootharray(double dx,double dy,double dt,GriMatrix<double> & z,GriMatrix<double> & zS,GriMatrix<bool> & legit,GriMatrix<bool> & legitS,int nx,int ny,int method)58 gr_smootharray(double dx,
59 double dy,
60 double dt,
61 GriMatrix<double> &z,
62 GriMatrix<double> &zS,
63 GriMatrix<bool> &legit,
64 GriMatrix<bool> &legitS,
65 int nx,
66 int ny,
67 int method)
68 {
69 register int i, j;
70 double sum;
71 int good; // number data used per gridpoint
72 int nx_1, ny_1;
73 // Test for errors
74 if (nx <= 0 || ny <= 0)
75 return false;
76 nx_1 = nx - 1;
77 ny_1 = ny - 1;
78 dx = dy = dt = 0.0; // Kill warning on non-use
79 switch (method) {
80 case 0:
81 for (i = 0; i < nx; i++) {
82 for (j = 0; j < ny; j++) {
83 good = 0;
84 sum = 0.0;
85 // Datum at centre
86 if (legit(i, j) == true) {
87 sum += 0.5 * z(i, j);
88 good++;
89 }
90 // Datum to left
91 if (i == 0) {
92 // Left edge: interpolate
93 if (legit(i + 1, j) == true
94 && legit(i, j) == true) {
95 sum += 0.125 * (2.0 * z(i, j) - z(i + 1, j));
96 good++;
97 }
98 } else {
99 // Interior
100 if (legit(i - 1, j) == true) {
101 sum += 0.125 * z(i - 1, j);
102 good++;
103 }
104 }
105 // Datum to right
106 if (i == nx_1) {
107 // Right edge: interpolate
108 if (legit(i - 1, j) == true
109 && legit(i, j) == true) {
110 sum += 0.125 * (2.0 * z(i, j) - z(i - 1, j));
111 good++;
112 }
113 } else {
114 // Interior
115 if (legit(i + 1, j) == true) {
116 sum += 0.125 * z(i + 1, j);
117 good++;
118 }
119 }
120 // Datum below
121 if (j == 0) {
122 // Bottom: interpolate
123 if (legit(i, j + 1) == true
124 && legit(i, j) == true) {
125 sum += 0.125 * (2.0 * z(i, j) - z(i, j + 1));
126 good++;
127 }
128 } else {
129 if (legit(i, j - 1) == true) {
130 sum += 0.125 * z(i, j - 1);
131 good++;
132 }
133 }
134 // Datum above
135 if (j == ny_1) {
136 // Top: interpolate
137 if (legit(i, j - 1) == true
138 && legit(i, j) == true) {
139 sum += 0.125 * (2.0 * z(i, j) - z(i, j - 1));
140 good++;
141 }
142 } else {
143 if (legit(i, j + 1) == true) {
144 sum += 0.125 * z(i, j + 1);
145 good++;
146 }
147 }
148 if (good == 5) {
149 zS(i, j) = sum;
150 legitS(i, j) = true;
151 } else {
152 zS(i, j) = z(i, j); // won't be used anyway
153 legitS(i, j) = false;
154 }
155 }
156 }
157 break;
158 case 1: // smooth across x
159 for (i = 0; i < nx; i++) {
160 for (j = 0; j < ny; j++) {
161 sum = 0.0;
162 good = 0;
163 // Datum at centre
164 if (legit(i, j) == true) {
165 sum += 0.5 * z(i, j);
166 good++;
167 }
168 // Datum to left
169 if (i == 0) {
170 // Left edge: interpolate
171 if (legit(i + 1, j) == true
172 && legit(i, j) == true) {
173 sum += 0.25 * (2.0 * z(i, j) - z(i + 1, j));
174 good++;
175 }
176 } else {
177 // Interior
178 if (legit(i - 1, j) == true) {
179 sum += 0.25 * z(i - 1, j);
180 good++;
181 }
182 }
183 // Datum to right
184 if (i == nx_1) {
185 // Right edge: interpolate
186 if (legit(i - 1, j) == true
187 && legit(i, j) == true) {
188 sum += 0.25 * (2.0 * z(i, j) - z(i - 1, j));
189 good++;
190 }
191 } else {
192 // Interior
193 if (legit(i + 1, j) == true) {
194 sum += 0.25 * z(i + 1, j);
195 good++;
196 }
197 }
198 if (good == 5) {
199 zS(i, j) = sum;
200 legitS(i, j) = true;
201 } else {
202 zS(i, j) = z(i, j); // won't be used anyway
203 legitS(i, j) = false;
204 }
205 }
206 }
207 break;
208 case 2: // smooth across y
209 for (i = 0; i < nx; i++) {
210 for (j = 0; j < ny; j++) {
211 sum = 0.0;
212 good = 0;
213 // Datum at centre
214 if (legit(i, j) == true) {
215 sum += 0.5 * z(i, j);
216 good++;
217 }
218 // Datum below
219 if (j == 0) {
220 // Bottom: interpolate
221 if (legit(i, j + 1) == true
222 && legit(i, j) == true) {
223 sum += 0.25 * (2.0 * z(i, j) - z(i, j + 1));
224 good++;
225 }
226 } else {
227 if (legit(i, j - 1) == true) {
228 sum += 0.25 * z(i, j - 1);
229 good++;
230 }
231 }
232 // Datum above
233 if (j == ny_1) {
234 // Top: interpolate
235 if (legit(i, j - 1) == true
236 && legit(i, j) == true) {
237 sum += 0.25 * (2.0 * z(i, j) - z(i, j - 1));
238 good++;
239 }
240 } else {
241 if (legit(i, j + 1) == true) {
242 sum += 0.25 * z(i, j + 1);
243 good++;
244 }
245 }
246 if (good == 5) {
247 zS(i, j) = sum;
248 legitS(i, j) = true;
249 } else {
250 zS(i, j) = z(i, j); // won't be used anyway
251 legitS(i, j) = false;
252 }
253 }
254 }
255 break;
256 default:
257 return false; // unknown
258 }
259 return true;
260 }
261