How to Fix Inconsistent Data Conversions in SFCGAL Expressions

Question

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 ... ???

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? — John Powell, Commented Apr 30, 2018 at 14:36
Could you possibly indent/format your SQL a bit. It is somewhat hard to read. — John Powell, Commented Apr 30, 2018 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, Commented May 1, 2018 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. — John Powell, Commented May 1, 2018 at 10:37

John Powell · Accepted Answer · 2018-05-01 12:01:20Z

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.

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, Commented May 2, 2018 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. — John Powell, Commented May 2, 2018 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. — John Powell, Commented May 2, 2018 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. — John Powell, Commented May 2, 2018 at 11:21

craigL · Accepted Answer · 2018-05-24 00:29:47Z

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.

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. — John Powell, Commented May 4, 2018 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, Commented May 6, 2018 at 0:01

Stack Exchange Network

How to Fix Inconsistent Data Conversions in SFCGAL Expressions

2 Answers 2

Your Answer

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

Hot Network Questions

How to Fix Inconsistent Data Conversions in SFCGAL Expressions

2 Answers 2

Your Answer

Sign up or log in

Post as a guest

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

Related

Hot Network Questions

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