3

It is written in docs that ST_ApproximateMedialAxis doesn't drop z-index or third coordinate, but in practice it does drop z-index and returns 2D geometry.

I pass 3D geometry into ST_ApproximateMedialAxis and expect to get 3D geometry as a result.

What I tried:

SELECT ST_AsText(ST_ApproximateMedialAxis(ST_GeomFromText('POLYGON((0 0 10, 0 50 10, 10 50 10, 10 0 10, 0 0 10))')))

+-----------------------------+
| st_astext                   |
|-----------------------------|
| MULTILINESTRING((5 5,5 45)) |
+-----------------------------+


SELECT ST_AsText(ST_ApproximateMedialAxis(ST_GeomFromText('POLYGON Z ((0 0 10, 0 50 10, 10 50 10, 10 0 10, 0 0 10))')))

+-----------------------------+
| st_astext                   |
|-----------------------------|
| MULTILINESTRING((5 5,5 45)) |
+-----------------------------+


SELECT ST_AsEWKT(ST_ApproximateMedialAxis(GeomFromEWKT('POLYGON((0 0 10, 0 50 10, 10 50 10, 10 0 10, 0 0 10))')))

+-----------------------------+
| st_asewkt                   |
|-----------------------------|
| MULTILINESTRING((5 5,5 45)) |
+-----------------------------+


SELECT ST_AsEWKT(ST_ApproximateMedialAxis(GeomFromEWKT('POLYGON Z((0 0 10, 0 50 10, 10 50 10, 10 0 10, 0 0 10))')))

+-----------------------------+
| st_asewkt                   |
|-----------------------------|
| MULTILINESTRING((5 5,5 45)) |
+-----------------------------+

What am I doing wrong?

  • 2
    I think you are correct, I get the same answers, and looking at the source code I see no evidence of 3D versions of this function. It is possibly worth asking on Postgis mail group or irc. – John Powell Feb 11 at 12:11
  • 1
    @JohnPowell indeed, straightSkeleton transforms to and validates 2D polygons only. – ThingumaBob Feb 11 at 12:57
  • @ThingumaBob. Which implies that the docs are wrong, no? – John Powell Feb 11 at 13:04
  • 1
    One possible solution to this would be to rotate your polygon so that it is on the x-y plane, calcuate the medial axis or straight skeleton, add a default z coordinate, with ST_Force3D, and then do the inverse rotation. As the polygons should be planar, this ought to give the correct result. In your case, as the z coordinates are always 10, you simply need to do SELECT ST_Translate(ST_Force3D(ST_ApproximateMedialAxis(geom)), 0, 0, 10); – John Powell Feb 11 at 13:11
  • @JohnPowell I'd say so, yes. I just had that with another function (could safely run with 3D coordinates but Z index is dropped) but can't remember which one... – ThingumaBob Feb 11 at 13:16
4

So, looking at the source code for the ST_ApproximateMedialAxis and ST_StraightSkeleton indicates that both do in fact drop the z-axis, ie, they work, but ignore the z values.

However, there is a workaround. If you consider that a medial axis is a LineString calculated from a planar polygon, all you need to do is rotate your POLYGONZ so that it lies in the x/y plane, calculute the medial axis and then do the inverse rotation. There is a function, ST_Affine that allows for simultaneous rotations and translations, though you need to calculate all the parameters.

The equation for the plane is:

Ax + By + Cz + D = 0

which can be found by using a cross product and once you have that, you can plug these numbers into the equation found in this answer from Math Stack Exchange, which will allow you to translate and rotate the polygon so it lies in the x/y plane. The equations are pretty hideous to look at, but are basically just plug and play.

Once you have calculated the coefficients for rotation and translation, a, b, c, d, e, f, g, h, i, xoff, yoff, zoff, you problem can be reduced to:

 SELECT 
     ST_Translate(-xoff, -yoff, -zoff(
        ST_Affine(a, d, g, b, e, h, c, f, i, 0, 0, 0,
           ST_Force3D(
              ST_ApproximateMedialAxis(
                 ST_Translate(xoff, yoff, zoff(
                    ST_Affine(a, b, c, d, e, f, g, h, i, 0, 0, 0, geom)
                 )
              )
           )
       )
 );

Note, you need the use ST_Force3D in order to add a z coordinate to the geometry representing the medial axis in the x/y plane -- this will default to 0, but get its correct value once the geometry is rotated and translated back to its original location. To reverse the transfrom, you need to do the inverse of the rotation, which conveniently is the transpose of the rotation matrix, and then do the reverse translation.

I leave it to you to calculate a, b, c, xoff, etc :D.

In your case, as the polygon has z = 10 in all cases, you could just run:

SELECT 
  ST_AsText( 
        ST_Translate(
              ST_Force3D(
                 ST_ApproximateMedialAxis(
                     ST_GeomFromText('POLYGON((0 0 10, 0 50 10, 10 50 10, 10 0 10, 0 0 10))')
                )
             )
       ,0,0,10)
   );

which gives:

MULTILINESTRING Z ((5 5 10,5 45 10))

  • Thank you John! It sounds reasonable. Will check it out! – Andrey Semakin Feb 11 at 14:06
  • 1
    No problem, sorry for being lazy, I couldn't be bothered to actually calculate the plane and rotation matrices :-) Quite a few 3D problems can actually be reduced to 2D problems via projections, rotations and the like. – John Powell Feb 11 at 14:15
  • +1, but I guess this is getting ridiculously intense for arbitrary polygons ,) one thing to the links: the SFCGAL approximateMedialAxis wraps around CGAL's Straight_skeleton_2 function, with some extra dangle removal. you are linking the automated unit test there. and cascading through the functions (e.g. here and here) it is quite obvious that the docs might need to get changed... – ThingumaBob Feb 11 at 14:19
  • @ThingumaBob. Lol. I have no control over the docs. Ridiculously intense, you say, would you like to write a plpgsql function that will rotate an arbitrary planar PolygonZ onto the x/y plane, thereby rendering it somewhat less intense :D. – John Powell Feb 11 at 14:27
  • 1
    Andrey, sorry, I forgot to to reverse a translation and a rotation, you have to separate them, and, the reverse rotation matrix is just the transpose. So, I have now separated this in the answer. – John Powell Feb 11 at 16:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.