1 /************************************************/
2 /* Test RSPL in 3D */
3 /************************************************/
4
5 /* Author: Graeme Gill
6 * Date: 22/4/96
7 * Derived from cmatch.c
8 * Copyright 1995 Graeme W. Gill
9 *
10 * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
11 * see the License.txt file for licencing details.
12 */
13
14
15 #define DEBUG
16 #define DETAILED
17
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <fcntl.h>
21 #include <math.h>
22 #include "rspl.h"
23 #include "tiffio.h"
24 #include "plot.h"
25 #include "ui.h"
26
27 #ifdef NEVER
28 #define INTERP spline_interp
29 #else
30 #define INTERP interp
31 #endif
32
33 #ifdef NEVER
34 FILE *verbose_out = stdout;
35 int verbose_level = 6; /* Current verbosity level */
36 /* 0 = none */
37 /* !0 = diagnostics */
38 #endif /* NEVER */
39
40 #define PLOTRES 256
41 #define WIDTH 400 /* Raster size */
42 #define HEIGHT 400
43
44 #define MAX_ITS 500
45 #define IT_TOL 0.0005
46 #define GRES0 33 /* Default rspl resolutions */
47 #define GRES1 33
48 #define GRES2 33
49 #undef NEVER
50 #define ALWAYS
51
52 //double t1xa[PNTS] = { 0.2, 0.25, 0.30, 0.35, 0.40, 0.44, 0.48, 0.51, 0.64, 0.75 };
53 //double t1ya[PNTS] = { 0.3, 0.35, 0.4, 0.41, 0.42, 0.46, 0.5, 0.575, 0.48, 0.75 };
54
55 /* 1D test function repeated 3 times along x = y = 0.5 */
56 co test_points1[] = {
57 {{ 0.4,0.4,0.20 },{ 0.30 }}, /* 0 */
58 {{ 0.4,0.4,0.25 },{ 0.35 }}, /* 1 */
59 {{ 0.4,0.4,0.30 },{ 0.40 }}, /* 2 */
60 {{ 0.4,0.4,0.35 },{ 0.41 }}, /* 3 */
61 {{ 0.4,0.4,0.40 },{ 0.42 }}, /* 4 */
62 {{ 0.4,0.4,0.44 },{ 0.46 }}, /* 5 */
63 {{ 0.4,0.4,0.48 },{ 0.50 }}, /* 6 */
64 {{ 0.4,0.4,0.51 },{ 0.575 }}, /* 7 */
65 {{ 0.4,0.4,0.64 },{ 0.48 }}, /* 8 */
66 {{ 0.4,0.4,0.75 },{ 0.75 }}, /* 9 */
67
68 {{ 0.5,0.4,0.20 },{ 0.30 }}, /* 0 */
69 {{ 0.5,0.4,0.25 },{ 0.35 }}, /* 1 */
70 {{ 0.5,0.4,0.30 },{ 0.40 }}, /* 2 */
71 {{ 0.5,0.4,0.35 },{ 0.41 }}, /* 3 */
72 {{ 0.5,0.4,0.40 },{ 0.42 }}, /* 4 */
73 {{ 0.5,0.4,0.44 },{ 0.46 }}, /* 5 */
74 {{ 0.5,0.4,0.48 },{ 0.50 }}, /* 6 */
75 {{ 0.5,0.4,0.51 },{ 0.575 }}, /* 7 */
76 {{ 0.5,0.4,0.64 },{ 0.48 }}, /* 8 */
77 {{ 0.5,0.4,0.75 },{ 0.75 }}, /* 9 */
78
79 {{ 0.6,0.4,0.20 },{ 0.30 }}, /* 0 */
80 {{ 0.6,0.4,0.25 },{ 0.35 }}, /* 1 */
81 {{ 0.6,0.4,0.30 },{ 0.40 }}, /* 2 */
82 {{ 0.6,0.4,0.35 },{ 0.41 }}, /* 3 */
83 {{ 0.6,0.4,0.40 },{ 0.42 }}, /* 4 */
84 {{ 0.6,0.4,0.44 },{ 0.46 }}, /* 5 */
85 {{ 0.6,0.4,0.48 },{ 0.50 }}, /* 6 */
86 {{ 0.6,0.4,0.51 },{ 0.575 }}, /* 7 */
87 {{ 0.6,0.4,0.64 },{ 0.48 }}, /* 8 */
88 {{ 0.6,0.4,0.75 },{ 0.75 }}, /* 9 */
89
90
91 {{ 0.4,0.5,0.20 },{ 0.30 }}, /* 0 */
92 {{ 0.4,0.5,0.25 },{ 0.35 }}, /* 1 */
93 {{ 0.4,0.5,0.30 },{ 0.40 }}, /* 2 */
94 {{ 0.4,0.5,0.35 },{ 0.41 }}, /* 3 */
95 {{ 0.4,0.5,0.40 },{ 0.42 }}, /* 4 */
96 {{ 0.4,0.5,0.44 },{ 0.46 }}, /* 5 */
97 {{ 0.4,0.5,0.48 },{ 0.50 }}, /* 6 */
98 {{ 0.4,0.5,0.51 },{ 0.575 }}, /* 7 */
99 {{ 0.4,0.5,0.64 },{ 0.48 }}, /* 8 */
100 {{ 0.4,0.5,0.75 },{ 0.75 }}, /* 9 */
101
102 {{ 0.5,0.5,0.20 },{ 0.30 }}, /* 0 */
103 {{ 0.5,0.5,0.25 },{ 0.35 }}, /* 1 */
104 {{ 0.5,0.5,0.30 },{ 0.40 }}, /* 2 */
105 {{ 0.5,0.5,0.35 },{ 0.41 }}, /* 3 */
106 {{ 0.5,0.5,0.40 },{ 0.42 }}, /* 4 */
107 {{ 0.5,0.5,0.44 },{ 0.46 }}, /* 5 */
108 {{ 0.5,0.5,0.48 },{ 0.50 }}, /* 6 */
109 {{ 0.5,0.5,0.51 },{ 0.575 }}, /* 7 */
110 {{ 0.5,0.5,0.64 },{ 0.48 }}, /* 8 */
111 {{ 0.5,0.5,0.75 },{ 0.75 }}, /* 9 */
112
113 {{ 0.6,0.5,0.20 },{ 0.30 }}, /* 0 */
114 {{ 0.6,0.5,0.25 },{ 0.35 }}, /* 1 */
115 {{ 0.6,0.5,0.30 },{ 0.40 }}, /* 2 */
116 {{ 0.6,0.5,0.35 },{ 0.41 }}, /* 3 */
117 {{ 0.6,0.5,0.40 },{ 0.42 }}, /* 4 */
118 {{ 0.6,0.5,0.44 },{ 0.46 }}, /* 5 */
119 {{ 0.6,0.5,0.48 },{ 0.50 }}, /* 6 */
120 {{ 0.6,0.5,0.51 },{ 0.575 }}, /* 7 */
121 {{ 0.6,0.5,0.64 },{ 0.48 }}, /* 8 */
122 {{ 0.6,0.5,0.75 },{ 0.75 }}, /* 9 */
123
124
125 {{ 0.4,0.6,0.20 },{ 0.30 }}, /* 0 */
126 {{ 0.4,0.6,0.25 },{ 0.35 }}, /* 1 */
127 {{ 0.4,0.6,0.30 },{ 0.40 }}, /* 2 */
128 {{ 0.4,0.6,0.35 },{ 0.41 }}, /* 3 */
129 {{ 0.4,0.6,0.40 },{ 0.42 }}, /* 4 */
130 {{ 0.4,0.6,0.44 },{ 0.46 }}, /* 5 */
131 {{ 0.4,0.6,0.48 },{ 0.50 }}, /* 6 */
132 {{ 0.4,0.6,0.51 },{ 0.575 }}, /* 7 */
133 {{ 0.4,0.6,0.64 },{ 0.48 }}, /* 8 */
134 {{ 0.4,0.6,0.75 },{ 0.75 }}, /* 9 */
135
136 {{ 0.5,0.6,0.20 },{ 0.30 }}, /* 0 */
137 {{ 0.5,0.6,0.25 },{ 0.35 }}, /* 1 */
138 {{ 0.5,0.6,0.30 },{ 0.40 }}, /* 2 */
139 {{ 0.5,0.6,0.35 },{ 0.41 }}, /* 3 */
140 {{ 0.5,0.6,0.40 },{ 0.42 }}, /* 4 */
141 {{ 0.5,0.6,0.44 },{ 0.46 }}, /* 5 */
142 {{ 0.5,0.6,0.48 },{ 0.50 }}, /* 6 */
143 {{ 0.5,0.6,0.51 },{ 0.575 }}, /* 7 */
144 {{ 0.5,0.6,0.64 },{ 0.48 }}, /* 8 */
145 {{ 0.5,0.6,0.75 },{ 0.75 }}, /* 9 */
146
147 {{ 0.6,0.6,0.20 },{ 0.30 }}, /* 0 */
148 {{ 0.6,0.6,0.25 },{ 0.35 }}, /* 1 */
149 {{ 0.6,0.6,0.30 },{ 0.40 }}, /* 2 */
150 {{ 0.6,0.6,0.35 },{ 0.41 }}, /* 3 */
151 {{ 0.6,0.6,0.40 },{ 0.42 }}, /* 4 */
152 {{ 0.6,0.6,0.44 },{ 0.46 }}, /* 5 */
153 {{ 0.6,0.6,0.48 },{ 0.50 }}, /* 6 */
154 {{ 0.6,0.6,0.51 },{ 0.575 }}, /* 7 */
155 {{ 0.6,0.6,0.64 },{ 0.48 }}, /* 8 */
156 {{ 0.6,0.6,0.75 },{ 0.75 }} /* 9 */
157 };
158
159 /* function */
160 co test_points2[] = {
161 {{ 0.50915, 0.50936, 0.57048 },{ 0.36209169635 }},
162 {{ 0.85943, 0.84331, 0.81487 },{ 0.91571493013 }},
163 {{ 0.11381, 0.80378, 0.82951 },{ 0.16052707023 }},
164 {{ 0.79087, 0.11157, 0.83913 },{ 0.30641388193 }},
165 {{ 0.16297, 0.23090, 0.96417 },{ 0.15005479047 }},
166 {{ 0.78181, 0.80097, 0.00192 },{ 0.32179402798 }},
167 {{ 0.16141, 0.84321, 0.16561 },{ 0.14446082013 }},
168 {{ 0.79859, 0.11111, 0.15547 },{ -0.1308162293 }},
169 {{ 0.12959, 0.16184, 0.21825 },{ 0.03520247555 }},
170 {{ 1.00000, 0.46395, 0.84399 },{ 0.63914200000 }},
171 {{ 0.48724, 0.76024, 1.00000 },{ 0.64150992880 }},
172 {{ 0.40744, 0.00000, 0.65533 },{ 0.13060244736 }},
173 {{ 0.10931, 0.43216, 0.25195 },{ 0.06329759515 }},
174 {{ 0.69401, 0.99412, 0.55335 },{ 0.75632862795 }},
175 {{ 0.51759, 0.30372, 0.13622 },{ 0.07965761859 }},
176 {{ 0.98628, 0.40857, 0.47688 },{ 0.29286006552 }},
177 {{ 0.44474, 0.26024, 1.00000 },{ 0.37263430380 }},
178 {{ 0.15229, 0.52694, 0.75101 },{ 0.16014862087 }},
179 {{ 0.50968, 0.74522, 0.14058 },{ 0.30725752992 }},
180 {{ 0.82291, 0.10614, 0.49956 },{ 0.07762756903 }},
181 {{ 0.15674, 0.78167, 0.50766 },{ 0.17389174472 }},
182 {{ 0.12961, 0.17234, 0.65990 },{ 0.08236132255 }},
183 {{ 1.00000, 0.83202, 0.35206 },{ 0.61366800000 }},
184 {{ 0.36633, 0.89555, 0.76014 },{ 0.48373766601 }},
185 {{ 0.85641, 0.39194, 0.18691 },{ 0.09699956583 }},
186 {{ 0.50918, 0.23923, 0.42132 },{ 0.16380116928 }},
187 {{ 0.65513, 0.40585, 0.83832 },{ 0.49065371733 }},
188 {{ 0.66726, 0.75472, 0.25419 },{ 0.41666516892 }},
189 {{ 0.22193, 0.58555, 0.13920 },{ 0.13003877385 }},
190 {{ 0.41507, 0.89581, 0.25539 },{ 0.37048608609 }},
191 {{ 0.43231, 0.12362, 0.17297 },{ 0.01981752271 }},
192 {{ 0.78191, 0.63297, 0.72639 },{ 0.64361123257 }}
193 };
194
195 co test_points3[] = {
196 {{ 0.50915, 0.50936, 0.57048 },{ -0.00010 }},
197 {{ 0.85943, 0.84331, 0.81487 },{ 0.46573 }},
198 {{ 0.11381, 0.80378, 0.82951 },{ 0.56085 }},
199 {{ 0.79087, 0.11157, 0.83913 },{ 0.33378 }},
200 {{ 0.16297, 0.23090, 0.96417 },{ 0.38872 }},
201 {{ 0.78181, 0.80097, 0.00192 },{ 0.28654 }},
202 {{ 0.16141, 0.84321, 0.16561 },{ 0.43410 }},
203 {{ 0.79859, 0.11111, 0.15547 },{ 0.35140 }},
204 {{ 0.12959, 0.16184, 0.21825 },{ 0.46711 }},
205 {{ 1.00000, 0.46395, 0.84399 },{ 0.94213 }},
206 {{ 0.48724, 0.76024, 1.00000 },{ 0.00058 }},
207 {{ 0.40744, 0.00000, 0.65533 },{ 0.00859 }},
208 {{ 0.10931, 0.43216, 0.25195 },{ 0.54679 }},
209 {{ 0.69401, 0.99412, 0.55335 },{ 0.12494 }},
210 {{ 0.51759, 0.30372, 0.13622 },{ -0.00126 }},
211 {{ 0.98628, 0.40857, 0.47688 },{ 0.89573 }},
212 {{ 0.44474, 0.26024, 1.00000 },{ 0.00280 }},
213 {{ 0.15229, 0.52694, 0.75101 },{ 0.43807 }},
214 {{ 0.50968, 0.74522, 0.14058 },{ 0.00004 }},
215 {{ 0.82291, 0.10614, 0.49956 },{ 0.40990 }},
216 {{ 0.15674, 0.78167, 0.50766 },{ 0.44271 }},
217 {{ 0.12961, 0.17234, 0.65990 },{ 0.46804 }},
218 {{ 1.00000, 0.83202, 0.35206 },{ 0.90938 }},
219 {{ 0.36633, 0.89555, 0.76014 },{ 0.07020 }},
220 {{ 0.85641, 0.39194, 0.18691 },{ 0.48559 }},
221 {{ 0.50918, 0.23923, 0.42132 },{ -0.00391 }},
222 {{ 0.65513, 0.40585, 0.83832 },{ 0.09449 }},
223 {{ 0.66726, 0.75472, 0.25419 },{ 0.10065 }},
224 {{ 0.22193, 0.58555, 0.13920 },{ 0.28186 }},
225 {{ 0.41507, 0.89581, 0.25539 },{ 0.02877 }},
226 {{ 0.43231, 0.12362, 0.17297 },{ 0.00237 }},
227 {{ 0.78191, 0.63297, 0.72639 },{ 0.29573 }}
228 };
229
230
231 /* x + y^2 + z^1/3 function with one non-monotonic point */
232 co test_points4[] = {
233 {{ 0.1,0.1,0.5 },{ 0.11 }}, /* 0 */
234 {{ 0.2,0.7,0.1 },{ 0.69 }}, /* 1 */
235 {{ 0.8,0.8,0.8 },{ 1.44 }}, /* 2 */
236 {{ 0.5,0.6,0.4 },{ 0.86 }}, /* 3 */
237 {{ 0.2,0.5,0.2 },{ 0.45 }}, /* 4 */
238 {{ 0.3,0.7,0.2 },{ 0.35 }}, /* nm 5 */
239 {{ 0.5,0.4,0.9 },{ 0.66 }}, /* 6 */
240 {{ 0.1,0.9,0.7 },{ 0.91 }}, /* 7 */
241 {{ 0.7,0.2,0.1 },{ 0.74 }}, /* 8 */
242 {{ 0.8,0.4,0.3 },{ 0.96 }}, /* 9 */
243 {{ 0.3,0.3,0.4 },{ 0.39 }} /* 10 */
244 };
245
246 /* doubled up x + y^2 + z^1/3 function with one non-monotonic point */
247 co test_points5[] = {
248 {{ 0.1,0.1,0.5 },{ 0.11 }}, /* 0 */
249 {{ 0.101,0.101,0.501 },{ 0.11 }}, /* 0d */
250 {{ 0.2,0.7,0.1 },{ 0.69 }}, /* 1 */
251 {{ 0.201,0.701,0.101 },{ 0.69 }}, /* 1d */
252 {{ 0.8,0.8,0.8 },{ 1.44 }}, /* 2 */
253 {{ 0.801,0.801,0.801 },{ 1.44 }}, /* 2d */
254 {{ 0.5,0.6,0.4 },{ 0.86 }}, /* 3 */
255 {{ 0.501,0.601,0.401 },{ 0.86 }}, /* 3d */
256 {{ 0.2,0.5,0.2 },{ 0.45 }}, /* 4 */
257 {{ 0.201,0.501,0.201 },{ 0.45 }}, /* 4d */
258 {{ 0.3,0.7,0.2 },{ 0.35 }}, /* nm 5 */
259 {{ 0.301,0.701,0.201 },{ 0.35 }}, /* nm 5d */
260 {{ 0.5,0.4,0.9 },{ 0.66 }}, /* 6 */
261 {{ 0.501,0.401,0.901 },{ 0.66 }}, /* 6d */
262 {{ 0.1,0.9,0.7 },{ 0.91 }}, /* 7 */
263 {{ 0.101,0.901,0.701 },{ 0.91 }}, /* 7d */
264 {{ 0.7,0.2,0.1 },{ 0.74 }}, /* 8 */
265 {{ 0.701,0.201,0.101 },{ 0.74 }}, /* 8d */
266 {{ 0.8,0.4,0.3 },{ 0.96 }}, /* 9 */
267 {{ 0.801,0.401,0.301 },{ 0.96 }}, /* 9d */
268 {{ 0.3,0.3,0.4 },{ 0.39 }}, /* 10 */
269 {{ 0.301,0.301,0.401 },{ 0.39 }} /* 10d */
270 };
271
272 co test_points6[] = {
273 {{ 0.0069, 0.0071, 0.0061 },{ 0.0726 }},
274 {{ 0.0068, 0.0071, 0.0060 },{ 0.0704 }},
275 {{ 0.0069, 0.0072, 0.0062 },{ 0.0720 }},
276 {{ 0.0069, 0.0072, 0.0061 },{ 0.0734 }},
277 {{ 0.0069, 0.0072, 0.0063 },{ 0.0750 }},
278 {{ 0.0070, 0.0072, 0.0062 },{ 0.0779 }},
279 {{ 0.0070, 0.0072, 0.0063 },{ 0.0741 }},
280 {{ 0.0069, 0.0072, 0.0061 },{ 0.0745 }},
281 {{ 0.0069, 0.0072, 0.0061 },{ 0.0747 }},
282 {{ 0.0071, 0.0073, 0.0063 },{ 0.0760 }},
283 {{ 0.0070, 0.0073, 0.0063 },{ 0.0751 }},
284 {{ 0.0070, 0.0073, 0.0062 },{ 0.0759 }},
285 {{ 0.0071, 0.0074, 0.0062 },{ 0.0693 }},
286 {{ 0.0071, 0.0074, 0.0064 },{ 0.0740 }},
287 {{ 0.0072, 0.0075, 0.0064 },{ 0.0741 }},
288 {{ 0.0199, 0.0209, 0.0184 },{ 0.1019 }},
289 {{ 0.0296, 0.0306, 0.0257 },{ 0.1213 }},
290 {{ 0.0627, 0.0651, 0.0548 },{ 0.1779 }},
291 {{ 0.0831, 0.0863, 0.0718 },{ 0.2095 }},
292 {{ 0.1091, 0.1134, 0.0946 },{ 0.2487 }},
293 {{ 0.1442, 0.1497, 0.1227 },{ 0.2949 }},
294 {{ 0.1745, 0.1814, 0.1495 },{ 0.3360 }},
295 {{ 0.1747, 0.1816, 0.1498 },{ 0.3367 }},
296 {{ 0.1747, 0.1816, 0.1496 },{ 0.3364 }},
297 {{ 0.1748, 0.1816, 0.1497 },{ 0.3355 }},
298 {{ 0.1749, 0.1817, 0.1497 },{ 0.3344 }},
299 {{ 0.1748, 0.1817, 0.1498 },{ 0.3356 }},
300 {{ 0.1748, 0.1817, 0.1498 },{ 0.3354 }},
301 {{ 0.1749, 0.1817, 0.1496 },{ 0.3361 }},
302 {{ 0.1749, 0.1818, 0.1498 },{ 0.3368 }},
303 {{ 0.1749, 0.1818, 0.1498 },{ 0.3335 }},
304 {{ 0.1750, 0.1818, 0.1499 },{ 0.3367 }},
305 {{ 0.1750, 0.1819, 0.1500 },{ 0.3362 }},
306 {{ 0.1750, 0.1819, 0.1498 },{ 0.3359 }},
307 {{ 0.1751, 0.1820, 0.1500 },{ 0.3354 }},
308 {{ 0.1752, 0.1821, 0.1501 },{ 0.3355 }},
309 {{ 0.1754, 0.1823, 0.1502 },{ 0.3369 }},
310 {{ 0.1756, 0.1824, 0.1504 },{ 0.3360 }},
311 {{ 0.2743, 0.2842, 0.2367 },{ 0.4381 }},
312 {{ 0.3289, 0.3411, 0.2834 },{ 0.4922 }},
313 {{ 0.4036, 0.4184, 0.3475 },{ 0.5617 }},
314 {{ 0.4689, 0.4854, 0.4020 },{ 0.6147 }},
315 {{ 0.5379, 0.5567, 0.4606 },{ 0.6709 }},
316 {{ 0.7137, 0.7420, 0.6169 },{ 0.8045 }},
317 {{ 0.8730, 0.9105, 0.7433 },{ 0.9150 }},
318 {{ 0.8738, 0.9113, 0.7435 },{ 0.9141 }},
319 {{ 0.8741, 0.9116, 0.7445 },{ 0.9120 }},
320 {{ 0.8744, 0.9118, 0.7443 },{ 0.9173 }},
321 {{ 0.8748, 0.9123, 0.7457 },{ 0.9219 }},
322 {{ 0.8748, 0.9123, 0.7450 },{ 0.9133 }},
323 {{ 0.8748, 0.9124, 0.7445 },{ 0.9210 }},
324 {{ 0.8751, 0.9127, 0.7462 },{ 0.9207 }},
325 {{ 0.8751, 0.9127, 0.7457 },{ 0.9225 }},
326 {{ 0.8754, 0.9130, 0.7454 },{ 0.9137 }},
327 {{ 0.8757, 0.9133, 0.7456 },{ 0.9219 }},
328 {{ 0.8759, 0.9135, 0.7470 },{ 0.9166 }},
329 {{ 0.8761, 0.9137, 0.7469 },{ 0.9162 }},
330 {{ 0.8759, 0.9137, 0.7469 },{ 0.9151 }},
331 {{ 0.8765, 0.9141, 0.7470 },{ 0.9167 }},
332 };
333
334
335 #ifdef NEVER
336 #ifdef __STDC__
337 #include <stdarg.h>
338 void error(char *fmt, ...), warning(char *fmt, ...), verbose(int level, char *fmt, ...);
339 #else
340 #include <varargs.h>
341 void error(), warning(), verbose();
342 #endif
343 #endif /* NEVER */
344
345 void write_rgb_tiff(char *name, int width, int height, unsigned char *data);
346
usage(void)347 void usage(void) {
348 fprintf(stderr,"Test 3D rspl interpolation\n");
349 fprintf(stderr,"Author: Graeme W. Gill\n");
350 fprintf(stderr,"usage: t2d [options]\n");
351 fprintf(stderr," -t n Test set:\n");
352 fprintf(stderr," * 1 = 1D test set along x = y = 0.5\n");
353 fprintf(stderr," 2 = Test set 1\n");
354 fprintf(stderr," 3 = Test set 2\n");
355 fprintf(stderr," 4 = x + y^2 + z^1/3 with nonmon point\n");
356 fprintf(stderr," 5 = doubled up x + y^2 + z^1/3 with nonmon point\n");
357 fprintf(stderr," 6 = neutral axis extrapolation\n");
358 fprintf(stderr," -r resx,resy,resz Set grid resolutions (def %d %d %d)\n",GRES0,GRES1,GRES2);
359 fprintf(stderr," -h Test half scale resolution too\n");
360 fprintf(stderr," -q Test quarter scale resolution too\n");
361 fprintf(stderr," -x Use auto smoothing\n");
362 fprintf(stderr," -s Test symetric smoothness\n");
363 fprintf(stderr," -p plot 4 slices, xy = 0.5, yz = 0.5, xz = 0.5, x=y=z\n");
364 fprintf(stderr," -P x1:y1:z1:x2:y2:z2 plot slice from x1,y1,z1,x2,y2,z2\n");
365 fprintf(stderr," -S factor smoothing factor (default 1.0)\n");
366 exit(1);
367 }
368
main(int argc,char * argv[])369 int main(int argc, char *argv[]) {
370 int fa,nfa; /* argument we're looking at */
371 rspl *rss; /* Regularized spline structure */
372 rspl *rss2 = NULL; /* Regularized spline structure at half/quarter resolution */
373 datai low,high;
374 int gres[MXDI];
375 int gres2[MXDI];
376 double avgdev[MXDO];
377 co *test_points = test_points1;
378 int npoints = sizeof(test_points1)/sizeof(co);
379 int dosym = 0;
380 int autosm = 0;
381 int doplot = 0;
382 double plotpts[2][3]; /* doplot == 2 start/end points */
383 int doh = 0; /* half scale */
384 int doq = 0;
385 int rsv;
386 double smoothf = 1.0;
387 int flags = RSPL_NOFLAGS;
388
389 low[0] = 0.0;
390 low[1] = 0.0;
391 low[2] = 0.0;
392 high[0] = 1.0;
393 high[1] = 1.0;
394 high[2] = 1.0;
395 gres[0] = GRES0;
396 gres[1] = GRES1;
397 gres[2] = GRES2;
398 avgdev[0] = 0.0;
399 avgdev[1] = 0.0;
400 avgdev[2] = 0.0;
401
402 /* Process the arguments */
403 for(fa = 1;fa < argc;fa++) {
404 nfa = fa; /* skip to nfa if next argument is used */
405 if (argv[fa][0] == '-') { /* Look for any flags */
406 char *na = NULL; /* next argument after flag, null if none */
407
408 if (argv[fa][2] != '\000')
409 na = &argv[fa][2]; /* next is directly after flag */
410 else {
411 if ((fa+1) < argc) {
412 if (argv[fa+1][0] != '-') {
413 nfa = fa + 1;
414 na = argv[nfa]; /* next is seperate non-flag argument */
415 }
416 }
417 }
418
419 if (argv[fa][1] == '?')
420 usage();
421
422 /* test set */
423 else if (argv[fa][1] == 't' || argv[fa][1] == 'T') {
424 int ix;
425 fa = nfa;
426 if (na == NULL) usage();
427 ix = atoi(na);
428 switch (ix) {
429 case 1:
430 test_points = test_points1;
431 npoints = sizeof(test_points1)/sizeof(co);
432 break;
433 case 2:
434 test_points = test_points2;
435 npoints = sizeof(test_points2)/sizeof(co);
436 break;
437 case 3:
438 test_points = test_points3;
439 npoints = sizeof(test_points3)/sizeof(co);
440 break;
441 case 4:
442 test_points = test_points4;
443 npoints = sizeof(test_points4)/sizeof(co);
444 break;
445 case 5:
446 test_points = test_points5;
447 npoints = sizeof(test_points5)/sizeof(co);
448 break;
449 case 6:
450 test_points = test_points6;
451 npoints = sizeof(test_points6)/sizeof(co);
452 break;
453 default:
454 usage();
455 }
456
457 } else if (argv[fa][1] == 'r' || argv[fa][1] == 'R') {
458 fa = nfa;
459 if (na == NULL) usage();
460 if (sscanf(na, " %d,%d,%d ", &gres[0], &gres[1], &gres[2]) != 3)
461 usage();
462
463 } else if (argv[fa][1] == 'h' || argv[fa][1] == 'H') {
464 doh = 1;
465
466 } else if (argv[fa][1] == 'q' || argv[fa][1] == 'Q') {
467 doh = 1;
468 doq = 1;
469
470 } else if (argv[fa][1] == 'p') {
471 doplot = 1;
472
473 } else if (argv[fa][1] == 'P') {
474 doplot = 2;
475 fa = nfa;
476 if (na == NULL) usage();
477 if (sscanf(na,"%lf:%lf:%lf:%lf:%lf:%lf",&plotpts[0][0],&plotpts[0][1],&plotpts[0][2],&plotpts[1][0],&plotpts[1][1],&plotpts[1][2]) != 6) {
478 usage();
479 }
480
481 } else if (argv[fa][1] == 's') {
482 dosym = 1;
483
484 } else if (argv[fa][1] == 'x') {
485 autosm = 1;
486
487 /* smoothing factor */
488 } else if (argv[fa][1] == 'S') {
489 int ix;
490 fa = nfa;
491 if (na == NULL) usage();
492 smoothf = atof(na);
493
494 } else
495 usage();
496 } else
497 break;
498 }
499
500
501 if (autosm)
502 flags |= RSPL_AUTOSMOOTH;
503
504 if (dosym)
505 flags |= RSPL_SYMDOMAIN;
506
507
508 /* Create the object */
509 rss = new_rspl(RSPL_NOFLAGS, 3, 1);
510
511 /* Fit to scattered data */
512 rss->fit_rspl(rss,
513 flags, /* Non-mon and clip flags */
514 test_points, /* Test points */
515 npoints, /* Number of test points */
516 low, high, gres, /* Low, high, resolution of grid */
517 NULL, NULL, /* Default data scale */
518 smoothf, /* Smoothing */
519 avgdev, /* Average deviation */
520 NULL); /* iwidth */
521 if (doh) {
522
523 if (doq) {
524 gres2[0] = gres[0]/4;
525 gres2[1] = gres[1]/4;
526 gres2[2] = gres[2]/4;
527 } else {
528 gres2[0] = gres[0]/2;
529 gres2[1] = gres[1]/2;
530 gres2[2] = gres[2]/2;
531 }
532
533 rss2 = new_rspl(RSPL_NOFLAGS, 3, 1);
534
535 /* Fit to scattered data */
536 rss2->fit_rspl(rss2,
537 flags, /* Non-mon and clip flags */
538 test_points, /* Test points */
539 npoints, /* Number of test points */
540 low, high, gres2, /* Low, high, resolution of grid */
541 NULL, NULL, /* Default data scale */
542 smoothf, /* Smoothing */
543 avgdev, /* Average deviation */
544 NULL); /* iwidth */
545 }
546
547 /* Test the interpolation with a slice in 2D */
548 for (rsv = 0; rsv <= doh; rsv++) {
549 double z[2][2] = { { 0.1, 0.5 } , { 0.5, 0.9 } };
550 double x1 = -0.2;
551 double x2 = 1.2;
552 double y1 = -0.2;
553 double y2 = 1.2;
554 double min = -0.0;
555 double max = 1.0;
556 rspl *rs;
557 unsigned char pa[HEIGHT][WIDTH][3];
558 co tco; /* Test point */
559 double sx,sy;
560 int i,j,k;
561
562 if (rsv == 0)
563 rs = rss;
564 else
565 rs = rss2;
566
567 sx = (x2 - x1)/(double)WIDTH;
568 sy = (y2 - y1)/(double)HEIGHT;
569
570 for (j=0; j < HEIGHT; j++) {
571 double jj = j/(HEIGHT-1.0);
572 tco.p[1] = (double)((HEIGHT-1) - j) * sy + y1;
573 for (i = 0; i < WIDTH; i++) {
574 double ii = j/(HEIGHT-1.0);
575 tco.p[0] = (double)i * sx + x1;
576
577 tco.p[2] = (1.0-ii) * (1.0-jj) * z[0][0]
578 + (1.0-ii) * jj * z[0][1]
579 + ii * (1.0-jj) * z[1][0]
580 + ii * jj * z[1][1];
581
582 if (rs->INTERP(rs, &tco)) {
583 pa[j][i][0] = 0; /* Out of bounds in green */
584 pa[j][i][1] = 100;
585 pa[j][i][2] = 0;
586 } else {
587 int m;
588 /* printf("%d %d, %f %f returned %f\n",i,j,tco.p[0],tco.p[1],tco.v[0]); */
589 m = (int)((255.0 * (tco.v[0] - min)/(max - min)) + 0.5);
590 if (m < 0) {
591 pa[j][i][0] = 20; /* Dark blue */
592 pa[j][i][1] = 20;
593 pa[j][i][2] = 50;
594 } else if (m > 255) {
595 pa[j][i][0] = 230; /* Light blue */
596 pa[j][i][1] = 230;
597 pa[j][i][2] = 255;
598 } else {
599 pa[j][i][0] = m; /* Level in grey */
600 pa[j][i][1] = m;
601 pa[j][i][2] = m;
602 }
603 }
604 }
605 }
606
607 /* Mark verticies in red */
608 for(k = 0; k < npoints; k++) {
609 j = (int)((HEIGHT * (y2 - test_points[k].p[1])/(y2 - y1)) + 0.5);
610 i = (int)((WIDTH * (test_points[k].p[0] - x1)/(x2 - x1)) + 0.5);
611 pa[j][i][0] = 255;
612 pa[j][i][1] = 0;
613 pa[j][i][2] = 0;
614 }
615 write_rgb_tiff(rsv == 0 ? "t3d.tif" : "t3dh.tif" ,WIDTH,HEIGHT,(unsigned char *)pa);
616 }
617
618 /* Plot out 4 slices */
619 if (doplot == 1) {
620 int slice;
621
622 for (slice = 0; slice < 4; slice++) {
623 co tp; /* Test point */
624 double x[PLOTRES];
625 double ya[PLOTRES];
626 double yb[PLOTRES];
627 double xx,yy,zz;
628 double x1,x2,y1,y2,z1,z2;
629 double sx,sy,sz;
630 int i,n;
631
632 /* Set up slice to plot */
633 if (slice == 0) {
634 x1 = 0.5; y1 = 0.5, z1 = 0.0;
635 x2 = 0.5; y2 = 0.5, z2 = 1.0;
636 printf("Plot along z at x = y = 0.5\n");
637 n = PLOTRES;
638 } else if (slice == 1) {
639 x1 = 0.0; y1 = 0.5, z1 = 0.5;
640 x2 = 1.0; y2 = 0.5, z2 = 0.5;
641 printf("Plot along x at y = z = 0.5\n");
642 n = PLOTRES;
643 } else if (slice == 2) {
644 x1 = 0.5; y1 = 0.0, z1 = 0.5;
645 x2 = 0.5; y2 = 1.0, z2 = 0.5;
646 printf("Plot along y at x = z = 0.5\n");
647 n = PLOTRES;
648 } else {
649 x1 = 0.0; y1 = 0.0, z1 = 0.0;
650 x2 = 1.0; y2 = 1.0, z2 = 1.0;
651 printf("Plot along x = y = z\n");
652 n = PLOTRES;
653 }
654
655 sx = (x2 - x1)/n;
656 sy = (y2 - y1)/n;
657 sz = (z2 - z1)/n;
658
659 xx = x1;
660 yy = y1;
661 zz = z1;
662 for (i = 0; i < n; i++) {
663 double vv = i/(n-1.0);
664 x[i] = vv;
665 tp.p[0] = xx;
666 tp.p[1] = yy;
667 tp.p[2] = zz;
668
669 if (rss->INTERP(rss, &tp))
670 tp.v[0] = -0.1;
671 ya[i] = tp.v[0];
672
673 if (doh) {
674 if (rss2->INTERP(rss2, &tp))
675 tp.v[0] = -0.1;
676 yb[i] = tp.v[0];
677 }
678
679 xx += sx;
680 yy += sy;
681 zz += sz;
682 }
683
684 /* Plot the result */
685 if (doh)
686 do_plot(x,ya,yb,NULL,n);
687 else
688 do_plot(x,ya,NULL,NULL,n);
689 }
690 } else if (doplot == 2) {
691 co tp; /* Test point */
692 double x[PLOTRES];
693 double ya[PLOTRES];
694 double yb[PLOTRES];
695 double xx,yy,zz;
696 double x1,x2,y1,y2,z1,z2;
697 double sx,sy,sz;
698 int i,n;
699
700
701 x1 = plotpts[0][0];
702 y1 = plotpts[0][1];
703 z1 = plotpts[0][2];
704 x2 = plotpts[1][0];
705 y2 = plotpts[1][1];
706 z2 = plotpts[1][2];
707
708 printf("Plot along z at x = y = 0.5\n");
709 n = PLOTRES;
710
711 sx = (x2 - x1)/n;
712 sy = (y2 - y1)/n;
713 sz = (z2 - z1)/n;
714
715 xx = x1;
716 yy = y1;
717 zz = z1;
718 for (i = 0; i < n; i++) {
719 double vv = i/(n-1.0);
720 x[i] = vv;
721 tp.p[0] = xx;
722 tp.p[1] = yy;
723 tp.p[2] = zz;
724
725 if (rss->INTERP(rss, &tp))
726 tp.v[0] = -0.1;
727 ya[i] = tp.v[0];
728
729 if (doh) {
730 if (rss2->INTERP(rss2, &tp))
731 tp.v[0] = -0.1;
732 yb[i] = tp.v[0];
733 }
734
735 xx += sx;
736 yy += sy;
737 zz += sz;
738 }
739
740 /* Plot the result */
741 if (doh)
742 do_plot(x,ya,yb,NULL,n);
743 else
744 do_plot(x,ya,NULL,NULL,n);
745 }
746
747 /* Report the fit */
748 {
749 co tco; /* Test point */
750 int k;
751 double avg = 0;
752 double max = 0.0;
753
754 for(k = 0; k < npoints; k++) {
755 double err;
756 tco.p[0] = test_points[k].p[0];
757 tco.p[1] = test_points[k].p[1];
758 tco.p[2] = test_points[k].p[2];
759 rss->INTERP(rss, &tco);
760
761 err = tco.v[0] - test_points[k].v[0];
762 err = fabs(err);
763
764 avg += err;
765 if (err > max)
766 max = err;
767 }
768 avg /= (double)npoints;
769 printf("Max error %f%%, average %f%%\n",100.0 * max, 100.0 * avg);
770 }
771 return 0;
772 }
773
774 /* ---------------------- */
775 /* Tiff diagnostic output */
776
777 void
write_rgb_tiff(char * name,int width,int height,unsigned char * data)778 write_rgb_tiff(
779 char *name,
780 int width,
781 int height,
782 unsigned char *data
783 ) {
784 int y;
785 unsigned char *dp;
786 TIFF *tif;
787
788 if ((tif = TIFFOpen(name, "w")) == NULL) {
789 fprintf(stderr,"Failed to open output TIFF file '%s'\n",name);
790 exit (-1);
791 }
792
793 TIFFSetField(tif, TIFFTAG_IMAGEWIDTH, width);
794 TIFFSetField(tif, TIFFTAG_IMAGELENGTH, height);
795 TIFFSetField(tif, TIFFTAG_ORIENTATION, ORIENTATION_TOPLEFT);
796 TIFFSetField(tif, TIFFTAG_SAMPLESPERPIXEL, 3);
797 TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE, 8);
798 TIFFSetField(tif, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
799 TIFFSetField(tif, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
800 TIFFSetField(tif, TIFFTAG_COMPRESSION, COMPRESSION_NONE);
801
802 for (dp = data, y = 0; y < height; y++, dp += 3 * width) {
803 if (TIFFWriteScanline(tif, (tdata_t)dp, y, 0) < 0) {
804 fprintf(stderr,"WriteScanline Failed at line %d\n",y);
805 exit (-1);
806 }
807 }
808 (void) TIFFClose(tif);
809 }
810
811 #ifdef NEVER
812 /******************************************************************/
813 /* Error/debug output routines */
814 /******************************************************************/
815
816 /* Basic printf type error() and warning() routines */
817
818 #ifdef __STDC__
819 void
error(char * fmt,...)820 error(char *fmt, ...)
821 #else
822 void
823 error(va_alist)
824 va_dcl
825 #endif
826 {
827 va_list args;
828 #ifndef __STDC__
829 char *fmt;
830 #endif
831
832 fprintf(stderr,"cmatch: Error - ");
833 #ifdef __STDC__
834 va_start(args, fmt);
835 #else
836 va_start(args);
837 fmt = va_arg(args, char *);
838 #endif
839 vfprintf(stderr, fmt, args);
840 va_end(args);
841 fprintf(stderr, "\n");
842 fflush(stdout);
843 exit (-1);
844 }
845
846 #ifdef __STDC__
847 void
warning(char * fmt,...)848 warning(char *fmt, ...)
849 #else
850 void
851 warning(va_alist)
852 va_dcl
853 #endif
854 {
855 va_list args;
856 #ifndef __STDC__
857 char *fmt;
858 #endif
859
860 fprintf(stderr,"cmatch: Warning - ");
861 #ifdef __STDC__
862 va_start(args, fmt);
863 #else
864 va_start(args);
865 fmt = va_arg(args, char *);
866 #endif
867 vfprintf(stderr, fmt, args);
868 va_end(args);
869 fprintf(stderr, "\n");
870 }
871
872 #ifdef __STDC__
873 void
verbose(int level,char * fmt,...)874 verbose(int level, char *fmt, ...)
875 {
876 va_list args;
877 va_start(args, fmt);
878 #else
879 verbose(va_alist)
880 va_dcl
881 {
882 va_list args;
883 int level;
884 char *fmt;
885 va_start(args);
886 level = va_arg(args, int);
887 fmt = va_arg(args, char *);
888 #endif
889 if (verbose_level >= level)
890 {
891 fprintf(verbose_out,"cmatch: ");
892 vfprintf(verbose_out, fmt, args);
893 fprintf(verbose_out, "\n");
894 fflush(verbose_out);
895 }
896 va_end(args);
897 }
898 #endif /* NEVER */
899