Add +proj=stack +push=idx,... / +pop=idx,... / +swap operation

[rouault]

Implements https://lists.osgeo.org/pipermail/proj/2024-February/011251.html

Push/pop/swap operation operating on a stack of coordinate tuples, attached to a pipeline.

Superseds existing push/pop

https://github.com/rouault/PROJ/blob/stack/docs/source/operations/conversions/stack.rst for more details

The below example demonstrates a pipeline transforming coordinates in ETRS 89 (longitude, latitude, ellipsoidal height) to (longitude, latitude) in the S-JTSK/05 datum and orthometric height in the Baltic 1957 datum.

::

echo 15 50 100 | cct -d 10 +proj=pipeline \
            +step +proj=stack +push=3 +omit_inv \           # (1)
            +step +proj=vgridshift +grids=CR2005.tif \      # (2)
            +step +proj=stack +push=3 \                     # (3)
            +step +proj=stack +swap \                       # (4)
            +step +proj=stack +pop=3 \                      # (5)
            +step +proj=cart +ellps=GRS80 \                 # (6)
            +step +inv +proj=helmert +x=572.213 +y=85.334 +z=461.94 \
                +rx=-4.9732 +ry=-1.529 +rz=-5.2484 +s=3.5378 +convention=coordinate_frame \
            +step +inv +proj=cart +ellps=bessel \
            +step +proj=stack +pop=3                        # (7)

     15.0011680291   50.0007534747  55.0384863419

Let's examine step by step, when executing the pipeline in the forward direction:

1. Save the ETRS89 ellipsoidal height on the stack
2. Apply the geoid model to transform the ETRS89 ellipsoidal height into a Baltic 1957 orthometric height
3. Save the Baltic 1957 on the stack ("above" the ETRS89 ellipsoidal height)
4. Swap the top 2 tuples of the stack, that is now the ETRS89 height will be on top of the Baltic 1957 one.
5. Pop the ETRS89 height from the stack as the active Z value.
6. Apply a 3D Helmert transformation to go from ETRS89 to S-JTSK/05
7. Pop the Baltic 1957 height from the stack as the active Z value.

When run in the inverse direction, the steps are interpreted as:

7. Push the Baltic 1957 height on the stack
8. Apply the inverse 3D Helmert transformation to go from S-JTSK/05 to ETRS89
9. Push the ETRS89 height on the stack
10. Swap the top 2 tuples of the stack, that is now the Baltic 1957 height will be on top of the ETRS89 one.
11. Pop the Baltic 1957 height from the stack as the active Z value.
12. Apply the inverse geoid model to transform the Baltic 1957 orthometric height into a ETRS89 ellipsoidal one.
13. Do not apply this step in the reverse direction ! We got what we want.

[busstoptaktik]

If I understand the code correctly, this is still an unnecessarily complex and restricted implementation. What I tried to describe in our mailing list thread was just a stack of double. Something like this:

1. The pipeline-operator includes a std:vector<double> thestack in its opaque section (or stack-allocated at the start of each pipeline invocation, for additional thread-robustness).

2. There's a stack_items object for each stack step. It contains a std::vector<size_t> of push/pop args, initialized from the push=n,m,... at construction time

3. Prior to executing the first step in the pipeline, pipeline_forward_4d calls pipeline->thestack.clear()

4. The stack push=n,m,... and stack pop=n,m,... operators can take any number of arguments, as long as they are all in the 4D range 1-4.

5. stack push does (pseudocode):

   for arg in args
       pipeline->thestack.push(coord[arg-1]);
   

6. stack pop does:

   if args.len() > thestack.len()
       coord[0] = coord[1] = coord[2] = coord[3] = DBL_MAX,
   else for arg in args
       coord[arg-1] = pipeline->thestack.pop();
   

   i.e. if there's too few elements on the stack, pop stomps
   on the operand.

7. stack swap does:

   if thestack.len() > 2 {
     auto first = thestack.pop;
     auto second = thestack.pop;
     thestack.push(first);
     thestack.push(second);
   }
   

   if there's less than 2 elements on the stack, swap is a no-op

8. If you want to swap more than two elements, also implement the Postscript roll operator: roll=4,2 rotates the upper 4 elements of the stack by two positions, keeping the internal order (i.e. swaps the upper pairs of 2D coordinates)

9. pipeline_inverse_4d works identically, but switches the meaning of push and pop, and reads their arguments in reverse order.

10. Balancing pushs and pops is not necessary, due to item 3 above. Although an unbalanced pipeline is invalid in one of the directions, where the stack will underflow, and hence stomp on the operands.

I have implemented this in Rust Geodesy here. See especially the tests and explanations from line 264 of the file.

And note that due to Geodesy working directly on full arrays of coordinates, this implementation needs allocation on its way, and is way more complicated than what can be achieved in PROJ.

(although the complication is a consecquence of the Geodesy overall goal of demonstrating a cleaner data flow architecture for PROJ. But you cannot win everywhere: Here it comes at the cost of inflating to 22 lines the single line PROJ pseudocode for push above).

[rouault]

What I tried to describe in our mailing list thread was just a stack of double

hum, so "+push=1,2" would need to be paired with "+pop=2,1" to restore arguments to their initial values?

While it may have some sense from a pure low-level computer science point of view, I don't like it as it is completely counter-intuitive for someone just wanting to save their (x,y) coordinates and restore them.

I believe I'm going to entirely give up on that topic, as I can't find anything I like enough.

[busstoptaktik]

hum, so "+push=1,2" would need to be paired with "+pop=2,1" to restore arguments to their initial values?

@rouault That is just a matter of selecting the syntax as you prefer:

You could just by definition choose that pop (or push) arguments are to be read backwards - then push=1,2/pop=1,2 does what you prefer.

I think this is an extremely good idea, because it provides the highly readable situation that

• if writing the stack contents sequentially from the bottom, going from left to right,
• then the picture of the stack fits exactly with both
  • the order of the arguments as written in the push=m,n,... case AND
  • the order of the arguments as written in the pop=m,n,...

Implementation-internally, for the inverse case, the code would still need to swap both the meaning of push/pop and the order arguments are read (i.e. push reads backwards, pop forwards) for consistency: Otherwise you get the wrong item at the top of the stack.

But from a UX perspective this is probably just a matter of "trusting that PROJ just does the right thing": I'm quite sure that I would probably only check my own pipelines in the forward direction, and rest assured that as long as pushs and pops are balanced, things will just work. And if not, I did something stupid, and DBL_MAX will tell me.

[busstoptaktik]

@rouault sorry for temporally messing up this thread by, by mistake, editing your message above. I have now found your original text in a github notification, and cleaned up as good as I could

[rouault]

I'm closing that as I have no plan to work on that in the foreseeable future