diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 5530026d..827cb3bf 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -11,6 +11,7 @@ Unreleased
 
 **Added**
 * Documentation updates
+* rdp libigl function
 
 **Changed**
 
diff --git a/docs/examples/01_planar_slicing_simple.rst b/docs/examples/01_planar_slicing_simple.rst
index b71f64b5..577b3912 100644
--- a/docs/examples/01_planar_slicing_simple.rst
+++ b/docs/examples/01_planar_slicing_simple.rst
@@ -26,7 +26,7 @@ The first step is to import the required functions:
     from compas_slicer.slicers import PlanarSlicer
     from compas_slicer.post_processing import generate_brim
     from compas_slicer.post_processing import generate_raft
-    from compas_slicer.post_processing import simplify_paths_rdp
+    from compas_slicer.post_processing import simplify_paths_rdp_igl
     from compas_slicer.post_processing import seams_smooth
     from compas_slicer.print_organization import PlanarPrintOrganizer
     from compas_slicer.print_organization import set_extruder_toggle
@@ -106,14 +106,14 @@ that are offset from the bottom layer, to improve adhesion to the build plate
                   raft_layers=1)
 
 Depending on the amount of faces that your input mesh has, a very large amount of 
-points can be generated. ``simplify_paths_rdp`` is a function that removes points
+points can be generated. ``simplify_paths_rdp`` or ``simplify_paths_rdp_igl`` are functions that remove points
 that do not have a high impact on the final shape of the polyline. Increase the
 threshold value to remove more points, decrease it to remove less. For more 
 information on how the algorithm works see: `Ramer–Douglas–Peucker algorithm <https://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm>`_
 
 .. code-block:: python
 
-    simplify_paths_rdp(slicer, threshold=0.6)
+    simplify_paths_rdp_igl(slicer, threshold=0.6)
 
 Currently the 'seam' between different layers of our shape is a 'hard seam',
 the printer would move up almost vertically to move to the next layer. 
@@ -213,7 +213,7 @@ The completed final script can be found below:
     from compas_slicer.slicers import PlanarSlicer
     from compas_slicer.post_processing import generate_brim
     from compas_slicer.post_processing import generate_raft
-    from compas_slicer.post_processing import simplify_paths_rdp
+    from compas_slicer.post_processing import simplify_paths_rdp_igl
     from compas_slicer.post_processing import seams_smooth
     from compas_slicer.print_organization import PlanarPrintOrganizer
     from compas_slicer.print_organization import set_extruder_toggle
@@ -277,7 +277,7 @@ The completed final script can be found below:
     # Simplify the paths by removing points with a certain threshold
     # change the threshold value to remove more or less points
     # ==========================================================================
-    simplify_paths_rdp(slicer, threshold=0.6)
+    simplify_paths_rdp_igl(slicer, threshold=0.6)
 
     # ==========================================================================
     # Smooth the seams between layers
diff --git a/docs/examples/02_curved_slicing_simple.rst b/docs/examples/02_curved_slicing_simple.rst
index 96e73337..70d8f065 100644
--- a/docs/examples/02_curved_slicing_simple.rst
+++ b/docs/examples/02_curved_slicing_simple.rst
@@ -29,7 +29,7 @@ Imports and initialization
     import logging
     import compas_slicer.utilities as utils
     from compas_slicer.slicers import InterpolationSlicer
-    from compas_slicer.post_processing import simplify_paths_rdp
+    from compas_slicer.post_processing import simplify_paths_rdp_igl
     from compas_slicer.pre_processing import InterpolationSlicingPreprocessor
     from compas_slicer.print_organization import set_extruder_toggle, set_linear_velocity_by_range
     from compas_slicer.print_organization import add_safety_printpoints
@@ -104,7 +104,7 @@ options are available for all slicers.
     # post processing
     generate_brim(slicer, layer_width=3.0, number_of_brim_offsets=5)
     seams_smooth(slicer, smooth_distance=10)
-    simplify_paths_rdp(slicer, threshold=1.0)
+    simplify_paths_rdp_igl(slicer, threshold=1.0)
     slicer.printout_info()
     utils.save_to_json(slicer.to_data(), OUTPUT_PATH, 'curved_slicer.json')
 
@@ -166,7 +166,7 @@ The completed final script can be found below:
     import logging
     import compas_slicer.utilities as utils
     from compas_slicer.slicers import InterpolationSlicer
-    from compas_slicer.post_processing import simplify_paths_rdp
+    from compas_slicer.post_processing import simplify_paths_rdp_igl
     from compas_slicer.pre_processing import InterpolationSlicingPreprocessor
     from compas_slicer.print_organization import set_extruder_toggle, set_linear_velocity_by_range
     from compas_slicer.print_organization import add_safety_printpoints
@@ -217,7 +217,7 @@ The completed final script can be found below:
     generate_brim(slicer, layer_width=3.0, number_of_brim_offsets=5)
     seams_smooth(slicer, smooth_distance=10)
 
-    simplify_paths_rdp(slicer, threshold=1.0)
+    simplify_paths_rdp_igl(slicer, threshold=0.5)
     slicer.printout_info()
     utils.save_to_json(slicer.to_data(), OUTPUT_PATH, 'curved_slicer.json')
 
diff --git a/docs/examples/03_planar_slicing_vertical_sorting.rst b/docs/examples/03_planar_slicing_vertical_sorting.rst
index e79357f7..211ab5f7 100644
--- a/docs/examples/03_planar_slicing_vertical_sorting.rst
+++ b/docs/examples/03_planar_slicing_vertical_sorting.rst
@@ -28,7 +28,7 @@ the default method. The example below demonstrates how planar paths can be sorte
     from compas_slicer.pre_processing import move_mesh_to_point
     from compas_slicer.slicers import PlanarSlicer
     from compas_slicer.post_processing import generate_brim
-    from compas_slicer.post_processing import simplify_paths_rdp
+    from compas_slicer.post_processing import simplify_paths_rdp_igl
     from compas_slicer.post_processing import sort_into_vertical_layers
     from compas_slicer.post_processing import reorder_vertical_layers
     from compas_slicer.post_processing import seams_smooth
@@ -69,7 +69,7 @@ the default method. The example below demonstrates how planar paths can be sorte
 
         # Post-processing
         generate_brim(slicer, layer_width=3.0, number_of_brim_offsets=5)
-        simplify_paths_rdp(slicer, threshold=0.7)
+        simplify_paths_rdp_igl(slicer, threshold=0.7)
         seams_smooth(slicer, smooth_distance=10)
         slicer.printout_info()
         save_to_json(slicer.to_data(), OUTPUT_DIR, 'slicer_data.json')
diff --git a/docs/examples/04_gcode_generation.rst b/docs/examples/04_gcode_generation.rst
index 7bac2b85..d9f40bb0 100644
--- a/docs/examples/04_gcode_generation.rst
+++ b/docs/examples/04_gcode_generation.rst
@@ -16,7 +16,7 @@ to materialize them in a typical desktop 3D printer. The gcode generation is sti
     from compas_slicer.pre_processing import move_mesh_to_point
     from compas_slicer.slicers import PlanarSlicer
     from compas_slicer.post_processing import generate_brim
-    from compas_slicer.post_processing import simplify_paths_rdp
+    from compas_slicer.post_processing import simplify_paths_rdp_igl
     from compas_slicer.post_processing import seams_smooth
     from compas_slicer.print_organization import PlanarPrintOrganizer
     from compas_slicer.print_organization import set_extruder_toggle
@@ -49,7 +49,7 @@ to materialize them in a typical desktop 3D printer. The gcode generation is sti
         slicer = PlanarSlicer(compas_mesh, slicer_type="cgal", layer_height=4.5)
         slicer.slice_model()
         generate_brim(slicer, layer_width=3.0, number_of_brim_offsets=4)
-        simplify_paths_rdp(slicer, threshold=0.6)
+        simplify_paths_rdp_igl(slicer, threshold=0.6)
         seams_smooth(slicer, smooth_distance=10)
         slicer.printout_info()
         save_to_json(slicer.to_data(), OUTPUT_DIR, 'slicer_data.json')
diff --git a/docs/examples/05_non_planar_slicing_on_custom_base.rst b/docs/examples/05_non_planar_slicing_on_custom_base.rst
index 89aeaf41..81c1fc7b 100644
--- a/docs/examples/05_non_planar_slicing_on_custom_base.rst
+++ b/docs/examples/05_non_planar_slicing_on_custom_base.rst
@@ -22,7 +22,7 @@ vertex of the mesh. In this case we create a scalar field with the distance of e
     from compas.datastructures import Mesh
     import os
     import compas_slicer.utilities as slicer_utils
-    from compas_slicer.post_processing import simplify_paths_rdp
+    from compas_slicer.post_processing import simplify_paths_rdp_igl
     from compas_slicer.slicers import ScalarFieldSlicer
     import compas_slicer.utilities as utils
     from compas_slicer.print_organization import ScalarFieldPrintOrganizer
@@ -57,11 +57,10 @@ vertex of the mesh. In this case we create a scalar field with the distance of e
         # --- Slice model by generating contours of scalar field
         slicer = ScalarFieldSlicer(mesh, u, no_of_isocurves=50)
         slicer.slice_model()
-        # simplify_paths_rdp(slicer, threshold=0.3)
         slicer_utils.save_to_json(slicer.to_data(), OUTPUT_PATH, 'isocontours.json')  # save results to json
 
         # --- Print organization calculations (i.e. generation of printpoints with fabrication-related information)
-        simplify_paths_rdp(slicer, threshold=0.3)
+        simplify_paths_rdp_igl(slicer, threshold=0.3)
         print_organizer = ScalarFieldPrintOrganizer(slicer, parameters={}, DATA_PATH=DATA_PATH)
         print_organizer.create_printpoints()
 
diff --git a/docs/examples/06_attributes_transfer.rst b/docs/examples/06_attributes_transfer.rst
index 52280406..70d9d67d 100644
--- a/docs/examples/06_attributes_transfer.rst
+++ b/docs/examples/06_attributes_transfer.rst
@@ -25,7 +25,7 @@ barycentric coordinates.
     from compas.geometry import Point, Vector, distance_point_plane, normalize_vector
     from compas.datastructures import Mesh
     import compas_slicer.utilities as slicer_utils
-    from compas_slicer.post_processing import simplify_paths_rdp
+    from compas_slicer.post_processing import simplify_paths_rdp_igl
     from compas_slicer.slicers import PlanarSlicer
     import compas_slicer.utilities.utils as utils
     from compas_slicer.utilities.attributes_transfer import transfer_mesh_attributes_to_printpoints
@@ -77,7 +77,7 @@ barycentric coordinates.
         # --------------- Slice mesh
         slicer = PlanarSlicer(mesh, slicer_type="default", layer_height=5.0)
         slicer.slice_model()
-        simplify_paths_rdp(slicer, threshold=1.0)
+        simplify_paths_rdp_igl(slicer, threshold=1.0)
         slicer_utils.save_to_json(slicer.to_data(), OUTPUT_PATH, 'slicer_data.json')
 
         # --------------- Create printpoints
diff --git a/docs/installation.rst b/docs/installation.rst
index 2620c9e9..af2bc40c 100644
--- a/docs/installation.rst
+++ b/docs/installation.rst
@@ -107,3 +107,17 @@ Make sure you are in the correct environment and type:
 .. code-block:: bash
 
     pip install numpy==1.19.3
+
+Fractions error
+-----------
+.. code-block:: bash
+
+    ImportError: cannot import name 'gcd' from 'fractions' (C:\ProgramData\Anaconda3\envs\compas_slicer\lib\fractions.py)
+
+This issue can be solved, as explained here:  https://stackoverflow.com/questions/66174862/import-error-cant-import-name-gcd-from-fractions
+by typing the following command (make sure you are in the correct environment)
+
+.. code-block:: bash
+
+    conda install -c conda-forge networkx=2.5
+
diff --git a/examples/1_planar_slicing_simple/example_1_planar_slicing_simple.py b/examples/1_planar_slicing_simple/example_1_planar_slicing_simple.py
index 10d425f9..757a27ec 100644
--- a/examples/1_planar_slicing_simple/example_1_planar_slicing_simple.py
+++ b/examples/1_planar_slicing_simple/example_1_planar_slicing_simple.py
@@ -7,7 +7,7 @@
 from compas_slicer.slicers import PlanarSlicer
 from compas_slicer.post_processing import generate_brim
 from compas_slicer.post_processing import generate_raft
-from compas_slicer.post_processing import simplify_paths_rdp
+from compas_slicer.post_processing import simplify_paths_rdp_igl
 from compas_slicer.post_processing import seams_smooth
 from compas_slicer.print_organization import PlanarPrintOrganizer
 from compas_slicer.print_organization import set_extruder_toggle
@@ -72,7 +72,7 @@ def main():
     # Simplify the paths by removing points with a certain threshold
     # change the threshold value to remove more or less points
     # ==========================================================================
-    simplify_paths_rdp(slicer, threshold=0.6)
+    simplify_paths_rdp_igl(slicer, threshold=0.6)
 
     # ==========================================================================
     # Smooth the seams between layers
diff --git a/examples/2_curved_slicing/ex2_curved_slicing_advanced.py b/examples/2_curved_slicing/ex2_curved_slicing_advanced.py
index 5bc18cf3..252f98c6 100644
--- a/examples/2_curved_slicing/ex2_curved_slicing_advanced.py
+++ b/examples/2_curved_slicing/ex2_curved_slicing_advanced.py
@@ -3,7 +3,7 @@
 import logging
 import compas_slicer.utilities as utils
 from compas_slicer.slicers import InterpolationSlicer
-from compas_slicer.post_processing import simplify_paths_rdp
+from compas_slicer.post_processing import simplify_paths_rdp_igl
 from compas_slicer.pre_processing import InterpolationSlicingPreprocessor
 from compas_slicer.pre_processing import create_mesh_boundary_attributes
 from compas_slicer.print_organization import InterpolationPrintOrganizer
@@ -86,7 +86,7 @@ def main():
                 generate_brim(slicer, layer_width=3.0, number_of_brim_offsets=5)
 
             seams_smooth(slicer, smooth_distance=0.1)
-            simplify_paths_rdp(slicer, threshold=0.25)
+            simplify_paths_rdp_igl(slicer, threshold=0.25)
             utils.save_to_json(slicer.to_data(), OUTPUT_PATH, 'curved_slicer_%d.json' % i)
             slicers.append(slicer)
 
diff --git a/examples/2_curved_slicing/ex2_curved_slicing_basic.py b/examples/2_curved_slicing/ex2_curved_slicing_basic.py
index a5b96642..02868cac 100644
--- a/examples/2_curved_slicing/ex2_curved_slicing_basic.py
+++ b/examples/2_curved_slicing/ex2_curved_slicing_basic.py
@@ -3,7 +3,7 @@
 import logging
 import compas_slicer.utilities as utils
 from compas_slicer.slicers import InterpolationSlicer
-from compas_slicer.post_processing import simplify_paths_rdp
+from compas_slicer.post_processing import simplify_paths_rdp_igl
 from compas_slicer.pre_processing import InterpolationSlicingPreprocessor
 from compas_slicer.print_organization import set_extruder_toggle, set_linear_velocity_by_range
 from compas_slicer.print_organization import add_safety_printpoints
@@ -55,7 +55,7 @@ def main():
     generate_brim(slicer, layer_width=3.0, number_of_brim_offsets=5)
     seams_smooth(slicer, smooth_distance=10)
 
-    simplify_paths_rdp(slicer, threshold=1.0)
+    simplify_paths_rdp_igl(slicer, threshold=0.5)
     slicer.printout_info()
     utils.save_to_json(slicer.to_data(), OUTPUT_PATH, 'curved_slicer.json')
 
diff --git a/examples/3_planar_slicing_vertical_sorting/example_3_planar_vertical_sorting.py b/examples/3_planar_slicing_vertical_sorting/example_3_planar_vertical_sorting.py
index 7da981a4..119b95c6 100644
--- a/examples/3_planar_slicing_vertical_sorting/example_3_planar_vertical_sorting.py
+++ b/examples/3_planar_slicing_vertical_sorting/example_3_planar_vertical_sorting.py
@@ -5,7 +5,7 @@
 from compas_slicer.pre_processing import move_mesh_to_point
 from compas_slicer.slicers import PlanarSlicer
 from compas_slicer.post_processing import generate_brim
-from compas_slicer.post_processing import simplify_paths_rdp
+from compas_slicer.post_processing import simplify_paths_rdp_igl
 from compas_slicer.post_processing import sort_into_vertical_layers
 from compas_slicer.post_processing import reorder_vertical_layers
 from compas_slicer.post_processing import seams_smooth
@@ -46,7 +46,7 @@ def main():
 
     # Post-processing
     generate_brim(slicer, layer_width=3.0, number_of_brim_offsets=5)
-    simplify_paths_rdp(slicer, threshold=0.7)
+    simplify_paths_rdp_igl(slicer, threshold=0.7)
     seams_smooth(slicer, smooth_distance=10)
     slicer.printout_info()
     save_to_json(slicer.to_data(), OUTPUT_DIR, 'slicer_data.json')
diff --git a/examples/4_gcode_generation/example_4_gcode.py b/examples/4_gcode_generation/example_4_gcode.py
index 75127d8c..9583e236 100644
--- a/examples/4_gcode_generation/example_4_gcode.py
+++ b/examples/4_gcode_generation/example_4_gcode.py
@@ -4,7 +4,7 @@
 from compas_slicer.pre_processing import move_mesh_to_point
 from compas_slicer.slicers import PlanarSlicer
 from compas_slicer.post_processing import generate_brim
-from compas_slicer.post_processing import simplify_paths_rdp
+from compas_slicer.post_processing import simplify_paths_rdp_igl
 from compas_slicer.post_processing import seams_smooth
 from compas_slicer.print_organization import PlanarPrintOrganizer
 from compas_slicer.print_organization import set_extruder_toggle
@@ -37,7 +37,7 @@ def main():
     slicer = PlanarSlicer(compas_mesh, slicer_type="cgal", layer_height=4.5)
     slicer.slice_model()
     generate_brim(slicer, layer_width=3.0, number_of_brim_offsets=4)
-    simplify_paths_rdp(slicer, threshold=0.6)
+    simplify_paths_rdp_igl(slicer, threshold=0.6)
     seams_smooth(slicer, smooth_distance=10)
     slicer.printout_info()
     save_to_json(slicer.to_data(), OUTPUT_DIR, 'slicer_data.json')
diff --git a/examples/5_non_planar_slicing_on_custom_base/scalar_field_slicing.py b/examples/5_non_planar_slicing_on_custom_base/scalar_field_slicing.py
index 1d93168e..ffd1b45f 100644
--- a/examples/5_non_planar_slicing_on_custom_base/scalar_field_slicing.py
+++ b/examples/5_non_planar_slicing_on_custom_base/scalar_field_slicing.py
@@ -3,7 +3,7 @@
 from compas.datastructures import Mesh
 import os
 import compas_slicer.utilities as slicer_utils
-from compas_slicer.post_processing import simplify_paths_rdp
+from compas_slicer.post_processing import simplify_paths_rdp_igl
 from compas_slicer.slicers import ScalarFieldSlicer
 import compas_slicer.utilities as utils
 from compas_slicer.print_organization import ScalarFieldPrintOrganizer
@@ -42,7 +42,7 @@
     slicer_utils.save_to_json(slicer.to_data(), OUTPUT_PATH, 'isocontours.json')  # save results to json
 
     # --- Print organization calculations (i.e. generation of printpoints with fabrication-related information)
-    simplify_paths_rdp(slicer, threshold=0.3)
+    simplify_paths_rdp_igl(slicer, threshold=0.3)
     print_organizer = ScalarFieldPrintOrganizer(slicer, parameters={}, DATA_PATH=DATA_PATH)
     print_organizer.create_printpoints()
 
diff --git a/examples/6_attributes_transfer/example_6_attributes_transfer.py b/examples/6_attributes_transfer/example_6_attributes_transfer.py
index 40b91005..10a7646d 100644
--- a/examples/6_attributes_transfer/example_6_attributes_transfer.py
+++ b/examples/6_attributes_transfer/example_6_attributes_transfer.py
@@ -3,7 +3,7 @@
 from compas.geometry import Point, Vector, distance_point_plane, normalize_vector
 from compas.datastructures import Mesh
 import compas_slicer.utilities as slicer_utils
-from compas_slicer.post_processing import simplify_paths_rdp
+from compas_slicer.post_processing import simplify_paths_rdp_igl
 from compas_slicer.slicers import PlanarSlicer
 import compas_slicer.utilities.utils as utils
 from compas_slicer.utilities.attributes_transfer import transfer_mesh_attributes_to_printpoints
@@ -55,7 +55,7 @@
     # --------------- Slice mesh
     slicer = PlanarSlicer(mesh, slicer_type="default", layer_height=5.0)
     slicer.slice_model()
-    simplify_paths_rdp(slicer, threshold=1.0)
+    simplify_paths_rdp_igl(slicer, threshold=1.0)
     slicer_utils.save_to_json(slicer.to_data(), OUTPUT_PATH, 'slicer_data.json')
 
     # --------------- Create printpoints
diff --git a/scripts/planar_slicing_igl.py b/scripts/planar_slicing_igl.py
index fb0d4f47..50a8fa4e 100644
--- a/scripts/planar_slicing_igl.py
+++ b/scripts/planar_slicing_igl.py
@@ -1,77 +1,161 @@
-import compas
-from compas.datastructures import Mesh
-from compas.geometry import Point, distance_point_point
-from compas_slicer.geometry import Path
+import itertools
+from compas.geometry import Point, distance_point_point_sqrd
 from compas_slicer.geometry import Layer
-from compas_slicer.utilities import save_to_json
-import logging
+from compas_slicer.geometry import Path
 import numpy as np
-import scipy
+import copy
 import networkx as nx
-
-import igl
-
-logger = logging.getLogger('logger')
+import logging
+import compas_slicer.utilities as utils
+from compas_slicer.slicers.slice_utilities import sort_graph_connected_components
+import progressbar
 
 __all__ = ['create_planar_paths_igl']
 
+logger = logging.getLogger('logger')
 
-def create_planar_paths_igl(mesh, min_z, max_z, planes):
-
-    vs, fi = mesh.to_vertices_and_faces()
-    z = [v[2] for v in vs]  # height value
-
-    isoV, isoE = igl.isolines(np.array(vs), np.array(fi), np.array(z), len(planes))
 
-    # create adjacency matrix A
-    G = nx.Graph()
+def try_to_create_connection(G, isoV, ei, ej, i, j, side_i, side_j, connections_found, tol):
+    vi = isoV[ei[side_i]]
+    vj = isoV[ej[side_j]]
+    if distance_point_point_sqrd(vi, vj) < tol:
+        G.add_edge(i, j)
+        connections_found[i][side_i] = True
+        connections_found[j][side_j] = True
+        return True
+    return False
+
+def create_planar_paths_igl(mesh, n):
+    """
+    Creates planar contours using the libigl function: https://libigl.github.io/libigl-python-bindings/igl_docs/#isolines
+    TODO: ??? It is currently the only method that can return identify OPEN versus CLOSED paths.
+
+    Parameters
+    ----------
+    mesh: :class: 'compas.datastructures.Mesh'
+        The mesh to be sliced
+    n: number of contours
+    """
+    utils.check_package_is_installed('igl')
+    import igl
+    v, f = mesh.to_vertices_and_faces()
+
+    ds = np.array([vertex[2] for vertex in v])
+    # save distances of each vertex from the floor: this will be the scalar field for the slicing operation.
+
+    # --- generate disconnected segments of isolines using the libigl function
+    v = np.array(v)
+    f = np.array(f)
+    isoV, isoE = igl.isolines(v, f, ds, n)
+
+    # --- group resulting segments per level
+    print('Grouping segments per level')
+    segments_per_level = {}
+    tol = 1e-6
     for e in isoE:
-        if e[0] != e[1]:
-            G.add_edge(*e)
-    A = nx.to_scipy_sparse_matrix(G)
-    C = igl.connected_components(A)
-
-    print(C[0])
-    # print(type(C[0]))
-    # print(type(C[1]))
-    # print(type(C[2]))
-
-    iso = {}
-    iso['isoV'] = [[v[0], v[1], v[2]] for v in isoV]
-    iso['isoE'] = [[int(v[0]), int(v[1])] for v in isoE]
-
-    save_to_json(data=iso, filepath='/examples/WIP/3_curved_slicing/data/',
-                 name='iso.json')
-    print(nx.number_connected_components(G))
-
-    es = {}
-    for i, e in enumerate(nx.dfs_edges(G, 3000)):
-        print(e)
-        es[i] = [int(e[0]), int(e[1])]
-
-    save_to_json(data=es, filepath='/examples/WIP/3_curved_slicing/data/',
-                 name='es.json')
-
-    layers = {}
-    for j, cp in enumerate(nx.connected_components(G)):
-        if j==0:
-            print (cp)
-        layer = []
-        g = nx.Graph();
-        for point_index in cp:
-            p = isoV[point_index]
-            layer.append([p[0], p[1], p[2]])
-
-        layers[j] = layer
-
-    print(layers[0])
-
-    save_to_json(data=layers, filepath ='/examples/WIP/3_curved_slicing/data/', name='isolines.json')
-
-    raise ValueError
-
-
-
-
-if __name__ == "__main__":
-    pass
+        v0 = isoV[e[0]]
+        v1 = isoV[e[1]]
+        assert(abs(v0[2] - v1[2]) < tol)
+        h = v0[2]
+
+        found = False
+        for key in segments_per_level:
+            if abs(key - h) < tol:
+                if e[0] != e[1]:  # do not add null segments
+                    segments_per_level[key].append(e)
+                    found = True
+
+        if not found:
+            segments_per_level[h] = [e]
+
+    #assert(len(segments_per_level) == n)
+
+    utils.save_to_json(utils.point_list_to_dict(isoV), "/examples/1_planar_slicing_simple/data/output", 'isoV.json')
+
+    sorted_keys = [key for key in segments_per_level]
+    sorted(sorted_keys)
+
+    layers = []
+
+    # --- Create connectivity graph G for every group, and sort edges based on G
+    with progressbar.ProgressBar(max_value=len(segments_per_level)) as bar:
+        for index, key in enumerate(sorted_keys):
+            #print(' Creating connectivity for level z = ', key)
+
+            es = segments_per_level[key]  # list of the current edges, which whose indices we are working below
+
+            G = nx.Graph()
+
+            connections_found = [[False, False] for _ in es]
+
+            for i, e in enumerate(es):
+                G.add_node(i)  # node, attribute
+
+            for i, ei in enumerate(es):  # ei here is the edge, not the index!
+                for side_i in range(2):
+                    if not connections_found[i][side_i]:
+                        for jj, ej in enumerate(es[i+1:]):  # ej here is the edge, not the index!
+                            j = i + jj + 1
+                            for side_j in range(2):
+                                if not connections_found[j][side_j]:
+                                    vi = isoV[ei[side_i]]
+                                    vj = isoV[ej[side_j]]
+                                    if distance_point_point_sqrd(vi, vj) < tol:
+                                        G.add_edge(i, j)
+                                        connections_found[i][side_i] = True
+                                        connections_found[j][side_j] = True
+
+            # sort connected components of G
+            sorted_indices_dict = sort_graph_connected_components(G)
+
+            # for i, s_key in enumerate(sorted_indices_dict):
+            #     sorted_eis = sorted_indices_dict[s_key]
+            #     this_segment_pts = []
+            #
+            #     for j, ei in enumerate(sorted_eis):
+            #         e = es[ei]
+            #         # segment pts
+            #         for k in range(2):
+            #             this_segment_pts.append(isoV[e[k]])
+            #
+            #     utils.save_to_json(utils.point_list_to_dict(this_segment_pts),
+            #                        "C:/dev/compas_slicer/examples/1_planar_slicing_simple/data/output", 'isoV.json')
+            #     print('saved')
+
+            # get points list from sorted edge indices (for each connected component)
+            paths_pts = [[] for _ in sorted_indices_dict]
+            are_closed = [True for _ in sorted_indices_dict]
+            for i, s_key in enumerate(sorted_indices_dict):
+                sorted_indices = sorted_indices_dict[s_key]
+
+                for j, s_ei in enumerate(sorted_indices):
+                    v0 = isoV[es[s_ei][0]]
+                    v1 = isoV[es[s_ei][1]]
+                    if j == 0:
+                        s_ei_next = sorted_indices[j+1]
+                        v0_next = isoV[es[s_ei_next][0]]
+                        v1_next = isoV[es[s_ei_next][1]]
+                        v_mid_next = (Point(*v0_next) + Point(*v1_next)) * 0.5
+                        if distance_point_point_sqrd(v0, v_mid_next) < distance_point_point_sqrd(v1, v_mid_next):
+                            paths_pts[i].extend([Point(*v1), Point(*v0)])
+                        else:
+                            paths_pts[i].extend([Point(*v0), Point(*v1)])
+
+                    else:
+                        v_prev = paths_pts[i][-1]
+                        if distance_point_point_sqrd(v0, v_prev) < distance_point_point_sqrd(v1, v_prev):
+                            paths_pts[i].append(Point(*v1))
+                        else:
+                            paths_pts[i].append(Point(*v0))
+
+                if distance_point_point_sqrd(paths_pts[i][0], paths_pts[i][1]) > tol:
+                    are_closed[i] = False
+
+            paths = []
+            for i in range(len(sorted_indices_dict)):
+                paths.append(Path(points=paths_pts[i], is_closed=are_closed[i]))
+
+            layers.append(Layer(paths))
+            bar.update(index)
+
+    return layers
diff --git a/src/compas_slicer/post_processing/simplify_paths_rdp.py b/src/compas_slicer/post_processing/simplify_paths_rdp.py
index ff52e121..2cc8d0d3 100644
--- a/src/compas_slicer/post_processing/simplify_paths_rdp.py
+++ b/src/compas_slicer/post_processing/simplify_paths_rdp.py
@@ -2,14 +2,18 @@
 import numpy as np
 import logging
 import progressbar
-
 from compas.geometry import Point
+import compas_slicer.utilities as utils
+from compas.plugins import PluginNotInstalledError
 
-# from compas_slicer.geometry import PrintPoint, Contour
+packages = utils.TerminalCommand('conda list').get_split_output_strings()
+if 'igl' in packages:
+    import igl
 
 logger = logging.getLogger('logger')
 
-__all__ = ['simplify_paths_rdp']
+__all__ = ['simplify_paths_rdp',
+           'simplify_paths_rdp_igl']
 
 
 def simplify_paths_rdp(slicer, threshold):
@@ -23,7 +27,6 @@ def simplify_paths_rdp(slicer, threshold):
     threshold: float
         Controls the degree of polyline simplification.
         Low threshold removes few points, high threshold removes many points.
-
     """
 
     logger.info("Paths simplification rdp")
@@ -40,5 +43,34 @@ def simplify_paths_rdp(slicer, threshold):
         logger.info('%d Points remaining after rdp simplification' % remaining_pts_num)
 
 
+def simplify_paths_rdp_igl(slicer, threshold):
+    """
+    https://libigl.github.io/libigl-python-bindings/igl_docs/#ramer_douglas_peucker
+    Parameters
+    ----------
+    slicer: :class:`compas_slicer.slicers.BaseSlicer`
+        An instance of one of the compas_slicer.slicers classes.
+    threshold: float
+        Controls the degree of polyline simplification.
+        Low threshold removes few points, high threshold removes many points.
+    """
+    try:
+        utils.check_package_is_installed('igl')
+        logger.info("Paths simplification rdp - igl")
+        remaining_pts_num = 0
+
+        for i, layer in enumerate(slicer.layers):
+            if not layer.is_raft:  # no simplification necessary for raft layer
+                for path in layer.paths:
+                    pts = np.array([[pt[0], pt[1], pt[2]] for pt in path.points])
+                    S, J, Q = igl.ramer_douglas_peucker(pts, threshold)
+                    path.points = [Point(pt[0], pt[1], pt[2]) for pt in S]
+        logger.info('%d Points remaining after rdp simplification' % remaining_pts_num)
+
+    except PluginNotInstalledError:
+        logger.info("Libigl is not installed. Falling back to python rdp function")
+        simplify_paths_rdp(slicer, threshold)
+
+
 if __name__ == "__main__":
     pass
diff --git a/src/compas_slicer/pre_processing/preprocessing_utils/geodesics.py b/src/compas_slicer/pre_processing/preprocessing_utils/geodesics.py
index 53f38ca5..18d553e9 100644
--- a/src/compas_slicer/pre_processing/preprocessing_utils/geodesics.py
+++ b/src/compas_slicer/pre_processing/preprocessing_utils/geodesics.py
@@ -1,30 +1,17 @@
 import numpy as np
 import logging
 import compas_slicer.utilities as utils
-from compas.plugins import PluginNotInstalledError
 from compas_slicer.pre_processing.preprocessing_utils.gradient import get_scalar_field_from_gradient, \
     get_face_gradient_from_scalar_field, normalize_gradient
 import scipy
 import math
 
-packages = utils.TerminalCommand('conda list').get_split_output_strings()
-if 'igl' in packages:
-    import igl
-
 logger = logging.getLogger('logger')
 
 __all__ = ['get_igl_EXACT_geodesic_distances',
            'get_custom_HEAT_geodesic_distances']
 
 
-def check_igl_is_installed():
-    """ Throws an error if igl python bindings are not installed in the current environment. """
-    if 'igl' not in packages:
-        raise PluginNotInstalledError("--------ATTENTION! ----------- \
-                        Libigl library is missing! \
-                        Install it with 'conda install -c conda-forge igl'")
-
-
 def get_igl_EXACT_geodesic_distances(mesh, vertices_start):
     """
     Calculate geodesic distances using libigl.
@@ -34,7 +21,8 @@ def get_igl_EXACT_geodesic_distances(mesh, vertices_start):
     mesh: :class: 'compas.datastructures.Mesh'
     vertices_start: list, int
     """
-    check_igl_is_installed()
+    utils.check_package_is_installed('igl')
+    import igl
 
     v, f = mesh.to_vertices_and_faces()
     v = np.array(v)
@@ -47,8 +35,6 @@ def get_igl_EXACT_geodesic_distances(mesh, vertices_start):
 
 def get_custom_HEAT_geodesic_distances(mesh, vi_sources, OUTPUT_PATH, v_equalize=None, anisotropic_scaling=False):
     """ Calculate geodesic distances using the heat method. """
-    check_igl_is_installed()
-
     geodesics_solver = GeodesicsSolver(mesh, OUTPUT_PATH)
     u = geodesics_solver.diffuse_heat(vi_sources, v_equalize, method='simulation')
     geodesic_dist = geodesics_solver.get_geodesic_distances(u, vi_sources, v_equalize)
@@ -75,6 +61,9 @@ class GeodesicsSolver:
     """
 
     def __init__(self, mesh, OUTPUT_PATH):
+        utils.check_package_is_installed('igl')
+        import igl
+
         logger.info('GeodesicsSolver')
         self.mesh = mesh
         self.OUTPUT_PATH = OUTPUT_PATH
diff --git a/src/compas_slicer/print_organization/base_print_organizer.py b/src/compas_slicer/print_organization/base_print_organizer.py
index c8c64343..be595271 100644
--- a/src/compas_slicer/print_organization/base_print_organizer.py
+++ b/src/compas_slicer/print_organization/base_print_organizer.py
@@ -133,7 +133,7 @@ def number_of_paths_on_layer(self, layer_index):
     ######################
 
     def remove_duplicate_points_in_path(self, layer_key, path_key, tolerance=0.0001):
-        """Remove subsequent points that are within a certain tolerance.
+        """Remove subsequent points that are within a certain threshold.
 
         Parameters
         ----------
diff --git a/src/compas_slicer/slicers/planar_slicing/planar_slicing_cgal.py b/src/compas_slicer/slicers/planar_slicing/planar_slicing_cgal.py
index 53c60464..a6ed1c0a 100644
--- a/src/compas_slicer/slicers/planar_slicing/planar_slicing_cgal.py
+++ b/src/compas_slicer/slicers/planar_slicing/planar_slicing_cgal.py
@@ -7,11 +7,6 @@
 import compas_slicer.utilities as utils
 from compas.plugins import PluginNotInstalledError
 
-packages = utils.TerminalCommand('conda list').get_split_output_strings()
-
-if 'compas-cgal' in packages or 'compas_cgal' in packages:
-    from compas_cgal.slicer import slice_mesh
-
 logger = logging.getLogger('logger')
 
 __all__ = ['create_planar_paths_cgal']
@@ -28,7 +23,11 @@ def create_planar_paths_cgal(mesh, planes):
         A compas mesh.
     planes: list, :class: 'compas.geometry.Plane'
     """
-    if 'compas-cgal' not in packages and 'compas_cgal' not in packages:
+    packages = utils.TerminalCommand('conda list').get_split_output_strings()
+
+    if 'compas-cgal' in packages or 'compas_cgal' in packages:
+        from compas_cgal.slicer import slice_mesh
+    else:
         raise PluginNotInstalledError("--------ATTENTION! ----------- \
                         Compas_cgal library is missing! \
                         You can't use this planar slicing method without it. \
diff --git a/src/compas_slicer/utilities/__init__.py b/src/compas_slicer/utilities/__init__.py
index 090464e1..3b39559a 100644
--- a/src/compas_slicer/utilities/__init__.py
+++ b/src/compas_slicer/utilities/__init__.py
@@ -34,8 +34,8 @@
 from __future__ import division
 from __future__ import print_function
 
-from .utils import *  # noqa: F401 E402 F403
 from .terminal_command import *  # noqa: F401 E402 F403
+from .utils import *  # noqa: F401 E402 F403
 from .attributes_transfer import *  # noqa: F401 E402 F403
 
 __all__ = [name for name in dir() if not name.startswith('_')]
diff --git a/src/compas_slicer/utilities/utils.py b/src/compas_slicer/utilities/utils.py
index aca8817b..eaf8a9b0 100644
--- a/src/compas_slicer/utilities/utils.py
+++ b/src/compas_slicer/utilities/utils.py
@@ -7,6 +7,8 @@
 import networkx as nx
 import numpy as np
 import scipy
+from compas.plugins import PluginNotInstalledError
+from compas_slicer.utilities import TerminalCommand
 
 logger = logging.getLogger('logger')
 
@@ -36,7 +38,8 @@
            'smooth_vectors',
            'get_normal_of_path_on_xy_plane',
            'get_all_files_with_name',
-           'get_closest_mesh_normal_to_pt']
+           'get_closest_mesh_normal_to_pt',
+           'check_package_is_installed']
 
 
 def remap(input_val, in_from, in_to, out_from, out_to):
@@ -399,6 +402,7 @@ def get_mesh_cotmatrix_igl(mesh, fix_boundaries=True):
     :class: 'scipy.sparse.csr_matrix'
         sparse matrix (dimensions: #V x #V), laplace operator, each row i corresponding to v(i, :)
     """
+    check_package_is_installed('igl')
     import igl
     v, f = mesh.to_vertices_and_faces()
     C = igl.cotmatrix(np.array(v), np.array(f))
@@ -427,6 +431,7 @@ def get_mesh_cotans_igl(mesh):
     :class: 'np.array'
         Dimensions: F by 3 list of 1/2*cotangents corresponding angles
     """
+    check_package_is_installed('igl')
     import igl
     v, f = mesh.to_vertices_and_faces()
     return igl.cotmatrix_entries(np.array(v), np.array(f))
@@ -600,5 +605,17 @@ def get_all_files_with_name(startswith, endswith, DATA_PATH):
     return files
 
 
+#######################################
+# check installation
+
+
+def check_package_is_installed(package_name):
+    """ Throws an error if igl python bindings are not installed in the current environment. """
+    packages = TerminalCommand('conda list').get_split_output_strings()
+    if package_name not in packages:
+        raise PluginNotInstalledError(" ATTENTION! Package : " + package_name +
+                                      " is missing! Please follow installation guide to install it.")
+
+
 if __name__ == "__main__":
     pass