2

I need to compound several sfcgal expressions to create the overall effect of:

SELECT 
   ST_Volume(
       ST_MakeValid(
            ST_MakeSolid(
                 ST_3DDifference(
                   ST_3DUnion(PolyhedralSurface1, PolyhedralSurface2), 
                        PolyhedralSurface1)))) ...;

However, the nested functions fail either due to incompatible i/o types, or invalid geometric structures.

    SELECT ST_AsText(ST_3DUnion(p1,p2))
     FROM   
        (SELECT 
            (ST_GeomFromText('POLYHEDRALSURFACEZ( 
            ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0)),
            ((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0)), 
            ((0 0 0, 1 0 0, 1 0 1, 0 0 1, 0 0 0)),
            ((1 1 0, 1 1 1, 1 0 1, 1 0 0, 1 1 0)),
            ((0 1 0, 0 1 1, 1 1 1, 1 1 0, 0 1 0)), 
            ((0 0 1, 1 0 1, 1 1 1, 0 1 1, 0 0 1)) )')) AS p1,
            (ST_GeomFromText('POLYHEDRALSURFACEZ( 
            ((0.5 0 0, 0.5 0 1, 0.5 1 1, 0.5 1 0, 0.5 0 0)),
            ((0.5 0 0, 0.5 1 0, 1.5 1 0, 1.5 0 0, 0.5 0 0)), 
            ((0.5 0 0, 1.5 0 0, 1.5 0 1, 0.5 0 1, 0.5 0 0)),
            ((1.5 1 0, 1.5 1 1, 1.5 0 1, 1.5 0 0, 1.5 1 0)),
            ((0.5 1 0, 0.5 1 1, 1.5 1 1, 1.5 1 0, 0.5 1 0)), 
            ((0.5 0 1, 1.5 0 1, 1.5 1 1, 0.5 1 1, 0.5 0 1)) )'))  AS p2) AS t:

returns:

    "GEOMETRYCOLLECTION Z (TIN Z (
    ((-0 -0 -0,-0 -0 1,-0 1 1,-0 -0 -0)),
    ((-0 1 -0,-0 -0 -0,-0 1 1,-0 1 -0)),
    ((0.5 0.5 -0,-0 1 -0,0.5 1 -0,0.5 0.5 -0)),
    ((0.5 0.5 -0,-0 -0 -0,-0 1 -0,0.5 0.5 -0)),
    ((1 1 -0,0.5 0.5 -0,0.5 1 -0,1 1 -0)),
    ((0.5 -0 -0,-0 -0 -0,0.5 0.5 - (...)"

Then:

    SELECT ST_3DDifference(ST_3DUnion(p1,p2),p1)

    FROM    (SELECT (ST_GeomFromText('POLYHEDRALSURFACEZ( 
    ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0)),
    ((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0)), 
    ((0 0 0, 1 0 0, 1 0 1, 0 0 1, 0 0 0)),
    ((1 1 0, 1 1 1, 1 0 1, 1 0 0, 1 1 0)),
    ((0 1 0, 0 1 1, 1 1 1, 1 1 0, 0 1 0)), 
    ((0 0 1, 1 0 1, 1 1 1, 0 1 1, 0 0 1)) )')) AS p1,

    (ST_GeomFromText('POLYHEDRALSURFACEZ( 
    ((0.5 0 0, 0.5 0 1, 0.5 1 1, 0.5 1 0, 0.5 0 0)),
    ((0.5 0 0, 0.5 1 0, 1.5 1 0, 1.5 0 0, 0.5 0 0)), 
    ((0.5 0 0, 1.5 0 0, 1.5 0 1, 0.5 0 1, 0.5 0 0)),
    ((1.5 1 0, 1.5 1 1, 1.5 0 1, 1.5 0 0, 1.5 1 0)),
    ((0.5 1 0, 0.5 1 1, 1.5 1 1, 1.5 1 0, 0.5 1 0)), 
    ((0.5 0 1, 1.5 0 1, 1.5 1 1, 0.5 1 1, 0.5 0 1)) )'))  AS p2) AS t;

creates the error:

    NOTICE:  During difference_3d(A,B) :

    NOTICE:    with A: GEOMETRYCOLLECTION(TIN(
    ((0/1 0/1 0/1,0/1 0/1 1/1,0/1 1/1 1/1,0/1 0/1 0/1)),
    ((0/1 1/1 0/1,0/1 0/1 0/1,0/1 1/1 1/1,0/1 1/1 0/1)),
    ((1/2 1/2 0/1,0/1 1/1 0/1,1/2 1/1 0/1,1/2 1/2 0/1)),
    ((1/2 1/2 0/1,0/1 0/1 0/1,0/1 1/1 0/1,1/2 1/2 0/1)),
    ((1/1 1/1 0/1,

    NOTICE:     and B: POLYHEDRALSURFACE(
    ((0/1 0/1 0/1,0/1 0/1 1/1,0/1 1/1 1/1,0/1 1/1 0/1,0/1 0/1 0/1)),
    ((0/1 0/1 0/1,0/1 1/1 0/1,1/1 1/1 0/1,1/1 0/1 0/1,0/1 0/1 0/1)),
    ((0/1 0/1 0/1,1/1 0/1 0/1,1/1 0/1 1/1,0/1 0/1 1/1,0/1 0/1 0/1)),
    ((1/1 1/1 0/1,1/1 1/1 1/1,1/1 0/1 1...


    ERROR:  GeometryCollection is invalid : TriangulatedSurface 0 is invalid: self intersects : GEOMETRYCOLLECTION(TIN
    (((0/1 0/1 0/1,0/1 0/1 1/1,0/1 1/1 1/1,0/1 0/1 0/1)),
    ((0/1 1/1 0/1,0/1 0/1 0/1,0/1 1/1 1/1,0/1 1/1 0/1)),
    ((1/2 1/2 0/1,0/1 1/1 0/1,1/2 1/1 0/1,1/2 1/...

    ********** Error **********

    ERROR: GeometryCollection is invalid : TriangulatedSurface 0 is invalid: self intersects : GEOMETRYCOLLECTION(TIN(
    ((0/1 0/1 0/1,0/1 0/1 1/1,0/1 1/1 1/1,0/1 0/1 0/1)),
    ((0/1 1/1 0/1,0/1 0/1 0/1,0/1 1/1 1/1,0/1 1/1 0/1)),
    ((1/2 1/2 0/1,0/1 1/1 0/1,1/2 1/1 0/1,1/2 1/
    SQL state: XX000

Hence the ST_Union command appears to return a valid result, although this is typed first as a TIN and then as a geometrycollection.

But then these as a direct input to ST_3DDifference is invalid. This seems to imply that the outputs of the inner sfcgal functions are not compatible as inputs to the outer functions ... ???

4
  • I don't know what shapes your polyhedral surfaces are, but if I run a query like: WITH cube1 (geom) AS (SELECT ST_Extrude(ST_Force3D(ST_Expand(ST_Makepoint(0, 0), 1, 1)), 0,0,2) ) , cube2 (geom) AS (SELECT ST_Translate(geom, 0.5, 0.5, 0.5) FROM cube1) select ST_AsText(ST_3dIntersection(c1.geom, c2.geom)) FROM cube1 c1, cube2 c2; ie, the intersection of two cubes, one offset from the other, I get a PolyhedralSurface Z back. I don't know why you are getting TINs. It is because they are coincident in certain places? Apr 30 '18 at 14:36
  • Could you possibly indent/format your SQL a bit. It is somewhat hard to read. Apr 30 '18 at 14:46
  • The shapes are largely overlapping and hence somewhat coincident in volume. p1 and p2 are both 3D, 1 x 1 x 1 cubes aligned with the x, y, z axes. p2 is translated by 0.5 units along the x-axis. Being coincident suits the sfcgal operators and is consistent with the postgis documentation examples for the individual functions. ST_3DDifference should therefore return a block that is half the length in x, while the ST_Volume function should return 0.5 cubic units.
    – craigL
    May 1 '18 at 1:27
  • Thanks for the formatting, I can now read it :D. As you can see, I have put a mega comment as an answer. I am fairly sure it is to do with the ordering of your sides, as I can replicate what you are doing, without errors or TINs. May 1 '18 at 10:37
1

I am posting this as an answer, even though it is not a full answer, as it lacks a complete explanation of why what you are doing isn't working, but this approach does work (and because it is way too long for a comment).

I suspect one of the issues is around the cube construction (I am guessing you hand crafted your coordinates).

If you run the following, which makes a cube with volume 1, centred on (0.5, 0.5, 0.5), by using ST_Expand, ST_Force3D and ST_Extrude, as with yours:

WITH 
  cube1 (geom) AS (
     SELECT 
         ST_Extrude(
             ST_Force3D(
                 ST_Expand(
                     ST_Makepoint(0.5, 0.5), 
                  0.5, 0.5)), 
         0,0,1)
       )
      SELECT ST_AsText(geom) 
        FROM cube1;

you get:

 POLYHEDRALSURFACE Z (
     ((0 0 0,0 1 0,1 1 0,1 0 0,0 0 0)),
     ((0 0 1,1 0 1,1 1 1,0 1 1,0 0 1)),
     ((0 0 0,0 0 1,0 1 1,0 1 0,0 0 0)),
     ((0 1 0,0 1 1,1 1 1,1 1 0,0 1 0)),
     ((1 1 0,1 1 1,1 0 1,1 0 0,1 1 0)),
     ((1 0 0,1 0 1,0 0 1,0 0 0,1 0 0)))

which has a different ordering of sides.

If we now offset this 0.5 along the x-axis, by using ST_Translate, and do the 3D union, we get:

WITH 
  cube1 (geom) AS (
     SELECT 
         ST_Extrude(
             ST_Force3D(
                 ST_Expand(
                     ST_Makepoint(0.5, 0.5), 
                  0.5, 0.5)), 
         0,0,1)
       ),
  cube2 (geom) AS (
       SELECT 
             ST_Translate(c1.geom, 0.5, 0, 0) 
         FROM cube1 c1)
 SELECT 
     ST_AsText(ST_3DUnion(c1.geom, c2.geom)) 
   FROM cube1 c1, cube2 c2;

which yields a Polyhedral Surface Z, and not a TIN. The first few sides of which are:

POLYHEDRALSURFACE Z (
      ((1 0 0,0.5 -0 0.5,0.5 0 0,1 0 0)),
      ((0.5 -0 0.5,1 0 0,1 -0 0.5,0.5 -0 0.5)),
      ((0.5 0 0,0.5 -0 0.5,0 0 0,0.5 0 0)),
      ((1 0 0,0.5 0 0,1 0.5 -0,1 0 0)),........

And the volume is 1.5, as you would expect.

Finally, if we run the 3D difference of the 3D union:

WITH 
  cube1 (geom) AS (
     SELECT 
         ST_Extrude(
             ST_Force3D(
                 ST_Expand(
                     ST_Makepoint(0.5, 0.5), 
                  0.5, 0.5)), 
         0,0,1)
       ),
  cube2 (geom) AS (
       SELECT 
             ST_Translate(c1.geom, 0.5, 0, 0) 
         FROM cube1 c1)
      SELECT 
           ST_AsText(ST_3DDifference(ST_3DUnion(c1.geom, c2.geom), c1.geom))  
        FROM cube1 c1, cube2 c2;

we get another polyhedral surface:

POLYHEDRALSURFACE Z (
       ((1 0 0,1 -0 0.5,1 0.5 -0,1 0 0)),
       ((1 0.5 -0,1.5 0 0,1 0 0,1 0.5 -0)),
       ((1 -0.5,1 0 1,1 0.5 -0,1 -0 0.5)),....

So, I realize that this doesn't exactly answer you question, but it does suggest that if you use sfcgal and Postgis functions to contruct the cubes, the resulting intersections, unions, differences behave as expected, ie, there is no error, and Polyhedral Surfaces instead of TINs are returned. I have had a few of the errors you describe, but have not yet managed to get a clear view on how the ordering of sides needs to be in a Polyhedral Surface, though this appears to be key.

4
  • Thanks for this, it's a plausible explanation that the problem is in the ordering of sides. The difficulty though is that the cube is only a simple example to try to get it working. The larger objective is to conduct the spatial algebra on point clouds converted to meshes, so the geometries are potentially highly complex and don't have simple geometrical shapes.
    – craigL
    May 2 '18 at 3:45
  • Yes, I am working on similar problems, and have also found some of the issues you describe. It might be worth posting to Postgis or sfcgal list, as there aren't many people on here with 3d Postgis experience yet. May 2 '18 at 6:30
  • There is one other possible approach to your problem with is ST_SnapToGrid. In Geos, there is a similar error about self-intersections and the like that is caused by arbitrary and spurious precision (16 decimal places, kind of thing) and can be solved using this function. I haven't investigated too far with the 3D functionality yet, but it is conceivable it might help. I have seen errors like yours when doing intersections of cylinders, and extruded buildings for line of sight purposes. Please let me know if you solve this problem. May 2 '18 at 9:08
  • Ignore me. I have checked ST_SnapToGrid and it works on PolygonZ, but not on Polyhedral Surfaces or TINs. So, you would have to dump all the sides, snap them, and reassemble. There must be a better way. May 2 '18 at 11:21
1

The answer is to send the POLYGONZ text strings through ST_MakeSolid before performing any other operations on them. E.g.

 SELECT ST_Volume(ST_3DDifference(stopeDesign,ST_3DIntersection(s1,s2)))

 FROM   (SELECT (ST_MakeSolid(

 'POLYHEDRALSURFACEZ( 
                            ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0)),
                            ((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0)), 
                            ((0 0 0, 1 0 0, 1 0 1, 0 0 1, 0 0 0)),
                            ((1 1 0, 1 1 1, 1 0 1, 1 0 0, 1 1 0)),
                            ((0 1 0, 0 1 1, 1 1 1, 1 1 0, 0 1 0)), 
                            ((0 0 1, 1 0 1, 1 1 1, 0 1 1, 0 0 1)) )')) AS s1,

 (ST_MakeSolid('POLYHEDRALSURFACEZ( 
                            ((0.5 0 0, 0.5 0 1, 0.5 1 1, 0.5 1 0, 0.5 0 0)),
                            ((0.5 0 0, 0.5 1 0, 1.5 1 0, 1.5 0 0, 0.5 0 0)), 
                            ((0.5 0 0, 1.5 0 0, 1.5 0 1, 0.5 0 1, 0.5 0 0)),
                            ((1.5 1 0, 1.5 1 1, 1.5 0 1, 1.5 0 0, 1.5 1 0)),
                            ((0.5 1 0, 0.5 1 1, 1.5 1 1, 1.5 1 0, 0.5 1 0)), 
                            ((0.5 0 1, 1.5 0 1, 1.5 1 1, 0.5 1 1, 0.5 0 1)) )'))  AS s2 ) AS t;

This yields the correct answer of volume = 0.5.

Incidentally, I used a visualisation tool for POLYHEDRALSURFACEZ and TINZ, which made an enormous difference in being able to see what was happening: in the original form, the ST_Union performed on two geometries defined via ST_GeomFromText('POLYHEDRALSURFACEZ(... was creating a unified result but not eliminating interior surfaces from the unified form, which thereby violated the simple geometry criteria.

3
  • Which viz tool did you use, as a matter of interest? May 4 '18 at 7:01
  • I have used ST_MakeSolid prior to ST_Volume, which seems reasonable, but the docs are not at all clear about using it for unions and intersections in 3D. There is an implication that it is automatic. Anyway, good find. May 4 '18 at 7:54
  • The viz tool was made by one of our programmers using Unity3D. Very simple but extremely useful for debugging postgis/sfcgal expressions. I'll see if we can put it online somewhere.
    – craigL
    May 6 '18 at 0:01

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.