Merge branch 'master' of https://github.com/mwoussen/lightvegemanager

src/lightvegemanager/CARIBUinputs.py

@@ -313,9 +313,9 @@ def create_caribu_legume_sensors(dxyz,
     # generate the sensors in plantGL format among the grid
     ID_capt = start_id
     dico_translat = {}
-    for ix in range(nxyz[0]):
+    for iz in range(nxyz[2] - skylayer):
         for iy in range(nxyz[1]):
-            for iz in range(nxyz[2] - skylayer):
+            for ix in range(nxyz[0]):
                 # translation vector
                 tx = ix * dxyz[0]
                 ty = iy * dxyz[1]

src/lightvegemanager/GLTF.py

@@ -0,0 +1,256 @@
+'''
+    GLTF
+    ***
+
+    Writes GLTF files from LightVegeManager geometry.
+    Used for visualisation
+'''
+import itertools
+import numpy
+import json
+import base64
+
+def GLTFtriangles(trimesh, var=[], varname=[], filename="output.gltf"):
+    """Writes GLTF 2.0 files from a triangulation mesh.
+
+    :param trimesh: triangles mesh aggregated by indice elements
+
+        .. code-block:: 
+
+            { id : [triangle1, triangle2, ...]}
+
+    :type trimesh: dict of list or list
+    :param var: list of physical values associated to each triangle.
+    :type var: list of list
+    :param varname: list of variable names
+    :type varname: lits of string
+    :param filename: name and path for the output file
+    :type filename: string
+    """    
+    if isinstance(trimesh, dict) :
+        list_triangles = list(itertools.chain(*[v for v in trimesh.values()]))
+    elif isinstance(trimesh, list) : 
+        list_triangles = trimesh
+
+    # Flatten vertices
+    points = []
+    for tr in list_triangles:
+        for p in tr:
+            points.extend(p)
+
+    if not points:
+        return
+
+    # Convert to numpy array (float32)
+    np_points = numpy.array(points, dtype=numpy.float32)
+
+    # Calculate min/max for accessor (required by GLTF spec for POSITION)
+    np_points_vec3 = np_points.reshape(-1, 3)
+    min_pos = np_points_vec3.min(axis=0).tolist()
+    max_pos = np_points_vec3.max(axis=0).tolist()
+
+    # Prepare buffers and accessors
+    buffers_list = []
+    buffer_views = []
+    accessors = []
+    attributes = {}
+
+    current_byte_offset = 0
+
+    # 1. POSITION
+    blob_pos = np_points.tobytes()
+    len_pos = len(blob_pos)
+    buffers_list.append(blob_pos)
+
+    buffer_views.append({
+        "buffer": 0,
+        "byteOffset": current_byte_offset,
+        "byteLength": len_pos,
+        "target": 34962 # ARRAY_BUFFER
+    })
+
+    accessors.append({
+        "bufferView": 0,
+        "byteOffset": 0,
+        "componentType": 5126, # FLOAT
+        "count": len(np_points_vec3),
+        "type": "VEC3",
+        "max": max_pos,
+        "min": min_pos
+    })
+    attributes["POSITION"] = 0
+    current_byte_offset += len_pos
+
+    # 2. Variables (Colors and Custom Attributes)
+    if var and len(var) > 0:
+        for i, v_list in enumerate(var):
+            if len(v_list) != len(list_triangles):
+                continue
+
+            values = numpy.array(v_list, dtype=numpy.float32)
+            # Expand to vertices (flat shading: same value for 3 vertices of triangle)
+            val_per_vertex = numpy.repeat(values, 3)
+
+            # A. Map first variable to COLOR_0 for visualization
+            if i == 0:
+                v_min, v_max = values.min(), values.max()
+                if v_max - v_min > 1e-9:
+                    norm_values = (val_per_vertex - v_min) / (v_max - v_min)
+                else:
+                    norm_values = numpy.zeros_like(val_per_vertex)
+
+                # Simple Colormap: Blue (low) -> Green -> Red (high)
+                colors = numpy.zeros((len(val_per_vertex), 3), dtype=numpy.float32)
+
+                # Blue to Green (0.0 to 0.5)
+                mask1 = norm_values <= 0.5
+                colors[mask1, 2] = 1.0 - 2.0 * norm_values[mask1] # Blue component
+                colors[mask1, 1] = 2.0 * norm_values[mask1]       # Green component
+
+                # Green to Red (0.5 to 1.0)
+                mask2 = norm_values > 0.5
+                colors[mask2, 1] = 1.0 - 2.0 * (norm_values[mask2] - 0.5) # Green component
+                colors[mask2, 0] = 2.0 * (norm_values[mask2] - 0.5)       # Red component
+
+                blob_col = colors.tobytes()
+                len_col = len(blob_col)
+                buffers_list.append(blob_col)
+
+                buffer_views.append({"buffer": 0, "byteOffset": current_byte_offset, "byteLength": len_col, "target": 34962})
+                accessors.append({"bufferView": len(buffer_views)-1, "byteOffset": 0, "componentType": 5126, "count": len(colors), "type": "VEC3", "max": colors.max(axis=0).tolist(), "min": colors.min(axis=0).tolist()})
+                attributes["COLOR_0"] = len(accessors) - 1
+                current_byte_offset += len_col
+
+            # B. Save raw data as Custom Attribute _VARNAME
+            if i < len(varname):
+                # Sanitize name (uppercase, alphanumeric)
+                safe_name = "".join(c if c.isalnum() else "_" for c in varname[i].upper())
+                attr_name = "_" + safe_name
+            else:
+                attr_name = f"_VAR_{i}"
+
+            blob_scalar = val_per_vertex.tobytes()
+            len_scalar = len(blob_scalar)
+            buffers_list.append(blob_scalar)
+
+            buffer_views.append({"buffer": 0, "byteOffset": current_byte_offset, "byteLength": len_scalar, "target": 34962})
+            accessors.append({"bufferView": len(buffer_views)-1, "byteOffset": 0, "componentType": 5126, "count": len(val_per_vertex), "type": "SCALAR", "max": [float(values.max())], "min": [float(values.min())]})
+            attributes[attr_name] = len(accessors) - 1
+            current_byte_offset += len_scalar
+
+            # C. Save as TEXCOORD_{i} for easy access in Game Engines (UE5, Unity)
+            # We pack the scalar value into the X component of a VEC2.
+            zeros = numpy.zeros_like(val_per_vertex)
+            vec2_values = numpy.column_stack((val_per_vertex, zeros)).astype(numpy.float32)
+
+            blob_uv = vec2_values.tobytes()
+            len_uv = len(blob_uv)
+            buffers_list.append(blob_uv)
+
+            buffer_views.append({"buffer": 0, "byteOffset": current_byte_offset, "byteLength": len_uv, "target": 34962})
+            accessors.append({
+                "bufferView": len(buffer_views)-1, 
+                "byteOffset": 0, 
+                "componentType": 5126, 
+                "count": len(vec2_values), 
+                "type": "VEC2", 
+                "max": vec2_values.max(axis=0).tolist(), 
+                "min": vec2_values.min(axis=0).tolist()
+            })
+            attributes[f"TEXCOORD_{i}"] = len(accessors) - 1
+            current_byte_offset += len_uv
+
+    # Combine all blobs
+    full_blob = b"".join(buffers_list)
+    uri = "data:application/octet-stream;base64," + base64.b64encode(full_blob).decode('utf-8')
+
+    gltf = {
+        "asset": {"version": "2.0", "generator": "LightVegeManager"},
+        "scene": 0,
+        "scenes": [{"nodes": [0]}],
+        "nodes": [{"mesh": 0}],
+        "meshes": [{
+            "primitives": [{
+                "attributes": attributes,
+                "mode": 4 # TRIANGLES
+            }]
+        }],
+        "buffers": [{"byteLength": len(full_blob), "uri": uri}],
+        "bufferViews": buffer_views,
+        "accessors": accessors
+    }
+
+    with open(filename, 'w') as f:
+        json.dump(gltf, f, indent=2)
+
+def VTKtoGLTF(vtk_filename, gltf_filename):
+    """Converts a VTK file (generated by VTKtriangles) to a GLTF file.
+
+    :param vtk_filename: path to the input VTK file
+    :type vtk_filename: string
+    :param gltf_filename: path to the output GLTF file
+    :type gltf_filename: string
+    """
+    with open(vtk_filename, 'r') as f:
+        lines = f.readlines()
+
+    iterator = iter(lines)
+
+    points = []
+    var = []
+    varname = []
+
+    try:
+        while True:
+            line = next(iterator).strip()
+            if line.startswith("POINTS"):
+                parts = line.split()
+                nb_points = int(parts[1])
+                for _ in range(nb_points):
+                    p_line = next(iterator).strip().split()
+                    points.append((float(p_line[0]), float(p_line[1]), float(p_line[2])))
+
+            elif line.startswith("CELLS"):
+                parts = line.split()
+                nb_cells = int(parts[1])
+                for _ in range(nb_cells):
+                    next(iterator)
+
+            elif line.startswith("CELL_TYPES"):
+                parts = line.split()
+                n_types = int(parts[1])
+                for _ in range(n_types):
+                    next(iterator)
+
+            elif line.startswith("FIELD"):
+                parts = line.split()
+                num_arrays = int(parts[2])
+
+                for _ in range(num_arrays):
+                    # Skip empty lines to find next header
+                    header = next(iterator).strip()
+                    while not header:
+                        header = next(iterator).strip()
+
+                    h_parts = header.split()
+                    name = h_parts[0]
+                    num_tuples = int(h_parts[2])
+
+                    values = []
+                    for _ in range(num_tuples):
+                        val_line = next(iterator).strip()
+                        values.append(float(val_line))
+
+                    varname.append(name)
+                    var.append(values)
+
+    except StopIteration:
+        pass
+
+    # Reconstruct trimesh from points (assuming 3 points per triangle as per VTKtriangles)
+    trimesh = []
+    for i in range(0, len(points), 3):
+        if i + 2 < len(points):
+            trimesh.append([points[i], points[i+1], points[i+2]])
+
+    GLTFtriangles(trimesh, var, varname, gltf_filename)

src/lightvegemanager/LVM.py

@@ -96,6 +96,7 @@
     compute_area_max,
     compute_minmax_coord,
     compute_trilenght_max,
+    compute_total_area
 )
 from lightvegemanager.VTK import VTKtriangles
 from lightvegemanager.defaultvalues import default_LightVegeManager_inputs
@@ -245,6 +246,7 @@ def build(self, geometry={}, global_scene_tesselate_level=0):
             )
 
         self.__areamax = compute_area_max(self.__complete_trimesh)
+        self.__totalarea = compute_total_area(self.__complete_trimesh)
 
         # global tesselation of triangulation
         if self.__matching_ids and global_scene_tesselate_level > 0:
@@ -294,17 +296,34 @@ def build(self, geometry={}, global_scene_tesselate_level=0):
                     arg = (dxyz, nxyz, orig, self.__pmax, self.__complete_trimesh, self.__matching_ids, None, True)
                     sensors_caribu, sensors_plantgl, Pmax_capt = create_caribu_legume_sensors(*arg)
 
-                    id = [-1]
+
                     self.__sensors_plantgl = sensors_plantgl
                     self.__sensors_caribu = sensors_caribu
-                    self.__nb0 = reduce_layers_from_trimesh(
-                        self.__complete_trimesh,
-                        self.__pmax,
-                        self.__lightmodel_parameters["sensors"][1],
-                        self.__lightmodel_parameters["sensors"][2],
-                        self.__matching_ids,
-                        id,
-                    )
+
+
+                    ###je ne comprends pas la logique derrière ce calcul. 
+                    # comme id = -1 n'est pas dans matching_id, on écrase toutes les couches 
+                    # id = [-1]
+                    # self.__nb0 = reduce_layers_from_trimesh(
+                    #     self.__complete_trimesh,
+                    #     self.__pmax,
+                    #     self.__lightmodel_parameters["sensors"][1],
+                    #     self.__lightmodel_parameters["sensors"][2],
+                    #     self.__matching_ids,
+                    #     id,
+                    # )
+
+                    #je remplace de cette manière : 
+                    zmax_sensors = Pmax_capt[2]
+
+                    skylayer = zmax_sensors // dxyz[2]
+
+                    if skylayer < nxyz[2] and zmax_sensors > 0:
+                        self.__nb0 = int(nxyz[2] - 1 - skylayer)
+                    else:
+                        self.__nb0 = int(nxyz[2] - 1)
+
+
                 elif isinstance(self.__lightmodel_parameters["sensors"], dict):
                     start_id = 0
                     sensors_plantgl = pgl.Scene()
@@ -463,6 +482,8 @@ def run(self, energy=0.0, day=0, hour=0, parunit="micromol.m-2.s-1", truesolarti
         """
         self.__energy = energy
 
+        epsilon = 1e-14
+
         ## RATP ##
         if self.__lightmodel == "ratp":
             from lightvegemanager.RATPinputs import RATP_vegetation, RATP_meteo
@@ -594,9 +615,23 @@ def run(self, energy=0.0, day=0, hour=0, parunit="micromol.m-2.s-1", truesolarti
                         scene_unit=self.__main_unit,
                         pattern=self.__geometry["domain"],
                         soil_mesh=self.__lightmodel_parameters["soil mesh"],
+                        #z_soil=-epsilon,
                         debug=debug,
                     )
 
+
+                    ##### test de relever la couche inférieure de epsilon
+                    epsilon = 0.001 #0.000001
+
+                    if sensors_caribu is not None :
+                        for cle, polygones in sensors_caribu.items():
+                            sensors_caribu[cle] = [
+                                [(x, y, epsilon if z == 0.0 else z) for x, y, z in polygone]
+                                for polygone in polygones
+                                ]
+                        self.__sensors_caribu = sensors_caribu
+                    #####
+
                     # Runs CARIBU
                     arg = [
                         c_scene,
@@ -1081,7 +1116,7 @@ def to_VTK(
         :raises AttributeError: it needs to have a grid of virtual sensors
         """
         if lighting:
-            # deactivate lighting if if no lighting results
+            # deactivate lighting if no lighting results
             if self.__lightmodel == "ratp" and (not hasattr(self, "_LightVegeManager__voxels_outputs")):
                 print("--- VTK:  No light data, run the simulation")
                 lighting = False
@@ -1362,6 +1397,20 @@ def maxtrianglearea(self):
             return self.__areamax
         except AttributeError:
             return 0.0
+
+    @property
+    def totaltrianglearea(self):
+        """Returns the total area of all triangles of triangles mesh
+        Computed in :meth:`build`
+
+        :return: area of all triangles
+        :rtype: float
+        """
+        try:
+            return self.__totalarea
+        except AttributeError:
+            return 0.0
+
 
     @property
     def legume_empty_layers(self):