diff --git a/.gitignore b/.gitignore
index 08b71877..d6262489 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,3 @@
 **/__pycache__
 tmp
-main.log
\ No newline at end of file
+**.log
\ No newline at end of file
diff --git a/configs/configs_lidro.yaml b/configs/configs_lidro.yaml
index 117e4114..2ad66f58 100644
--- a/configs/configs_lidro.yaml
+++ b/configs/configs_lidro.yaml
@@ -8,6 +8,8 @@ work_dir: ${hydra:runtime.cwd}
 
 io:
   input_filename: null
+  input_mask_hydro: null
+  input_skeleton: null
   input_dir: null
   output_dir: null
   srid: 2154
@@ -17,24 +19,33 @@ io:
   no_data_value: -9999
   tile_size: 1000
 
-raster:
-  # size for dilatation
-  dilation_size: 3
+mask_generation:
+  raster:
+    # size for dilatation
+    dilation_size: 3
+  filter:
+    # Classes to be considered as "non-water"
+    keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 64, 65, 66, 67] # All classes
+  vector:
+    # Filter water's area (m²)
+    min_water_area: 150
+    # Parameters for buffer
+    buffer_positive: 1
+    buffer_negative: -1.5 # negative buffer should be bigger than positive buffer to prevent protruding over the banks
+    # Tolerance from Douglas-Peucker
+    tolerance: 1
 
-vector:
-  # Filter water's area (m²)
-  min_water_area: 150
-  # Parameters for buffer
-  buffer_positive: 1
-  buffer_negative: -1.5 # negative buffer should be bigger than positive buffer to prevent protruding over the banks
-  # Tolerance from Douglas-Peucker
-  tolerance: 1
-  # distance in meters between 2 consecutive points 
-  distance_meters: 2
-
-filter:
-  # Classes to be considered as "non-water"
-  keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 64, 65, 66, 67] # All classes
+virtual_point:
+  filter:
+    # Keep ground and water pointclouds between Hydro Mask and Hydro Mask buffer
+    keep_neighbors_classes: [2, 9]
+  vector:
+    # distance in meters between 2 consecutive points 
+    distance_meters: 2
+    # buffer for searching the points classification (default. "2") of the input LAS/LAZ file
+    buffer: 2
+    # The number of nearest neighbors to find with KNeighbors
+    k: 20
 
 hydra:
   output_subdir: null
diff --git a/environment.yml b/environment.yml
index 3d43c141..a44b1992 100644
--- a/environment.yml
+++ b/environment.yml
@@ -11,7 +11,7 @@ dependencies:
   - scipy
   - geojson
   - rasterio
-  - geopandas
+  - geopandas==0.*
   - pyproj
   - pdal>=2.6
   - python-pdal>=3.2.1
@@ -24,4 +24,6 @@ dependencies:
   - flake8
   - isort
   - pre-commit
+  - pip:
+    - -r requirements.txt
 
diff --git a/lidro/create_virtual_point/pointcloud/crop_las.py b/lidro/create_virtual_point/pointcloud/crop_las.py
new file mode 100644
index 00000000..424e0081
--- /dev/null
+++ b/lidro/create_virtual_point/pointcloud/crop_las.py
@@ -0,0 +1,35 @@
+# -*- coding: utf-8 -*-
+""" Filter pointcloud to keep only the ground class """
+from typing import List
+
+import numpy as np
+import pdal
+from shapely.geometry import MultiPolygon
+
+
+def read_filter_and_crop_pointcloud(input_points: str, geom: MultiPolygon, classes: List[int]) -> np.array:
+    """Filter pointcloud to keep only the ground class
+
+    Args:
+        input_points (str): Path to the input LAS/LAZ file
+        geom (MultiPolygon):  An array of WKT or GeoJSON 2D MultiPolygon (Mask Hydro with buffer)
+        classes (List[int]): Classification values to keep for input points
+
+    Returns:
+        np.array : Numpy array containing point coordinates (X, Y, Z) after filtering and croping
+    """
+    # Crop pointcloud by point
+    pipeline = (
+        pdal.Reader.las(filename=input_points, nosrs=True)
+        | pdal.Filter.range(
+            limits=",".join(f"Classification[{v}:{v}]" for v in classes),
+        )
+        | pdal.Filter.crop(
+            polygon=geom,
+        )
+    )
+    pipeline.execute()
+    # extract points
+    cropped_points = pipeline.arrays[0]
+
+    return np.array([cropped_points["X"], cropped_points["Y"], cropped_points["Z"]]).T
diff --git a/lidro/create_virtual_point/stats/calculate_stat.py b/lidro/create_virtual_point/stats/calculate_stat.py
new file mode 100644
index 00000000..d781f2ef
--- /dev/null
+++ b/lidro/create_virtual_point/stats/calculate_stat.py
@@ -0,0 +1,37 @@
+# -*- coding: utf-8 -*-
+""" Calculates statistics from np.array (ex. pointcloud) """
+import numpy as np
+
+
+def calculate_percentile(points: np.array, q: int) -> float:
+    """
+    Calculates the quartile's value of Z coordinates
+
+    Parameters:
+        - points (numpy.ndarray): An array of 3D points [x, y, z]
+        - q (int): Percentage for the percentile to compute.
+                   Values must be between 0 and 100 inclusive.
+
+    Returns:
+        - float: The quartile of Z coordinates
+    """
+    altitudes = points[:, 2]  # Extract the altitude column
+    n_quartile = np.round(np.percentile(altitudes, q), 2)
+
+    return n_quartile
+
+
+def calculate_median(points: np.array) -> float:
+    """
+    Calculates the median's value of Z coordinates
+
+    Parameters:
+        - points (numpy.ndarray): An array of 3D points [x, y, z]
+
+    Returns:
+        - float: The median of Z coordinates
+    """
+    altitudes = points[:, 2]  # Extract the altitude column
+    n_median = np.round(np.median(altitudes), 2)
+
+    return n_median
diff --git a/lidro/create_virtual_point/stats/knn_distance.py b/lidro/create_virtual_point/stats/knn_distance.py
new file mode 100644
index 00000000..7362bcbe
--- /dev/null
+++ b/lidro/create_virtual_point/stats/knn_distance.py
@@ -0,0 +1,37 @@
+# -*- coding: utf-8 -*-
+""" KNN """
+import numpy as np
+from sklearn.neighbors import NearestNeighbors
+
+
+def find_k_nearest_neighbors(points_skeleton: np.array, points_ground_3d: np.array, k: int) -> np.array:
+    """Finds the K nearest neighbors of a given point from a list of 3D points
+
+    Args:
+        points_skeleton (np.array): a single point representing the skeleton Hydro
+                                    an array of 2D point [x, y]
+        points_ground_3d (np.array): ground's points from LIDAR.
+                                    an array of 3D points [x, y, z]
+        k (int): the number of nearest neighbors to find
+
+    Returns:
+        numpy.ndarray: an array of the K nearest 3D points [x, y, z]
+    """
+    # Extract 2D coordinates (x, y) from the 3D points
+    points_ground_2d = points_ground_3d[:, :2]
+
+    # Convert point_2d to nump.array
+    points_skeleton_array = np.array(points_skeleton).reshape(1, -1)
+
+    # Initialize and fit the NearestNeighbors model
+    nbrs = NearestNeighbors(n_neighbors=k, algorithm="auto").fit(points_ground_2d)
+
+    # Find the distances and indices of the nearest neighbors
+    indices = nbrs.kneighbors(points_skeleton_array, return_distance=False)
+
+    # Retrieve the points corresponding to the indices
+    # Use indices  to retrieve the closest 3D points
+    k_nearest_points = points_ground_3d[indices[0]]
+    # = points_3d[indices.flatten()]
+
+    return k_nearest_points
diff --git a/lidro/create_virtual_point/vectors/clip_points_with_bounding_box.py b/lidro/create_virtual_point/vectors/clip_points_with_bounding_box.py
new file mode 100644
index 00000000..77293469
--- /dev/null
+++ b/lidro/create_virtual_point/vectors/clip_points_with_bounding_box.py
@@ -0,0 +1,30 @@
+# -*- coding: utf-8 -*-
+""" Clip Skeleton points by tile (pointcloud)
+"""
+import geopandas as gpd
+from shapely.geometry import box
+
+
+def clip_points_with_box(gdf_points: gpd.GeoDataFrame, bbox: tuple) -> gpd.GeoDataFrame:
+    """Clip skeleton points by tile (bounding box)
+
+    Args:
+        gdf_points (gpd.GeoDataFrame): Points every 2 meters (by default) along skeleton hydro
+        bbox (tuple): bounding box from tile (pointcloud)
+
+    Returns:
+       gpd.GeoDataframe : Points every 2 meters (by default) along skeleton hydro by tile
+    """
+    # Extract the bounding box limits
+    xmin = bbox[0][0]
+    xmax = bbox[0][1]
+    ymin = bbox[1][0]
+    ymax = bbox[1][1]
+
+    # Create a polygon representing the bounding box
+    bounding_box = box(xmin, ymin, xmax, ymax)
+
+    # Clip points using the bounding box
+    clipped_points = gdf_points[gdf_points.within(bounding_box)]
+
+    return clipped_points
diff --git a/lidro/create_virtual_point/vectors/extract_points_around_skeleton.py b/lidro/create_virtual_point/vectors/extract_points_around_skeleton.py
new file mode 100644
index 00000000..a7a5cbda
--- /dev/null
+++ b/lidro/create_virtual_point/vectors/extract_points_around_skeleton.py
@@ -0,0 +1,63 @@
+# -*- coding: utf-8 -*-
+""" Extract points around skeleton by tile
+"""
+import logging
+import os
+from typing import List
+
+import geopandas as gpd
+from pdaltools.las_info import las_get_xy_bounds
+from shapely.geometry import Point
+
+from lidro.create_virtual_point.pointcloud.crop_las import (
+    read_filter_and_crop_pointcloud,
+)
+from lidro.create_virtual_point.vectors.clip_points_with_bounding_box import (
+    clip_points_with_box,
+)
+from lidro.create_virtual_point.vectors.las_around_point import filter_las_around_point
+
+
+def extract_points_around_skeleton_points_one_tile(
+    filename: str,
+    input_dir: str,
+    input_mask_hydro_buffer: gpd.GeoDataFrame,
+    points_skeleton_gdf: gpd.GeoDataFrame,
+    classes: List[int:int],
+    k: int,
+):
+    """Severals steps :
+        1-  Crop filtered pointcloud by Mask Hydro with buffer
+        2-  Extract a Z elevation value along the hydrographic skeleton
+
+    Args:
+        filename (str): filename to the LAS file
+        input_dir (str): input folder
+        input_mask_hydro_buffer (gpd.GeoDataFrame): hydro mask with buffer
+        points_skeleton_gdf (gpd.GeoDataFrame): Points every 2 meters (by default) along skeleton hydro
+        classes (List):  List of classes to use for the filtering
+        k (int): the number of nearest neighbors to find
+
+    Returns:
+        points_clip (np.array) : Numpy array containing point coordinates (X, Y, Z) after filtering and croping
+    """
+    # Step 1 : Crop filtered pointcloud by Mask Hydro with buffer
+    input_dir_points = os.path.join(input_dir, "pointcloud")
+    tilename = os.path.splitext(filename)[0]  # filename to the LAS file
+    input_pointcloud = os.path.join(input_dir_points, filename)  # path to the LAS file
+    logging.info(f"\nCroping filtered pointcloud by Mask Hydro with buffer for tile : {tilename}")
+    # Croping filtered pointcloud by Mask Hydro with buffer for tile
+    points_clip = read_filter_and_crop_pointcloud(input_pointcloud, str(input_mask_hydro_buffer), classes)
+    logging.info(f"\nCropping skeleton points for tile: {tilename}")
+    # Extract bounding box for clipping points by tile
+    bbox = las_get_xy_bounds(input_pointcloud)
+    points_skeleton_clip = clip_points_with_box(points_skeleton_gdf, bbox)
+    # Create list with Skeleton's points with Z for step 4
+    points_skeleton_with_z_clip = [
+        ([geom.x, geom.y]) for geom in points_skeleton_clip.geometry if isinstance(geom, Point)
+    ]
+    # Step 2 : Extract a Z elevation value along the hydrographic skeleton
+    logging.info(f"\nExtract a Z elevation value along the hydrographic skeleton for tile : {tilename}")
+    result = filter_las_around_point(points_skeleton_with_z_clip, points_clip, k)
+
+    return result
diff --git a/lidro/create_virtual_point/vectors/las_around_point.py b/lidro/create_virtual_point/vectors/las_around_point.py
new file mode 100644
index 00000000..ece0dec5
--- /dev/null
+++ b/lidro/create_virtual_point/vectors/las_around_point.py
@@ -0,0 +1,37 @@
+# -*- coding: utf-8 -*-
+""" Extract a Z elevation value every N meters along the hydrographic skeleton
+"""
+from typing import List
+
+import numpy as np
+from shapely.geometry import Point
+
+from lidro.create_virtual_point.stats.calculate_stat import calculate_percentile
+from lidro.create_virtual_point.stats.knn_distance import find_k_nearest_neighbors
+
+
+def filter_las_around_point(points_skeleton: List, points_clip: np.array, k: int) -> List:
+    """Extract a Z elevation value every N meters along the hydrographic skeleton
+
+    Args:
+        points_skeleton (list) : points every N meters (by default: 2) along skeleton Hydro
+        points_clip (np.array): Numpy array containing point coordinates (X, Y, Z) after filtering and croping
+        k (int): The number of nearest neighbors to find
+
+    Returns:
+        List : Result {'geometry': Point 3D, 'z_q1': points KNN}
+    """
+    # Finds the K nearest neighbors of a given point from a list of 3D points
+    points_knn_list = [
+        ({"geometry": Point(point), "points_knn": find_k_nearest_neighbors(point, points_clip, k)})
+        for point in points_skeleton
+    ]
+
+    # Calcule Z "Q1" for each points every N meters along skeleton hydro
+    results = [
+        ({"geometry": p["geometry"], "z_q1": calculate_percentile(p["points_knn"], 10)})
+        for p in points_knn_list
+        if not np.all(np.isnan(p["points_knn"]))
+    ]
+
+    return results
diff --git a/lidro/create_virtual_point/vectors/mask_hydro_with_buffer.py b/lidro/create_virtual_point/vectors/mask_hydro_with_buffer.py
new file mode 100644
index 00000000..c8ab157b
--- /dev/null
+++ b/lidro/create_virtual_point/vectors/mask_hydro_with_buffer.py
@@ -0,0 +1,28 @@
+# -*- coding: utf-8 -*-
+""" Extract a Z elevation value every N meters along the hydrographic skeleton
+"""
+import geopandas as gpd
+from shapely.geometry import CAP_STYLE
+
+
+def import_mask_hydro_with_buffer(file_mask: str, buffer: float, crs: str | int) -> gpd.GeoDataFrame:
+    """Apply buffer (2 meters by default) from Mask Hydro
+
+    Args:
+        file_mask (str): Path from Mask Hydro
+        buffer (float): buffer to apply to the input geometry
+        crs (str | int): Make a CRS from an EPSG code : CRS WKT string
+
+    Returns:
+        gpd.GeoDataFrame : geometry columns are encoded to WKT
+    """
+    # Read Mask Hydro merged
+    gdf = gpd.read_file(file_mask, crs=crs).unary_union
+
+    # Apply buffer (2 meters by default) from Mask Hydro
+    gdf_buffer = gdf.buffer(buffer, cap_style=CAP_STYLE.square)
+
+    # Return a polygon representing the limit of the bank with a buffer of N meters
+    limit_bank = gdf_buffer.difference(gdf)
+
+    return limit_bank
diff --git a/lidro/main_create_mask.py b/lidro/main_create_mask.py
index 3e01fb1c..6c4927f5 100644
--- a/lidro/main_create_mask.py
+++ b/lidro/main_create_mask.py
@@ -45,8 +45,8 @@ def main(config: DictConfig):
     pixel_size = config.io.pixel_size
     tile_size = config.io.tile_size
     crs = CRS.from_user_input(config.io.srid)
-    classe = config.filter.keep_classes
-    dilation_size = config.raster.dilation_size
+    classe = config.mask_generation.filter.keep_classes
+    dilation_size = config.mask_generation.raster.dilation_size
 
     def main_on_one_tile(filename):
         """Lauch main.py on one tile
diff --git a/lidro/main_create_virtual_point.py b/lidro/main_create_virtual_point.py
new file mode 100644
index 00000000..fc610166
--- /dev/null
+++ b/lidro/main_create_virtual_point.py
@@ -0,0 +1,103 @@
+""" Main script for calculate Mask HYDRO 1
+"""
+
+import logging
+import os
+
+import geopandas as gpd
+import hydra
+import pandas as pd
+from omegaconf import DictConfig
+from pyproj import CRS
+
+from lidro.create_virtual_point.vectors.extract_points_around_skeleton import (
+    extract_points_around_skeleton_points_one_tile,
+)
+from lidro.create_virtual_point.vectors.mask_hydro_with_buffer import (
+    import_mask_hydro_with_buffer,
+)
+from lidro.create_virtual_point.vectors.points_along_skeleton import (
+    generate_points_along_skeleton,
+)
+
+
+@hydra.main(config_path="../configs/", config_name="configs_lidro.yaml", version_base="1.2")
+def main(config: DictConfig):
+    """Create a virtual point along hydro surfaces from the points classification of
+    the input LAS/LAZ file and the Hyro Skeleton (GeoJSON) and save it as LAS file.
+
+    It can run either on a single file, or on each file of a folder
+
+    Args:
+        config (DictConfig): hydra configuration (configs/configs_lidro.yaml by default)
+        It contains the algorithm parameters and the input/output parameters
+    """
+    logging.basicConfig(level=logging.INFO)
+
+    # Check input/output files and folders
+    input_dir = config.io.input_dir
+    if input_dir is None:
+        raise ValueError("""config.io.input_dir is empty, please provide an input directory in the configuration""")
+
+    if not os.path.isdir(input_dir):
+        raise FileNotFoundError(f"""The input directory ({input_dir}) doesn't exist.""")
+
+    output_dir = config.io.output_dir
+    if output_dir is None:
+        raise ValueError("""config.io.output_dir is empty, please provide an input directory in the configuration""")
+
+    os.makedirs(output_dir, exist_ok=True)
+
+    # If input filename is not provided, lidro runs on the whole input_dir directory
+    initial_las_filename = config.io.input_filename
+
+    input_mask_hydro = config.io.input_mask_hydro
+    input_skeleton = config.io.input_skeleton
+
+    # Parameters for creating virtual point
+    distance_meters = config.virtual_point.vector.distance_meters
+    buffer = config.virtual_point.vector.buffer
+    crs = CRS.from_user_input(config.io.srid)
+    classes = config.virtual_point.filter.keep_neighbors_classes
+    k = config.virtual_point.vector.k
+
+    # Step 1 : Import Mask Hydro, then apply a buffer
+    # Return GeoDataframe
+    input_mask_hydro_buffer = import_mask_hydro_with_buffer(input_mask_hydro, buffer, crs).wkt
+
+    # Step 2 : Create several points every 2 meters (by default) along skeleton Hydro
+    # Return GeoDataframe
+    points_skeleton_gdf = generate_points_along_skeleton(input_skeleton, distance_meters, crs)
+
+    # Step 3 : Extract points around skeleton by tile
+    if initial_las_filename:
+        # Lauch croping filtered pointcloud by Mask Hydro with buffer by one tile:
+        points_clip = extract_points_around_skeleton_points_one_tile(
+            initial_las_filename, input_dir, input_mask_hydro_buffer, points_skeleton_gdf, classes, k
+        )
+
+    else:
+        # Lauch  croping filtered pointcloud by Mask Hydro with buffer tile by tile
+        input_dir_points = os.path.join(input_dir, "pointcloud")
+        points_clip_list = [
+            extract_points_around_skeleton_points_one_tile(
+                file, input_dir, input_mask_hydro_buffer, points_skeleton_gdf, classes, k
+            )
+            for file in os.listdir(input_dir_points)
+        ]
+        # Flatten the list of lists into a single list of dictionaries
+        points_clip = [item for sublist in points_clip_list for item in sublist]
+
+    # Create a pandas DataFrame from the flattened list
+    df = pd.DataFrame(points_clip)
+    # Create a GeoDataFrame from the pandas DataFrame
+    result_gdf = gpd.GeoDataFrame(df, geometry="geometry")
+    result_gdf.set_crs(crs, inplace=True)
+    # path to the Points along Skeleton Hydro file
+    output_file = os.path.join(output_dir, "Points_Skeleton.GeoJSON")
+    # Save to GeoJSON
+    result_gdf.to_file(output_file, driver="GeoJSON")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/lidro/main_merge_mask.py b/lidro/main_merge_mask.py
index 350625de..49854413 100644
--- a/lidro/main_merge_mask.py
+++ b/lidro/main_merge_mask.py
@@ -43,10 +43,10 @@ def main(config: DictConfig):
           input_dir, 
           output_dir, 
           CRS.from_user_input(config.io.srid),
-          config.vector.min_water_area, 
-          config.vector.buffer_positive, 
-          config.vector.buffer_negative, 
-          config.vector.tolerance)
+          config.mask_generation.vector.min_water_area, 
+          config.mask_generation.vector.buffer_positive, 
+          config.mask_generation.vector.buffer_negative, 
+          config.mask_generation.vector.tolerance)
 
 
 if __name__ == "__main__":
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 00000000..691ef763
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,2 @@
+ign-pdal-tools==1.5.*
+
diff --git a/example_create_mask_by_tile.sh b/scripts/example_create_mask_by_tile.sh
similarity index 100%
rename from example_create_mask_by_tile.sh
rename to scripts/example_create_mask_by_tile.sh
diff --git a/example_create_mask_default.sh b/scripts/example_create_mask_default.sh
similarity index 100%
rename from example_create_mask_default.sh
rename to scripts/example_create_mask_default.sh
diff --git a/scripts/example_create_virtual_point_by_tile.sh b/scripts/example_create_virtual_point_by_tile.sh
new file mode 100644
index 00000000..6a50c135
--- /dev/null
+++ b/scripts/example_create_virtual_point_by_tile.sh
@@ -0,0 +1,10 @@
+# For lauching create virtual points
+python -m lidro.main_create_virtual_point \
+io.input_dir=./data/ \
+io.input_filename=Semis_2021_0830_6291_LA93_IGN69.laz \
+io.input_mask_hydro=./data/merge_mask_hydro/MaskHydro_merge.geosjon \
+io.input_skeleton=./data/skeleton_hydro/skeleton_hydro.geojson \
+io.output_dir=./tmp/ \
+
+
+
diff --git a/scripts/example_create_virtual_point_default.sh b/scripts/example_create_virtual_point_default.sh
new file mode 100644
index 00000000..1ce57db3
--- /dev/null
+++ b/scripts/example_create_virtual_point_default.sh
@@ -0,0 +1,9 @@
+# Launch hydro mask merging
+python -m lidro.main_create_virtual_point \
+io.input_dir=./data/ \
+io.input_mask_hydro=./data/merge_mask_hydro/MaskHydro_merge.geosjon \
+io.input_skeleton=./data/skeleton_hydro/skeleton_hydro.geojson \
+io.output_dir=./tmp/ \
+
+
+
diff --git a/example_merge_mask_default.sh b/scripts/example_merge_mask_default.sh
similarity index 77%
rename from example_merge_mask_default.sh
rename to scripts/example_merge_mask_default.sh
index a2d6e65f..7e162fc7 100644
--- a/example_merge_mask_default.sh
+++ b/scripts/example_merge_mask_default.sh
@@ -1,4 +1,4 @@
-# For lauching Mask Hydro
+# For lauching Mask Hydro merged
 python -m lidro.main_merge_mask \
 io.input_dir=./data/mask_hydro/ \
 io.output_dir=./tmp/merge_mask_hydro/ \
diff --git a/test/pointcloud/test_crop_las.py b/test/pointcloud/test_crop_las.py
new file mode 100644
index 00000000..9f1fac5b
--- /dev/null
+++ b/test/pointcloud/test_crop_las.py
@@ -0,0 +1,39 @@
+import os
+import shutil
+from pathlib import Path
+
+import numpy as np
+
+from lidro.create_virtual_point.pointcloud.crop_las import (
+    read_filter_and_crop_pointcloud,
+)
+
+TMP_PATH = Path("./tmp/create_virtual_point/pointcloud/crop_las")
+LAS_FILE = "./data/pointcloud/Semis_2021_0830_6291_LA93_IGN69.laz"
+
+
+def setup_module(module):
+    if TMP_PATH.is_dir():
+        shutil.rmtree(TMP_PATH)
+    os.makedirs(TMP_PATH)
+
+
+def test_crop_pointcloud_default():
+    classes = [1, 2]
+
+    geom = str(
+        "MULTIPOLYGON (((830873.1249999998 6290475.625000002, \
+            830610.0480769228 6290335.43269231, \
+            830546.0096153843 6290405.528846156, \
+            830803.0288461535 6290559.567307694, \
+        830873.1249999998 6290475.625000002)), \
+            ((830210.5288461539 6290861.298076924, \
+            830205.0480769231 6290760.336538463, \
+            830101.2019230769 6290759.471153848, \
+            830103.2211538461 6290861.586538463, \
+        830210.5288461539 6290861.298076924)))"
+    )
+    output = read_filter_and_crop_pointcloud(LAS_FILE, geom, classes)
+
+    assert isinstance(output, np.ndarray) is True
+    assert isinstance((output == 0).all(), bool) is False
diff --git a/test/stats/test_calculate_stats.py b/test/stats/test_calculate_stats.py
new file mode 100644
index 00000000..3d43af73
--- /dev/null
+++ b/test/stats/test_calculate_stats.py
@@ -0,0 +1,44 @@
+import numpy as np
+import pytest
+
+from lidro.create_virtual_point.stats.calculate_stat import (
+    calculate_median,
+    calculate_percentile,
+)
+
+
+def test_calculate_percentile_25():
+    points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]])
+    assert calculate_percentile(points, 25) == 6.0
+
+
+def test_calculate_percentile_50():
+    points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]])
+    assert calculate_percentile(points, 50) == 9.0
+
+
+def test_calculate_percentile_75():
+    points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]])
+    assert calculate_percentile(points, 75) == 12.0
+
+
+def test_calculate_percentile_100():
+    points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]])
+    assert calculate_percentile(points, 100) == 15.0
+
+
+def test_calculate_percentile_invalid_q():
+    with pytest.raises(ValueError):
+        points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]])
+        calculate_percentile(points, 110)
+
+
+def test_calculate_percentile_negative_q():
+    with pytest.raises(ValueError):
+        points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]])
+        calculate_percentile(points, -10)
+
+
+def test_calculate_median():
+    points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]])
+    assert calculate_median(points) == 9.0
diff --git a/test/stats/test_knn_distance.py b/test/stats/test_knn_distance.py
new file mode 100644
index 00000000..122368c8
--- /dev/null
+++ b/test/stats/test_knn_distance.py
@@ -0,0 +1,21 @@
+import numpy as np
+
+from lidro.create_virtual_point.stats.knn_distance import find_k_nearest_neighbors
+
+
+def test_find_k_nearest_neighbors_default():
+    points_array = np.array(
+        [
+            [830438.91, 6290854.32, 2.56],
+            [830721.84, 6290447.79, 2.23],
+            [830861.04, 6290242.06, 2.78],
+            [830867.61, 6290202.62, 2.89],
+        ]
+    )
+    point = np.array([[830574.89, 6290648.53]])
+    k = 1
+
+    result = find_k_nearest_neighbors(point, points_array, k)
+
+    assert isinstance(result, np.ndarray) is True
+    assert np.array_equal(result, np.array([[830438.91, 6290854.32, 2.56]]))
diff --git a/test/test_main_create_virtual_point.py b/test/test_main_create_virtual_point.py
new file mode 100644
index 00000000..27401864
--- /dev/null
+++ b/test/test_main_create_virtual_point.py
@@ -0,0 +1,217 @@
+import os
+import subprocess as sp
+from pathlib import Path
+
+import pytest
+from hydra import compose, initialize
+
+from lidro.main_create_virtual_point import main
+
+INPUT_DIR = Path("data")
+OUTPUT_DIR = Path("tmp") / "create_virtual_point/main"
+
+
+def setup_module(module):
+    os.makedirs("tmp/create_virtual_point/main", exist_ok=True)
+
+
+def test_main_run_okay():
+    repo_dir = Path.cwd().parent
+    cmd = f"""python -m lidro.main_create_virtual_point \
+        io.input_dir="{repo_dir}/lidro/data/"\
+        io.input_filename=Semis_2021_0830_6291_LA93_IGN69.laz \
+        io.input_mask_hydro="{repo_dir}/lidro/data/merge_mask_hydro/MaskHydro_merge.geojson"\
+        io.input_skeleton="{repo_dir}/lidro/data/skeleton_hydro/Skeleton_Hydro.geojson"\
+        io.output_dir="{repo_dir}/lidro/tmp/create_virtual_point/main/"
+        """
+    sp.run(cmd, shell=True, check=True)
+
+
+def test_main_lidro_default():
+    input_dir = INPUT_DIR
+    output_dir = OUTPUT_DIR / "main_lidro_default"
+    input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson"
+    input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro.geojson"
+    distances_meters = 10
+    buffer = 10
+    srid = 2154
+    k = 20
+
+    with initialize(version_base="1.2", config_path="../configs"):
+        # config is relative to a module
+        cfg = compose(
+            config_name="configs_lidro",
+            overrides=[
+                f"io.input_dir={input_dir}",
+                f"io.output_dir={output_dir}",
+                f"io.input_mask_hydro={input_mask_hydro}",
+                f"io.input_skeleton={input_skeleton}",
+                f"virtual_point.vector.distance_meters={distances_meters}",
+                f"virtual_point.vector.buffer={buffer}",
+                f"virtual_point.vector.k={k}",
+                f"io.srid={srid}",
+            ],
+        )
+    main(cfg)
+    assert (Path(output_dir) / "Points_Skeleton.GeoJSON").is_file()
+
+
+def test_main_lidro_input_file():
+    input_dir = INPUT_DIR
+    output_dir = OUTPUT_DIR / "main_lidro_input_file"
+    input_filename = "Semis_2021_0830_6291_LA93_IGN69.laz"
+    input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson"
+    input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro.geojson"
+    distances_meters = 5
+    buffer = 2
+    srid = 2154
+    k = 10
+    with initialize(version_base="1.2", config_path="../configs"):
+        # config is relative to a module
+        cfg = compose(
+            config_name="configs_lidro",
+            overrides=[
+                f"io.input_filename={input_filename}",
+                f"io.input_dir={input_dir}",
+                f"io.output_dir={output_dir}",
+                f"io.input_mask_hydro={input_mask_hydro}",
+                f"io.input_skeleton={input_skeleton}",
+                f"virtual_point.vector.distance_meters={distances_meters}",
+                f"virtual_point.vector.buffer={buffer}",
+                f"virtual_point.vector.k={k}",
+                f"io.srid={srid}",
+            ],
+        )
+    main(cfg)
+    assert (Path(output_dir) / "Points_Skeleton.GeoJSON").is_file()
+
+
+def test_main_lidro_fail_no_input_dir():
+    output_dir = OUTPUT_DIR / "main_no_input_dir"
+    input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson"
+    input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro.geojson"
+    distances_meters = 5
+    buffer = 2
+    srid = 2154
+    k = 10
+    with initialize(version_base="1.2", config_path="../configs"):
+        # config is relative to a module
+        cfg = compose(
+            config_name="configs_lidro",
+            overrides=[
+                f"io.input_mask_hydro={input_mask_hydro}",
+                f"io.input_skeleton={input_skeleton}",
+                f"io.output_dir={output_dir}",
+                f"virtual_point.vector.distance_meters={distances_meters}",
+                f"virtual_point.vector.buffer={buffer}",
+                f"virtual_point.vector.k={k}",
+                f"io.srid={srid}",
+            ],
+        )
+    with pytest.raises(ValueError):
+        main(cfg)
+
+
+def test_main_lidro_fail_wrong_input_dir():
+    output_dir = OUTPUT_DIR / "main_wrong_input_dir"
+    input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson"
+    input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro.geojson"
+    distances_meters = 5
+    buffer = 2
+    srid = 2154
+    k = 10
+    with initialize(version_base="1.2", config_path="../configs"):
+        # config is relative to a module
+        cfg = compose(
+            config_name="configs_lidro",
+            overrides=[
+                "io.input_dir=some_random_input_dir",
+                f"io.input_mask_hydro={input_mask_hydro}",
+                f"io.input_skeleton={input_skeleton}",
+                f"io.output_dir={output_dir}",
+                f"virtual_point.vector.distance_meters={distances_meters}",
+                f"virtual_point.vector.buffer={buffer}",
+                f"virtual_point.vector.k={k}",
+                f"io.srid={srid}",
+            ],
+        )
+    with pytest.raises(FileNotFoundError):
+        main(cfg)
+
+
+def test_main_lidro_fail_no_output_dir():
+    input_dir = INPUT_DIR
+    input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson"
+    input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro.geojson"
+    distances_meters = 5
+    buffer = 2
+    srid = 2154
+    k = 10
+    with initialize(version_base="1.2", config_path="../configs"):
+        # config is relative to a module
+        cfg = compose(
+            config_name="configs_lidro",
+            overrides=[
+                f"io.input_dir={input_dir}",
+                f"io.input_mask_hydro={input_mask_hydro}",
+                f"io.input_skeleton={input_skeleton}",
+                f"virtual_point.vector.distance_meters={distances_meters}",
+                f"virtual_point.vector.buffer={buffer}",
+                f"virtual_point.vector.k={k}",
+                f"io.srid={srid}",
+            ],
+        )
+    with pytest.raises(ValueError):
+        main(cfg)
+
+
+def test_main_lidro_fail_no_input_mask_hydro():
+    input_dir = INPUT_DIR
+    output_dir = OUTPUT_DIR / "main_no_input_dir"
+    input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro.geojson"
+    distances_meters = 5
+    buffer = 2
+    srid = 2154
+    k = 10
+    with initialize(version_base="1.2", config_path="../configs"):
+        cfg = compose(
+            config_name="configs_lidro",
+            overrides=[
+                f"io.input_dir={input_dir}",
+                f"io.output_dir={output_dir}",
+                "io.input_mask_hydro=some_random_input_mask_hydro",
+                f"io.input_skeleton={input_skeleton}",
+                f"virtual_point.vector.distance_meters={distances_meters}",
+                f"virtual_point.vector.buffer={buffer}",
+                f"virtual_point.vector.k={k}",
+                f"io.srid={srid}",
+            ],
+        )
+    with pytest.raises(ValueError):
+        main(cfg)
+
+
+def test_main_lidro_fail_no_input_skeleton():
+    input_dir = INPUT_DIR
+    output_dir = OUTPUT_DIR / "main_no_input_dir"
+    input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson"
+    distances_meters = 5
+    buffer = 2
+    srid = 2154
+    k = 10
+    with initialize(version_base="1.2", config_path="../configs"):
+        cfg = compose(
+            config_name="configs_lidro",
+            overrides=[
+                f"io.input_dir={input_dir}",
+                f"io.output_dir={output_dir}",
+                f"io.input_mask_hydro={input_mask_hydro}",
+                "io.input_skeleton=some_random_input_skeleton",
+                f"virtual_point.vector.distance_meters={distances_meters}",
+                f"virtual_point.vector.buffer={buffer}",
+                f"virtual_point.vector.k={k}",
+                f"io.srid={srid}",
+            ],
+        )
+    with pytest.raises(ValueError):
+        main(cfg)
diff --git a/test/test_main_merge_mask.py b/test/test_main_merge_mask.py
index a4d8beb3..e7750a09 100644
--- a/test/test_main_merge_mask.py
+++ b/test/test_main_merge_mask.py
@@ -39,10 +39,10 @@ def test_main_lidro_fail_no_input_dir():
             overrides=[
                 f"io.output_dir={output_dir}",
                 f"io.pixel_size={pixel_size}",
-                f"vector.min_water_area={min_water_area}",
-                f"vector.buffer_positive={buffer_positive}",
-                f"vector.buffer_negative={buffer_negative}",
-                f"vector.tolerance={tolerance}",
+                f"mask_generation.vector.min_water_area={min_water_area}",
+                f"mask_generation.vector.buffer_positive={buffer_positive}",
+                f"mask_generation.vector.buffer_negative={buffer_negative}",
+                f"mask_generation.vector.tolerance={tolerance}",
                 f"io.srid={srid}",
             ],
         )
@@ -66,10 +66,10 @@ def test_main_lidro_fail_wrong_input_dir():
                 "io.input_dir=some_random_input_dir",
                 f"io.output_dir={output_dir}",
                 f"io.pixel_size={pixel_size}",
-                f"vector.min_water_area={min_water_area}",
-                f"vector.buffer_positive={buffer_positive}",
-                f"vector.buffer_negative={buffer_negative}",
-                f"vector.tolerance={tolerance}",
+                f"mask_generation.vector.min_water_area={min_water_area}",
+                f"mask_generation.vector.buffer_positive={buffer_positive}",
+                f"mask_generation.vector.buffer_negative={buffer_negative}",
+                f"mask_generation.vector.tolerance={tolerance}",
                 f"io.srid={srid}",
             ],
         )
@@ -92,10 +92,10 @@ def test_main_lidro_fail_no_output_dir():
             overrides=[
                 f"io.input_dir={input_dir}",
                 f"io.pixel_size={pixel_size}",
-                f"vector.min_water_area={min_water_area}",
-                f"vector.buffer_positive={buffer_positive}",
-                f"vector.buffer_negative={buffer_negative}",
-                f"vector.tolerance={tolerance}",
+                f"mask_generation.vector.min_water_area={min_water_area}",
+                f"mask_generation.vector.buffer_positive={buffer_positive}",
+                f"mask_generation.vector.buffer_negative={buffer_negative}",
+                f"mask_generation.vector.tolerance={tolerance}",
                 f"io.srid={srid}",
             ],
         )
diff --git a/test/vectors/test_clip_points_with_bounding_box.py b/test/vectors/test_clip_points_with_bounding_box.py
new file mode 100644
index 00000000..d60b42e0
--- /dev/null
+++ b/test/vectors/test_clip_points_with_bounding_box.py
@@ -0,0 +1,39 @@
+import geopandas as gpd
+from shapely.geometry import Point
+
+from lidro.create_virtual_point.vectors.clip_points_with_bounding_box import (
+    clip_points_with_box,
+)
+
+
+def test_create_hydro_vector_mask_default():
+    # Define test points (both inside and outside the bounding box)
+    points = [
+        (830378.719, 6290999.42),  # Inside
+        (830382.382, 6290990.115),  # Inside
+        (830386.045, 6290980.81),  # Inside
+        (830389.709, 6290971.505),  # Inside
+        (830216.893, 6291850.548),  # Outside
+        (830217.273, 6291840.555),  # Outside
+        (830375.000, 6291001.000),  # Outside
+        (830400.000, 6291000.000),  # Outside
+    ]
+    # Convert points to GeoDataFrame
+    points_gdf = gpd.GeoDataFrame(geometry=[Point(point) for point in points])
+
+    # Define a bounding box that should clip out some points
+    bbox = ((830375, 830390), (6290970, 6291000))
+
+    # Expected points within the bounding box
+    expected_points = [
+        (830378.719, 6290999.42),
+        (830382.382, 6290990.115),
+        (830386.045, 6290980.81),
+        (830389.709, 6290971.505),
+    ]
+    # Convert expected points to GeoDataFrame for comparison
+    expected_gdf = gpd.GeoDataFrame(geometry=[Point(point) for point in expected_points])
+
+    result_gdf = clip_points_with_box(points_gdf, bbox)
+
+    assert (result_gdf.equals(expected_gdf)) is True
diff --git a/test/vectors/test_create_mask_hydro_buffer.py b/test/vectors/test_create_mask_hydro_buffer.py
new file mode 100644
index 00000000..6591425a
--- /dev/null
+++ b/test/vectors/test_create_mask_hydro_buffer.py
@@ -0,0 +1,44 @@
+import os
+import shutil
+from pathlib import Path
+
+import geopandas as gpd
+from pyproj import CRS
+from shapely.geometry import MultiPolygon
+
+from lidro.create_virtual_point.vectors.mask_hydro_with_buffer import (
+    import_mask_hydro_with_buffer,
+)
+
+TMP_PATH = Path("./tmp/create_virtual_point/vectors/create_mask_hydro_buffer/")
+
+file_mask = "./data/merge_mask_hydro/MaskHydro_merge.geojson"
+output = Path("./tmp/create_virtual_point/vectors/create_mask_hydro_buffer/MaskHydro_merge_buffer.geojson")
+
+
+def setup_module(module):
+    if TMP_PATH.is_dir():
+        shutil.rmtree(TMP_PATH)
+    os.makedirs(TMP_PATH)
+
+
+def test_create_mask_hydro_buffer_default():
+    # Parameters
+    buffer = 2
+    crs = CRS.from_epsg(2154)
+
+    result = import_mask_hydro_with_buffer(file_mask, buffer, crs)
+
+    # Save the result to GeoJSON
+    d = {"nom": "mask_hydro_merge_buffer", "geometry": [result]}
+    gdf = gpd.GeoDataFrame(d, crs=crs)
+    gdf.to_file(output, driver="GeoJSON", crs=crs)
+
+    assert Path(output).exists()
+
+    gdf = gpd.read_file(output)
+
+    assert not gdf.empty  # GeoDataFrame shouldn't empty
+    assert gdf.crs.to_string() == crs  # CRS is identical
+
+    assert all(isinstance(geom, MultiPolygon) for geom in gdf.geometry)  # All geometry should MultiPolygone
diff --git a/test/vectors/test_las_around_point.py b/test/vectors/test_las_around_point.py
new file mode 100644
index 00000000..5cdfb756
--- /dev/null
+++ b/test/vectors/test_las_around_point.py
@@ -0,0 +1,53 @@
+import os
+import shutil
+from pathlib import Path
+
+import geopandas as gpd
+import numpy as np
+from pyproj import CRS
+from shapely.geometry import Point
+
+from lidro.create_virtual_point.vectors.las_around_point import filter_las_around_point
+
+TMP_PATH = Path("./tmp/create_virtual_point/vectors/las_around_point")
+
+output = Path("./tmp/create_virtual_point/vectors/las_around_point/Points.geojson")
+
+
+def setup_module(module):
+    if TMP_PATH.is_dir():
+        shutil.rmtree(TMP_PATH)
+    os.makedirs(TMP_PATH)
+
+
+def test_las_around_point_default():
+    # Parameters
+    points_along_skeleton = [[830864.5181373736, 6290217.943739296], [830866.5957826116, 6290208.162525126]]
+
+    points_clip = np.array(
+        [
+            [8.30822700e05, 6.29000133e06, 2.59000000e00],
+            [8.30836950e05, 6.29000254e06, 2.47000000e00],
+            [8.30837730e05, 6.29000379e06, 2.58000000e00],
+            [8.30042740e05, 6.29041476e06, 3.09000000e00],
+            [8.30033610e05, 6.29041739e06, 3.87000000e00],
+            [8.30042180e05, 6.29041447e06, 3.10000000e00],
+        ]
+    )
+
+    k = 6
+    crs = CRS.from_epsg(2154)
+
+    result = filter_las_around_point(points_along_skeleton, points_clip, k)
+
+    # Convert results to GeoDataFrame
+    result_gdf = gpd.GeoDataFrame(result)
+    result_gdf.set_crs(crs, inplace=True)
+    # Save to GeoJSON
+    result_gdf.to_file(output, driver="GeoJSON")
+
+    gdf = gpd.read_file(output)
+
+    assert not gdf.empty  # GeoDataFrame shouldn't empty
+    assert gdf.crs.to_string() == crs  # CRS is identical
+    assert all(isinstance(geom, Point) for geom in gdf.geometry)  # All geometry should Points
diff --git a/test/vectors/test_points_along_skeleton.py b/test/vectors/test_points_along_skeleton.py
index 97c3329c..19e09be7 100644
--- a/test/vectors/test_points_along_skeleton.py
+++ b/test/vectors/test_points_along_skeleton.py
@@ -11,10 +11,10 @@
     generate_points_along_skeleton,
 )
 
-TMP_PATH = Path("./tmp/create_point_virtual/vectors/points_along_skeleton")
+TMP_PATH = Path("./tmp/create_virtual_point/vectors/points_along_skeleton")
 
 file = "./data/skeleton_hydro/Skeleton_Hydro.geojson"
-output = Path("./tmp/create_point_virtual/vectors/points_along_skeleton/Points.geojson")
+output = Path("./tmp/create_virtual_point/vectors/points_along_skeleton/Points.geojson")
 
 
 def setup_module(module):