converting ply to postgis, script works but volume is invalid

Question

I'm trying to convert a ply to a polyhedron stored in postgis. In this example I have a cube from 0,0,0 to 1,1,1 the script below brings in a geometry to postgis, it is there, no problem, however I go to calculate the volume and it gives the following error:

query:

SELECT ST_Volume(st_geomfromtext) FROM public.test3

error:

ERROR:  PolyhedralSurface is invalid : inconsistent orientation of PolyhedralSurface detected at edge 1 (4-3) of polygon 11 : POLYHEDRALSURFACE(((1/1 1/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 0/1,1/1 0/1 1/1,0/1 0/1 0/1,1/1 0/1 0/1)),((1/1 0/1 0/1,1
SQL state: XX000

The object is a cube 0,0,0 to 1,1,1, here it is as an array of vertices and triangles.

vertices:

array([[ 1.,  1., -0.],
       [ 1.,  0., -0.],
       [ 0.,  0., -0.],
       [ 0.,  1., -0.],
       [ 1.,  1.,  1.],
       [ 0.,  1.,  1.],
       [ 0.,  0.,  1.],
       [ 1.,  0.,  1.],
       [ 1.,  1., -0.],
       [ 1.,  1.,  1.],
       [ 1.,  0.,  1.],
       [ 1.,  0., -0.],
       [ 1.,  0., -0.],
       [ 1.,  0.,  1.],
       [ 0.,  0.,  1.],
       [ 0.,  0., -0.],
       [ 0.,  0., -0.],
       [ 0.,  0.,  1.],
       [ 0.,  1.,  1.],
       [ 0.,  1., -0.],
       [ 1.,  1.,  1.],
       [ 1.,  1., -0.],
       [ 0.,  1., -0.],
       [ 0.,  1.,  1.]])

The triangles are defined as:

array([[ 0,  1,  2],
       [ 0,  2,  3],
       [ 4,  5,  6],
       [ 4,  6,  7],
       [ 8,  9, 10],
       [ 8, 10, 11],
       [12, 13, 14],
       [12, 14, 15],
       [16, 17, 18],
       [16, 18, 19],
       [20, 21, 22],
       [20, 22, 23]], dtype=int32)

I put together this conversion script to take a ply and add it to postgis:

import numpy as np
from open3d import *
import psycopg2

import dbconfig

def connect_db():
    global connection
    connection = psycopg2.connect(host=dbconfig.DATABASE_CONFIG['host'],
                                  user=dbconfig.DATABASE_CONFIG['user'],
                                  password=dbconfig.DATABASE_CONFIG['password'],
                                  dbname=dbconfig.DATABASE_CONFIG['dbname'])
    return connection
# mesh = read_triangle_mesh('C:/Users/garyn/PycharmProjects/pointcloudprocessor/tmp/contexts/99.526/396.ply')

mesh = read_triangle_mesh('C:/Users/garyn/PycharmProjects/pointcloudprocessor/tmp/contexts/cube3.ply')
verts = mesh.vertices
verts = np.asarray(verts)
tri = mesh.triangles
tri = np.asarray(tri)

data = ''
header = ("'POLYHEDRALSURFACE(")
for i in range(len(tri)):
# for i in range(0,2):
    x1 = (tri[i][0]) # 3
    y1 = (tri[i][1]) # 44
    z1 = (tri[i][2]) # 1
    x_coords1 = verts[x1][0]
    y_coords1 = verts[y1][0]
    z_coords1 = verts[z1][0]
    x_coords2 = verts[x1][1]
    y_coords2 = verts[y1][1]
    z_coords2 = verts[z1][1]
    x_coords3 = verts[x1][2]
    y_coords3 = verts[y1][2]
    z_coords3 = verts[z1][2]

    data += "((%s %s %s, %s %s %s, %s %s %s, %s %s %s))," % \
        (x_coords1, y_coords1, z_coords1, \
        x_coords2, y_coords2, z_coords2, \
        x_coords3, y_coords3, z_coords3, \
        x_coords1, y_coords1, z_coords1)
data = data[:-1]

projection = ")',32635)"
create_stmt = "CREATE TABLE test3 AS "
select_stmt = "SELECT ST_GeomFromText("
polyhedron = header + data + projection
query = create_stmt + select_stmt + polyhedron

conn = connect_db()
cur = conn.cursor()
cur.execute(query)
conn.commit()
cur.close()
conn.close()

Which results in:

CREATE TABLE test3 AS 
   SELECT ST_GeomFromText('POLYHEDRALSURFACE(((1.0 1.0 0.0, 1.0 0.0 0.0, -0.0 -0.0 -0.0, 1.0 1.0 0.0)),((1.0 0.0 0.0, 1.0 0.0 1.0, -0.0 -0.0 -0.0, 1.0 0.0 0.0)),((1.0 0.0 0.0, 1.0 1.0 0.0, 1.0 1.0 1.0, 1.0 0.0 0.0)),((1.0 0.0 1.0, 1.0 0.0 0.0, 1.0 1.0 1.0, 1.0 0.0 1.0)),((1.0 1.0 1.0, 1.0 1.0 0.0, -0.0 1.0 1.0, 1.0 1.0 1.0)),((1.0 1.0 1.0, 1.0 0.0 0.0, -0.0 1.0 -0.0, 1.0 1.0 1.0)),((1.0 1.0 0.0, 0.0 0.0 0.0, -0.0 1.0 1.0, 1.0 1.0 0.0)),((1.0 0.0 0.0, 0.0 0.0 0.0, -0.0 1.0 -0.0, 1.0 0.0 0.0)),((0.0 0.0 0.0, 0.0 0.0 1.0, -0.0 1.0 1.0, 0.0 0.0 0.0)),((0.0 0.0 0.0, 0.0 1.0 1.0, -0.0 1.0 -0.0, 0.0 0.0 0.0)),((1.0 1.0 0.0, 1.0 1.0 1.0, 1.0 -0.0 -0.0, 1.0 1.0 0.0)),((1.0 0.0 0.0, 1.0 1.0 1.0, 1.0 -0.0 1.0, 1.0 0.0 0.0)))',32635)

It looks ok to me and postgis accepts it as a polyhedral surface, but how do I make sure the cube is constructed correctly? P.S. and postgis 3D viewers out there, I'm doing this blind.

So I reconstructed the cube in blender, the 12 faces are correct, but some overlap, overs transect the cube internally. So it is clear I have mis-understood the relationship between the vertices and the triangles. — Spatial Digger, Mar 22 '19 at 12:09

score 0 · Answer 1 · answered Mar 22 '19 at 14:16

I solved it, I had the ordering wrong, here is the correct ordering.

data = ''
header = ("'POLYHEDRALSURFACE(")
for i in range(len(tri)):
# for i in range(0,2):
    x1 = (tri[i][0]) # 3
    y1 = (tri[i][1]) # 44
    z1 = (tri[i][2]) # 1
    coords1 = verts[x1][0]
    coords2 = verts[x1][1]
    coords3 = verts[x1][2]

    coords4 = verts[y1][0]
    coords5 = verts[y1][1]
    coords6 = verts[y1][2]

    coords7 = verts[z1][0]
    coords8 = verts[z1][1]
    coords9 = verts[z1][2]

    data += "((%s %s %s, %s %s %s, %s %s %s, %s %s %s))," % \
        (
        coords1, coords2, coords3,
        coords4, coords5, coords6,
        coords7, coords8, coords9,
        coords1, coords2, coords3
        )
data = data[:-1]

converting ply to postgis, script works but volume is invalid

1 Answers1