diff --git a/nibio_preprocessing/distance_filtering_dem_based.py b/nibio_preprocessing/distance_filtering_dem_based.py
new file mode 100644
index 0000000000000000000000000000000000000000..224300bc00699a782dc246079e0fb367b6fefc68
--- /dev/null
+++ b/nibio_preprocessing/distance_filtering_dem_based.py
@@ -0,0 +1,153 @@
+
+import argparse
+
+import laspy
+import numpy as np
+import pandas as pd
+
+from scipy.interpolate import griddata
+from scipy.spatial import KDTree
+
+
+class DistanceFilteringDemBased(object):
+ # MLS data : 1 - ground, 2 - vegetation, 3 - CWD, 4 - trunk
+ GROUND_CLASS = 1
+ TARGET_CLASS = 2
+
+ def __init__(self, las_file, distance, verbose=False):
+ # compute distance filtering based on DEM
+ # remove points which are smaller than distance from the DEM
+ # low vegetation - 0.01 m
+
+ self.point_cloud_file = las_file
+ self.distance = distance
+ self.verbose = verbose
+
+ def read_las_and_put_to_pandas(self, las_file):
+ # Read the file
+ file_content = laspy.read(las_file)
+
+ # Put x, y, z, label into a numpy array
+ points = np.vstack((file_content.x, file_content.y, file_content.z, file_content.label, file_content.treeID)).T
+
+ # Put x, y, z, label into a pandas dataframe
+ points_df = pd.DataFrame(points, columns=['x', 'y', 'z', 'label', 'treeID'])
+
+ # get points which belong to the ground class
+ ground_points = points_df[points_df['label'] == self.GROUND_CLASS]
+
+ # get points which belong to the target class
+ target_points = points_df[points_df['label'] == self.TARGET_CLASS]
+
+ # get just x, y, z columns from ground_points
+ ground_points = ground_points[['x', 'y', 'z']]
+
+ # get just x, y, z columns from target_points
+ target_points = target_points[['x', 'y', 'z']]
+
+ return points_df, ground_points, target_points
+
+ def compute_dem_for_ground(self, ground_points):
+ # DEM - Digital Elevation Model
+
+ # create digital elevation model (DEM) from label 1 points
+ x = ground_points['x']
+ y = ground_points['y']
+ z = ground_points['z']
+
+ # create a grid of points
+ xi = np.linspace(x.min(), x.max(), 1000)
+ yi = np.linspace(y.min(), y.max(), 1000)
+ xi, yi = np.meshgrid(xi, yi)
+
+ # interpolate
+ zi = griddata((x, y), z, (xi, yi), method='linear')
+
+ # fill the NaN values and missing and numertical values with 0
+ zi = np.nan_to_num(zi, copy=False)
+
+ # replace non-numerical values with 0
+ zi[zi == np.inf] = 0
+
+ # reduce precision of the DEM to 2 decimal places
+ zi = np.around(zi, decimals=2)
+
+ # put DEM into pandas dataframe
+ dem = pd.DataFrame({'x': xi.flatten(), 'y': yi.flatten(), 'z': zi.flatten(), })
+
+ return dem
+
+ def compute_distance_between_dem_and_target(self, dem, target_points):
+ original_target_points = target_points.copy()
+ target_points = target_points.values
+
+ # Create the KD-tree using the DEM data
+ dem_tree = KDTree(dem[:, :2])
+
+ # Compute distances for all target points
+ _, id = dem_tree.query(target_points[:, :2], workers=-1)
+
+ # get distance in z direction between target points and dem
+ original_target_points['z_distance'] = original_target_points['z'] - dem[id, 2]
+
+ # add the distance in z direction to the original dataframe as the column 'z_distance'
+ original_target_points['z_distance'] = original_target_points['z_distance'].values
+
+ return original_target_points
+
+
+ def filter_points(self, target_points_with_distances):
+ # filter points
+ target_points_with_distances = target_points_with_distances[target_points_with_distances['z_distance'] < self.distance]
+
+ return target_points_with_distances
+
+
+ def update_las_file(self, points_df, target_points_with_distances):
+ # remove all the points which belong to the target class and are in the target_points_with_distances dataframe
+ points_df = points_df[~points_df.isin(target_points_with_distances)].dropna()
+
+ return points_df
+
+ def save_las_file(self):
+ # save las file
+ pass
+
+ def run(self):
+ points_df, ground_points, target_points = self.read_las_and_put_to_pandas(las_file=self.point_cloud_file)
+ dem = self.compute_dem_for_ground(ground_points=ground_points)
+
+ # save pandas dataframe to csv
+ dem.to_csv('maciek/dem_from_class.csv', index=False, header=True, sep=',')
+
+ target_points_with_distances = self.compute_distance_between_dem_and_target(dem=dem.values, target_points=target_points)
+
+ # save pandas dataframe to csv
+ target_points_with_distances.to_csv('maciek/target_points_with_distances.csv', index=False, header=True, sep=',')
+
+ filtered_points = self.filter_points(target_points_with_distances=target_points_with_distances)
+
+ # save pandas dataframe to csv
+ filtered_points.to_csv('maciek/filtered_points.csv', index=False, header=True, sep=',')
+
+ points_df = self.update_las_file(points_df=points_df, target_points_with_distances=filtered_points)
+
+ # save pandas dataframe to csv
+ points_df.to_csv('maciek/updated_points.csv', index=False, header=True, sep=',')
+
+
+if __name__ == '__main__':
+ # parse the arguments
+ parser = argparse.ArgumentParser(description='Distance filtering based on DEM')
+ parser.add_argument('-i', '--input', help='Input file', required=True)
+ parser.add_argument('-d', '--distance', help='Distance', default=0.01, required=False, type=float)
+ parser.add_argument('-v', '--verbose', help='Verbose', required=False)
+ args = vars(parser.parse_args())
+
+ # get the arguments
+ LAS_FILE = args['input']
+ DISTANCE = args['distance']
+ VERBOSE = args['verbose']
+
+ # run the script
+ DistanceFilteringDemBased(LAS_FILE, DISTANCE, VERBOSE).run()