CloudCompare / CloudComPy

Python wrapper for CloudCompare
Other
283 stars 40 forks source link

Subsampling a pointcloud #48

Open s-du opened 2 years ago

s-du commented 2 years ago

I am struggling to subsample a point cloud using CloudComPy. My problem seems to be this 'GenericIndexedCloudPersist' class that is needed as input. Can anyone indicates me on how to convert a ccPointCloud into a genericIndexedCloudPersist?

Many thanks!

prascle commented 2 years ago

Hello, genericIndexedCloudPersist is a base class for ccPointCloud, so you don't need to convert it. You can find an example of usage in the tests: test019.py. Can you provide a small example of data o a script that does not work? Paul

s-du commented 2 years ago

Thanks a lot for the reply! test019.py helped me indeed! Here is the code that worked for me:

params = cc.SFModulationParams() sub_cloud = cc.CloudSamplingTools.resampleCloudSpatially(cloud, 0.05, params) (spatialCloud, res) = cloud.partialClone(sub_cloud)

Kind regards

ajdesalvio commented 8 months ago

I know the above situation is for a cloud but I'm working with a 3D mesh and don't know how to proceed. I'd like to use the cut2d feature to make several 3D subsamples of a large mesh. I have a shapefile that I imported via geopandas gpd.read_file('shapefile.shp') and for this example, I've saved one of the coordinate sets (call it polygon_coords) as a list of 5 tuples of size 2 (to make a rectangle). An example of what polygon_coords looks like is attached as an image. But I am struggling with formatting this list of vertices in a way that is accepted by cc.Polyline. I'm receiving this error:

poly = cc.Polyline(polygon_coords)

TypeError: init(): incompatible constructor arguments. The following argument types are supported:

  1. _cloudComPy.Polyline(arg0: _cloudComPy.GenericIndexedCloudPersist)

I cannot figure out how to put the vertices in the format of GenericIndexedCloudPersist. Do you have any advice? Thank you very much for your time.

Best regards, Aaron polygon_coords

prascle commented 8 months ago

Hello, It's not straightforward to build a polyline from a set of coordinates. You have to create a point cloud, using 0 for the z coordinates, then you are able to create the polyline from the cloud. Add nodes in order to the cloud. Don't repeat the first node at the end of the list, but use the setClosed method:

coord2 = ((0., 0.), (0., 1.), (1., 1.), (1., 0.))

cloud = cc.ccPointCloud('temp')
cloud.reserve(len(coord2))
for coord in coord2:
    cloud.addPoint((coord[0], coord[1], 0.))

poly = cc.ccPolyline(cloud)
poly.addChild(cloud)
poly.addPointIndex(0, cloud.size())
poly.setClosed(True)

Maybe I can wrap this in a new polyline constructor taking in argument a list of coordinates. Regards, Paul

ajdesalvio commented 8 months ago

Hi Paul,

Thank you very much for the quick response. Your method worked seamlessly for creating a polyline. I have been able to define a polyline from the first bounding box within my shapefile and I'm including my script thus far below. I'm not getting any errors, but the cut_mesh object (which is supposed to be one plant within the field I'm cutting up) is not a ccMesh, it's being labeled as NoneType for some reason. I was hoping to extract a small mesh for every boundary within the shapefile via a loop I attempted to write below. If I'm interpreting correctly, I'm defining the integer value within crop2D as 2 for the Z axis, since the XY plane of the mesh is the ground and I'm using the polyline to cut through it in the Z plane? Please correct me if this is wrong. Thank you again for your time, I'll be sure to cite this repository in future publications and my research lab will benefit greatly from this project.

Best regards, Aaron

import cloudComPy as cc # import the CloudComPy module
cc.initCC() # to do once before dealing with plugins
import geopandas as gpd
import os
from gendata import getSampleCloud, getSamplePoly, dataDir, isCoordEqual, createSymbolicLinks
createSymbolicLinks() # required for tests on build, before cc.initCC.init

# Load mesh
mesh = cc.loadMesh('3Dmesh.obj')
mesh.setName('mesh')
cc.computeNormals([mesh], computePerVertexNormals=False)

# Read shapefile
shapefile = gpd.read_file('shapefile.shp')

# Processing polylines in the shapefile
for index, row in shapefile.iterrows():
    plant_id = row['id']  # Get the plant ID
    geometry = row['geometry']  # Get the geometry for cutting
    polygon_coords = [(x, y) for x, y in geometry.exterior.coords]

    coord2 = tuple(polygon_coords[0:4])
    cloud = cc.ccPointCloud('temp')
    cloud.reserve(len(coord2))
    for coord in coord2:
        cloud.addPoint((coord[0], coord[1], 0.))

    poly = cc.ccPolyline(cloud)
    poly.addChild(cloud)
    poly.addPointIndex(0, cloud.size())
    poly.setClosed(True)

    cut_mesh = mesh.crop2D(poly, 2, False)
    output_filename = f"{plant_id}_cut_mesh.obj"
    cc.SaveMesh(cut_mesh, output_filename)

    break
prascle commented 8 months ago

Hi Aaron, thank you for your encouragement!

The first possible explanation is that you are using unshifted coordinates for your polygons, whereas your mesh is internally shifted during import. See the doc: Load without explicit shift parameters and getting the shift value. If you have large coordinates, you'll need to offset your data to ensure accuracy: this is the default setting when you load the mesh. Retrieve the offset value, and apply it to the coordinates of your polygons to use the same coordinate system. Let me know if it works. Best regards, Paul

ajdesalvio commented 8 months ago

Hi Paul,

Absolutely, and thanks again for your help. As you predicted, the 3D mesh was shifted during import. I accounted for this and updated the script. In case any other users are trying to perform a similar iterative function involving cutting several portions out of a mesh (or cloud), I'll include the functional script below. Thank you again and I am looking forward to continuing to use this wrapper.

Best regards, Aaron

import cloudComPy as cc # import the CloudComPy module
cc.initCC() # to do once before dealing with plugins
import geopandas as gpd
import os
from gendata import getSampleCloud, getSamplePoly, dataDir, isCoordEqual, createSymbolicLinks
createSymbolicLinks() # required for tests on build, before cc.initCC.init

# Load mesh
mesh = cc.loadMesh('3Dmesh.obj')
mesh.setName('mesh')
cc.computeNormals([mesh], computePerVertexNormals=False)
shift = mesh.getGlobalShift() # Retrieves shift value for imported mesh
x_shift = shift[0]
y_shift = shift[1]

# Read shapefile
shapefile = gpd.read_file('shapefile.shp')

# Processing polylines in the shapefile
for index, row in shapefile.iterrows():
    plant_id = row['id']  # Get the plant ID
    geometry = row['geometry']  # Get the geometry for cutting
    polygon_coords = [(x, y) for x, y in geometry.exterior.coords]

    coord2 = tuple(polygon_coords[0:4])
    coord2 = [(x + x_shift, y + y_shift) for x, y in coord2]

    cloud = cc.ccPointCloud('temp')
    cloud.reserve(len(coord2))
    for coord in coord2:
        cloud.addPoint((coord[0], coord[1], 0.))

    poly = cc.ccPolyline(cloud)
    poly.addChild(cloud)
    poly.addPointIndex(0, cloud.size())
    poly.setClosed(True)

    cut_mesh = mesh.crop2D(poly, 2, True)
    output_filename = f"{plant_id}_cut_mesh.obj"
    cc.SaveMesh(cut_mesh, output_filename)
ajdesalvio commented 8 months ago

I thought of something else I wanted to run by you. The large OBJ file contains an associated texture file (in the form of a JPEG or TIFF). Currently, the script produces an .mtl file with each cut segment. But is there a way for it to cut the texture file for each segment too? I'd like to use these segments for deep learning and ideally would have both shape and texture (color) to train the model. Thanks again. Best regards, Aaron

prascle commented 8 months ago

Hi Aaron, I don't know of a tool for cutting out the texture, but I think there may be a way to recreate an image of the cut-out texture, if the top view of the cut-out object is sufficient: the idea is to produce a raster (image.tif of the cloud associated with the mesh). The texture will then be projected onto the cloud points. The first step is to resample the cloud to obtain enough points to define the image (say 1000000). Next, the raster process produces a regular grid based on the points, and the color (texture) is interpolated onto the grid. The raster image is saved as a .tif file. I've tried this with CloudCompare's GUI, and the same can be done with CloudComPy. There are a few parameters to set, depending on your problem. If you need help writing the script, can you post a sample of a cut object with its texture? Best, Paul

ajdesalvio commented 8 months ago

Hi Paul, Thank you for the detailed explanation! That sounds like a viable option. Here is a link to several of the individual .obj files (the plants) and their associated .mtl files, along with the texture file corresponding to the entire field (a .tif). I have found that the .obj files have the correct texture as long as the large texture file is in the same directory. If there are any other files you need, please let me know and I'll upload them. If you have time to assist with the script I would greatly appreciate it, not sure how I can repay you but thank you nonetheless! Best regards, Aaron

prascle commented 8 months ago

Hi Aaron, I wrote the following script, which creates a .tif file colored with texture for each .obj file:

import os
import sys
import math
import fnmatch

os.environ["_CCTRACE_"]="ON" # only if you want C++ debug traces

import cloudComPy as cc

for file in os.listdir('.'):
    if fnmatch.fnmatch(file, '*.obj'):
        name = os.path.splitext(os.path.basename(file))[0]
        print(name)
        mesh = cc.loadMesh(file)
        cloud = mesh.samplePoints(True, 5000000)
        cloud.setName(name)
        cc.RasterizeGeoTiffOnly(cloud,
                                gridStep=0.001,
                                outputRasterRGB =True,
                                pathToImages = '.',
                                emptyCellFillStrategy = cc.EmptyCellFillOption.INTERPOLATE_DELAUNAY,
                                customHeight = 0.)

I'm not sure it's useful in your case, it's only a top view and may be missing a lot of useful information.

Best, Paul

ajdesalvio commented 8 months ago

Hi Paul,

Thank you for your time and effort. I will give this a shot and reply here if I have further questions.

Best regards, Aaron