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.setName("poly1")
poly1.addChild(cloud1)
poly1.addPointIndex(0, cloud1.size())
poly1.setClosed(True)
The above code snippets are from test026.py
.
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.
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
.
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
.
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.
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 named Approx. distances
,
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 method cloudComPy.DistanceComputationTools.computeCloud2MeshDistances()
generates a new scalar field in the cloud,
with accurate distances to the mesh. It’s name begins with C2M absolute distances
or C2M signed distances
following the boolean parameter signedDistances
value. If the parameter maxSearchDist
is used, the scalar field name is suffixed
with [val]
where val
is the maxSearchDist
value.
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
.
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 method cloudComPy.DistanceComputationTools.computeCloud2CloudDistances()
generates a new scalar field in the cloud,
with accurate distances to the cloud. It’s name begins with C2C absolute distances
.
If the parameter maxSearchDist
is used, the scalar field name is suffixed
with [val]
where val
is the maxSearchDist
value.
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
.
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
.
4. Cloud registration
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.
Note: the ‘aligned’ member of the result structure is the same object as the original cloud to align: clone it before applying
the transformation if you want to keep the original.
res=cc.ICP(data=cloud2, model=cloud1, finalOverlapRatio=0.1)
tr2 = res.transMat
print(tr2.toString())
cloud3 = res.aligned.cloneThis() # note: res.aligned is clone2!
cloud3.applyRigidTransformation(tr2)
cloud3.setName("cloud2_transformed_afterICP")
The above code snippets are from test010.py
.
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
.
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. Fitting ccPlane
The method Fit()
allows to adjust a plane primitive on a cloud.
The method getEquation()
gives the 4 coefficients of the plane ([a, b, c, d] as ax+by+cz=d):
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.
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.
8. Cloud sampling, noise filter, filter by scalar field
The class cloudComPy.CloudSamplingTools
offers several cloud sampling filters.
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() < 7470 or refCloud.size() > 7570:
raise RuntimeError
if res != 0:
raise RuntimeError
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()
refCloud = cc.CloudSamplingTools.resampleCloudSpatially(cloud, 0.05)
(spatialCloud, res) = cloud.partialClone(refCloud)
spatialCloud.setName("spatialCloud")
if refCloud.size() != 55465:
raise RuntimeError
if res != 0:
raise RuntimeError
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
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
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
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
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
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
.
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
.
9. Generate histograms
To draw plots or to generate csv files from scalar field data, we use numpy, matplotlib and csv.
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
.
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
.
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 |
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,
projectionType=cc.PROJ_MINIMUM_VALUE,
groundEmptyCellFillStrategy=cc.INTERPOLATE_DELAUNAY,
groundMaxEdgeLength=0,
ceilEmptyCellFillStrategy=cc.INTERPOLATE_DELAUNAY,
ceilMaxEdgeLength=0)
check the results:
if not isOk:
raise RuntimeError
if not math.isclose(report.volume, 299.07, rel_tol=1e-03):
raise RuntimeError
if not math.isclose(report.surface, 133.152, rel_tol=1e-03):
raise RuntimeError
if not math.isclose(report.addedVolume, 301.305, rel_tol=1e-03):
raise RuntimeError
if not math.isclose(report.removedVolume, 2.237, 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, rel_tol=1e-03):
raise RuntimeError
if not math.isclose(report.groundNonMatchingPercent, 0., rel_tol=1e-03):
raise RuntimeError
if not math.isclose(report.averageNeighborsPerCell, 8.0, rel_tol=1e-03):
raise RuntimeError
The above code snippets are from test023.py
.
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)
rcloud.setName("raster_0")
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)
rcloud1.setName("raster_1")
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_DELAUNAY,
export_perCellCount = True,
export_perCellAvgHeight = True,
export_perCellMedian = True,
export_perCellPercentile = True)
rcloud2.setName("raster_2")
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)
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
.
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
.
14. Working with quaternions for rotations
Quaternions are a convenient way to express rotations: see for instance Quaternions and spatial rotation. With numpy, the quaternion package provides quaternion objects very useful to work with rotations.
A quaternion is based on an array of four floats (q1, q2, q3, q4). A normalized quaternion such as |q| = sqrt(q1*q1 + q2*q2 + q3*q3 + q4*q4) = 1 can also be expressed as:
q = {cos(θ/2), ux*sin(θ/2), uy*sin(θ/2), uz*sin(θ/2)}
where (ux, uy, uz) represent the components of a unit vector u, and θ is an angle. The quaternion defines a rotation around the vector u, of angle θ.
To use the quaternion package, it must be imported:
import numpy as np
import quaternion
The following exemple function gives a quaternion with a vector and angle as arguments.
def normalized_quaternion(theta, ux, uy, uz):
norm = math.sqrt(ux*ux + uy*uy +uz*uz)
q1 = math.cos(theta/2.)
q2 = ux*math.sin(theta/2.)/norm
q3 = uy*math.sin(theta/2.)/norm
q4 = uz*math.sin(theta/2.)/norm
return np.quaternion(q1, q2, q3 ,q4)
A rotation in a cloudComPy transformation tr1 can be expressed as a quaternion q1
with the toQuaternion()
method
from ccGLMatrixd
or ccGLMatrix
.
The quaternion can also be expressed as a rotation matrix rot1.
tr1 = cc.ccGLMatrixd()
tr1.initFromParameters(math.pi/3., (1., 1., 1.), (0,0,0)) # 60°
q1 = quaternion.as_quat_array(tr1.toQuaternion())
rot1 = quaternion.as_rotation_matrix(q1)
A cloudCompy transformation can be built from a quaternion and an optional translation vector.
with the FromQuaternionAndTranslation()
method
from ccGLMatrixd
or ccGLMatrix
.
trb = cc.ccGLMatrixd.FromQuaternionAndTranslation(qb.components)
datb = trb.getParameters1()
angle = datb.alpha_rad*180/math.pi
axis = datb.axis3D
The above code snippets are from test058.py
.
15. 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
.
16. 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
.
18. Boolean operations on meshes with plugin MeshBoolean or plugin Cork
The MeshBoolean plugin is described in the CloudCompare Plugins wiki - Mesh_Bolean.
The Cork plugin is described in the CloudCompare Plugins wiki - Cork.
These plugins can be used to perform Boolean operations on meshes (also called CSG = Constructive Solid Geometry). The MeshBoolean plugin is slower but supposedly more robust.
- In the first 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 not cc.isPluginMeshBoolean():
print("Test skipped")
sys.exit()
import cloudComPy.MeshBoolean
mesh = cc.MeshBoolean.computeMeshBoolean(sphere, cylinder,
cc.MeshBoolean.CSG_OPERATION.INTERSECT)
The above code snippet is from test034.py
.
- In the second example, we try an difference of two primitives, a cylinder and a sphere,
with
compute()
.
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 not cc.isPluginCork():
print("Test skipped")
sys.exit()
import cloudComPy.Cork
mesh = cc.Cork.Cork.compute(cylinder, sphere, cc.Cork.DIFF)
The above code snippet is from test056.py
.
19. 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)
For each type of primitive (plane, sphere, cylinder,..) it is possible to get the parameters of the equation. For instance, a plane is identified with
- ::
mesh.isA((cc.CC_TYPES.PLANE)
For the plane primitive, the method getEquation()
returns the 4 coefficients of the plane equation: [a, b, c, d] as ax+by+cz=d
The above code snippets are from test035.py
.
20. 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”, plus optionally, the cloth mesh.
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, computeMesh=True)
The above code snippets are from test043.py
.
21. Classify a point cloud with Canupo plugin and a trained classifier
The concepts are presented in the CloudCompare wiki - CANUPO (plugin).
You have to import the The Canupo
plugin.
For the test, an example of cloud and a classifier are provided:
# example data available here: https://nicolas.brodu.net/common/recherche/canupo/benchmark.tar.gz
# classifier are available here : http://www.cloudcompare.org/forum/viewtopic.php?f=17&t=808&start=90#p11588
if not os.path.isfile(os.path.join(dataExtDir,"recombi_10.txt")):
if not os.path.exists(dataExtDir):
os.makedirs(dataExtDir)
url = "https://www.simulation.openfields.fr/phocadownload/recombi_10.txt"
with open(os.path.join(dataExtDir,"recombi_10.txt"), 'wb') as f:
for line in urllib.request.urlopen(url):
f.write(line)
if not os.path.isfile(os.path.join(dataExtDir,"vegetTidal.prm")):
if not os.path.exists(dataExtDir):
os.makedirs(dataExtDir)
url = "https://www.simulation.openfields.fr/phocadownload/vegetTidal.prm"
r = requests.get(url)
with open(os.path.join(dataExtDir,"vegetTidal.prm"), 'wb') as f:
f.write(r.content)
cloud=cc.loadPointCloud(os.path.join(dataExtDir,"recombi_10.txt"))
The Classify()
function takes in input a ccPointCloud, a classifier file and several optional parameters.
The function produce two scalar fields in the cloud: ‘CANUPO.class’ and ‘CANUPO.confidence’.
if cc.isPluginCanupo():
import cloudComPy.Canupo
res = cc.Canupo.Classify(cloud, os.path.join(dataExtDir,"vegetTidal.prm"))
dic = cloud.getScalarFieldDic()
sf1 = cloud.getScalarField(dic['CANUPO.class'])
if sf1 is None:
raise RuntimeError
sf2 = cloud.getScalarField(dic['CANUPO.confidence'])
if sf2 is None:
raise RuntimeError
res = cc.SaveEntities([cloud], os.path.join(dataDir, "cloudCanupo.bin"))
The above code snippets are from test046.py
.
22. Compute distance between a cloud and a surface of revolution, with SRA plugin
The concepts are presented in the CloudCompare wiki - SRA (plugin).
You have to import the The SRA
plugin.
The surface of revolution is defined by a profile.
See CloudCompare wiki - SRA (plugin)
for a descritpion of the text file format for the profile.
This profile can be loaded with the loadProfile()
function.
The doComputeRadialDists()
method of qSRA
class
computes the radial distances between the cloud and the surface of revolution
and store the result as a scalar field named ‘Radial distance’.
An example:
if cc.isPluginSRA():
import cloudComPy.SRA
poly = cc.SRA.loadProfile(os.path.join(dataDir, "profil.txt"), 2, False)
sra = cc.SRA.qSRA()
res=sra.doComputeRadialDists(cl, poly)
dic = cl.getScalarFieldDic()
sf = cl.getScalarField(dic['Radial distance'])
Distance map can be exported as a cloud or a mesh with the exportMapAsCloud()
and exportMapAsMesh()
functions:
clmap = cc.SRA.exportMapAsCloud(cl, poly, sf, 0.5, 0.01, 0., 10., baseRadius=2)
meshmap = cc.SRA.exportMapAsMesh(cl, poly, sf, 0.5, 0.01, 0., 10., colScale=cc.SRA.DEFAULT_SCALES.YELLOW_BROWN)
The above code snippets are from test045.py
.
23. Surface Mesh reconstruction, with Poisson Reconstruction plugin
The concepts are presented in the CloudCompare wiki - Poisson Surface Reconstruction (plugin).
You have to import the The PoissonRecon
plugin.
To use this plugin, you have to provide a cloud with normals.
The PoissonReconstruction()
static method of PR
class
computes the triangular mesh on the given cloud.
mesh = cc.PoissonRecon.PR.PoissonReconstruction(pc=cloud, threads=multiprocessing.cpu_count(), density=True)
The above code snippet is from test053.py
.
24. 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
.
25. Extract Sections (generate cloud slices and profiles)
The concepts are presented in the CloudCompare wiki - Extract Sections.
This set of tools allows the user to import poylines on top of a point cloud so as to extract sections and profiles.
For the test, we use a cloud representing the altimetry around a river, and a 2D polyline giving the river axis.
The generateOrthoSections()
method automatically generate orthogonal sections with a given polyline at regular intervals.
cloud=cc.loadPointCloud(os.path.join(dataExtDir, "garonne_L93.xyz"))
shift=cloud.getGlobalShift()
poly=cc.loadPolyline(os.path.join(dataExtDir, "ligne_eau.poly"), cc.CC_SHIFT_MODE.XYZ, 0, shift[0], shift[1], shift[2] )
orthoPolys=poly.generateOrthoSections(1000, 600, 2)
With clouds and polylines, the extractPointsAlongSections()
function allows to extract either cloud slices along the polylines,
or to define vertical profiles (polylines) built on the clouds.
The unfoldPointsAlongPolylines()
function unfolds a polyline and the clouds points near the polyline (with a given thickness).
The result is a new cloud.
sections = cc.extractPointsAlongSections([cloud], [poly])
orthoSections = cc.extractPointsAlongSections([cloud], orthoPolys)
cloudSections = cc.extractPointsAlongSections([cloud], [poly], extractSectionsAsClouds=True, extractSectionsAsEnvelopes=False)
unfoldedClouds = cc.unfoldPointsAlongPolylines([cloud], [poly], 1000)
The above code snippets are from test052.py
.
26. Extract or Label Connected Components
The tools used here are described in CloudCompare wiki - Label Connected Components.
There are two tools relative to connected components: ExtractConnectedComponents()
and LabelConnectedComponents()
.
ExtractConnectedComponents()
Creates separate clouds and
LabelConnectedComponents()
creates a scalar field.
These tools segment the 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]
Here, we take the ExtractConnectedComponents()
function.
It uses 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, list of residual components). The “residual component” corresponds to all the remaining nodes of an input cloud: see for instance:
cloud = cc.loadPointCloud(os.path.join(dataExtDir, "pumpARowColumnIndexNoInvalidPoints.e57"))
res = cc.ExtractConnectedComponents(clouds=[cloud], randomColors=True)
cloud01 = cc.MergeEntities(res[1], createSFcloudIndex=True)
cloud01.setName("pump_Extract_Components")
cloud02 = cc.MergeEntities(res[2], createSFcloudIndex=True)
cloud02.setName("pump_residual_Components")
res2 = res[1] + res[2] # connected components plus regrouped residual components
cloud2 = cc.MergeEntities(res2, createSFcloudIndex=True)
cloud2.setName("pump_extract_Residual_Components")
cloud3 = cc.MergeEntities(res2, deleteOriginalClouds=True)
cloud3.setName("pump_extract_Residual_Components")
res = None
The above code snippets are from test037.py
and test048.py
.
27. Merge several entities (clouds or meshes)
The MergeEntities()
function allows to merge several clouds, or several meshes.
Merging a list of clouds:
cloud = cc.loadPointCloud(os.path.join(dataExtDir, "pumpARowColumnIndexNoInvalidPoints.e57"))
res = cc.ExtractConnectedComponents(clouds=[cloud], randomColors=True)
cloud01 = cc.MergeEntities(res[1], createSFcloudIndex=True)
cloud01.setName("pump_Extract_Components")
cloud02 = cc.MergeEntities(res[2], createSFcloudIndex=True)
cloud02.setName("pump_residual_Components")
res2 = res[1] + res[2] # connected components plus regrouped residual components
cloud2 = cc.MergeEntities(res2, createSFcloudIndex=True)
cloud2.setName("pump_extract_Residual_Components")
cloud3 = cc.MergeEntities(res2, deleteOriginalClouds=True)
cloud3.setName("pump_extract_Residual_Components")
res = None
Merging a list of meshes:
mesh1 = cc.MergeEntities([box, cone, cylinder])
mesh2 = cc.MergeEntities([box, cone, cylinder], createSubMeshes=True)
The above code snippets are from test048.py
.
28. 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
.
29. Render a 3D scene to an image file
In order to create a 3D view of one or more entities and write a 2D file image (png, jpeg,…),
the entities can be added to the scene with the cloudComPy.addToRenderScene()
function.
The cloudComPy.render()
function renders and writes the image.
cloud1 = cc.loadPointCloud(os.path.join(dataDir, "boule.bin"))
cloud1.setCurrentScalarField(0)
cloud1.setCurrentDisplayedScalarField(0)
cc.addToRenderScene(cloud1)
struct = cc.importFile(os.path.join(dataExtDir, "manitou.e57"))
for cloud in struct[1]:
cc.computeNormals([cloud])
cloud.showNormals(True)
cc.addToRenderScene(cloud)
cc.setGlobalZoom()
cc.render(os.path.join(dataDir, "rendera.png"), 2000,1500)
The cloudComPy.removeFromRenderScene()
function allows to remove an entity from the scene, but it deletes the C++ object.
Thus, the Python object becomes invalid.
cc.removeFromRenderScene(cloud1)
cloud1=None
Several functions allow to define a standard point of view for the scene:
Different modes of perspective are available:
To define a zoom on the whole scene or on selected entities, use the
cloudComPy.setGlobalZoom()
andcloudComPy.zoomOnSelectedEntity()
functions.
The cloudComPy.setCustomView()
and cloudComPy.setCameraPos()
functions are used to define a custom point of view.
The size of cloud points can be changed with setPointSize()
method on a cloud.
Several functions and methods can be used to set the background color, display the scalar bar…
cloud = cc.loadPointCloud(getSampleCloud(5.0))
cloud.exportCoordToSF(False, False, True)
cloud.setCurrentScalarField(0)
cloud.setCurrentDisplayedScalarField(0)
cloud.showSFColorsScale(True)
cc.addToRenderScene(cloud)
cc.setOrthoView()
cc.setGlobalZoom()
cc.setBackgroundColor(False, 255, 255, 255)
cc.setTextDefaultCol(0, 0, 0)
cc.setColorScaleShowHistogram(True)
cc.render(os.path.join(dataDir, "render0.png"), 2000,1500)
cc.setCameraPos((0., 0., 20.))
cc.setCustomView((0., 1., 0.), (0., 0., 1.))
cc.render(os.path.join(dataDir, "renderRef.png"), 2000,1500)
cloud.setPointSize(4)
cc.render(os.path.join(dataDir, "renderRef4.png"), 2000,1500)
cloud.setPointSize(0)
cc.setCenteredPerspectiveView()
cc.render(os.path.join(dataDir, "renderRefCenterPers.png"), 2000,1500)
dist=20.
alphaRad=math.pi/20.
cc.setCenteredPerspectiveView()
for i in range(15):
cc.setCustomView((0., math.cos(i*alphaRad), -math.sin(i*alphaRad)), (0., math.sin(i*alphaRad), math.cos(i*alphaRad)))
cc.setGlobalZoom()
cc.render(os.path.join(dataDir, "renderangle_%d.png"%i), 2000,1500, False)
By default, the cloudComPy.render()
function returns immediately after rendering the image.
It is possible to select an interactive mode by setting the last optional parameter (bool isInteractive) to True.
In this mode, the Python script is suspended, it is possible to adjust the 3D view. When the vue is correct,
you have to select the action “resume Python script” in the “options” menu to save the render and resume the Python script execution.
The above code snippets are from test051.py
.
30. choosing a specific color scale for a scalar field rendering
The default color scale for a scalar field is not always convenient for your rendering.
You can choose a color scale from the predefined list of color scales. In the following example, we select the VIRIDIS color scale for one scalar field:
cloud = cc.loadPointCloud(getSampleCloud(5.0))
res = cloud.exportCoordToSF(True, True, True)
sfx=cloud.getScalarField(0)
sfy=cloud.getScalarField(1)
sfz=cloud.getScalarField(2)
colorScalesManager = cc.ccColorScalesManager.GetUniqueInstance()
scale=colorScalesManager.getDefaultScale(cc.DEFAULT_SCALES.VIRIDIS)
sfx.setColorScale(scale)
It is also possible to build a specific color scale by defining several steps (relative position, associated color), the relative positions 0 and 1 corresponding to the start and end of the scale must always be defined. You can add several intermediate steps.
newScale=cc.ccColorScale.Create("myNewScale")
newScale.insert(0, cc.QColor("blue"))
newScale.insert(1, cc.QColor("red"))
newScale.insert(0.5, cc.QColor("green"))
sfy.setColorScale(newScale)
The above code snippets are from test057.py
.