How to filter mesh after Poisson reconstruction in CloudComPy?

[joshua19]

I am trying to automate the creation of meshes from point clouds so I can do it quickly and in a consistent way for many point clouds in my project. I am using the Python connector CloudComPy to do it (see code below), but I cannot figure out how to filter the "fake" points using the density scalar field that is output by the Poisson reconstruction. In the GUI, I would typically decide how much of the mesh I want to include by adjusting the range in the SF display params window and then filter by value on the density (edit -> scalar fields -> filter by value). I can't seem to figure out how to do this via the Python connector, though. Is it possible?

I can filter the mesh's associated point cloud, but I cannot assign the filtered point cloud back to the mesh, so I just get a copy of the filtered point cloud. I can make a new mesh from that point cloud via triangulation, but that messes up the colors really badly, which is problematic for my application. I also tried reconstructing the mesh's associated point cloud from the filtered point cloud using the fromNParray, but this just creates malformed mesh files. For now I am just calling the filter via a system call (the workaround solution at the end of the Python code below) to the CloudCompare command line tool, which does work but can cause problems when automating on a larger scale.

Code: Select all

#!/usr/bin/env python3

import argparse
import multiprocessing
import numpy
import os
import subprocess

import cloudComPy as cc
import cloudComPy.PoissonRecon
from gendata import getSampleCloud2, getSamplePoly, dataDir

# Initialize CloudCompare
cc.initCC()

parser = argparse.ArgumentParser(description="Perform Poisson Surface Reconstruction on a point cloud")
parser.add_argument("input_file", help="Absolute path to the input point cloud file")
parser.add_argument("--output_file", default=None, help="Absolute path to the output .ply file")
args = parser.parse_args()

# Convert input file path to absolute path
input_file_abs = os.path.abspath(args.input_file)

# Generate output file name
base, ext = os.path.splitext(input_file_abs)

output_file_abs = "".join(input_file_abs.split(".")[:-1]) + ".ply"
if( args.output_file is not None ):
    output_file = os.path.abspath(args.output_file)

print(f"{input_file_abs} --> {output_file_abs}")

print(f"loading input file {input_file_abs}...")
cloud = cc.loadPointCloud(input_file_abs)

print(f"shifting cloud to center around (x=0, y=0) and min(z) = 0")
boundingBox = cloud.getOwnBB()
cmin = numpy.array(boundingBox.minCorner())
cmax = numpy.array(boundingBox.maxCorner())
center = 0.5 * (cmin + cmax)
translation = (-center[0], -center[1], -cmin[0]) 
cloud.translate(translation)

radius = 0.1
print(f"calculating verticality, {radius=}")
cc.computeFeature(cc.GeomFeature.Verticality, radius, [cloud])

radius = 1
print(f"calculating planarity, {radius=}")
cc.computeFeature(cc.GeomFeature.Planarity, radius, [cloud])
print(f"calculating omnivariance, {radius=}")
cc.computeFeature(cc.GeomFeature.Omnivariance, radius, [cloud])
print(f"calculating surface_variance, {radius=}")
cc.computeFeature(cc.GeomFeature.SurfaceVariation, radius, [cloud])

print("computing normals...")
cc.computeNormals([cloud], defaultRadius=0.1, model=cc.LOCAL_MODEL_TYPES.QUADRIC, orientNormals=True, preferredOrientation=cc.Orientation.PLUS_Z)

print("orienting normals...")
cloud.orientNormalsWithMST(octreeLevel=10)

print("doing Poisson reconstruction...")
mesh = cc.PoissonRecon.PR.PoissonReconstruction(pc=cloud, threads=multiprocessing.cpu_count(), density=True, withColors=True, depth=10, samplesPerNode=2)

#print("filtering reconstructed mesh by point density...")
#density_scalar_field_index = 0 # we assume there are no other scalar fields yet 
#density_scalar_field_index = next(i for i, key in enumerate(mesh.getAssociatedCloud().getScalarFieldDic().keys()) if 'density' in key.lower())
#density_scalar_field = mesh.getAssociatedCloud().getScalarField(density_scalar_field_index)
#mesh.getAssociatedCloud().setCurrentOutScalarField(density_scalar_field_index)
#filtered_cloud = cc.filterBySFValue( 8, density_scalar_field.getMax(), mesh.getAssociatedCloud())
#filtered_mesh = cc.ccMesh.triangulate(filtered_cloud, cc.TRIANGULATION_TYPES.DELAUNAY_2D_AXIS_ALIGNED) # this works but messes up the colors
#filtered_mesh.getAssociatedCloud().interpolateColorsFrom(filtered_cloud) # doesn't help
#ret = cc.SaveMesh(filtered_mesh, "filtered_test.ply")

print("interpolating all scalar fields from the pointcloud to the mesh...")
parameters = cc.interpolatorParameters()
parameters.method = cc.INTERPOL_METHOD.RADIUS
parameters.algos = cc.INTERPOL_ALGO.NORMAL_DIST
parameters.radius = 0.1
parameters.sigma = 0.04
ret = cc.interpolateScalarFieldsFrom(mesh.getAssociatedCloud(), cloud, list(range(len(cloud.getScalarFieldDic()))), parameters)

print("setting density as active scalar field for later filtering...")
density_scalar_field_index = 0 # we assume there are no other scalar fields yet 
density_scalar_field_index = next(i for i, key in enumerate(mesh.getAssociatedCloud().getScalarFieldDic().keys()) if 'density' in key.lower())
density_scalar_field = mesh.getAssociatedCloud().getScalarField(density_scalar_field_index)
mesh.getAssociatedCloud().setCurrentOutScalarField(density_scalar_field_index)
mesh.showColors(True)

print(f"saving mesh to {output_file_abs}")
ret = cc.SaveMesh(mesh, output_file_abs)
if ret:
    print(f"Mesh saved successfully to {output_file_abs}")
else:
    print(f"Failed to save mesh to {output_file_abs}")

# I hate to do this but it's the easiest solution for now
command = f"CloudCompare -SILENT -O {output_file_abs} -AUTO_SAVE OFF -FILTER_SF 8 MAX -M_EXPORT_FMT PLY -PLY_EXPORT_FMT ASCII -SAVE_MESHES"
command = command.split(" ")
result = subprocess.run(command, capture_output=True, text=True)

#import code
#code.interact(local=locals())

[daniel]

Questions about CloudCompy should be asked directly on the project's github 'Issues' tracker.

In CloudCompare's code, there's a dedicated method to filter a mesh based on the scalar field (you first hide the vertices based on their scalar value with ccPointCloud::hidePointsByScalarValue, and then you create the new mesh instance with ccMesh::createNewMeshFromSelection).