sparsification implemented

nibio_preprocessing/density_filtering.py

@@ -6,6 +6,16 @@ import glob
 import os
 from pathlib import Path
 
+
+"""
+The pdal density command calculates the point density for each point in the point cloud data. 
+Points in regions of the point cloud with density below the specified min_density threshold
+are flagged or marked. The script then generates a density map and stores it as a shapefile (.shp) on disk.
+
+The density map is later used to filter the point cloud, essentially
+removing regions with a density below the specified min_density threshold. 
+"""
+
 class DensityFiltering:
     def __init__(self, path_data, min_density=1, count_threshold=15000, buffer_size=0.01, verbose=False):
         self.path_data = Path(path_data)

nibio_preprocessing/shrink_point_clound.py

@@ -6,6 +6,12 @@ import laspy # requires laspy==2.1.2 (check if doesn't work)
 import numpy as np
 from tqdm import tqdm
 
+"""
+ The script essentially removes points that are close to the boundaries of the point cloud
+   in each axis (x, y, and optionally z). Specifically, it filters points based on a threshold which is 
+   calculated from the minimum and maximum coordinates along each axis. 
+"""
+
 
 class ShrinkPointCloud():
     def __init__(self, input_folder, output_folder, treshold=0.9, verbose=False):

nibio_preprocessing/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 SparsifyPointCloudToDensity:
+    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.verbose = verbose
+        if self.verbose:
+            logging.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):
+        density = self.calculate_density_convex_hull(point_cloud)
+        if self.verbose:
+            logging.info(f"Point cloud density: {density} points per square meter.")
+        x = point_cloud.x
+        keep_count = int(len(x) * (self.target_density / density))
+        sampled_indices = random.sample(range(len(x)), keep_count)
+        filtered_point_cloud = point_cloud.points[sampled_indices]
+        if self.verbose:
+            logging.info(f"Reduced point cloud size from {len(x)} to {len(filtered_point_cloud)} points.")
+            logging.info(f"Reduced point cloud by {(1 - len(filtered_point_cloud) / len(x)) * 100}%.")
+            density = self.calculate_density_convex_hull(filtered_point_cloud)
+            density = round(density, 2)
+            logging.info(f"Filtered point cloud density: {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)
+        if self.verbose:
+            logging.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
+        # save the point cloud to the output folder with the same name as the input file but suffixed with _sparse and desired density
+        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__":
+    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 = SparsifyPointCloudToDensity(args.input_file, args.output_folder, args.target_density, args.verbose)
+    sparsifier.process()

nibio_preprocessing/sparsify_las_file.py

+import argparse
+import os
+import random
+import laspy
+
+class SparsifyPointCloud:
+    def __init__(self, input_file, output_folder=None, keep_percentage=50.0):
+        self.input_file = input_file
+        if output_folder is None:
+            self.output_folder = os.path.dirname(input_file)
+        self.keep_percentage = keep_percentage
+
+    def sparsify(self, point_cloud):
+        original_count = len(point_cloud.points)
+        keep_count = int(original_count * (self.keep_percentage / 100.0))
+
+        # Randomly sample points
+        sampled_indices = random.sample(range(original_count), keep_count)
+        sparse_point_cloud = point_cloud.points[sampled_indices]
+
+        print(f"Reduced point cloud from {original_count} to {len(sparse_point_cloud)} points.")
+        print(f"Reduced point cloud by {(1 - len(sparse_point_cloud) / original_count) * 100}%.")
+
+        return sparse_point_cloud
+
+    def process(self):
+        # Create output folder if it doesn't exist
+        if not os.path.exists(self.output_folder):
+            os.makedirs(self.output_folder)
+
+        # Read the input .las file
+        inFile = laspy.read(self.input_file)
+
+        # Sparsify the point cloud
+        sparse_point_cloud = self.sparsify(inFile)
+
+        # Create an empty las file and copy the header from the original las file
+        outFile = laspy.create(point_format=inFile.header.point_format, file_version=inFile.header.version)
+        outFile.header = inFile.header
+
+        # Update point count in the header
+        outFile.header.point_count = len(sparse_point_cloud)
+
+        # Assign the points to the output file
+        outFile.points = sparse_point_cloud
+
+        # Write the sparse point cloud to a new .las file
+        output_file_path = os.path.join(self.output_folder, os.path.basename(self.input_file).replace(".las", "_sparse.las"))
+        outFile.write(output_file_path)
+
+if __name__ == "__main__":
+    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", "--keep_percentage", help="The percentage of points to keep.", default=50.0, type=float)
+
+    args = parser.parse_args()
+
+    sparsifier = SparsifyPointCloud(args.input_file, args.output_folder, args.keep_percentage)
+    sparsifier.process()