/src/PROJ/src/pipeline.cpp
Line | Count | Source (jump to first uncovered line) |
1 | | /******************************************************************************* |
2 | | |
3 | | Transformation pipeline manager |
4 | | |
5 | | Thomas Knudsen, 2016-05-20/2016-11-20 |
6 | | |
7 | | ******************************************************************************** |
8 | | |
9 | | Geodetic transformations are typically organized in a number of |
10 | | steps. For example, a datum shift could be carried out through |
11 | | these steps: |
12 | | |
13 | | 1. Convert (latitude, longitude, ellipsoidal height) to |
14 | | 3D geocentric cartesian coordinates (X, Y, Z) |
15 | | 2. Transform the (X, Y, Z) coordinates to the new datum, using a |
16 | | 7 parameter Helmert transformation. |
17 | | 3. Convert (X, Y, Z) back to (latitude, longitude, ellipsoidal height) |
18 | | |
19 | | If the height system used is orthometric, rather than ellipsoidal, |
20 | | another step is needed at each end of the process: |
21 | | |
22 | | 1. Add the local geoid undulation (N) to the orthometric height |
23 | | to obtain the ellipsoidal (i.e. geometric) height. |
24 | | 2. Convert (latitude, longitude, ellipsoidal height) to |
25 | | 3D geocentric cartesian coordinates (X, Y, Z) |
26 | | 3. Transform the (X, Y, Z) coordinates to the new datum, using a |
27 | | 7 parameter Helmert transformation. |
28 | | 4. Convert (X, Y, Z) back to (latitude, longitude, ellipsoidal height) |
29 | | 5. Subtract the local geoid undulation (N) from the ellipsoidal height |
30 | | to obtain the orthometric height. |
31 | | |
32 | | Additional steps can be added for e.g. change of vertical datum, so the |
33 | | list can grow fairly long. None of the steps are, however, particularly |
34 | | complex, and data flow is strictly from top to bottom. |
35 | | |
36 | | Hence, in principle, the first example above could be implemented using |
37 | | Unix pipelines: |
38 | | |
39 | | cat my_coordinates | geographic_to_xyz | helmert | xyz_to_geographic > |
40 | | my_transformed_coordinates |
41 | | |
42 | | in the grand tradition of Software Tools [1]. |
43 | | |
44 | | The proj pipeline driver implements a similar concept: Stringing together |
45 | | a number of steps, feeding the output of one step to the input of the next. |
46 | | |
47 | | It is a very powerful concept, that increases the range of relevance of the |
48 | | proj.4 system substantially. It is, however, not a particularly intrusive |
49 | | addition to the PROJ.4 code base: The implementation is by and large |
50 | | completed by adding an extra projection called "pipeline" (i.e. this file), |
51 | | which handles all business, and a small amount of added functionality in the |
52 | | pj_init code, implementing support for multilevel, embedded pipelines. |
53 | | |
54 | | Syntactically, the pipeline system introduces the "+step" keyword (which |
55 | | indicates the start of each transformation step), and reintroduces the +inv |
56 | | keyword (indicating that a given transformation step should run in reverse, |
57 | | i.e. forward, when the pipeline is executed in inverse direction, and vice |
58 | | versa). |
59 | | |
60 | | Hence, the first transformation example above, can be implemented as: |
61 | | |
62 | | +proj=pipeline +step proj=cart +step proj=helmert <ARGS> +step proj=cart |
63 | | +inv |
64 | | |
65 | | Where <ARGS> indicate the Helmert arguments: 3 translations (+x=..., +y=..., |
66 | | +z=...), 3 rotations (+rx=..., +ry=..., +rz=...) and a scale factor |
67 | | (+s=...). Following geodetic conventions, the rotations are given in arcseconds, |
68 | | and the scale factor is given as parts-per-million. |
69 | | |
70 | | [1] B. W. Kernighan & P. J. Plauger: Software tools. |
71 | | Reading, Massachusetts, Addison-Wesley, 1976, 338 pp. |
72 | | |
73 | | ******************************************************************************** |
74 | | |
75 | | Thomas Knudsen, thokn@sdfe.dk, 2016-05-20 |
76 | | |
77 | | ******************************************************************************** |
78 | | * Copyright (c) 2016, 2017, 2018 Thomas Knudsen / SDFE |
79 | | * |
80 | | * Permission is hereby granted, free of charge, to any person obtaining a |
81 | | * copy of this software and associated documentation files (the "Software"), |
82 | | * to deal in the Software without restriction, including without limitation |
83 | | * the rights to use, copy, modify, merge, publish, distribute, sublicense, |
84 | | * and/or sell copies of the Software, and to permit persons to whom the |
85 | | * Software is furnished to do so, subject to the following conditions: |
86 | | * |
87 | | * The above copyright notice and this permission notice shall be included |
88 | | * in all copies or substantial portions of the Software. |
89 | | * |
90 | | * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS |
91 | | * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
92 | | * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL |
93 | | * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
94 | | * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING |
95 | | * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER |
96 | | * DEALINGS IN THE SOFTWARE. |
97 | | * |
98 | | ********************************************************************************/ |
99 | | |
100 | | #include <math.h> |
101 | | #include <stack> |
102 | | #include <stddef.h> |
103 | | #include <string.h> |
104 | | #include <vector> |
105 | | |
106 | | #include "geodesic.h" |
107 | | #include "proj.h" |
108 | | #include "proj_internal.h" |
109 | | |
110 | | PROJ_HEAD(pipeline, "Transformation pipeline manager"); |
111 | | PROJ_HEAD(pop, "Retrieve coordinate value from pipeline stack"); |
112 | | PROJ_HEAD(push, "Save coordinate value on pipeline stack"); |
113 | | |
114 | | /* Projection specific elements for the PJ object */ |
115 | | namespace { // anonymous namespace |
116 | | |
117 | | struct Step { |
118 | | PJ *pj = nullptr; |
119 | | bool omit_fwd = false; |
120 | | bool omit_inv = false; |
121 | | |
122 | | Step(PJ *pjIn, bool omitFwdIn, bool omitInvIn) |
123 | 103k | : pj(pjIn), omit_fwd(omitFwdIn), omit_inv(omitInvIn) {} |
124 | | Step(Step &&other) |
125 | 105k | : pj(std::move(other.pj)), omit_fwd(other.omit_fwd), |
126 | 105k | omit_inv(other.omit_inv) { |
127 | 105k | other.pj = nullptr; |
128 | 105k | } |
129 | | Step(const Step &) = delete; |
130 | | Step &operator=(const Step &) = delete; |
131 | | |
132 | 208k | ~Step() { proj_destroy(pj); } |
133 | | }; |
134 | | |
135 | | struct Pipeline { |
136 | | char **argv = nullptr; |
137 | | char **current_argv = nullptr; |
138 | | std::vector<Step> steps{}; |
139 | | std::stack<double> stack[4]; |
140 | | }; |
141 | | |
142 | | struct PushPop { |
143 | | bool v1; |
144 | | bool v2; |
145 | | bool v3; |
146 | | bool v4; |
147 | | }; |
148 | | } // anonymous namespace |
149 | | |
150 | | static void pipeline_forward_4d(PJ_COORD &point, PJ *P); |
151 | | static void pipeline_reverse_4d(PJ_COORD &point, PJ *P); |
152 | | static PJ_XYZ pipeline_forward_3d(PJ_LPZ lpz, PJ *P); |
153 | | static PJ_LPZ pipeline_reverse_3d(PJ_XYZ xyz, PJ *P); |
154 | | static PJ_XY pipeline_forward(PJ_LP lp, PJ *P); |
155 | | static PJ_LP pipeline_reverse(PJ_XY xy, PJ *P); |
156 | | |
157 | 0 | static void pipeline_reassign_context(PJ *P, PJ_CONTEXT *ctx) { |
158 | 0 | auto pipeline = static_cast<struct Pipeline *>(P->opaque); |
159 | 0 | for (auto &step : pipeline->steps) |
160 | 0 | proj_assign_context(step.pj, ctx); |
161 | 0 | } |
162 | | |
163 | 1.39M | static void pipeline_forward_4d(PJ_COORD &point, PJ *P) { |
164 | 1.39M | auto pipeline = static_cast<struct Pipeline *>(P->opaque); |
165 | 3.36M | for (auto &step : pipeline->steps) { |
166 | 3.36M | if (!step.omit_fwd) { |
167 | 3.36M | if (!step.pj->inverted) |
168 | 3.15M | pj_fwd4d(point, step.pj); |
169 | 208k | else |
170 | 208k | pj_inv4d(point, step.pj); |
171 | 3.36M | if (point.xyzt.x == HUGE_VAL) { |
172 | 14.5k | break; |
173 | 14.5k | } |
174 | 3.36M | } |
175 | 3.36M | } |
176 | 1.39M | } |
177 | | |
178 | 0 | static void pipeline_reverse_4d(PJ_COORD &point, PJ *P) { |
179 | 0 | auto pipeline = static_cast<struct Pipeline *>(P->opaque); |
180 | 0 | for (auto iterStep = pipeline->steps.rbegin(); |
181 | 0 | iterStep != pipeline->steps.rend(); ++iterStep) { |
182 | 0 | const auto &step = *iterStep; |
183 | 0 | if (!step.omit_inv) { |
184 | 0 | if (step.pj->inverted) |
185 | 0 | pj_fwd4d(point, step.pj); |
186 | 0 | else |
187 | 0 | pj_inv4d(point, step.pj); |
188 | 0 | if (point.xyzt.x == HUGE_VAL) { |
189 | 0 | break; |
190 | 0 | } |
191 | 0 | } |
192 | 0 | } |
193 | 0 | } |
194 | | |
195 | 0 | static PJ_XYZ pipeline_forward_3d(PJ_LPZ lpz, PJ *P) { |
196 | 0 | PJ_COORD point = {{0, 0, 0, 0}}; |
197 | 0 | point.lpz = lpz; |
198 | 0 | auto pipeline = static_cast<struct Pipeline *>(P->opaque); |
199 | 0 | for (auto &step : pipeline->steps) { |
200 | 0 | if (!step.omit_fwd) { |
201 | 0 | point = pj_approx_3D_trans(step.pj, PJ_FWD, point); |
202 | 0 | if (point.xyzt.x == HUGE_VAL) { |
203 | 0 | break; |
204 | 0 | } |
205 | 0 | } |
206 | 0 | } |
207 | |
|
208 | 0 | return point.xyz; |
209 | 0 | } |
210 | | |
211 | 0 | static PJ_LPZ pipeline_reverse_3d(PJ_XYZ xyz, PJ *P) { |
212 | 0 | PJ_COORD point = {{0, 0, 0, 0}}; |
213 | 0 | point.xyz = xyz; |
214 | 0 | auto pipeline = static_cast<struct Pipeline *>(P->opaque); |
215 | 0 | for (auto iterStep = pipeline->steps.rbegin(); |
216 | 0 | iterStep != pipeline->steps.rend(); ++iterStep) { |
217 | 0 | const auto &step = *iterStep; |
218 | 0 | if (!step.omit_inv) { |
219 | 0 | point = proj_trans(step.pj, PJ_INV, point); |
220 | 0 | if (point.xyzt.x == HUGE_VAL) { |
221 | 0 | break; |
222 | 0 | } |
223 | 0 | } |
224 | 0 | } |
225 | |
|
226 | 0 | return point.lpz; |
227 | 0 | } |
228 | | |
229 | 0 | static PJ_XY pipeline_forward(PJ_LP lp, PJ *P) { |
230 | 0 | PJ_COORD point = {{0, 0, 0, 0}}; |
231 | 0 | point.lp = lp; |
232 | 0 | auto pipeline = static_cast<struct Pipeline *>(P->opaque); |
233 | 0 | for (auto &step : pipeline->steps) { |
234 | 0 | if (!step.omit_fwd) { |
235 | 0 | point = pj_approx_2D_trans(step.pj, PJ_FWD, point); |
236 | 0 | if (point.xyzt.x == HUGE_VAL) { |
237 | 0 | break; |
238 | 0 | } |
239 | 0 | } |
240 | 0 | } |
241 | |
|
242 | 0 | return point.xy; |
243 | 0 | } |
244 | | |
245 | 0 | static PJ_LP pipeline_reverse(PJ_XY xy, PJ *P) { |
246 | 0 | PJ_COORD point = {{0, 0, 0, 0}}; |
247 | 0 | point.xy = xy; |
248 | 0 | auto pipeline = static_cast<struct Pipeline *>(P->opaque); |
249 | 0 | for (auto iterStep = pipeline->steps.rbegin(); |
250 | 0 | iterStep != pipeline->steps.rend(); ++iterStep) { |
251 | 0 | const auto &step = *iterStep; |
252 | 0 | if (!step.omit_inv) { |
253 | 0 | point = pj_approx_2D_trans(step.pj, PJ_INV, point); |
254 | 0 | if (point.xyzt.x == HUGE_VAL) { |
255 | 0 | break; |
256 | 0 | } |
257 | 0 | } |
258 | 0 | } |
259 | |
|
260 | 0 | return point.lp; |
261 | 0 | } |
262 | | |
263 | 25.6k | static PJ *destructor(PJ *P, int errlev) { |
264 | 25.6k | if (nullptr == P) |
265 | 0 | return nullptr; |
266 | | |
267 | 25.6k | if (nullptr == P->opaque) |
268 | 0 | return pj_default_destructor(P, errlev); |
269 | | |
270 | 25.6k | auto pipeline = static_cast<struct Pipeline *>(P->opaque); |
271 | | |
272 | 25.6k | free(pipeline->argv); |
273 | 25.6k | free(pipeline->current_argv); |
274 | | |
275 | 25.6k | delete pipeline; |
276 | 25.6k | P->opaque = nullptr; |
277 | | |
278 | 25.6k | return pj_default_destructor(P, errlev); |
279 | 25.6k | } |
280 | | |
281 | | /* count the number of args in pipeline definition, and mark all args as used */ |
282 | 25.6k | static size_t argc_params(paralist *params) { |
283 | 25.6k | size_t argc = 0; |
284 | 1.25M | for (; params != nullptr; params = params->next) { |
285 | 1.23M | argc++; |
286 | 1.23M | params->used = 1; |
287 | 1.23M | } |
288 | 25.6k | return ++argc; /* one extra for the sentinel */ |
289 | 25.6k | } |
290 | | |
291 | | /* Sentinel for argument list */ |
292 | | static const char *argv_sentinel = "step"; |
293 | | |
294 | | /* turn paralist into argc/argv style argument list */ |
295 | 25.6k | static char **argv_params(paralist *params, size_t argc) { |
296 | 25.6k | char **argv; |
297 | 25.6k | size_t i = 0; |
298 | 25.6k | argv = static_cast<char **>(calloc(argc, sizeof(char *))); |
299 | 25.6k | if (nullptr == argv) |
300 | 0 | return nullptr; |
301 | 1.25M | for (; params != nullptr; params = params->next) |
302 | 1.23M | argv[i++] = params->param; |
303 | 25.6k | argv[i++] = const_cast<char *>(argv_sentinel); |
304 | 25.6k | return argv; |
305 | 25.6k | } |
306 | | |
307 | | /* Being the special operator that the pipeline is, we have to handle the */ |
308 | | /* ellipsoid differently than usual. In general, the pipeline operation does */ |
309 | | /* not need an ellipsoid, but in some cases it is beneficial nonetheless. */ |
310 | | /* Unfortunately we can't use the normal ellipsoid setter in pj_init, since */ |
311 | | /* it adds a +ellps parameter to the global args if nothing else is specified*/ |
312 | | /* This is problematic since that ellipsoid spec is then passed on to the */ |
313 | | /* pipeline children. This is rarely what we want, so here we implement our */ |
314 | | /* own logic instead. If an ellipsoid is set in the global args, it is used */ |
315 | | /* as the pipeline ellipsoid. Otherwise we use GRS80 parameters as default. */ |
316 | | /* At last we calculate the rest of the ellipsoid parameters and */ |
317 | | /* re-initialize P->geod. */ |
318 | 23.9k | static void set_ellipsoid(PJ *P) { |
319 | 23.9k | paralist *cur, *attachment; |
320 | 23.9k | int err = proj_errno_reset(P); |
321 | | |
322 | | /* Break the linked list after the global args */ |
323 | 23.9k | attachment = nullptr; |
324 | 87.5k | for (cur = P->params; cur != nullptr; cur = cur->next) |
325 | | /* cur->next will always be non 0 given argv_sentinel presence, */ |
326 | | /* but this is far from being obvious for a static analyzer */ |
327 | 87.5k | if (cur->next != nullptr && |
328 | 87.5k | strcmp(argv_sentinel, cur->next->param) == 0) { |
329 | 23.9k | attachment = cur->next; |
330 | 23.9k | cur->next = nullptr; |
331 | 23.9k | break; |
332 | 23.9k | } |
333 | | |
334 | | /* Check if there's any ellipsoid specification in the global params. */ |
335 | | /* If not, use GRS80 as default */ |
336 | 23.9k | if (0 != pj_ellipsoid(P)) { |
337 | 116 | P->a = 6378137.0; |
338 | 116 | P->f = 1.0 / 298.257222101; |
339 | 116 | P->es = 2 * P->f - P->f * P->f; |
340 | | |
341 | | /* reset an "unerror": In this special use case, the errno is */ |
342 | | /* not an error signal, but just a reply from pj_ellipsoid, */ |
343 | | /* telling us that "No - there was no ellipsoid definition in */ |
344 | | /* the PJ you provided". */ |
345 | 116 | proj_errno_reset(P); |
346 | 116 | } |
347 | 23.9k | P->a_orig = P->a; |
348 | 23.9k | P->es_orig = P->es; |
349 | | |
350 | 23.9k | if (pj_calc_ellipsoid_params(P, P->a, P->es) == 0) |
351 | 23.9k | geod_init(P->geod, P->a, P->f); |
352 | | |
353 | | /* Re-attach the dangling list */ |
354 | | /* Note: cur will always be non 0 given argv_sentinel presence, */ |
355 | | /* but this is far from being obvious for a static analyzer */ |
356 | 23.9k | if (cur != nullptr) |
357 | 23.9k | cur->next = attachment; |
358 | 23.9k | proj_errno_restore(P, err); |
359 | 23.9k | } |
360 | | |
361 | 25.6k | PJ *OPERATION(pipeline, 0) { |
362 | 25.6k | int i, nsteps = 0, argc; |
363 | 25.6k | int i_pipeline = -1, i_first_step = -1, i_current_step; |
364 | 25.6k | char **argv, **current_argv; |
365 | | |
366 | 25.6k | if (P->ctx->pipelineInitRecursiongCounter == 5) { |
367 | | // Can happen for a string like: |
368 | | // proj=pipeline step "x="""," u=" proj=pipeline step ste=""[" u=" |
369 | | // proj=pipeline step ste="[" u=" proj=pipeline step ste="[" u=" |
370 | | // proj=pipeline step ste="[" u=" proj=pipeline step ste="[" u=" |
371 | | // proj=pipeline step ste="[" u=" proj=pipeline step ste="[" u=" |
372 | | // proj=pipeline step ste="[" u=" proj=pipeline p step ste="[" u=" |
373 | | // proj=pipeline step ste="[" u=" proj=pipeline step ste="[" u=" |
374 | | // proj=pipeline step ste="[" u=" proj=pipeline step ""x=""""""""""" |
375 | | // Probably an issue with the quoting handling code |
376 | | // But doesn't hurt to add an extra safety check |
377 | 0 | proj_log_error(P, _("Pipeline: too deep recursion")); |
378 | 0 | return destructor( |
379 | 0 | P, PROJ_ERR_INVALID_OP_WRONG_SYNTAX); /* ERROR: nested pipelines */ |
380 | 0 | } |
381 | | |
382 | 25.6k | P->fwd4d = pipeline_forward_4d; |
383 | 25.6k | P->inv4d = pipeline_reverse_4d; |
384 | 25.6k | P->fwd3d = pipeline_forward_3d; |
385 | 25.6k | P->inv3d = pipeline_reverse_3d; |
386 | 25.6k | P->fwd = pipeline_forward; |
387 | 25.6k | P->inv = pipeline_reverse; |
388 | 25.6k | P->destructor = destructor; |
389 | 25.6k | P->reassign_context = pipeline_reassign_context; |
390 | | |
391 | | /* Currently, the pipeline driver is a raw bit mover, enabling other |
392 | | * operations */ |
393 | | /* to collaborate efficiently. All prep/fin stuff is done at the step |
394 | | * levels. |
395 | | */ |
396 | 25.6k | P->skip_fwd_prepare = 1; |
397 | 25.6k | P->skip_fwd_finalize = 1; |
398 | 25.6k | P->skip_inv_prepare = 1; |
399 | 25.6k | P->skip_inv_finalize = 1; |
400 | | |
401 | 25.6k | P->opaque = new (std::nothrow) Pipeline(); |
402 | 25.6k | if (nullptr == P->opaque) |
403 | 0 | return destructor(P, PROJ_ERR_INVALID_OP /* ENOMEM */); |
404 | | |
405 | 25.6k | argc = (int)argc_params(P->params); |
406 | 25.6k | auto pipeline = static_cast<struct Pipeline *>(P->opaque); |
407 | 25.6k | pipeline->argv = argv = argv_params(P->params, argc); |
408 | 25.6k | if (nullptr == argv) |
409 | 0 | return destructor(P, PROJ_ERR_INVALID_OP /* ENOMEM */); |
410 | | |
411 | 25.6k | pipeline->current_argv = current_argv = |
412 | 25.6k | static_cast<char **>(calloc(argc, sizeof(char *))); |
413 | 25.6k | if (nullptr == current_argv) |
414 | 0 | return destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); |
415 | | |
416 | | /* Do some syntactical sanity checking */ |
417 | 1.16M | for (i = 0; i < argc && argv[i] != nullptr; i++) { |
418 | 1.13M | if (0 == strcmp(argv_sentinel, argv[i])) { |
419 | 229k | if (-1 == i_pipeline) { |
420 | 0 | proj_log_error(P, _("Pipeline: +step before +proj=pipeline")); |
421 | 0 | return destructor(P, PROJ_ERR_INVALID_OP_WRONG_SYNTAX); |
422 | 0 | } |
423 | 229k | if (0 == nsteps) |
424 | 24.0k | i_first_step = i; |
425 | 229k | nsteps++; |
426 | 229k | continue; |
427 | 229k | } |
428 | | |
429 | 909k | if (0 == strcmp("proj=pipeline", argv[i])) { |
430 | 25.6k | if (-1 != i_pipeline) { |
431 | 0 | proj_log_error(P, _("Pipeline: Nesting only allowed when child " |
432 | 0 | "pipelines are wrapped in '+init's")); |
433 | 0 | return destructor( |
434 | 0 | P, PROJ_ERR_INVALID_OP_WRONG_SYNTAX); /* ERROR: nested |
435 | | pipelines */ |
436 | 0 | } |
437 | 25.6k | i_pipeline = i; |
438 | 884k | } else if (0 == nsteps && 0 == strncmp(argv[i], "proj=", 5)) { |
439 | | // Non-sensical to have proj= in the general pipeline parameters. |
440 | | // Would not be a big issue in itself, but this makes bad |
441 | | // performance in parsing hostile pipelines more likely, such as the |
442 | | // one of |
443 | | // https://bugs.chromium.org/p/oss-fuzz/issues/detail?id=41290 |
444 | 1.61k | proj_log_error( |
445 | 1.61k | P, _("Pipeline: proj= operator before first step not allowed")); |
446 | 1.61k | return destructor(P, PROJ_ERR_INVALID_OP_WRONG_SYNTAX); |
447 | 882k | } else if (0 == nsteps && 0 == strncmp(argv[i], "o_proj=", 7)) { |
448 | | // Same as above. |
449 | 14 | proj_log_error( |
450 | 14 | P, |
451 | 14 | _("Pipeline: o_proj= operator before first step not allowed")); |
452 | 14 | return destructor(P, PROJ_ERR_INVALID_OP_WRONG_SYNTAX); |
453 | 14 | } |
454 | 909k | } |
455 | 24.0k | nsteps--; /* Last instance of +step is just a sentinel */ |
456 | | |
457 | 24.0k | if (-1 == i_pipeline) |
458 | 0 | return destructor( |
459 | 0 | P, PROJ_ERR_INVALID_OP_WRONG_SYNTAX); /* ERROR: no pipeline def */ |
460 | | |
461 | 24.0k | if (0 == nsteps) |
462 | 60 | return destructor( |
463 | 60 | P, PROJ_ERR_INVALID_OP_WRONG_SYNTAX); /* ERROR: no pipeline def */ |
464 | | |
465 | 23.9k | set_ellipsoid(P); |
466 | | |
467 | | /* Now loop over all steps, building a new set of arguments for each init */ |
468 | 23.9k | i_current_step = i_first_step; |
469 | 127k | for (i = 0; i < nsteps; i++) { |
470 | 117k | int j; |
471 | 117k | int current_argc = 0; |
472 | 117k | int err; |
473 | 117k | PJ *next_step = nullptr; |
474 | | |
475 | | /* Build a set of setup args for the current step */ |
476 | 117k | proj_log_trace(P, "Pipeline: Building arg list for step no. %d", i); |
477 | | |
478 | | /* First add the step specific args */ |
479 | 575k | for (j = i_current_step + 1; 0 != strcmp("step", argv[j]); j++) |
480 | 458k | current_argv[current_argc++] = argv[j]; |
481 | | |
482 | 117k | i_current_step = j; |
483 | | |
484 | | /* Then add the global args */ |
485 | 713k | for (j = i_pipeline + 1; 0 != strcmp("step", argv[j]); j++) |
486 | 596k | current_argv[current_argc++] = argv[j]; |
487 | | |
488 | 117k | proj_log_trace(P, "Pipeline: init - %s, %d", current_argv[0], |
489 | 117k | current_argc); |
490 | 1.05M | for (j = 1; j < current_argc; j++) |
491 | 937k | proj_log_trace(P, " %s", current_argv[j]); |
492 | | |
493 | 117k | err = proj_errno_reset(P); |
494 | | |
495 | 117k | P->ctx->pipelineInitRecursiongCounter++; |
496 | 117k | next_step = pj_create_argv_internal(P->ctx, current_argc, current_argv); |
497 | 117k | P->ctx->pipelineInitRecursiongCounter--; |
498 | 117k | proj_log_trace(P, "Pipeline: Step %d (%s) at %p", i, current_argv[0], |
499 | 117k | next_step); |
500 | | |
501 | 117k | if (nullptr == next_step) { |
502 | | /* The step init failed, but possibly without setting errno. If so, |
503 | | * we say "malformed" */ |
504 | 14.1k | int err_to_report = proj_errno(P); |
505 | 14.1k | if (0 == err_to_report) |
506 | 0 | err_to_report = PROJ_ERR_INVALID_OP_WRONG_SYNTAX; |
507 | 14.1k | proj_log_error(P, _("Pipeline: Bad step definition: %s (%s)"), |
508 | 14.1k | current_argv[0], |
509 | 14.1k | proj_context_errno_string(P->ctx, err_to_report)); |
510 | 14.1k | return destructor(P, err_to_report); /* ERROR: bad pipeline def */ |
511 | 14.1k | } |
512 | 103k | next_step->parent = P; |
513 | | |
514 | 103k | proj_errno_restore(P, err); |
515 | | |
516 | | /* Is this step inverted? */ |
517 | 1.06M | for (j = 0; j < current_argc; j++) { |
518 | 961k | if (0 == strcmp("inv", current_argv[j])) { |
519 | | /* if +inv exists in both global and local args the forward |
520 | | * operation should be used */ |
521 | 33.7k | next_step->inverted = next_step->inverted == 0 ? 1 : 0; |
522 | 33.7k | } |
523 | 961k | } |
524 | | |
525 | 103k | bool omit_fwd = pj_param(P->ctx, next_step->params, "bomit_fwd").i != 0; |
526 | 103k | bool omit_inv = pj_param(P->ctx, next_step->params, "bomit_inv").i != 0; |
527 | 103k | pipeline->steps.emplace_back(next_step, omit_fwd, omit_inv); |
528 | | |
529 | 103k | proj_log_trace(P, "Pipeline at [%p]: step at [%p] (%s) done", P, |
530 | 103k | next_step, current_argv[0]); |
531 | 103k | } |
532 | | |
533 | | /* Require a forward path through the pipeline */ |
534 | 50.4k | for (auto &step : pipeline->steps) { |
535 | 50.4k | PJ *Q = step.pj; |
536 | 50.4k | if (step.omit_fwd) { |
537 | 79 | continue; |
538 | 79 | } |
539 | 50.3k | if (Q->inverted) { |
540 | 15.0k | if (Q->inv || Q->inv3d || Q->inv4d) { |
541 | 14.7k | continue; |
542 | 14.7k | } |
543 | 216 | proj_log_error( |
544 | 216 | P, _("Pipeline: Inverse operation for %s is not available"), |
545 | 216 | Q->short_name); |
546 | 216 | return destructor(P, PROJ_ERR_OTHER_NO_INVERSE_OP); |
547 | 35.3k | } else { |
548 | 35.3k | if (Q->fwd || Q->fwd3d || Q->fwd4d) { |
549 | 35.3k | continue; |
550 | 35.3k | } |
551 | 1 | proj_log_error( |
552 | 1 | P, _("Pipeline: Forward operation for %s is not available"), |
553 | 1 | Q->short_name); |
554 | 1 | return destructor(P, PROJ_ERR_INVALID_OP_WRONG_SYNTAX); |
555 | 35.3k | } |
556 | 50.3k | } |
557 | | |
558 | | /* determine if an inverse operation is possible */ |
559 | 49.7k | for (auto &step : pipeline->steps) { |
560 | 49.7k | PJ *Q = step.pj; |
561 | 49.7k | if (step.omit_inv || pj_has_inverse(Q)) { |
562 | 49.3k | continue; |
563 | 49.3k | } else { |
564 | 407 | P->inv = nullptr; |
565 | 407 | P->inv3d = nullptr; |
566 | 407 | P->inv4d = nullptr; |
567 | 407 | break; |
568 | 407 | } |
569 | 49.7k | } |
570 | | |
571 | | /* Replace PJ_IO_UNITS_WHATEVER with input/output units of neighbouring |
572 | | * steps where */ |
573 | | /* it make sense. It does in most cases but not always, for instance */ |
574 | | /* proj=pipeline step proj=unitconvert xy_in=deg xy_out=rad step ... */ |
575 | | /* where the left-hand side units of the first step shouldn't be changed to |
576 | | * RADIANS */ |
577 | | /* as it will result in deg->rad conversions in cs2cs and other |
578 | | * applications. |
579 | | */ |
580 | | |
581 | 49.9k | for (i = nsteps - 2; i >= 0; --i) { |
582 | 40.3k | auto pj = pipeline->steps[i].pj; |
583 | 40.3k | if (pj_left(pj) == PJ_IO_UNITS_WHATEVER && |
584 | 40.3k | pj_right(pj) == PJ_IO_UNITS_WHATEVER) { |
585 | 9.94k | const auto right_pj = pipeline->steps[i + 1].pj; |
586 | 9.94k | const auto right_pj_left = pj_left(right_pj); |
587 | 9.94k | const auto right_pj_right = pj_right(right_pj); |
588 | 9.94k | if (right_pj_left != right_pj_right || |
589 | 9.94k | right_pj_left != PJ_IO_UNITS_WHATEVER) { |
590 | 8.99k | pj->left = right_pj_left; |
591 | 8.99k | pj->right = right_pj_left; |
592 | 8.99k | } |
593 | 9.94k | } |
594 | 40.3k | } |
595 | | |
596 | 49.9k | for (i = 1; i < nsteps; i++) { |
597 | 40.3k | auto pj = pipeline->steps[i].pj; |
598 | 40.3k | if (pj_left(pj) == PJ_IO_UNITS_WHATEVER && |
599 | 40.3k | pj_right(pj) == PJ_IO_UNITS_WHATEVER) { |
600 | 2.40k | const auto left_pj = pipeline->steps[i - 1].pj; |
601 | 2.40k | const auto left_pj_left = pj_left(left_pj); |
602 | 2.40k | const auto left_pj_right = pj_right(left_pj); |
603 | 2.40k | if (left_pj_left != left_pj_right || |
604 | 2.40k | left_pj_right != PJ_IO_UNITS_WHATEVER) { |
605 | 1.74k | pj->left = left_pj_right; |
606 | 1.74k | pj->right = left_pj_right; |
607 | 1.74k | } |
608 | 2.40k | } |
609 | 40.3k | } |
610 | | |
611 | | /* Check that units between each steps match each other, fail if they don't |
612 | | */ |
613 | 45.9k | for (i = 0; i + 1 < nsteps; i++) { |
614 | 36.8k | enum pj_io_units curr_step_output = pj_right(pipeline->steps[i].pj); |
615 | 36.8k | enum pj_io_units next_step_input = pj_left(pipeline->steps[i + 1].pj); |
616 | | |
617 | 36.8k | if (curr_step_output == PJ_IO_UNITS_WHATEVER || |
618 | 36.8k | next_step_input == PJ_IO_UNITS_WHATEVER) |
619 | 3.62k | continue; |
620 | | |
621 | 33.2k | if (curr_step_output != next_step_input) { |
622 | 493 | proj_log_error( |
623 | 493 | P, _("Pipeline: Mismatched units between step %d and %d"), |
624 | 493 | i + 1, i + 2); |
625 | 493 | return destructor(P, PROJ_ERR_INVALID_OP_WRONG_SYNTAX); |
626 | 493 | } |
627 | 33.2k | } |
628 | | |
629 | 9.08k | proj_log_trace( |
630 | 9.08k | P, "Pipeline: %d steps built. Determining i/o characteristics", nsteps); |
631 | | |
632 | | /* Determine forward input (= reverse output) data type */ |
633 | 9.08k | P->left = pj_left(pipeline->steps.front().pj); |
634 | | |
635 | | /* Now, correspondingly determine forward output (= reverse input) data type |
636 | | */ |
637 | 9.08k | P->right = pj_right(pipeline->steps.back().pj); |
638 | 9.08k | return P; |
639 | 9.58k | } |
640 | | |
641 | 9.74k | static void push(PJ_COORD &point, PJ *P) { |
642 | 9.74k | if (P->parent == nullptr) |
643 | 0 | return; |
644 | | |
645 | 9.74k | struct Pipeline *pipeline = |
646 | 9.74k | static_cast<struct Pipeline *>(P->parent->opaque); |
647 | 9.74k | struct PushPop *pushpop = static_cast<struct PushPop *>(P->opaque); |
648 | | |
649 | 9.74k | if (pushpop->v1) |
650 | 7.05k | pipeline->stack[0].push(point.v[0]); |
651 | 9.74k | if (pushpop->v2) |
652 | 2.01k | pipeline->stack[1].push(point.v[1]); |
653 | 9.74k | if (pushpop->v3) |
654 | 0 | pipeline->stack[2].push(point.v[2]); |
655 | 9.74k | if (pushpop->v4) |
656 | 0 | pipeline->stack[3].push(point.v[3]); |
657 | 9.74k | } |
658 | | |
659 | 119k | static void pop(PJ_COORD &point, PJ *P) { |
660 | 119k | if (P->parent == nullptr) |
661 | 0 | return; |
662 | | |
663 | 119k | struct Pipeline *pipeline = |
664 | 119k | static_cast<struct Pipeline *>(P->parent->opaque); |
665 | 119k | struct PushPop *pushpop = static_cast<struct PushPop *>(P->opaque); |
666 | | |
667 | 119k | if (pushpop->v1 && !pipeline->stack[0].empty()) { |
668 | 5.03k | point.v[0] = pipeline->stack[0].top(); |
669 | 5.03k | pipeline->stack[0].pop(); |
670 | 5.03k | } |
671 | | |
672 | 119k | if (pushpop->v2 && !pipeline->stack[1].empty()) { |
673 | 1.67k | point.v[1] = pipeline->stack[1].top(); |
674 | 1.67k | pipeline->stack[1].pop(); |
675 | 1.67k | } |
676 | | |
677 | 119k | if (pushpop->v3 && !pipeline->stack[2].empty()) { |
678 | 0 | point.v[2] = pipeline->stack[2].top(); |
679 | 0 | pipeline->stack[2].pop(); |
680 | 0 | } |
681 | | |
682 | 119k | if (pushpop->v4 && !pipeline->stack[3].empty()) { |
683 | 0 | point.v[3] = pipeline->stack[3].top(); |
684 | 0 | pipeline->stack[3].pop(); |
685 | 0 | } |
686 | 119k | } |
687 | | |
688 | 23.5k | static PJ *setup_pushpop(PJ *P) { |
689 | 23.5k | auto pushpop = |
690 | 23.5k | static_cast<struct PushPop *>(calloc(1, sizeof(struct PushPop))); |
691 | 23.5k | P->opaque = pushpop; |
692 | 23.5k | if (nullptr == P->opaque) |
693 | 0 | return destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); |
694 | | |
695 | 23.5k | if (pj_param_exists(P->params, "v_1")) |
696 | 859 | pushpop->v1 = true; |
697 | | |
698 | 23.5k | if (pj_param_exists(P->params, "v_2")) |
699 | 490 | pushpop->v2 = true; |
700 | | |
701 | 23.5k | if (pj_param_exists(P->params, "v_3")) |
702 | 19.5k | pushpop->v3 = true; |
703 | | |
704 | 23.5k | if (pj_param_exists(P->params, "v_4")) |
705 | 5 | pushpop->v4 = true; |
706 | | |
707 | 23.5k | P->left = PJ_IO_UNITS_WHATEVER; |
708 | 23.5k | P->right = PJ_IO_UNITS_WHATEVER; |
709 | | |
710 | 23.5k | return P; |
711 | 23.5k | } |
712 | | |
713 | 11.5k | PJ *OPERATION(push, 0) { |
714 | 11.5k | P->fwd4d = push; |
715 | 11.5k | P->inv4d = pop; |
716 | | |
717 | 11.5k | return setup_pushpop(P); |
718 | 11.5k | } |
719 | | |
720 | 11.9k | PJ *OPERATION(pop, 0) { |
721 | 11.9k | P->inv4d = push; |
722 | 11.9k | P->fwd4d = pop; |
723 | | |
724 | 11.9k | return setup_pushpop(P); |
725 | 11.9k | } |