6. Some use cases

6.1. Create a point cloud with numpy

The way to set coordinates of a cloud from a Numpy array have been described in Read, modify or create cloud coordinates with Numpy and the method to create scalar fields from Numpy arrays have been exposed in Read, modify or create a scalar field with Numpy.

A first example generates a kind of spherical cloud with sinusoidal altitude fluctuations. One scalar field is defined with the local altitude. Nodes are randomly generated.


# --- generate a set of coords and a scalar field

npts = 10000000
phi = 2*np.pi*np.random.random((npts))
theta = 2*np.pi*np.random.random((npts))
r = 5 + 0.3*np.sin(2*2*np.pi*phi + 3*2*np.pi*theta)
x = np.float32(r*np.sin(phi)*np.cos(theta))
y = np.float32(r*np.sin(phi)*np.sin(theta))
z = np.float32(r*np.cos(phi))
coords = np.column_stack((x,y,z))
dr = np.float32(np.sqrt(x*x + y*y + z*z) -5)

# --- create the pointCloud, add the scalar field

cl = cc.ccPointCloud("boule")
cl.coordsFromNPArray_copy(coords)
cl.addScalarField("delta")
sf = cl.getScalarField(0)
sf.fromNpArrayCopy(dr)

# --- save the point cloud

res = cc.SavePointCloud(cl, os.path.join(dataDir, "boule.bin"))

The above code snippet is from test017.py.

Another similar example in which we create a kind of 2.5D wavelet and add the altitude scalar field after the cloud creation with exportCoordToSF(). Nodes are randomly generated.


# --- generate a set of coords and a scalar field

npts = 1000000
h = 2.
x = np.float32(-5. + 10.*np.random.random((npts)))
y = np.float32(-5. + 10.*np.random.random((npts)))
z = np.float32(np.sin(h * np.sqrt(x**2 + y**2)) / np.sqrt(x**2 + y**2))
coords = np.column_stack((x,y,z))

# --- create the pointCloud

cloud = cc.ccPointCloud("wave_%d"%h)
cloud.coordsFromNPArray_copy(coords)

# --- add the scalar field

res = cloud.exportCoordToSF(False, False, True)

The above code snippet is from test025.py.

A third example draw a polyline from a bounding box obtained, for instance, with getOwnBB().

Firstly, we create a Numpy array from the corner coordinates of the bounding box. The nodes are ordered to build a polyline:


cmin = boundingBox.minCorner()
cmax = boundingBox.maxCorner()

x = np.float32(np.array((cmin[0], cmax[0], cmax[0], cmin[0], cmin[0], cmax[0], cmax[0], cmin[0])))
y = np.float32(np.array((cmin[1], cmin[1], cmin[1], cmin[1], cmax[1], cmax[1], cmax[1], cmax[1])))
z = np.float32(np.array((cmin[2], cmin[2], cmax[2], cmax[2], cmax[2], cmax[2], cmin[2], cmin[2])))
coords = np.column_stack((x,y,z))

Then, we create the cloud from the Numpy array, apply a transformation, and build the polyline.


cloud1 = cc.ccPointCloud("boundingBox1")
cloud1.coordsFromNPArray_copy(coords)
cloud1.applyRigidTransformation(transform1)

poly1 = cc.ccPolyline(cloud1)
poly1.addChild(cloud1)
poly1.addPointIndex(0, cloud1.size())
poly1.setClosed(True)

The above code snippets are from test026.py.

6.2. Cut a cloud or a mesh with a polyline

This corresponds to the CloudCompare GUI Interactive Segmentation Tool when used with an imported polygon.

6.2.1. Cut a cloud with a polyline

The method crop2D() is used to cut a cloud by a 2D polyline. The normal to the plan of the polyline can only be one of the cardinal axes oX, oY or oZ. The original cloud is not modified, a new cloud is created with either the nodes inside or the nodes outside.

First step, get a cloud and the polyline tool:


cloud = cc.loadPointCloud(getSampleCloud(5.0))
poly = cc.loadPolyline(getSamplePoly("poly1"))

Second step, close the polyline:


poly.setClosed(True)

Third step, cut:


cloudCropZ = cloud.crop2D(poly, 2, True)

The above code snippets are from test007.py.

6.2.2. Cut a mesh with a polyline

As with the cloud, the method crop2D() is used to cut a mesh by a 2D polyline. The normal to the plan of the polyline can only be one of the cardinal axes oX, oY or oZ. The original mesh is not modified, a new mesh is created with either the triangles inside or the triangles outside.

In this example, scalar fields and normals are defined, to check they are kept in the result.


cloud1 = cc.loadPointCloud(getSampleCloud2(3.0,0, 0.1))
cloud1.setName("cloud1")
cloud1.exportCoordToSF(True, True, True)

mesh1 = cc.ccMesh.triangulate(cloud1, cc.TRIANGULATION_TYPES.DELAUNAY_2D_AXIS_ALIGNED)
mesh1.setName("mesh1")
cc.computeNormals([mesh1], computePerVertexNormals=False)

poly = cc.loadPolyline(getSamplePoly("poly1"))
meshCropZ = mesh1.crop2D(poly, 2, True)

The above code snippet is from test039.py.

6.3. Compute distances between clouds and meshes

The idea is to build a scalarField giving, for each node in a cloud, the distance between this node and another cloud or mesh. The process can be long, several tools help to find an optimal set of parameters.

The CloudCompare wiki provides a good introduction to the distances computation.

See also :ref: cloud_comparison_M3C2.

6.3.1. Compute distances between a cloud and a mesh (C2M)

The method cloudComPy.DistanceComputationTools.computeApproxCloud2MeshDistance() gives an approximate solution, without too much cpu time. You get an approximate distance scalar field, and statistics (min, max, mean, variance, error max).


stats = cc.DistanceComputationTools.computeApproxCloud2MeshDistance(cloud, cylinder)
print(stats) # min, max, mean, variance, error max

The structure cloudComPy.Cloud2MeshDistancesComputationParams() is used to define the optimal parameters for accurate distance calculation. For instance, if the cloud and mesh are close, filling the parameter maxSearchDist can greatly speed up the process, but setting it to a non-zero value force a single thread processing. The result is a new scalar field in the cloud, with accurate distances to the mesh.


nbCpu = psutil.cpu_count()
bestOctreeLevel = cc.DistanceComputationTools.determineBestOctreeLevel(cloud,cylinder)
params = cc.Cloud2MeshDistancesComputationParams()
params.maxThreadCount = nbCpu
params.octreeLevel = bestOctreeLevel
cc.DistanceComputationTools.computeCloud2MeshDistances(cloud, cylinder, params)

The above code snippet is from test009.py.

6.3.2. Compute distances between two clouds (C2C)

The method cloudComPy.DistanceComputationTools.computeApproxCloud2CloudDistance() gives an approximate solution, without too much cpu time. You get an approximate distance scalar field, and statistics (min, max, mean, variance, error max).


stats = cc.DistanceComputationTools.computeApproxCloud2CloudDistance(cloud2ref, cloud3)
print(stats) # min, max, mean, variance, error max

The structure cloudComPy.Cloud2CloudDistancesComputationParams() is used to define the optimal parameters for accurate distance calculation. For instance, if the two clouds are close, filling the parameter maxSearchDist can greatly speed up the process, but setting it to a non-zero value force a single thread processing. The result is a new scalar field in the cloud, with accurate distances to the cloud.


nbCpu = psutil.cpu_count()
bestOctreeLevel = cc.DistanceComputationTools.determineBestOctreeLevel(cloud2ref, None, cloud3)
params = cc.Cloud2CloudDistancesComputationParams()
params.maxThreadCount = nbCpu
params.octreeLevel = bestOctreeLevel
cc.DistanceComputationTools.computeCloud2CloudDistances(cloud2ref, cloud3, params)

The above code snippet is from test010.py.

6.3.3. Compute split distances between two clouds (C2C)

The use case is similar to Compute distances between two clouds (C2C) but we need to split the distance components according to the cardinal directions. Here, we set the parameter setSplitDistances to cloud.size() to create 3 scalar fields of the cloud size. They are not associated to the cloud, this will be done when the distances are computed.


stats = cc.DistanceComputationTools.computeApproxCloud2CloudDistance(cloud,
                                               cylinder.getAssociatedCloud())
print(stats) # min, max, mean, variance, error max

nbCpu = psutil.cpu_count()
bestOctreeLevel = cc.DistanceComputationTools.determineBestOctreeLevel(cloud,
                                          None, cylinder.getAssociatedCloud())
params = cc.Cloud2CloudDistancesComputationParams()
params.maxThreadCount = nbCpu
params.octreeLevel = bestOctreeLevel
params.setSplitDistances(cloud.size()) # creates 3 scalar fields of cloud.size()
                                       # associated to the cloud on compute

cc.DistanceComputationTools.computeCloud2CloudDistances(cloud, 
                                 cylinder.getAssociatedCloud(), params)

The cloud contains four new scalar fields:

  • C2C absolute distances

  • C2C absolute distances (X)

  • C2C absolute distances (Y)

  • C2C absolute distances (Z)

The above code snippet is from test009.py.

6.4. Cloud registration

6.4.1. Cloud registration with the CloudCompare implementation of ICP algorithm

The CloudCompare wiki provides a good introduction to the alignment and registration process.

To test the registration with ICP (iterative ClosestPoint) we get two clouds with an exact overlap of 10% and we apply a small rotation-translation on one cloud.


cloud1 = cc.loadPointCloud(getSampleCloud(5.0))
cloud1.setName("cloud1")

cloud2ref = cc.loadPointCloud(getSampleCloud(5.0, 9.0))
cloud2ref.setName("cloud2_reference")

cloud2 = cloud2ref.cloneThis()
tr1 = cc.ccGLMatrix()
# --------------------  z -- y -- x  (rotation x 0.1, translation z 0.3)
tr1.initFromParameters(0.0, 0.0, 0.1, (0.0, 0.0, 0.3))
cloud2.applyRigidTransformation(tr1)
cloud2.setName("cloud2_transformed")

cc.SaveEntities([cloud1, cloud2ref, cloud2], os.path.join(dataDir, "clouds2.bin"))

Then we apply the ICP algorithm cloudComPy.ICP() which gives a transformation to apply to get the precise overlap of the clouds.


res=cc.ICP(data=cloud2, model=cloud1, minRMSDecrease=1.e-5,
           maxIterationCount=20, randomSamplingLimit=50000, removeFarthestPoints=False,
           method=cc.CONVERGENCE_TYPE.MAX_ITER_CONVERGENCE,
           adjustScale=False, finalOverlapRatio=0.1)
tr2 = res.transMat
cloud3 = res.aligned
cloud3.applyRigidTransformation(tr2)
cloud3.setName("cloud2_transformed_afterICP")

The above code snippets are from test010.py.

6.4.2. Cloud registration with the PCL plugin

The PCL plugin of CloudCompare wraps several methods of the PCL Library. The Fast Global Registration is one of these methods.

To test this tool, we start by cloning a cloud and applying a rotation on it. The algorithm requires to compute normals on clouds. We keep clones of the clouds for testing purposes.


cloud1 = cc.loadPointCloud(getSampleCloud2(3.0, 0, 0.1))

cloud2 = cloud1.cloneThis()
tr1 = cc.ccGLMatrix()
tr1.initFromParameters(0.1, (0., 0.1, 0.9), (0.,0.,0.))
cloud2.applyRigidTransformation(tr1)

cc.computeNormals([cloud1, cloud2])

cloud2ref = cloud2.cloneThis()
cloud2bis = cloud2ref.cloneThis()

We need to import the PCL plugin:


if cc.isPluginPCL():
    import cloudComPy.PCL
    

The PCL plugin gives access to the class FastGlobalRegistrationFilter


    fpcl = cc.PCL.FastGlobalRegistrationFilter()
    fpcl.setParameters(cloud1, [cloud2])
    res=fpcl.compute()
    

The above code snippets are from test038.py.

6.5. Cloud triangulation: create a mesh

The cloudComPy.ccMesh.triangulate() methods is adapted to create a Delaunay 2.5D mesh from a point cloud that represents a kind of 2.5D elevated surface.

When the best fitting plane of the cloud is axis aligned (dim=2 means Z axis):


cloud1 = cc.loadPointCloud(getSampleCloud2(3.0, 0, 0.1))
cloud1.setName("cloud1")

mesh1 = cc.ccMesh.triangulate(cloud1, cc.TRIANGULATION_TYPES.DELAUNAY_2D_AXIS_ALIGNED, dim=2)
mesh1.setName("mesh1")

Otherwise we tell the algorithm to find the best fitting plane:


cloud2 = cc.loadPointCloud(getSampleCloud2(3.0, 0, 0.1))
cloud2.setName("cloud2")
tr1 = cc.ccGLMatrix()
tr1.initFromParameters(0.5, (0.7, 0.7, 0.7), (0.,0.,0.))
cloud2.applyRigidTransformation(tr1)

mesh4 = cc.ccMesh.triangulate(cloud2, cc.TRIANGULATION_TYPES.DELAUNAY_2D_BEST_LS_PLANE)
mesh4.setName("mesh4")

The above code snippets are from test011.py.

6.6. Fitting ccPlane

The method Fit() allows to adjust a plane primitive on a cloud:


cloud1 = cc.loadPointCloud(getSampleCloud2(3.0,0, 0.1))
cloud1.setName("cloud1")
tr1 = cc.ccGLMatrix()
tr1.initFromParameters(0.2*math.pi, (1., 1., 1.), (0.0, 0.0, 10.0))
cloud1.applyRigidTransformation(tr1)

plane = cc.ccPlane.Fit(cloud1)
equation = plane.getEquation()
tr2 = plane.getTransformation()

The above code snippet is from test012.py.

See Plane for the plane methods.

See also 2D polygons (facets) if you want the 2D polygon delimiting the projection of the cloud on the fitting plane.

6.7. Searching points in neighbourhood

To find cloud nodes in the neighbourhood of a given point, the octree offers several methods, to search in a sphere, a cylinder, a box. See octrees for examples of use.

6.8. Cloud sampling, noise filter, filter by scalar field

The class cloudComPy.CloudSamplingTools offers several cloud sampling filters.

6.8.1. Noise filter

The noiseFilter() method is based on the distance to the approximate local surface. This filter removes points based on their distance relatively to the best fit plane computed on their neighbors.

The mandatory parameters are the cloud, the neighbourhood radius and the number of sigmas:

The result is a selection (ReferenceCloud) that must be transformed in ccPointCloud with partialClone()


cloud = cc.loadPointCloud(getSampleCloud(5.0))

refCloud = cc.CloudSamplingTools.noiseFilter(cloud, 0.04, 1.0) # selection on cloud
origCloud = refCloud.getAssociatedCloud()            # the original cloud ~ cloud
(noiseCloud, res) = origCloud.partialClone(refCloud) # ccPointCloud from selection, status
noiseCloud.setName("noiseCloud")

if refCloud.__class__ != cc.ReferenceCloud:
    raise RuntimeError
if refCloud.size() != 7489:
    raise RuntimeError
if res != 0:
    raise RuntimeError

6.8.2. Resample cloud spatially

The resampleCloudSpatially() method resamples a point cloud (process based on inter point distance) The cloud is resampled so that there is no point nearer than a given distance to other points. It works by picking a reference point, removing all points which are to close to this point, and repeating these two steps until the result is reached.

The mandatory parameters are the cloud and the distance.

The result is a selection (ReferenceCloud) that must be transformed in ccPointCloud with partialClone()


params = cc.SFModulationParams()
refCloud = cc.CloudSamplingTools.resampleCloudSpatially(cloud, 0.05, params)
(spatialCloud, res) = cloud.partialClone(refCloud)
spatialCloud.setName("spatialCloud")

if refCloud.size() != 55465:
    raise RuntimeError
if res != 0:
    raise RuntimeError

6.8.3. Subsample cloud randomly

The subsampleCloudRandomly() method subsamples a point cloud (process based on random selections). This is a very simple subsampling algorithm that simply consists in selecting “n” different points, in a random way.

The mandatory parameters are the cloud and the total number of points to keep.

The result is a selection (ReferenceCloud) that must be transformed in ccPointCloud with partialClone()


refCloud = cc.CloudSamplingTools.subsampleCloudRandomly(cloud, 50000)
(randomCloud, res) = cloud.partialClone(refCloud)
randomCloud.setName("randomCloud")

if refCloud.size() != 50000:
    raise RuntimeError
if res != 0:
    raise RuntimeError

6.8.4. resample cloud with octree at level

The resampleCloudWithOctreeAtLevel() method is a resampling algorithm is applied inside each cell of the octree. The different resampling methods are represented as an enumerator (see RESAMPLING_CELL_METHOD) and consist in simple processes such as replacing all the points lying in a cell by the cell center or by the points gravity center.

The mandatory parameters are the cloud, the octree level and the resampling method.


resOctrAlCloud  = cc.CloudSamplingTools.resampleCloudWithOctreeAtLevel(randomCloud, octreeLevel=5,
                                         resamplingMethod=cc.RESAMPLING_CELL_METHOD.CELL_CENTER)
resOctrAlCloud.setName("resOctrAlCloud")

if resOctrAlCloud.size() < 2050 or resOctrAlCloud.size() > 2100:
    raise RuntimeError

6.8.5. resample cloud with octree

The resampleCloudWithOctree() method is the same as resampleCloudWithOctreeAtLevel() method, apart the fact that instead of giving a specific octree subdivision level as input parameter, one can specify an approximative number of points for the resulting cloud (the algorithm will automatically determine the corresponding octree level).

The mandatory parameters are the cloud, the approximate number of points and the resampling method.


resOctrCloud = cc.CloudSamplingTools.resampleCloudWithOctree(randomCloud, newNumberOfPoints=5000,
                                         resamplingMethod=cc.RESAMPLING_CELL_METHOD.CELL_CENTER)
resOctrCloud.setName("resOctrCloud")

if resOctrCloud.size() < 1900 or resOctrCloud.size() > 8000:
    raise RuntimeError

6.8.6. Statistical Outliers Removal (SOR) filter

The sorFilter() method removes points based on their mean distance to their distance (by comparing it to the average distance of all points to their neighbors). It is equivalent to PCL Library StatisticalOutlierRemoval filter (see Removing outliers using a StatisticalOutlierRemoval filter)

The only mandatory parameter is the cloud.

The result is a selection (ReferenceCloud) that must be transformed in ccPointCloud with partialClone()


refCloud = cc.CloudSamplingTools.sorFilter(randomCloud)
(sorCloud, res) = randomCloud.partialClone(refCloud)
sorCloud.setName("sorCloud")

if refCloud.size() < 43000 or refCloud.size() > 45000:
    raise RuntimeError
if res != 0:
    raise RuntimeError

6.8.7. subsample cloud with octree at level

The subsampleCloudWithOctreeAtLevel() method applies a subsampling algorithm inside each cell of the octree. The different subsampling methods are represented as an enumerator (see SUBSAMPLING_CELL_METHOD) and consist in simple processes such as choosing a random point, or the one closest to the cell center.

The mandatory parameters are the cloud, the octree level and the subsampling method.

The result is a selection (ReferenceCloud) that must be transformed in ccPointCloud with partialClone()


refCloud = cc.CloudSamplingTools.subsampleCloudWithOctreeAtLevel(cloud, 5,
                    cc.SUBSAMPLING_CELL_METHOD.NEAREST_POINT_TO_CELL_CENTER)
(subOctLevCloud, res) = cloud.partialClone(refCloud)
subOctLevCloud.setName("subOctLevCloud")

if refCloud.size() != 2262:
    raise RuntimeError
if res != 0:
    raise RuntimeError

6.8.8. subsample cloud with octree

The subsampleCloudWithOctree() method is the same as subsampleCloudWithOctreeAtLevel() method, apart the fact that instead of giving a specific octree subdivision level as input parameter, one can specify an approximative number of points for the resulting cloud (the algorithm will automatically determine the corresponding octree level).

The mandatory parameters are the cloud, the approximate number of points and the subsampling method.

The result is a selection (ReferenceCloud) that must be transformed in ccPointCloud with partialClone()


refCloud = cc.CloudSamplingTools.subsampleCloudWithOctree(cloud, 25000,
                                 cc.SUBSAMPLING_CELL_METHOD.RANDOM_POINT)
(subOctreeCloud, res) = cloud.partialClone(refCloud)
subOctreeCloud.setName("subOctreeCloud")

if refCloud.size() != 36777:
    raise RuntimeError
if res != 0:
    raise RuntimeError

All the above code snippets are from test019.py.

6.8.9. filter by scalar field values

The cloudComPy.filterBySFValue() method creates a new point cloud by filtering points using the current out ScalarField (see cloudComPy.ccPointCloud.setCurrentOutScalarField()). It keeps the points whose ScalarField value is between the min and max parameters.

Here, we create a scalar field based on curvature:


radius = 0.03
res = cc.computeCurvature(cc.CurvatureType.GAUSSIAN_CURV, radius, [cloud])
nsf = cloud.getNumberOfScalarFields()
sfc = cloud.getScalarField(nsf - 1)

if sfc.getName() != "Gaussian curvature (0.03)":
    raise RuntimeError

Then, we apply the filter:


cloud.setCurrentOutScalarField(nsf - 1)
fcloud = cc.filterBySFValue(0.01, sfc.getMax(), cloud)
filteredSize = fcloud.size()

The above code snippets are from test003.py.

6.9. Generate histograms

To draw plots or to generate csv files from scalar field data, we use numpy, matplotlib and csv.

6.9.1. histogram as a figure

Here, we want to obtain a histogram of distances as a figure from a scalar field. We get a numpy array from the scalar field.


sf = cloud.getScalarField(cloud.getScalarFieldDic()['C2M absolute distances'])
asf = sf.toNpArray()

We need matplotlib and numpy.


import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import colors

Building a png histogram from a Numpy array is a standard use of matplotlib.


matplotlib.use('agg') # png images

(n, bins, patches) = plt.hist(asf, bins=256, density=1) # histogram for matplotlib
fracs = bins / bins.max()
norm = colors.Normalize(fracs.min(), fracs.max())
for thisfrac, thispatch in zip(fracs, patches):
    color = plt.cm.rainbow(norm(thisfrac))
    thispatch.set_facecolor(color)

plt.savefig(os.path.join(dataDir, "histogram.png"))

The above code snippets are from test022.py.

6.9.2. histogram as a csv file

Here, we want to obtain a histogram of distances as a csv file from a scalar field. We get a numpy array from the scalar field.


sf = cloud.getScalarField(cloud.getScalarFieldDic()['C2M absolute distances'])
asf = sf.toNpArray()

We need matplotlib and csv.


import numpy as np
import csv

Building a csv histogram from a Numpy array is a standard use of csv.


(n, bins) = np.histogram(asf, bins=256) # numpy histogram (without graphics)
with open(os.path.join(dataDir, "histogram.csv"), 'w') as f:
    writer = csv.writer(f, delimiter=';')
    writer.writerow(("Class", "Value", "Class start", "Class end"))
    for i in range(n.size):
        writer.writerow((i, n[i], bins[i], bins[i+1]))

The above code snippets are from test022.py.

6.10. Compute a 2.5D volume

This corresponds to the CloudCompare GUI 2.5D Volume

The cloudComPy.ComputeVolume25D() method computes a 2.5D volume between a cloud and a ground plane,or two clouds, following a given direction (X, Y or Z). If only one cloud is given, the direction (X, Y or Z) defines the normal to the plane used for calculation. The calculation uses a cartesian grid, you provide the gridStep.

The cloudComPy.ReportInfoVol structure is completed by the calculation and the following information:

volume:

the resulting volume

addedVolume:

the positive part of the volume (ceil > floor)

removedVolume:

the negative part of the volume (ceil < floor)

surface:

the section of the point cloud along in the given direction

matchingPercent:

percentage of the section matching ceil and floor

ceilNonMatchingPercent:

percentage of the ceil section non matching floor

float groundNonMatchingPercent:

percentage of the floor section non matching ceil

int averageNeighborsPerCell:

average Neighbor number per cell

(see gridStep argument of ComputeVolume25D())

A first example with a flat floor:


cloud = cc.loadPointCloud(getSampleCloud(5.0))

report = cc.ReportInfoVol()
isOk = cc.ComputeVolume25D(report, ground=None, ceil=cloud, 
                           vertDim=2, gridStep=0.05, groundHeight=0, ceilHeight=0)

check the results:


if not isOk:
    raise RuntimeError
if not math.isclose(report.volume, 0.995, rel_tol=1e-03):
    raise RuntimeError
if not math.isclose(report.surface, 101.002, rel_tol=1e-03):
    raise RuntimeError
if not math.isclose(report.addedVolume, 11.726, rel_tol=1e-03):
    raise RuntimeError
if not math.isclose(report.removedVolume, 10.731, rel_tol=1e-03):
    raise RuntimeError
if not math.isclose(report.matchingPercent, 100., rel_tol=1e-03):
    raise RuntimeError
if not math.isclose(report.ceilNonMatchingPercent, 0., abs_tol=1e-03):
    raise RuntimeError
if not math.isclose(report.groundNonMatchingPercent, 0., abs_tol=1e-03):
    raise RuntimeError
if not math.isclose(report.averageNeighborsPerCell, 8., abs_tol=1e-03):
    raise RuntimeError

A second example with a cloud floor, translated to obtain non matching parts:


cloud2 = cc.loadPointCloud(getSampleCloud(2.0))
cloud2.translate((1,2, -3)) # creates a translated floor, 
                            # with a non matching part with the ceil 

report = cc.ReportInfoVol()
isOk = cc.ComputeVolume25D(report, ground=cloud2, ceil=cloud,
                           vertDim=2, gridStep=0.05, groundHeight=0, ceilHeight=0) 

check the results:


if not isOk:
    raise RuntimeError
if not math.isclose(report.volume, 216.609, rel_tol=1e-03):
    raise RuntimeError
if not math.isclose(report.surface, 72.852, rel_tol=1e-03):
    raise RuntimeError
if not math.isclose(report.addedVolume, 216.609, rel_tol=1e-03):
    raise RuntimeError
if not math.isclose(report.removedVolume, 0., rel_tol=1e-03):
    raise RuntimeError
if not math.isclose(report.matchingPercent, 56.4, rel_tol=1e-03):
    raise RuntimeError
if not math.isclose(report.ceilNonMatchingPercent, 21.8, rel_tol=1e-03):
    raise RuntimeError
if not math.isclose(report.groundNonMatchingPercent, 21.8, rel_tol=1e-03):
    raise RuntimeError
if not math.isclose(report.averageNeighborsPerCell, 7.93, rel_tol=1e-03):
    raise RuntimeError

The above code snippets are from test023.py.

6.11. Cloud rasterization

The main purpose of this tool is to ‘rasterize’ a point cloud (i.e. convert it to a 2.5D grid) and then export it as a new cloud or a raster image (GeoTiff) for instance. Concepts are introduced in the CloudCompare wiki Rasterize

Three functions are available for rasterization.

All functions have a lot of parameters, to produce or not GeoTiff files, with or without scalar fields or colors, to fill or not the empty cells, to export statistic scalar fields…

The gridStep parameter is mandatory in all rasterize methods.

The code snippets below are from test025.py. See the clouds and meshes generated by this test to understand the different options.

We start from a wave cloud created with a random generator, to obtain points with random (x,y) coordinates (not on a grid).


# --- generate a set of coords and a scalar field

npts = 1000000
h = 2.
x = np.float32(-5. + 10.*np.random.random((npts)))
y = np.float32(-5. + 10.*np.random.random((npts)))
z = np.float32(np.sin(h * np.sqrt(x**2 + y**2)) / np.sqrt(x**2 + y**2))
coords = np.column_stack((x,y,z))

# --- create the pointCloud

cloud = cc.ccPointCloud("wave_%d"%h)
cloud.coordsFromNPArray_copy(coords)

# --- add the scalar field

res = cloud.exportCoordToSF(False, False, True)

For a simple rasterisation to cloud:


rcloud = cc.RasterizeToCloud(cloud, 0.01)

With the default options, empty cells stays empty.

Below, an example where the empty cells are filled with a forced constant height, and with a GeoTiff file With raster altitude.


rcloud1 = cc.RasterizeToCloud(cloud,
                          gridStep=0.01, 
                          outputRasterZ = True,
                          pathToImages = dataDir,
                          emptyCellFillStrategy = cc.EmptyCellFillOption.FILL_CUSTOM_HEIGHT,
                          customHeight = 1.,
                          export_perCellCount = True)

Below, an example where the empty cells are filled with an iterpolated height, and with a GeoTiff file With raster altitude and scalar fields.


rcloud2 = cc.RasterizeToCloud(cloud,
                          gridStep=0.01, 
                          outputRasterZ = True,
                          outputRasterSFs = True,
                          pathToImages = dataDir,
                          emptyCellFillStrategy = cc.EmptyCellFillOption.INTERPOLATE,
                          export_perCellCount = True,
                          export_perCellAvgHeight = True)

Next, we generate only the GeoTiff file With raster altitude and scalar fields.


cloud.setName("wave_g2") # to distinguish the GeoTiff file from the previous one
cc.RasterizeGeoTiffOnly(cloud,
                        gridStep=0.01, 
                        outputRasterZ = True,
                        outputRasterSFs = True,
                        pathToImages = dataDir,
                        emptyCellFillStrategy = cc.EmptyCellFillOption.FILL_MAXIMUM_HEIGHT,
                        customHeight = 0.,
                        export_perCellCount = True)

Finally, a simple rasterization to mesh:


rmesh = cc.RasterizeToMesh(cloud, 0.03)

6.12. Interpolate scalar fields from one cloud to another

The concepts are introduced in the CloudCompare wiki interpolate scalar fields

In the following example, 3 scalar fields from a cloud are selected for interpolation on another cloud. the interpolatorParameters class is used to define the method, algorithm and values required for the interpolation. The cloudComPy.interpolateScalarFieldsFrom() function computes the interpolation.

Here, the distance between the clouds is locally greater than the radius, so the interpolated scalar fields are filed with NaN in some places.


cloud1 = cc.loadPointCloud(getSampleCloud(5.0))
cloud2 = cc.loadPointCloud(getSampleCloud(1.0))
cloud1.exportCoordToSF(True, True, True)
dic = cloud1.getScalarFieldDic()

params = cc.interpolatorParameters()
params.method = cc.INTERPOL_METHOD.RADIUS
params.algos = cc.INTERPOL_ALGO.NORMAL_DIST
params.radius = 0.15
params.sigma = 0.06

sfIndexes = [dic['Coord. X'], dic['Coord. Y'], dic['Coord. Z']]

ret = cc.interpolateScalarFieldsFrom(cloud2, cloud1, sfIndexes, params)

The above code snippet is from test044.py.

6.13. Finding an optimal bounding box

minimalBoundingBox is a pure Python plugin built with CloudComPy.

This tool is provided as an example of CloudComPy pure Python extension.

The findRotation() function finds a minimal bounding box for a cloud, not oriented along the main axes, and the corresponding rotation. It is not proven that the solution meets to the minimum, and the performance and robustness can certainly be improved.

We have to import minimalBoundingBox


from cloudComPy.minimalBoundingBox import findRotation

We create an ellipsoid, and apply a rotation on it, so that its bounding box is far from optimal.


sphere = cc.ccSphere(1.0)
cloud = sphere.samplePoints(False, 100000)
cloud.scale(1.0, 3.0, 9.0)

cloudBeforeRot = cloud.cloneThis()
cloudBeforeRot.setName("cloudBeforeRot")

transform1 = cc.ccGLMatrix()
transform1.initFromParameters(1., (1.5, 2.5, 2.), (0,0,0))
cloud.applyRigidTransformation(transform1)
cloud.setName("rotated object")

The findRotation() function returns several objects.

  • the optimized bounding box

  • the rotation (double precision) to go the object on axes to the original object

  • the polyline representing the optimal bounding box, aligned with the original object

  • the cloud associated to this polyline


boundingBox, rotinv, poly, clbbox = findRotation(cloud)

To check the results, we apply the rotation (inverse) to a copy of the cloud and polyline.


rotation = (cc.ccGLMatrix.fromDouble(rotinv)).inverse()
axisObj = cloud.cloneThis()
axisObj.applyRigidTransformation(rotation)
axisObj.setName("object on axes")

clb = clbbox.cloneThis()
clb.applyRigidTransformation(rotation)
axisPoly = cc.ccPolyline(clb)
axisPoly.addChild(clb)
axisPoly.addPointIndex(0, clb.size())
axisPoly.setClosed(True)

res = cc.SaveEntities([cloudBeforeRot, cloud, axisObj, poly, axisPoly], os.path.join(dataDir, "optimalBoundingBox.bin"))

The above code snippets are from test027.py.

6.14. Cloud comparison with plugin M3C2

The M3C2 plugin is introduced in the CloudCompare Plugins wiki - M3C2. The plugin computes signed distances between two clouds.

The computeM3C2() function relies on a large number of parameters grouped in a file. The CloudCompare GUI offers a ‘guess parameters’ option and a way to save the parameters to a file. The M3C2guessParamsToFile() function do the same job.

You have to import the The M3C2 plugin. The M3C2guessParamsToFile() function requires the two clouds, the params file, and a boolean to get a fast guess or not. The computeM3C2() function requires only the two clouds and the params file, and produces a cloud with new scalar field ‘M3C2 distance’:


if cc.isPluginM3C2():
    import cloudComPy.M3C2
    cloud = cc.loadPointCloud(getSampleCloud(5.0))
    cloud1 = cc.loadPointCloud(getSampleCloud(1.0))
    paramFilename = os.path.join(dataDir, "m3c2_params.txt")
    ret = cc.M3C2.M3C2guessParamsToFile([cloud,cloud1], paramFilename, True)
    cloud2 = cc.M3C2.computeM3C2([cloud,cloud1], paramFilename)

The params file has a simple syntax and can be created from scratch or edited. The M3C2 params file generated below corresponds to the ‘guess params’ option of the GUI for the clouds of this example:


import multiprocessing

m3c2_params_dic={}
m3c2_params_dic["ExportDensityAtProjScale"] = "false"
m3c2_params_dic["ExportStdDevInfo"] = "false"
m3c2_params_dic["M3C2VER"] = 1
m3c2_params_dic["MaxThreadCount"] = multiprocessing.cpu_count()
m3c2_params_dic["MinPoints4Stat"] = 5
m3c2_params_dic["NormalMaxScale"] = 0.283607
m3c2_params_dic["NormalMinScale"] = 0.070902
m3c2_params_dic["NormalMode"] = 0
m3c2_params_dic["NormalPreferedOri"] = 4
m3c2_params_dic["NormalScale"] = 0.141803
m3c2_params_dic["NormalStep"] = 0.070902
m3c2_params_dic["NormalUseCorePoints"] = "false"
m3c2_params_dic["PM1Scale"] = 1
m3c2_params_dic["PM2Scale"] = 1
m3c2_params_dic["PositiveSearchOnly"] = "false"
m3c2_params_dic["ProjDestIndex"] = 1
m3c2_params_dic["RegistrationError"] = 0
m3c2_params_dic["RegistrationErrorEnabled"] = "false"
m3c2_params_dic["SearchDepth"] = 0.709017
m3c2_params_dic["SearchScale"] = 0.141803
m3c2_params_dic["SubsampleEnabled"] = "true"
m3c2_params_dic["SubsampleRadius"] = 0.070902
m3c2_params_dic["UseMedian"] = "false"
m3c2_params_dic["UseMinPoints4Stat"] = "false"
m3c2_params_dic["UseOriginalCloud"] = "false"
m3c2_params_dic["UsePrecisionMaps"] = "false"
m3c2_params_dic["UseSinglePass4Depth"] = "false"

paramFilename1 = os.path.join(dataDir, "m3c2_params1.txt")
with open(paramFilename1, 'w') as f:
    f.write("[General]\n")
    for k,v in m3c2_params_dic.items():
        f.write("%s=%s\n"%(k,v))

To use the precision maps, you have to provide the six corresponding scalar fields in list (3 components for the first cloud and 3 components for the second cloud), and also two scales in a list. When these lists are provided, the computation uses the precision maps, and supersedes the corresponding option in the params file.


    sfs = []
    dic = cloud.getScalarFieldDic()
    sfs.append(cloud.getScalarField(dic["ux"]))
    sfs.append(cloud.getScalarField(dic["uy"]))
    sfs.append(cloud.getScalarField(dic["uz"]))
    dic1 = cloud1.getScalarFieldDic()
    sfs.append(cloud1.getScalarField(dic1["ux"]))
    sfs.append(cloud1.getScalarField(dic1["uy"]))
    sfs.append(cloud1.getScalarField(dic1["uz"]))
    scales = [1., 1.]
    paramFilename = os.path.join(dataDir, "m3c2_params2.txt")
    ret = cc.M3C2.M3C2guessParamsToFile([cloud,cloud1], paramFilename, True)
    cloud2 = cc.M3C2.computeM3C2([cloud,cloud1], paramFilename, sfs, scales)

The above code snippets are from test030.py.

6.15. ShadeVIS (ambiant occlusion) with Plugin PCV

The PCV plugin is documented in the CloudCompare Plugins wiki - PCV. It calculates the illumination of a point cloud (or vertices of a mesh) as if the light was coming from a theoretical hemisphere or sphere around the object (the user can also provide its own set of light directions - see below).

You have to import the The PCV plugin.


if cc.isPluginPCV():
    import cloudComPy.PCV

We first try to illuminate the cloud with a dish as light source, oriented to create unrealistic sharp shadows. The cloud used as light source must have normals.


    tr1 = cc.ccGLMatrix()
    tr1.initFromParameters(0.3*math.pi, (0., 1., 0.), (0.0, 0.0, 0.0))
    dish = cc.ccDish(2.0, 0.5, 0.0, tr1)
    cln =dish.getAssociatedCloud()
    cc.computeNormals([cln])

The computeShadeVIS() function takes a list of objects (clouds or meshes) as first parameter, and will add a scalar field illuminance (PCV) to these clouds or meshes. The second parameter is the optional light source. There are several other optional parameters.


    cc.PCV.computeShadeVIS([cloud],cln)
    dic = cloud.getScalarFieldDic()
    cloud.renameScalarField(dic["Illuminance (PCV)"], "IlluminanceDish")

The Scalar field is renamed to keep it and compute the illuminance with the default hemispheric light source.


    cc.PCV.computeShadeVIS([cloud])

The above code snippets are from test032.py.

6.16. Compute Hidden Point Removal with plugin HPR

The HPR plugin is described in the CloudCompare Plugins wiki - HPR.

The algorithm filters out the points of a cloud that wouldn’t be seen (by the current 3D camera) if the cloud was a closed surface. Therefore it tries to remove the points that should be hidden in the current viewport configuration.

Without GUI, in CloudComPy, we have to explitely provide the coordinates of the viewpoint in computeHPR().


if cc.isPluginHPR():
    import cloudComPy.HPR
    
    cloud = cc.loadPointCloud(getSampleCloud(5.0))
    
    cloudCut = cc.HPR.computeHPR(cloud, (0.,-15., 25.))
    

The above code snippet is from test033.py.

6.17. Boolean operations on meshes with plugin MeshBoolean

The MeshBoolean plugin is described in the CloudCompare Plugins wiki - Mesh_Bolean.

This plugin can be used to perform Boolean operations on meshes (also called CSG = Constructive Solid Geometry).

In the example, we try an intersection of two primitives, a sphere and a cylinder,

with computeMeshBoolean().


tr1 = cc.ccGLMatrix()
tr1.initFromParameters(0.0, (0., 0., 0.), (1.0, 0.0, 0.0))
sphere = cc.ccSphere(2, tr1, "aSphere")

cylinder = cc.ccCylinder(2.0, 5.0)

if cc.isPluginMeshBoolean():
    import cloudComPy.MeshBoolean
    mesh = cc.MeshBoolean.computeMeshBoolean(sphere, cylinder, 
                                             cc.MeshBoolean.CSG_OPERATION.INTERSECT)
    

The above code snippet is from test034.py.

6.18. Finding primitives on a cloud with plugin RANSAC-SD

The RANSAC-SD plugin is documented in the CloudCompare Plugins wiki - RANSAC Shape Detection.

The RANSAC-SD plugin (RANdom SAmple Consensus) is used to detect simple shapes (planes, spheres, cylinders, cones, torus) in a cloud.

You have to import the The RANSAC_SD plugin.

The computeRANSAC_SD() function requires a set of parameters: RansacParams.

The optimizeForCloud() method helps to get a correct set of parameters. By default, the algorithm searchs only planes, spheres and cylinders, you can enable or disable a primitive type with setPrimEnabled.

For the example we regroup 4 clouds (samples from 3 spheres and 1 cylinder) in a single cloud.


tr0 = cc.ccGLMatrix()
tr0.initFromParameters(math.pi/3., (0., 1., 0.), (3.0, 0.0, 4.0))
cylinder = cc.ccCylinder(0.5, 2.0, tr0)
c0 = cylinder.samplePoints(True, 1000)

tr1 = cc.ccGLMatrix()
tr1.initFromParameters(0.0, (0., 0., 0.), (-2.0, 5.0, 1.0))
sphere1 = cc.ccSphere(1.5, tr1)
c1 = sphere1.samplePoints(True, 1000)

tr2 = cc.ccGLMatrix()
tr2.initFromParameters(0.0, (0., 0., 0.), (6.0, -3.0, -2.0))
sphere2 = cc.ccSphere(2.0, tr2)
c2 = sphere2.samplePoints(True, 1000)

tr3 = cc.ccGLMatrix()
tr3.initFromParameters(0.0, (0., 0., 0.), (0.0, 1.0, 2.0))
sphere3 = cc.ccSphere(1.0, tr3)
c3 = sphere3.samplePoints(True, 1000)

cloud = c0.cloneThis()
cloud.fuse(c1)
cloud.fuse(c2)
cloud.fuse(c3)

The computeRANSAC_SD() function returns a list of meshes and clouds for the shapes found.


if cc.isPluginRANSAC_SD():
    import cloudComPy.RANSAC_SD
    params = cc.RANSAC_SD.RansacParams()
    params.optimizeForCloud(cloud)
    print(params.epsilon, params.bitmapEpsilon)
    meshes, clouds = cc.RANSAC_SD.computeRANSAC_SD(cloud, params)
    

The above code snippets are from test035.py.

6.19. compute Cloth Simulation Filter on a cloud with CSF plugin

The concepts are presented in the CloudCompare wiki - CSF (plugin).

You have to import the The CSF plugin.

The computeCSF() function takes a ccPointCloud in input and several optional parameters. The function produce two clouds: “ground points” and “off-ground points”.

A first example with default parameters:


if cc.isPluginCSF():
    import cloudComPy.CSF
    cloud = cc.loadPointCloud(getSampleCloud(5.0))
    clouds = cc.CSF.computeCSF(cloud)
    

And a second example:


    clouds2 = cc.CSF.computeCSF(cloud, csfRigidness=1, clothResolution=1.0, classThreshold=0.3)
    

The above code snippets are from test043.py.

6.20. Sclices and contours

The concepts are presented in the CloudCompare wiki - Cross Section.

Given a bounding box used as a cut tool, and a list of entities to sclice (clouds and meshes), the ExtractSlicesAndContours() function allows to:

  • segment the entities (slice)

  • extract the envelope of all the points inside the slice

  • extract one or several contours of the points inside each slice

  • repeat the segmentation or extraction processes above along one or several dimensions (to extract multiple ‘slices’ for instances)

For the example, we create a list of objects to slice, one cloud and two meshes.


tr1 = cc.ccGLMatrix()
tr1.initFromParameters(0.0, (0., 0., 0.), (-2.0, 0.0, 0.0))
sphere1 = cc.ccSphere(3.0, tr1)
c1 = sphere1.samplePoints(True, 1000)

tr2 = cc.ccGLMatrix()
tr2.initFromParameters(0.0, (0., 0., 0.), (+2.0, 0.0, 0.0))
sphere2 = cc.ccSphere(3.0, tr2)
c2 = sphere2.samplePoints(True, 1000)

cloud = c1.cloneThis()
cloud.fuse(c2)

toslice = [cloud, sphere1, sphere2] # one cloud, two meshes

We define a bounding box as a cutting tool:


bbox = cc.ccBBox((-5.0, -5.0, 0.), (5., 5., 0.5), True)

First, generate a single slice: the function returns a list containing a single slice


res=cc.ExtractSlicesAndContours(entities=toslice, bbox=bbox)

then, repeat the slices along Z with a gap: the function returns a list containing several slices


res=cc.ExtractSlicesAndContours(entities=toslice, bbox=bbox, singleSliceMode=False,
                                gap=0.5, generateRandomColors=True)

next, the same as above, plus envelopes: the function returns a tuple(list of slices, list of envelopes)


res=cc.ExtractSlicesAndContours(entities=toslice, bbox=bbox, singleSliceMode=False, 
                                gap=0.5, generateRandomColors=True,
                                extractEnvelopes=True, maxEdgeLength=0.1, envelopeType=2)

after, the same as above, with a rotation of the bounding box: the function returns a tuple(list of slices, list of envelopes)


tr0 = cc.ccGLMatrix()
tr0.initFromParameters(math.pi/12., (0., 1., 0.), (0.0, 0.0, 0.0))
res=cc.ExtractSlicesAndContours(entities=toslice, bbox=bbox, bboxTrans=tr0, 
                                singleSliceMode=False, gap=0.5, generateRandomColors=True,
                                extractEnvelopes=True, maxEdgeLength=0.1, envelopeType=2)

finally, the same as above, plus contours: the function returns a tuple(list of slices, list of envelopes, list of contours)


res=cc.ExtractSlicesAndContours(entities=toslice, bbox=bbox, bboxTrans=tr0, 
                                singleSliceMode=False, gap=0.5, generateRandomColors=True,
                                extractEnvelopes=True, maxEdgeLength=0.1, envelopeType=2,
                                extractLevelSet=True, levelSetGridStep=0.05, 
                                levelSetMinVertCount=100)

The above code snippets are from test036.py.

6.21. ExtractConnectedComponents

The tool used here is described in CloudCompare wiki - Label Connected Components.

This tool segments the selected cloud(s) in smaller parts separated by a minimum distance. Each part is a connected component (i.e. a set of ‘connected’ points).

For the test, we create a cloud with clear gaps, using a sclice operation:


cloud = cc.loadPointCloud(getSampleCloud(5.0))
bbox = cc.ccBBox((-5.0, -5.0, 0.), (5., 5., 1.), True)
res=cc.ExtractSlicesAndContours(entities=[cloud], bbox=bbox)
clouds = res[0] # result = [one cloud]

The ExtractConnectedComponents() function use the octree level to define the size of the gap between the components.


res2 = cc.ExtractConnectedComponents(clouds=clouds, octreeLevel=6, randomColors=True)
components = res2[1]
for comp in components:
    comp.showColors(True)

The function returns a tuple (number of clouds processed, list of components)

The above code snippets are from test037.py.

6.22. Avoid memory leaks in an iterative process

In a pure Python script, the garbage collector takes care of the objects that are no longer referenced. With CloudComPy, there is no shared reference count between a CloudComPy Python object and the CloudCompare C++ object wrapped. In other words, if a Python object is deleted, the wrapped C++ object remains. Conversely, if a CloudCompare C++ object is deleted, the corresponding Python object is not automatically deleted, but is no longer usable.

To avoid memory leaks, you should delete explicitely the objects that are no longer used with the cloudComPy.deleteEntity() function.

For instance, we create some data:


process = psutil.Process(os.getpid())
mem_start = process.memory_full_info().rss
cloud = cc.loadPointCloud(getSampleCloud(5.0))
bbox = cc.ccBBox((-5.0, -5.0, 0.), (5., 5., 1.), True)
res=cc.ExtractSlicesAndContours(entities=[cloud], bbox=bbox)
clouds = res[0]
print('input data memory usage:', process.memory_full_info().rss - mem_start)

Then, we iterate on ExtractConnectedComponents()


for i in range(10):
    mem_start = process.memory_full_info().rss
    res2 = cc.ExtractConnectedComponents(clouds=clouds, octreeLevel=6, randomColors=True)
    components = res2[1]
    for comp in components:
        comp.showColors(True)
    

I we do not delete explicitely the extracted components, they are kept forever in CloudCompare. To avoid the memory leak:


    for comp in components:
        cc.deleteEntity(comp)
    # be sure to have no more Python variable referencing the deleted items
    components = None
    res2 = None
    addedMemory = process.memory_full_info().rss - mem_start
    print(f'iteration {i}, ExtractConnectedComponents added memory: ', addedMemory)
    

The above code snippets are from test042.py.