Can't reconstruct surface from point cloud with constant x, y or z coord

[ldouteau]

Hi Marco. Thanks a lot for Vedo, great project !

I'm struggling at the moment with the function reconstruct_surface. I am unable to reconstruct surfaces from point cloud that are normal to either x, y, or z vector. Looks like somewhere in VTK a division by zero happens. Relevant error: ERROR: In vtkLookupTable.cxx, line 173 vtkLookupTable (0000016FCC469B90): Bad table range: [1e+299, -1e+299]

I wonder if there's something that can be done in Vedo to prevent that or if it's just the maths that VTK implement that create this issue.

By the way i'm also wondering why reconstructed surface show wider range than that of the point cloud. Do you have any insight on this ?

Here is a piece of code that reproduces the issue with divisions by zero (ran it on vedo 2023.4.6). Let me know if you need more information

Have a nice weekend

import numpy as np
from colorcet import bmy
from scipy.spatial.transform import Rotation
from vedo import Points, show


def add_reconst(n, i_color=2):
    # build mesh & color
    p = Points(np_pts)
    m = p.reconstruct_surface()
    m.cmap(input_cmap=bmy, input_array=m.points()[:, i_color]).add_scalarbar()
    # fill result tables
    names.append(n)
    pts.append(p)
    mesh.append(m)


names = []
pts = []
mesh = []

# x, y grid
x, y = np.meshgrid(np.linspace(0, 1), np.linspace(0, 1))

# grid with constant z=0
np_pts = np.column_stack((x.ravel(), y.ravel(), np.zeros_like(x.ravel())))
add_reconst("z=0")

# grid with constant z=1
np_pts[:, 2] = 1
add_reconst("z=1")

# grid with varying z
np_pts[:, 2] = np.sin(np_pts[:, 0])
add_reconst("sin z")

# constant y
np_pts[:, 2] = np_pts[:, 1]
np_pts[:, 1] = 0
add_reconst("y=0", 1)

# constant x
np_pts[:, 1] = np_pts[:, 0]
np_pts[:, 0] = 0
add_reconst("x=0", 0)


# rotated plan
r = Rotation.from_euler("xyz", [30, 40, 60])
np_pts = r.apply(np_pts)
add_reconst("tilted")

# original data
np_pts = np.column_stack((x.ravel(), y.ravel(), np.zeros_like(x.ravel())))
np_pts[0, 2] += np.finfo(float).eps
add_reconst("z=0 with eps")

show([t for t in zip(names, pts, mesh)], N=len(mesh), axes=1)

[marcomusy]

Hi @ldouteau that's because you have a zero range in x (or y or z). Consider this:

import numpy as np
from colorcet import bmy
from vedo import Points, Grid, show


def add_reconst(n, i_color=2):
    # build mesh & color
    p = Points(np_pts)

    bb = list(p.bounds())
    if bb[0] == bb[1]:
        bb[1] += 1
        bb[0] -= 1
    if bb[2] == bb[3]:
        bb[3] += 1
        bb[2] -= 1
    if bb[4] == bb[5]:
        bb[5] += 1
        bb[4] -= 1

    m = p.reconstruct_surface(bounds=bb)
    m.cmap(bmy, input_array=m.points()[:, i_color]).add_scalarbar()
    names.append(n)
    pts.append(p)
    mesh.append(m)


names = []
pts = []
mesh = []

grid = Grid(res=(20, 20))

# grid with constant z=0
np_pts = grid.clone().points()
add_reconst("z=0")

# grid with constant z=1
np_pts = grid.clone().z(1).points()
add_reconst("z=1")

# grid with varying z
np_pts = grid.clone().points()
np_pts[:, 2] = np.sin(np_pts[:, 0])
add_reconst("sin z")

# constant y
np_pts = grid.clone().rotate_x(90).points()
add_reconst("y=0", 1)

# constant x
np_pts = grid.clone().rotate_y(90).points()
add_reconst("x=0", 0)

# rotated plan
np_pts = grid.clone().rotate_x(90).rotate_y(40).rotate_z(60).points()
add_reconst("tilted")

show([t for t in zip(names, pts, mesh)], N=len(mesh), sharecam=False, axes=1)

BTW i love:

show([t for t in zip(names, pts, mesh)], N=len(mesh), sharecam=False, axes=1)

I never thought of doing it that way! :)

[marcomusy]

By the way i'm also wondering why reconstructed surface show wider range than that of the point cloud. Do you have any insight on this ?

that is because the pointcloud is embedded in a volume that is necessarily a box aligned to the cartesian axes.