sparsify_las_based_sq_m.py

import argparse
import os
import numpy as np
import random
import laspy
from scipy.spatial import ConvexHull
import logging


class SparsifyLasBasedSqM:
    def __init__(self, input_file, output_folder=None, target_density=10, verbose=False):
        self.input_file = input_file
        if output_folder is None:
            self.output_folder = os.path.dirname(input_file)
        self.target_density = target_density
        self.desity = None
        self.new_point_cloud_density = None
        self.verbose = verbose

        # Initialize logging
        self.logger = logging.getLogger(__name__)
        if self.verbose:
            self.logger.setLevel(logging.INFO)
        else:
            self.logger.setLevel(logging.WARNING)

        self.logger.info(f"Initialized with input file: {self.input_file}, target density: {self.target_density}")

    def calculate_density_convex_hull(self, las):
        points_3D = np.vstack((las.x, las.y, las.z)).transpose()
        points_2D = points_3D[:, :2]
        hull = ConvexHull(points_2D)
        area = hull.volume
        point_count = len(points_2D)
        density = point_count / area
        return density

    def sparsify(self, point_cloud):


        self.density = self.calculate_density_convex_hull(point_cloud)
        self.logger.info(f"Point cloud density: {self.density} points per square meter.")
        x = point_cloud.x
        keep_count = int(len(x) * (self.target_density / self.density))
        sampled_indices = random.sample(range(len(x)), keep_count)
        filtered_point_cloud = point_cloud.points[sampled_indices]
        self.new_point_cloud_density =self.calculate_density_convex_hull(filtered_point_cloud)
        self.logger.info(f"Reduced point cloud size from {len(x)} to {len(filtered_point_cloud)} points.")
        self.logger.info(f"Reduced point cloud by {(1 - len(filtered_point_cloud) / len(x)) * 100}%.")
        self.logger.info(f"New point cloud density: {self.new_point_cloud_density} points per square meter.")
        return filtered_point_cloud

    def process(self):
        if not os.path.exists(self.output_folder):
            os.makedirs(self.output_folder)
        inFile = laspy.read(self.input_file)
        filtered_points = self.sparsify(inFile)
        self.logger.info("Creating output laspy object")
        outFile = laspy.create(point_format=inFile.point_format, file_version=inFile.header.version)
        outFile.header = inFile.header
        outFile.points = filtered_points
        output_file_path = os.path.join(self.output_folder, os.path.basename(self.input_file).replace(".las", f"_sparse_{self.target_density}.las"))
        outFile.write(output_file_path)


if __name__ == "__main__":
    # Configure logging
    logging.basicConfig(level=logging.INFO)

    parser = argparse.ArgumentParser()
    parser.add_argument("-i", "--input_file", help="The .las file to sparsify.")
    parser.add_argument("--output_folder", default=None, help="The folder where the sparse point cloud will be saved.")
    parser.add_argument("-d", "--target_density", help="Target density of points per square meter.", default=10, type=int)
    parser.add_argument("-v", "--verbose", help="Enable verbose logging", action="store_true")
    args = parser.parse_args()

    sparsifier = SparsifyLasBasedSqM(args.input_file, args.output_folder, args.target_density, args.verbose)
    sparsifier.process()