diff --git a/src/openalea/caribu/CaribuScene.py b/src/openalea/caribu/CaribuScene.py
index faff400e..538b1652 100644
--- a/src/openalea/caribu/CaribuScene.py
+++ b/src/openalea/caribu/CaribuScene.py
@@ -29,7 +29,7 @@
 from openalea.caribu.display import jet_colors, generate_scene, nan_to_zero
 from openalea.caribu.caribu_shell import vperiodise, Path
 from functools import reduce
-from openalea.caribu.caributriangleset import AbstractCaribuTriangleSet, CaribuTriangleSet 
+from openalea.caribu.caribu_triangle_set import AbstractCaribuTriangleSet, CaribuTriangleSet
 
 import tempfile
 
diff --git a/src/openalea/caribu/caribu_shell.py b/src/openalea/caribu/caribu_shell.py
index 6a76ab53..cd43af34 100644
--- a/src/openalea/caribu/caribu_shell.py
+++ b/src/openalea/caribu/caribu_shell.py
@@ -9,26 +9,113 @@
 #       WebSite : https://github.com/openalea-incubator/caribu
 #
 # ==============================================================================
-""" Chaine les appel de s2v, mcsail et canestra, comme le fait cpfg via caribu
-    Syntaxe: caribu.csh Ds file.light file.can nz h file.8 file1.opt ... fileN.opt
-
-  Exemple /home/tld/chelle/QHS/Calgary/dev/Canestra/MC-SAIL/Test
-  caribu.csh 0.2 sky2.light Tout.can 3 18 Tout.8 test nir par
-
-  MC98
-  C. Pradal - 2009
-  MC09
-  Contact: chelle@grinon.inra.fr
-  INRA - INRIA - CIRAD
+""" Pythonic implementation of caribu.csh shell script (see original script at the end of this file)
 """
-from subprocess import Popen, STDOUT, PIPE
+
 import tempfile
-import platform
 from pathlib import Path
 import shutil
-import openalea.libcaribu.tools as lc
+from copy import deepcopy
+import openalea.libcaribu.algos as lcal
 import openalea.libcaribu.io as lcio
+import numpy as np
+# alias for set_scene
+from openalea.libcaribu.algos import set_scene
+
+
+
+
+def caribu(scene_path, bands=None, direct_only=True, toric=False, mixed_radiosity=0, layers=2, height=1,
+           screen_size=None, sensors=False, soil=False, artifacts=False, outdir=None, verbose=False):
+    """ Low level interface to caribu algorithms
+
+    Args:
+        - scene_path (Path) : a path to a dir containing scene data files with normalised names (scene.can, ...).
+          See openalea.caribu.caribu_shell.set_scene for instanciation
+        - bands (str or list of str): the name of the optical band to compute. If None (default), stem parts of *.opt files present in scene_path are used
+        - direct_only (bool): consider only first order illumination (default true)
+        - toric (bool): Consider a toric canopy. Needs a pattern file to be present in the scene_path to take effect
+        - mixed_radiosity (float) : control the mixed radiosity behavior. '-1' meeans 'no mixed radiosity' (pure radiosity, no sail),
+         '0' (default) means 'no radiosity' (sail only), any other number defines the diameter of the spherical boundary
+          between radiosity/non radiosity domain around primitives
+        - layers: number of layers to be considered for discretising the scene for sail
+        - height: height of the highest layer to use for sail
+        - screen_size: the size (pixel) of the diagonal of the projection screen
+        - soil: should soil computations be activated ? (default False). requires a scene.soil in the scene
+        - sensors: should virtual sensor be activated ? (default False). requires a scene.sensor file in the scene_path
+        - artifacts: should canestra debugging artifacts (B.dat, Bz.dat, ...) be generated (default False) ?
+        - outdir: path where output files are written. If None (default), no output files are generated
+
+    Returns:
+        A {band_0: (result, measures), ...} dict containing the incident and absorbed flux of energy for all primitives
+        """
+    if bands is None:
+        bands = [opt.stem for opt in scene_path.glob("*.opt")]
+
+    if isinstance(bands, str):
+        bands = [bands]
+    else:
+        bands = list(bands)
+
+    if outdir:
+        outdir = Path(outdir)
+        outdir.mkdir(exist_ok=True)
+
+    if not direct_only and mixed_radiosity >= 0:
+        toric = True
+        if not all([(scene_path / f'{band}.vec').exists() for band in bands]):
+            lcal.s2v(scene_path, bands=bands, layers=layers, height=height, verbose=verbose)
+
+    if toric and not (scene_path / 'motif.can').exists():
+        lcal.periodise(scene_path, verbose=verbose)
+
+    args = []
+    if screen_size:
+        args += ["-L", str(screen_size)]
+    if sensors:
+        args += ['-C', 'scene.sensor']
+    if verbose:
+        args += ['-v', '2']
+    if not artifacts:
+        args += ['-n']
+
+    res = {}
+    FF_path = None
+    for i, band in enumerate(bands):
+        more_args = []
+        more_args += args
+        if direct_only:
+            if i == 0:
+                if toric:
+                    res[band] = lcal.toric_raycasting(scene_path, band=band, soil=soil, more_args=more_args, verbose=verbose)
+                else:
+                    res[band] = lcal.raycasting(scene_path, band=band, soil=soil, more_args=more_args, verbose=verbose)
+            else:
+                opticals = lcio.read_opt(scene_path / f'{band}.opt')
+                r, s, m = deepcopy(res[bands[0]])
+                alpha = lcio.absorptance_from_labels(r['label'], opticals)
+                r['Eabs'] = alpha * r['Ei']
+                res[band] = r, s, m
+        else:
+            if i == 0:
+                FF_path = scene_path / 'FF'
+                FF_path.mkdir(exist_ok=True)
+                more_args += ['-t', str(FF_path),
+                              '-f', 'scene.FF']
+            else:
+                more_args += ['-t', str(FF_path),
+                              '-w', 'scene.FF']
+            if mixed_radiosity < 0:
+                res[band] = lcal.radiosity(scene_path, band=band, soil=soil, more_args=more_args, verbose=verbose)
+            else:
+                res[band] = lcal.mixed_radiosity(scene_path, band=band, soil=soil, sd=mixed_radiosity, more_args=more_args, verbose=verbose)
+
+        if outdir:
+            shutil.copy(scene_path / 'Etri.vec0', outdir / f'{band}.vec0')
+            if i == 0 and FF_path:
+                shutil.copy(FF_path / 'scene.FF', outdir)
 
+    return res
 
 def _safe_iter(obj, atomic_types=(str, int, float, complex)):
     """Equivalent to iter when obj is iterable and not defined as atomic.
diff --git a/src/openalea/caribu/caribu_triangle_set.py b/src/openalea/caribu/caribu_triangle_set.py
new file mode 100644
index 00000000..1151299c
--- /dev/null
+++ b/src/openalea/caribu/caribu_triangle_set.py
@@ -0,0 +1,169 @@
+from openalea.caribu.display import generate_scene
+from openalea.libcaribu.io import read_can, decode_labels
+
+import numpy
+from itertools import chain
+
+class AbstractCaribuTriangleSet:
+    def __init__(self):
+        self.bbox = None
+        pass
+
+    def getBoundingBox(self):
+        pass
+
+    def triangle_areas(self):
+        pass
+
+    def getZmin(self):
+        pass
+
+    def __getshape__(self, shapeid):
+        """ Return all triangles of a shape """
+        pass
+
+    def keys(self):
+        pass
+
+    def values(self):
+        pass
+
+    def shapes(self):
+        pass
+
+    def allvalues(self, copied=False):
+        pass
+
+    def allids(self):
+        pass
+
+    def getNumberOfTriangles(self, shapeid):
+        pass
+
+    def generate_scene(self, colorproperty):
+        pass
+
+    def __len__(self):
+        raise NotImplemented()
+
+class CaribuTriangleSet(AbstractCaribuTriangleSet):
+    def __init__(self, tri_dict, scales=None):
+        AbstractCaribuTriangleSet.__init__(self)
+        self.triangle_set = tri_dict
+        all_tris = chain.from_iterable(tri_dict.values())
+        self.all_triangles = numpy.fromiter(
+            (coord for tri in all_tris for pt in tri for coord in pt),
+            dtype=float
+        ).reshape(-1, 3, 3)
+        self.counts = numpy.fromiter((len(v) for v in tri_dict.values()), dtype=int)
+        self.scales = {'shape': list(tri_dict.keys())}
+        self.by_scale = {'shape': numpy.repeat(numpy.arange(len(self.counts)), self.counts)}
+        if scales is not None:
+            orgs, by_org = self.scales['shape'], self.by_scale['shape']
+            for scale, mapping in scales.items():
+                groups = [mapping[o] for o in orgs]
+                group_names, id_to_groupid = numpy.unique(groups, return_inverse=True)
+                self.scales[scale] = group_names
+                self.by_scale[scale]  = id_to_groupid[by_org]
+
+    @classmethod
+    def fromfile(cls, source):
+        triangles, labels = read_can(source)
+        tri_dict = {}
+        unique_ids, inv = numpy.unique(labels, return_inverse=True)
+        specie, plant, leaf, element = decode_labels(unique_ids)
+        for i, uid in enumerate(unique_ids):
+            tri_dict[uid] = triangles[inv == i]
+        plant_scale = dict(zip(unique_ids, plant))
+        elt_scale = dict(zip(unique_ids, element))
+        sp_scale = dict(zip(unique_ids, specie))
+        leaf_scale = dict(zip(unique_ids, leaf))
+        return cls(tri_dict, scales={'plant': plant_scale, 'specie': sp_scale, 'leaf': leaf_scale, 'elt': elt_scale})
+
+
+    def getBoundingBox(self):
+        if self.bbox is None:
+            x, y, z = self.all_triangles.reshape(-1, 3).T
+            self.bbox = (x.min(), y.min(), z.min()), (x.max(), y.max(), z.max()) 
+        return self.bbox     
+        
+    def triangle_areas(self):
+        """ compute area of elementary triangles in the scene """
+
+        A = self.all_triangles[:, 0, :]
+        B = self.all_triangles[:, 1, :]
+        C = self.all_triangles[:, 2, :]
+
+        return 0.5 * numpy.linalg.norm(numpy.cross(B - A, C - A), axis=1)
+
+
+    def sum_at_scale(self, arr, scale='shape'):
+            if isinstance(scale, (list, tuple)):
+                invs = [self.by_scale[s] for s in scale]
+                sizes = [len(self.scales[s]) for s in scale]
+
+                mult = numpy.cumprod([1] + sizes[:-1])
+                by = sum(i * m for i, m in zip(invs, mult))
+
+                # build combined keys
+                keys = [
+                    tuple(self.scales[s][(i // m) % size]
+                          for s, m, size in zip(scale, mult, sizes))
+                    for i in range(mult[-1] * sizes[-1])
+                ]
+
+            else:
+                keys = self.scales[scale]
+                by = self.by_scale[scale]
+
+            _sum = numpy.bincount(by, weights=arr, minlength=len(keys))
+
+            return dict(zip(keys, _sum))
+
+    def mean_at_scale(self, arr, scale=None):
+        return self.sum_at_scale(arr / self.counts, scale=scale)
+
+
+
+    def getZmin(self):
+        return self.getBoundingBox()[0][2]
+
+    def getZmax(self):
+        return self.getBoundingBox()[1][2]
+
+    def __getshape__(self, shapeid):
+        """ Return all triangles of a shape """
+        return self.triangle_set[shapeid]
+
+    def __len__(self):
+        return len(self.triangle_set)
+
+    def keys(self):
+        return self.triangle_set.keys()
+
+    def values(self):
+        return self.triangle_set.values()
+
+    def shapes(self):
+        return self.triangle_set.shapes()
+
+    def allvalues(self, copied=False):
+        from copy import copy
+        if copied : 
+            return copy(self.allpoints)
+        else:
+            return self.allpoints
+    
+    def allids(self):
+        return self.repeat_for_triangles(list(self.triangle_set.keys()))
+
+    def repeat_for_triangles(self, values):
+        return numpy.repeat(values, self.counts)
+
+
+    def getNumberOfTriangles(self, shapeid):
+        return len(self.triangle_set[shapeid])
+
+    def generate_scene(self, colorproperty):
+        return generate_scene(self.triangle_set, colorproperty)
+
diff --git a/src/openalea/caribu/caributriangleset.py b/src/openalea/caribu/caributriangleset.py
deleted file mode 100644
index 7d36b44c..00000000
--- a/src/openalea/caribu/caributriangleset.py
+++ /dev/null
@@ -1,110 +0,0 @@
-from openalea.caribu.display import generate_scene
-
-import numpy
-
-class AbstractCaribuTriangleSet:
-    def __init__(self):
-        pass
-
-    def getBoundingBox(self):
-        pass
-
-    def triangle_areas(self):
-        pass
-
-    def getZmin(self):
-        pass
-
-    def __getitem__(self, shapeid):
-        """ Return all triangles of a shape """
-        pass
-
-    def keys(self):
-        pass
-
-    def values(self):
-        pass
-
-    def items(self):
-        pass
-
-    def allvalues(self, copied=False):
-        pass
-
-    def allids(self):
-        pass
-
-    def getNumberOfTriangles(self, shapeid):
-        pass
-
-    def generate_scene(self, colorproperty):
-        pass
-
-    def __len__(self):
-        raise NotImplemented()
-
-class CaribuTriangleSet(AbstractCaribuTriangleSet):
-    def __init__(self, pointtuplelistdict):
-        AbstractCaribuTriangleSet.__init__(self)
-        self._values = pointtuplelistdict
-        import itertools
-        self.allpoints = list(itertools.chain(*self._values.values()))
-        self.bbox = None
-
-    def getBoundingBox(self):
-        if self.bbox is None:
-            x, y, z = map(numpy.array, zip(*map(lambda x: zip(*x), self.allpoints)))
-            self.bbox = (x.min(), y.min(), z.min()), (x.max(), y.max(), z.max()) 
-        return self.bbox     
-        
-    def triangle_areas(self):
-        """ compute mean area of elementary triangles in the scene """
-
-        def _surf(triangle):
-            a, b, c = tuple(map(numpy.array, triangle))
-            x, y, z = numpy.cross(b - a, c - a).tolist()
-            return numpy.sqrt(x ** 2 + y ** 2 + z ** 2) / 2.0
-
-        return numpy.array(list(map(_surf, self.allpoints)))
-
-    def getZmin(self):
-        return self.getBoundingBox()[0][2]
-
-    def getZmax(self):
-        return self.getBoundingBox()[1][2]
-
-    def __getitem__(self, shapeid):
-        """ Return all triangles of a shape """
-        return self._values[shapeid]
-
-    def __len__(self):
-        return len(self._values)
-
-    def keys(self):
-        return self._values.keys()
-
-    def values(self):
-        return self._values.values()
-
-    def items(self):
-        return self._values.items()
-
-    def allvalues(self, copied=False):
-        from copy import copy
-        if copied : 
-            return copy(self.allpoints)
-        else:
-            return self.allpoints
-    
-    def allids(self):
-        return self.repeat_for_triangles(self._values.keys())
-
-    def repeat_for_triangles(self, values):
-        return [v for v,nb in zip(values,[len(v) for v in self._values.values()]) for j in range(nb)]
-
-    def getNumberOfTriangles(self, shapeid):
-        return len(self._values[shapeid])
-
-    def generate_scene(self, colorproperty):
-        return generate_scene(self._values, colorproperty)
-
diff --git a/src/openalea/caribu/data/__init__.py b/src/openalea/caribu/data/__init__.py
index e69de29b..875879eb 100644
--- a/src/openalea/caribu/data/__init__.py
+++ b/src/openalea/caribu/data/__init__.py
@@ -0,0 +1,10 @@
+from importlib.resources import files
+
+
+def get_path():
+    """Return a pathlib.Path representing the data directory."""
+    return files(__package__)
+
+
+def get(filename: str):
+    return files(__package__) / filename
\ No newline at end of file
diff --git a/test/test_caribu_triangle_set.py b/test/test_caribu_triangle_set.py
new file mode 100644
index 00000000..1c62703c
--- /dev/null
+++ b/test/test_caribu_triangle_set.py
@@ -0,0 +1,11 @@
+import openalea.caribu.data as data
+from openalea.caribu.caribu_triangle_set import CaribuTriangleSet
+
+
+
+def test_instantiation_from_file():
+    cts = CaribuTriangleSet.fromfile(data.get('f331s1_100plantes.can'))
+    plants = cts.scales['plant']
+    assert len(plants) == max(plants) == 100
+    areas = cts.sum_at_scale(cts.triangle_areas(), scale='plant')
+    assert all([0.95 < a / 0.36 < 1.05 for a in areas.values()])
\ No newline at end of file
diff --git a/test/test_legacy_using_caribu.py b/test/test_legacy_using_caribu.py
new file mode 100644
index 00000000..1bb036fa
--- /dev/null
+++ b/test/test_legacy_using_caribu.py
@@ -0,0 +1,79 @@
+"""
+Test caribu legacy test suite, using libcaribu.algos
+"""
+import pytest
+from numpy.testing import assert_array_equal
+from importlib.resources import files
+import openalea.libcaribu.io as lcio
+from openalea.caribu.caribu_shell import set_scene, caribu
+
+data_dir = files('openalea.libcaribu.data')
+
+
+@pytest.fixture(scope="module")  # tmp_path is kept between tests
+def caribu_test_scene(tmp_path_factory):
+    tmp_path = tmp_path_factory.mktemp("caribu_test_scene")
+    set_scene(tmp_path,
+              canopy=data_dir / "filterT.can",
+              pattern=data_dir / "filter.8",
+              lights=data_dir / "zenith.light",
+              sensors=data_dir / "filterT.sensor",
+              opts=data_dir / "par.opt")
+    return tmp_path
+
+
+def test_projection_non_toric_scene(caribu_test_scene):
+    out = caribu(caribu_test_scene, direct_only=True, toric=False)
+    assert 'par' in out
+    res, _ = out['par']
+    expected = lcio.read_results(data_dir / 'projection_non_toric_scene.vec0')
+    for field in res:
+        assert_array_equal(res[field], expected[field])
+
+
+def test_sensor_non_toric_scene(caribu_test_scene):
+    out = caribu(caribu_test_scene, direct_only=True, toric=False, sensors=True)
+    assert 'par' in out
+    res, mes = out['par']
+    expected_res = lcio.read_results(data_dir / 'projection_non_toric_scene.vec0')
+    expected_mes = lcio.read_measures(data_dir / 'sensor_non_toric_scene.dat')
+    for field in res:
+        assert_array_equal(res[field], expected_res[field])
+    for field in mes:
+        assert_array_equal(mes[field], expected_mes[field])
+
+
+def test_projection_toric_scene(caribu_test_scene):
+    out = caribu(caribu_test_scene, direct_only=True, toric=True)
+    assert 'par' in out
+    res, _ = out['par']
+    expected = lcio.read_results(data_dir / 'projection_toric_scene.vec0')
+    for field in res:
+        assert_array_equal(res[field], expected[field])
+
+
+def test_radiosity_non_toric_scene(caribu_test_scene):
+    out = caribu(caribu_test_scene, direct_only=False, mixed_radiosity=-1)
+    assert 'par' in out
+    res, _ = out['par']
+    expected = lcio.read_results(data_dir / 'radiosity_non_toric_scene.vec0')
+    for field in res:
+        assert_array_equal(res[field], expected[field])
+
+
+def test_projection_sail_toric_scene(caribu_test_scene):
+    out = caribu(caribu_test_scene, direct_only=False, mixed_radiosity=0, layers=6, height=21)
+    assert 'par' in out
+    res, _ = out['par']
+    expected = lcio.read_results(data_dir / 'projection_sail_toric_scene.vec0')
+    for field in res:
+        assert_array_equal(res[field], expected[field])
+
+
+def test_nested_radiosity_toric_scene(caribu_test_scene):
+    out = caribu(caribu_test_scene, direct_only=False, mixed_radiosity=1, layers=6, height=21)
+    assert 'par' in out
+    res, _ = out['par']
+    expected = lcio.read_results(data_dir / 'nested_radiosity_toric_scene.vec0')
+    for field in res:
+        assert_array_equal(res[field], expected[field])