add new caribu function based on libcaribu

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
 

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: [email protected]
-  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.

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)
+