relief/STEP.py

# -*- coding: utf-8 -*-
"""
transformeSTEP.py

Lit un fichier .rcl (grille d'altitudes) et génère un solide au format STEP.

Format attendu du fichier .rcl :
    ligne 1 : ignorée (non utilisée)
    ligne 2 : X  -> largeur de la grille (nombre de colonnes)
    ligne 3 : Y  -> hauteur de la grille (nombre de lignes)
    suite   : X*Y valeurs flottantes (altitudes), écrites comme une liste
              Python "[v1, v2, v3, ...]" pouvant s'étendre sur plusieurs lignes.

Dépendances :
    pip install numpy
    pythonocc-core (conseillé via conda) :
        conda install -c conda-forge pythonocc-core

    Si vous n'avez pas conda, cherchez une roue (wheel) pythonocc-core
    compatible avec votre OS/version de Python, l'installation via pip seul
    ne fonctionne pas sur toutes les plateformes.
"""

import re
import numpy as np

from OCC.Core.gp import gp_Pnt
from OCC.Core.BRepBuilderAPI import (
    BRepBuilderAPI_MakePolygon,
    BRepBuilderAPI_MakeFace,
    BRepBuilderAPI_Sewing,
    BRepBuilderAPI_MakeSolid,
)
from OCC.Core.TopoDS import topods, TopoDS_Shell
from OCC.Core.BRep import BRep_Builder
from OCC.Core.TopExp import TopExp_Explorer
from OCC.Core.TopAbs import TopAbs_SHELL
from OCC.Core.STEPControl import STEPControl_Writer, STEPControl_AsIs
from OCC.Core.IFSelect import IFSelect_RetDone
from OCC.Core.BRepLib import breplib


# ----------------------------------------------------------------------
# 1) Lecture / extraction des informations du fichier .rcl
# ----------------------------------------------------------------------
def lire_et_extraire_infos(nomDonne):
    """
    Lit un fichier .rcl et retourne (X, Y, grille) où grille est un
    tableau numpy de forme (Y, X) contenant les altitudes.
    """
    with open(nomDonne, "r") as f:
        lignes = f.readlines()

    if len(lignes) < 4:
        raise ValueError(f"Fichier {nomDonne} trop court / mal formé.")

    # ligne 1 (index 0) : non utilisée
    X = int(lignes[1].strip())
    Y = int(lignes[2].strip())

    # Le reste du fichier contient la liste des altitudes, potentiellement
    # répartie sur plusieurs lignes -> on regroupe tout et on extrait
    # tous les nombres flottants (gère aussi la notation scientifique).
    reste = "".join(lignes[3:])
    valeurs_str = re.findall(r"-?\d+(?:\.\d+)?(?:[eE][-+]?\d+)?", reste)
    valeurs = [float(v) for v in valeurs_str]

    attendu = X * Y
    if len(valeurs) < attendu:
        raise ValueError(
            f"Nombre de valeurs insuffisant dans {nomDonne} : "
            f"{len(valeurs)} lues, {attendu} attendues (X={X}, Y={Y})."
        )
    valeurs = valeurs[:attendu]

    grille = np.array(valeurs, dtype=float).reshape(Y, X)
    return X, Y, grille


# ----------------------------------------------------------------------
# 2) Construction du solide (mur + toit + plancher) à partir de la grille
# ----------------------------------------------------------------------
def _face_triangle(p1, p2, p3):
    """Crée une face plane triangulaire à partir de 3 gp_Pnt."""
    poly = BRepBuilderAPI_MakePolygon()
    poly.Add(p1)
    poly.Add(p2)
    poly.Add(p3)
    poly.Close()
    return BRepBuilderAPI_MakeFace(poly.Wire()).Face()


def _face_quad(p1, p2, p3, p4):
    """Crée une face plane quadrangulaire à partir de 4 gp_Pnt (ordre = contour)."""
    poly = BRepBuilderAPI_MakePolygon()
    poly.Add(p1)
    poly.Add(p2)
    poly.Add(p3)
    poly.Add(p4)
    poly.Close()
    return BRepBuilderAPI_MakeFace(poly.Wire()).Face()


def construire_solide_depuis_grille(grille, echelle_xy=1.0, echelle_z=1.0,
                                     z_base=None, pas=1):
    """
    Construit un solide B-Rep (OpenCASCADE) à partir d'une grille d'altitudes.

    grille     : tableau numpy (Y, X) d'altitudes
    echelle_xy : facteur d'échelle appliqué aux coordonnées X/Y (pas de grille -> unités)
    echelle_z  : facteur d'échelle appliqué aux altitudes
    z_base     : altitude du plancher (par défaut : min(grille) - 10% de l'amplitude)
    pas        : sous-échantillonnage de la grille (1 = pleine résolution).
                 Pour de grandes grilles (ex. 1000x1000), un solide B-Rep
                 avec un million de facettes est très lourd : augmentez "pas"
                 (ex. 5, 10, 20 ...) pour réduire le nombre de faces.

    Retourne un TopoDS_Solid.
    """
    Y, X = grille.shape

    # Sous-échantillonnage éventuel pour garder un nombre de facettes raisonnable
    xs = list(range(0, X, pas))
    if xs[-1] != X - 1:
        xs.append(X - 1)
    ys = list(range(0, Y, pas))
    if ys[-1] != Y - 1:
        ys.append(Y - 1)

    sous_grille = grille[np.ix_(ys, xs)]
    ny, nx = sous_grille.shape

    if z_base is None:
        zmin, zmax = float(sous_grille.min()), float(sous_grille.max())
        amplitude = max(zmax - zmin, 1.0)
        z_base = zmin - 0.1 * amplitude

    # Points du toit (top[j][i]) et points du plancher directement en dessous
    top = [[gp_Pnt(i * echelle_xy, j * echelle_xy, sous_grille[j, i] * echelle_z)
            for i in range(nx)] for j in range(ny)]
    bas = [[gp_Pnt(i * echelle_xy, j * echelle_xy, z_base)
            for i in range(nx)] for j in range(ny)]

    faces = []

    # --- Toit : chaque case de la grille -> 2 triangles ---
    for j in range(ny - 1):
        for i in range(nx - 1):
            p1, p2 = top[j][i], top[j][i + 1]
            p3, p4 = top[j + 1][i + 1], top[j + 1][i]
            faces.append(_face_triangle(p1, p2, p3))
            faces.append(_face_triangle(p1, p3, p4))

    # --- Murs latéraux (4 côtés du rectangle) ---
    def mur(p_top_a, p_top_b, p_bas_a, p_bas_b):
        faces.append(_face_triangle(p_top_a, p_top_b, p_bas_b))
        faces.append(_face_triangle(p_top_a, p_bas_b, p_bas_a))

    # côté j = 0
    for i in range(nx - 1):
        mur(top[0][i], top[0][i + 1], bas[0][i], bas[0][i + 1])
    # côté j = ny-1
    for i in range(nx - 1):
        mur(top[ny - 1][i], top[ny - 1][i + 1], bas[ny - 1][i], bas[ny - 1][i + 1])
    # côté i = 0
    for j in range(ny - 1):
        mur(top[j][0], top[j + 1][0], bas[j][0], bas[j + 1][0])
    # côté i = nx-1
    for j in range(ny - 1):
        mur(top[j][nx - 1], top[j + 1][nx - 1], bas[j][nx - 1], bas[j + 1][nx - 1])

    # --- Plancher : un unique rectangle plan à z_base ---
    faces.append(_face_quad(bas[0][0], bas[0][nx - 1], bas[ny - 1][nx - 1], bas[ny - 1][0]))

    # --- Couture des faces en une coque fermée ---
    sewing = BRepBuilderAPI_Sewing(1e-6)
    for f in faces:
        sewing.Add(f)
    sewing.Perform()
    coque = sewing.SewedShape()

    # Récupère le TopoDS_Shell résultant de la couture
    explorer = TopExp_Explorer(coque, TopAbs_SHELL)
    if not explorer.More():
        raise RuntimeError("Échec de la couture : aucune coque produite.")
    shell = topods.Shell(explorer.Current())

    # --- Construction du solide à partir de la coque fermée ---
    make_solid = BRepBuilderAPI_MakeSolid(shell)
    if not make_solid.IsDone():
        raise RuntimeError("Échec de la création du solide à partir de la coque.")
    solide = make_solid.Solid()

    # S'assure que les normales sont correctement orientées vers l'extérieur
    breplib.OrientClosedSolid(solide)

    return solide


# ----------------------------------------------------------------------
# 3) Export STEP
# ----------------------------------------------------------------------
def exporter_step(solide, chemin_fichier):
    writer = STEPControl_Writer()
    writer.Transfer(solide, STEPControl_AsIs)
    status = writer.Write(chemin_fichier)
    if status != IFSelect_RetDone:
        raise RuntimeError(f"Échec de l'écriture du fichier STEP : {chemin_fichier}")


# ----------------------------------------------------------------------
# 4) Fonction principale, sur le modèle de transformeSTL
# ----------------------------------------------------------------------
def transformeSTEP(nom, echelle_xy=1.0, echelle_z=1.0, z_base=None, pas=1):
    """
    Lit nom + ".rcl", construit un solide à partir de la grille d'altitudes
    et écrit nom + ".step".

    pas : sous-échantillonnage de la grille (voir construire_solide_depuis_grille).
          Pour une grille 1000x1000, commencez par exemple avec pas=10 ou plus
          pour garder un temps de calcul raisonnable, puis diminuez si besoin.
    """
    nomDonne = nom + ".rcl"
    X, Y, donne = lire_et_extraire_infos(nomDonne)

    solide = construire_solide_depuis_grille(
        donne, echelle_xy=echelle_xy, echelle_z=echelle_z,
        z_base=z_base, pas=pas
    )

    nomSortie = nom + ".step"
    exporter_step(solide, nomSortie)
    print(f"Fichier STEP généré : {nomSortie}")
    return nomSortie