Initial commit

LICENSE

+Academic Free License ("AFL") v. 3.0
+
+This Academic Free License (the "License") applies to any original work of
+authorship (the "Original Work") whose owner (the "Licensor") has placed the
+following licensing notice adjacent to the copyright notice for the Original
+Work:
+
+  Licensed under the Academic Free License version 3.0
+
+1) Grant of Copyright License. Licensor grants You a worldwide, royalty-free,
+non-exclusive, sublicensable license, for the duration of the copyright, to do
+the following:
+
+  a) to reproduce the Original Work in copies, either alone or as part of a
+  collective work;
+
+  b) to translate, adapt, alter, transform, modify, or arrange the Original
+  Work, thereby creating derivative works ("Derivative Works") based upon the
+  Original Work;
+
+  c) to distribute or communicate copies of the Original Work and Derivative
+  Works to the public, under any license of your choice that does not
+  contradict the terms and conditions, including Licensor's reserved rights
+  and remedies, in this Academic Free License;
+
+  d) to perform the Original Work publicly; and
+
+  e) to display the Original Work publicly.
+
+2) Grant of Patent License. Licensor grants You a worldwide, royalty-free,
+non-exclusive, sublicensable license, under patent claims owned or controlled
+by the Licensor that are embodied in the Original Work as furnished by the
+Licensor, for the duration of the patents, to make, use, sell, offer for sale,
+have made, and import the Original Work and Derivative Works.
+
+3) Grant of Source Code License. The term "Source Code" means the preferred
+form of the Original Work for making modifications to it and all available
+documentation describing how to modify the Original Work. Licensor agrees to
+provide a machine-readable copy of the Source Code of the Original Work along
+with each copy of the Original Work that Licensor distributes. Licensor
+reserves the right to satisfy this obligation by placing a machine-readable
+copy of the Source Code in an information repository reasonably calculated to
+permit inexpensive and convenient access by You for as long as Licensor
+continues to distribute the Original Work.
+
+4) Exclusions From License Grant. Neither the names of Licensor, nor the names
+of any contributors to the Original Work, nor any of their trademarks or
+service marks, may be used to endorse or promote products derived from this
+Original Work without express prior permission of the Licensor. Except as
+expressly stated herein, nothing in this License grants any license to
+Licensor's trademarks, copyrights, patents, trade secrets or any other
+intellectual property. No patent license is granted to make, use, sell, offer
+for sale, have made, or import embodiments of any patent claims other than the
+licensed claims defined in Section 2. No license is granted to the trademarks
+of Licensor even if such marks are included in the Original Work. Nothing in
+this License shall be interpreted to prohibit Licensor from licensing under
+terms different from this License any Original Work that Licensor otherwise
+would have a right to license.
+
+5) External Deployment. The term "External Deployment" means the use,
+distribution, or communication of the Original Work or Derivative Works in any
+way such that the Original Work or Derivative Works may be used by anyone
+other than You, whether those works are distributed or communicated to those
+persons or made available as an application intended for use over a network.
+As an express condition for the grants of license hereunder, You must treat
+any External Deployment by You of the Original Work or a Derivative Work as a
+distribution under section 1(c).
+
+6) Attribution Rights. You must retain, in the Source Code of any Derivative
+Works that You create, all copyright, patent, or trademark notices from the
+Source Code of the Original Work, as well as any notices of licensing and any
+descriptive text identified therein as an "Attribution Notice." You must cause
+the Source Code for any Derivative Works that You create to carry a prominent
+Attribution Notice reasonably calculated to inform recipients that You have
+modified the Original Work.
+
+7) Warranty of Provenance and Disclaimer of Warranty. Licensor warrants that
+the copyright in and to the Original Work and the patent rights granted herein
+by Licensor are owned by the Licensor or are sublicensed to You under the
+terms of this License with the permission of the contributor(s) of those
+copyrights and patent rights. Except as expressly stated in the immediately
+preceding sentence, the Original Work is provided under this License on an "AS
+IS" BASIS and WITHOUT WARRANTY, either express or implied, including, without
+limitation, the warranties of non-infringement, merchantability or fitness for
+a particular purpose. THE ENTIRE RISK AS TO THE QUALITY OF THE ORIGINAL WORK
+IS WITH YOU. This DISCLAIMER OF WARRANTY constitutes an essential part of this
+License. No license to the Original Work is granted by this License except
+under this disclaimer.
+
+8) Limitation of Liability. Under no circumstances and under no legal theory,
+whether in tort (including negligence), contract, or otherwise, shall the
+Licensor be liable to anyone for any indirect, special, incidental, or
+consequential damages of any character arising as a result of this License or
+the use of the Original Work including, without limitation, damages for loss
+of goodwill, work stoppage, computer failure or malfunction, or any and all
+other commercial damages or losses. This limitation of liability shall not
+apply to the extent applicable law prohibits such limitation.
+
+9) Acceptance and Termination. If, at any time, You expressly assented to this
+License, that assent indicates your clear and irrevocable acceptance of this
+License and all of its terms and conditions. If You distribute or communicate
+copies of the Original Work or a Derivative Work, You must make a reasonable
+effort under the circumstances to obtain the express assent of recipients to
+the terms of this License. This License conditions your rights to undertake
+the activities listed in Section 1, including your right to create Derivative
+Works based upon the Original Work, and doing so without honoring these terms
+and conditions is prohibited by copyright law and international treaty.
+Nothing in this License is intended to affect copyright exceptions and
+limitations (including "fair use" or "fair dealing"). This License shall
+terminate immediately and You may no longer exercise any of the rights granted
+to You by this License upon your failure to honor the conditions in Section
+1(c).
+
+10) Termination for Patent Action. This License shall terminate automatically
+and You may no longer exercise any of the rights granted to You by this
+License as of the date You commence an action, including a cross-claim or
+counterclaim, against Licensor or any licensee alleging that the Original Work
+infringes a patent. This termination provision shall not apply for an action
+alleging patent infringement by combinations of the Original Work with other
+software or hardware.
+
+11) Jurisdiction, Venue and Governing Law. Any action or suit relating to this
+License may be brought only in the courts of a jurisdiction wherein the
+Licensor resides or in which Licensor conducts its primary business, and under
+the laws of that jurisdiction excluding its conflict-of-law provisions. The
+application of the United Nations Convention on Contracts for the
+International Sale of Goods is expressly excluded. Any use of the Original
+Work outside the scope of this License or after its termination shall be
+subject to the requirements and penalties of copyright or patent law in the
+appropriate jurisdiction. This section shall survive the termination of this
+License.
+
+12) Attorneys' Fees. In any action to enforce the terms of this License or
+seeking damages relating thereto, the prevailing party shall be entitled to
+recover its costs and expenses, including, without limitation, reasonable
+attorneys' fees and costs incurred in connection with such action, including
+any appeal of such action. This section shall survive the termination of this
+License.
+
+13) Miscellaneous. If any provision of this License is held to be
+unenforceable, such provision shall be reformed only to the extent necessary
+to make it enforceable.
+
+14) Definition of "You" in This License. "You" throughout this License,
+whether in upper or lower case, means an individual or a legal entity
+exercising rights under, and complying with all of the terms of, this License.
+For legal entities, "You" includes any entity that controls, is controlled by,
+or is under common control with you. For purposes of this definition,
+"control" means (i) the power, direct or indirect, to cause the direction or
+management of such entity, whether by contract or otherwise, or (ii) ownership
+of fifty percent (50%) or more of the outstanding shares, or (iii) beneficial
+ownership of such entity.
+
+15) Right to Use. You may use the Original Work in all ways not otherwise
+restricted or conditioned by this License or by law, and Licensor promises not
+to interfere with or be responsible for such uses by You.
+
+16) Modification of This License. This License is Copyright © 2005 Lawrence
+Rosen. Permission is granted to copy, distribute, or communicate this License
+without modification. Nothing in this License permits You to modify this
+License as applied to the Original Work or to Derivative Works. However, You
+may modify the text of this License and copy, distribute or communicate your
+modified version (the "Modified License") and apply it to other original works
+of authorship subject to the following conditions: (i) You may not indicate in
+any way that your Modified License is the "Academic Free License" or "AFL" and
+you may not use those names in the name of your Modified License; (ii) You
+must replace the notice specified in the first paragraph above with the notice
+"Licensed under <insert your license name here>" or with a notice of your own
+that is not confusingly similar to the notice in this License; and (iii) You
+may not claim that your original works are open source software unless your
+Modified License has been approved by Open Source Initiative (OSI) and You
+comply with its license review and certification process.

pyproject.toml

+[build-system]
+requires = ["setuptools>=61.0", "numpy", "matplotlib", "gmsh", "h5py"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "qosm-tools"
+version = "0.0.1"
+authors = [
+  { name="Gregory Gaudin", email="gregory.gaudin@imt-atlantique.fr" },
+]
+description = "Python tools for Quasi-Optical System Modelling"
+readme = "README.md"
+requires-python = ">=3.8"
+classifiers = [
+    "Programming Language :: Python :: 3",
+    "Operating System :: OS Independent",
+    "License :: OSI Approved :: Academic Free License v. 3.0",
+]
+
+[project.urls]
+Homepage = "https://gitlab.imt-atlantique.fr/quasi-optical-system-modelling/"
+Issues = "https://gitlab.imt-atlantique.fr/quasi-optical-system-modelling/qosm-core/-/issues"
\ No newline at end of file

src/qosm/__init__.py

+from .items import Mesh, MshMesh, ObjMesh, Object, SlabMesh, StepMesh, Triangle, utils
+from .propagation import SlabPropagation, utils
+from .sources import GaussianBeam, Horn, VirtualSource, utils
+from .utils import Field, Pose
+from .waves import VAGB
\ No newline at end of file

src/qosm/items/MshMesh.py

+import gmsh
+from numpy import pi, arcsin, cos
+
+from qosm.items.Mesh import Mesh
+from qosm.utils.Pose import Frame, Vector
+from qosm.sources import Source
+
+
+class MshMesh(Mesh):
+    def __init__(self, f_obj_ref: Frame = Frame(), source: Source = None):
+        super().__init__(f_obj_ref)
+
+    def load_mesh(self, element_size: float, scale: float = 1., view=False, show_vectors=True):
+        gmsh.option.setNumber("Geometry.Tolerance", 1e-8)
+        gmsh.option.setNumber("Mesh.Algorithm", 6)
+        gmsh.option.setNumber('Mesh.MeshSizeMin', element_size)
+        gmsh.option.setNumber('Mesh.MeshSizeMax', element_size)
+
+        gmsh.model.occ.removeAllDuplicates()
+        gmsh.model.occ.synchronize()
+        gmsh.model.mesh.generate(2)
+
+        gmsh.model.mesh.removeDuplicateNodes()
+
+        self.pre_build(scale, False, None, view, show_vectors, Vector(), False)

src/qosm/items/ObjMesh.py

+import numpy as np
+
+from qosm.items.Mesh import Mesh
+from qosm.utils.Pose import Frame, Vector, Quaternion
+
+ELEMENT_TYPES = {
+    "line": 2,
+    "line3": 3,
+    "triangle": 3,
+    "quadrangle": 4,
+    "tetrahedron": 4,
+    "hexahedron": 8,
+    "prism": 6,
+    "pyramid": 5,
+}
+
+
+def parse_vertices(vertices_strings):
+    v = np.zeros((len(vertices_strings), 3))
+    i = 0
+    for vertice_strings in vertices_strings:
+        v[i, :] = vertice_strings.split()
+        i += 1
+
+    return v
+
+
+def parse_faces(faces_strings):
+    faces_v = np.zeros((len(faces_strings), 3), dtype=int)
+    faces_n = np.zeros((len(faces_strings), 3), dtype=int)
+    i = 0
+    for face_string in faces_strings:
+        j = 0
+        for e in face_string.split():
+            if '//' in e:  # f v//vn
+                v_idx = int(e.split('//')[0]) - 1
+                n_idx = int(e.split('//')[1]) - 1
+            elif '/' in e:
+                if len(e.split('/')) == 2:  # f v/vt
+                    v_idx = int(e.split('/')[0]) - 1
+                    n_idx = int(e.split('/')[0]) - 1
+                else:  # f v/vt/vn
+                    v_idx = int(e.split('/')[0]) - 1
+                    n_idx = int(e.split('/')[2]) - 1
+            else:  # f v v v
+                v_idx = int(e.split()[0]) - 1
+                n_idx = int(e.split()[0]) - 1
+            faces_v[i, j] = v_idx
+            faces_n[i, j] = n_idx
+            j += 1
+        i += 1
+    return faces_v, faces_n
+
+
+def triangulate(faces: list):
+    new_faces = []
+    for face in faces:
+        elements = np.array(face.split())
+        if len(elements) == 3:
+            new_faces.append(face)
+            continue
+        new_faces.append(' '.join(elements[[0, 1, 2]]))
+        new_faces.append(' '.join(elements[[2, 3, 0]]))
+    return new_faces
+
+
+class ObjMesh(Mesh):
+    def __init__(self, f_obj_glo: Frame):
+        super(ObjMesh, self).__init__(f_obj_glo)
+        self.frame = f_obj_glo
+
+    def load(self, filepath: str, offset: Vector = Vector(), rot: Quaternion = Quaternion(), center: bool = False):
+        num_faces_3_edges = 0
+        num_faces_4_edges = 0
+
+        with open(filepath) as f:
+            lines = f.readlines()
+            vStrings = [x.replace('\n', '').strip('v') for x in lines if x.startswith('v ')]
+            vertices = parse_vertices(vStrings).view(Vector)
+
+            barycentre = Vector()
+
+            if center:
+                barycentre = np.mean(vertices, axis=0).reshape((-1, 3))
+            vertices -= barycentre
+
+            self.vertices = rot.rotate(vertices)
+
+            self.vertices += offset
+
+            vnStrings = [x.strip('vn') for x in lines if x.startswith('vn')]
+            """if not vnStrings:  # if There is no normal vectors in the obj file then compute them
+                normals = fillNormalsArray(len(vStrings))
+            else:"""
+            self.normals = parse_vertices(vnStrings)
+            self.normals = rot.rotate(self.normals.view(Vector))
+
+            faces = [x.replace('\n', '').strip('f') for x in lines if x.startswith('f')]
+
+            for face in faces:
+                if len(face.split()) == 3:
+                    num_faces_3_edges += 1
+                elif len(face.split()) == 4:
+                    num_faces_4_edges += 1
+                else:
+                    raise Exception('Faces with more than 4 edges are not supported')
+
+            """print("File:", filepath, "\nTotal number of faces:", len(faces),
+                  "\nNumber of faces with 3 vertices:", num_faces_3_edges,
+                  "\nNumber of faces with 4 vertices:", num_faces_4_edges,
+                  "\nTotal number of vertices:", len(self.vertices),
+                  "\nTotal number of normals:", len(self.normals))"""
+            if num_faces_4_edges > 0:
+                faces = triangulate(faces)
+
+            self.triangles, triangles_n = parse_faces(faces)
+
+        self.tri_vertices = self.vertices[self.triangles]
+        self.tri_normals = self.normals[triangles_n]
+        self.tri_tangents = self.tri_normals*0.
+        self.tri_curvatures = np.zeros((self.triangles.shape[0], 3))
+
+    def pre_build(self, scale=1e-3, create_obj=None, progress=None, view=False, show_vectors=True,
+                  offset: Vector = Vector()):
+        pass
\ No newline at end of file

src/qosm/items/Object.py

+from copy import deepcopy
+import numpy as np
+
+from qosm_core import Ray, Medium, Dioptre, Vec3, Frame
+
+class Item:
+    def __init__(self, f_obj_ref: Frame):
+        self.dioptre = None
+        self.frame = f_obj_ref
+
+    def copy(self):
+        return deepcopy(self)
+
+    def set_dioptre(self, material1: Medium, material2: Medium) -> None:
+        self.dioptre = Dioptre(material1, material2)
+
+
+class Object(Item):
+    def __init__(self, f_obj_glo: Frame):
+        super().__init__(f_obj_glo)
+        self.surfaces = []
+
+    def copy(self):
+        return deepcopy(self)
+
+    def add_surface(self, name: str, stype: int, pos: Vec3, params: dict) -> None:
+        if stype == 0:
+            surface = Disk(name, pos, params['radius'])
+        elif stype == 1:
+            surface = Sphere(name, pos, params['radius'], params['zlim'])
+        elif stype == 2:
+            surface = Cylinder(name, pos, params['radius'], params['h'])
+        elif stype == 3:
+            surface = Plane(name, pos, params['n'], params['w'], params['h'])
+        else:
+            return
+        self.surfaces.append(surface)
+
+    """def intersect(self, ray_src: BeamRay, intersection: Intersection) -> None:
+        # change frame GLO->LOC
+        ray_loc = ray_src.change_frame(self.frame)
+
+        for surf in self.surfaces:
+            (status, n, t, t0, N) = surf.intersect(ray_loc)
+            if status:
+                if Intersection.acne < t0 < intersection.t0:
+                    intersection.set(
+                        dist=t0,
+                        normal=self.frame.rot_to_ref(n),
+                        tangent=self.frame.rot_to_ref(t),
+                        ray=ray_src,
+                        curvature=surf.C,
+                        dioptre=self.dioptre
+                    )
+                    ray_src.z_range[1] = t0"""
+
+
+class Surface:
+    def __init__(self,
+                 surface_name: str,
+                 stype: int,
+                 pos: Vec3,
+                 curv_matrix=None):
+        self.name = surface_name
+        self.type = stype
+        self.r_surf_obj = pos
+
+        # curvature matrix
+        if curv_matrix is not None:
+            self.C = curv_matrix
+        else:
+            # flat surface
+            self.C = np.zeros((2, 2))
+
+    def intersect(self, ray: Ray) -> (bool, Vec3, Vec3, float):
+        return False, Vec3(), Vec3(), np.inf
+
+    def copy(self):
+        return deepcopy(self)
+
+
+class Disk(Surface):
+    def __init__(self, surface_name, pos, disk_radius):
+        super().__init__(surface_name, 0, pos)
+        self._disk_radius = disk_radius
+
+    def intersect(self, ray: Ray):
+        # only non-tilted plane so far
+
+        # normal vector
+        # local frame: can use the z-cpnt to check normal
+        t = Vec3(1, 0, 0)
+        n = Vec3(0, 0, -1)
+
+        # find t0 distance between ray origin and intersection point
+        a0 = ray.di * 1e-6
+        P = ray.ori + a0
+        L = self.r_surf_obj - P
+        t0 = (L.dot(n)) / (ray.dir.dot(n)) + 1e-6
+
+        I = ray.ori + ray.dir * t0
+
+        r = (I - self.r_surf_obj).norm()
+
+        status = t0 > 0 and r <= self._disk_radius
+
+        return status, n, t, t0[0], 1
+
+
+class Plane(Surface):
+    def __init__(self, surface_name, pos, n, w, h):
+        super().__init__(surface_name, 0, pos)
+        self._w = w
+        self._h = h
+        self._n = n
+
+    def intersect(self, ray: Ray):
+        # only non-tilted plane so far
+        # normal vector
+        # local frame: can use the z-cpnt to check normal
+        # n = self._n
+        t = Vec3(1, 0, 0)
+        n = Vec3(0, 0, -1)
+
+        # find t0 distance between ray origin and intersection point
+        a0 = ray.dir * 1e-6
+        P = ray.ori + a0
+        L = self.r_surf_obj - P
+        t0 = (L.dot(n)) / (ray.dir.dot(n)) + 1e-6
+
+        status = t0 > 0
+
+        return status, n, t, t0[0], 1
+
+
+class Cylinder(Surface):
+    def __init__(self, surface_name, pos, cyl_radius, h):
+        super().__init__(surface_name, 2, pos, curv_matrix=np.eye(2)/cyl_radius)
+        self._hd = h/2
+        self._r2 = np.abs(cyl_radius)**2
+
+    def intersect(self, ray: Ray):
+        a0 = ray.dir * 1e-6
+        P = ray.ori + a0
+        p = np.reshape(P[0, 0:2], (2, 1))
+        d = np.reshape(ray.dir[0, 0:2], (2, 1))
+        a = d[0]**2 + d[1]**2
+        b = p[0]*d[0] + p[1]*d[1]
+        c = p[0]**2 + p[1]**2 - self._r2
+        N = 0
+        delta = b**2 - a*c
+        status = False
+        if delta < 0:
+            t0 = np.nan
+            n = Vec3()
+        else:
+            t = np.array([(-b - np.sqrt(delta))/a, (-b + np.sqrt(delta))/a]) + 1e-6
+            t[t < 1e-6] = np.nan
+
+            I = P + ray.dir * t
+            t[np.abs(I[:, 2]) > self._hd] = np.nan
+
+            if np.isnan(t[0]) and np.isnan(t[1]):
+                t0 = np.nan
+            else:
+                N = np.nansum(t > 0)
+                t0 = np.nanmin(t)
+                status = True
+
+            I = ray.ori+a0 + ray.dir * t0
+            n = Vec3([I[0, 0], I[0, 1], 0])
+            n.normalise()
+
+        return status, n, Vec3(), t0, N
+
+
+class Sphere(Surface):
+    def __init__(self, surface_name, pos, sph_radius, zlim):
+        super().__init__(surface_name, 1, pos, curv_matrix=np.eye(2)/sph_radius)
+        self._zlim = zlim
+        self._sph_radius = np.abs(sph_radius)
+
+    def intersect(self, ray: Ray):
+        a0 = ray.dir * 1e-6
+
+        # find t0 distance between ray origin and intersection point
+        L = ray.ori + a0 - self.r_surf_obj
+        r2 = self._sph_radius ** 2
+        b = 2 * ray.dir.dot(L)
+        c = L.dot(L) - r2
+        D = b ** 2 - 4 * c
+        N = 0
+        if D >= 0:
+            sb = np.sign(b)
+            if sb == 0:
+                sb = 1
+            q = -.5 * (b + sb * np.sqrt(D))
+            t = np.array([q, c / q]) + 1e-6
+
+            t[t < 1e-6] = np.nan
+            pts = ray.ori + ray.dir * t
+
+            t[pts[:, 2] < self._zlim[0]] = np.nan
+
+            if np.isnan(t[0]) and np.isnan(t[1]):
+                status = False
+                t0 = np.nan
+            else:
+                N = np.nansum(t > 0)
+                t0 = np.nanmin(t)
+                status = True
+
+            pts = ray.ori + ray.dir * t0
+
+            # normal vector: always output
+            n = (pts - self.r_surf_obj)
+            n.normalise()
+        else:
+            t0 = np.inf
+            n = Vec3()
+            status = False
+
+        return status, n, Vec3(), t0, N

src/qosm/items/SlabMesh.py

+import gmsh
+from numpy import pi
+from qosm_core import Frame, Vec3
+
+from qosm.items.Mesh import Mesh
+
+ELEMENT_TYPES = {
+    "line": 2,
+    "line3": 3,
+    "triangle": 3,
+    "quadrangle": 4,
+    "tetrahedron": 4,
+    "hexahedron": 8,
+    "prism": 6,
+    "pyramid": 5,
+}
+
+
+class SlabMesh(Mesh):
+    def __init__(self, f_obj_glo: Frame):
+        super(SlabMesh, self).__init__(f_obj_glo)
+        self.frame = f_obj_glo
+
+    def load(self, element_size: float, shape: str, size: tuple, flip_normal: bool = False, obj_file=None,
+             offset: Vec3 = Vec3(), progress=None, view=False):
+        gmsh.clear()
+        gmsh.model.add("OBJECT")
+
+        if shape == 'sphere':
+            if len(size) == 4:
+                radius, angle1, angle2, angle3 = size
+            else:
+                radius = size[0]
+                angle1, angle2, angle3 = (-pi/2, pi/2, 2*pi)
+            gmsh.model.occ.addSphere(xc=0, yc=0, zc=0, radius=radius, angle1=angle1, angle2=angle2, angle3=angle3)
+        if shape == 'cylinder':
+            radius, length = size
+            gmsh.model.occ.addCylinder(x=0, y=0, z=-length/2, dx=0, dy=0, dz=length, r=radius)
+        if shape == 'disk':
+            if len(size) == 2:
+                L = 0.
+                ru, rv = size
+            else:
+                ru, rv, L = size
+            if ru > rv:
+                rx = ru
+                ry = rv
+                x_axis = (1, 0, 0)
+            else:
+                rx = rv
+                ry = ru
+                x_axis = (0, 1, 0)
+            L *= 1 - 2*float(flip_normal)
+            gmsh.model.occ.addDisk(xc=0, yc=0, zc=0, rx=rx, ry=ry,
+                                   zAxis=(0, 0, 1 - 2*float(flip_normal)), xAxis=x_axis)
+        elif shape == 'box':
+            width, height, length = size
+            gmsh.model.occ.addBox(x=-width/2, y=-height/2, z=-length/2, dx=width, dy=height, dz=length)
+
+        elif shape == 'rect':
+            width, height = size
+            tag = gmsh.model.occ.addRectangle(x=-width/2, y=-height/2, z=0, dx=width, dy=height)
+            if flip_normal:
+                gmsh.model.occ.rotate([(tag, 2)], 0, 0, 0, 1, 0, 0, pi)
+
+        # Mesh
+        gmsh.model.occ.synchronize()
+        gmsh.option.setNumber("Mesh.Algorithm", 6)
+        gmsh.option.setNumber('Mesh.MeshSizeMin', element_size)
+        gmsh.option.setNumber('Mesh.MeshSizeMax', element_size)
+
+        gmsh.model.mesh.generate(2)
+
+        create_obj = obj_file is not None
+        self.pre_build(1.0, create_obj, progress, view, offset=offset)
+
+        if create_obj:
+            with open('%s.obj' % obj_file, 'w') as f:
+                f.write('%s' % self._obj_txt)
+                f.close()

src/qosm/items/StepMesh.py

+import gmsh
+from qosm_core import Frame, Vec3
+
+from qosm.items.Mesh import Mesh
+
+
+class StepMesh(Mesh):
+    def __init__(self, f_obj_ref: Frame = Frame()):
+        super().__init__(f_obj_ref)
+
+    def load_step(self, filename, element_size: float, scale: float = 1e-3, create_obj=None, progress=None,
+                  view=False, show_vectors=True, offset: Vec3 = Vec3(), centre_shape=False):
+        gmsh.clear()
+        gmsh.model.add(filename)
+
+        # Mesh
+        gmsh.option.setNumber("Geometry.Tolerance", 1e-8)
+        gmsh.option.setNumber("Mesh.Algorithm", 6)
+        gmsh.option.setNumber('Mesh.MeshSizeMin', element_size)
+        gmsh.option.setNumber('Mesh.MeshSizeMax', element_size)
+
+        gmsh.option.setNumber('Mesh.SmoothNormals', 0)
+        gmsh.option.setNumber('Mesh.AngleSmoothNormals', 0)
+        gmsh.option.setNumber("Geometry.OCCImportLabels", 0)  # import colors from STEP
+
+        gmsh.model.occ.importShapes(filename, format="step")
+
+        gmsh.model.occ.removeAllDuplicates()
+        gmsh.model.occ.synchronize()
+
+        gmsh.model.mesh.generate(1)
+        gmsh.model.mesh.generate(2)
+        gmsh.model.mesh.removeDuplicateNodes()
+
+        self.pre_build(scale, create_obj, progress, view, show_vectors, offset, centre_shape)
+
+        if create_obj:
+            with open('%s.obj' % filename, 'w') as f:
+                f.write('%s' % self._obj_txt)
+                f.close()

src/qosm/items/Triangle.py

+import numpy as np
+
+from qosm.utils.Pose import Vector
+import qosm_core
+
+
+def is_outward_normal(vertex_normal, vertex_position, centroid):
+    # Compute vector from the vertex to the centroid
+    vector_to_centroid = centroid - vertex_position
+    vector_to_centroid /= np.sqrt(np.sum(vector_to_centroid**2))
+
+    # Compute dot product of the vertex normal and the vector
+    dot_product = np.dot(vertex_normal, vector_to_centroid)
+
+    # If dot product is negative, normal is outward (opposite direction)
+    return dot_product < 0
+
+
+class Triangle:
+    def __init__(self,
+                 vertices,
+                 normals,
+                 tangents,
+                 curvatures,
+                 dioptre,
+                 obj_centroid=None,
+                 id_tri: int = -1,
+                 bounces=None):
+        if bounces is None:
+            bounces = [-1, 1]
+        self.dioptre = dioptre
+        self.id_tri = id_tri
+        self.bounces = bounces
+        v0, v1, v2 = vertices
+
+        _edge1, _edge2, _edge3 = (v1 - v0, v2 - v0, v2 - v1)
+        n0, n1, n2 = normals
+        t0, t1, t2 = tangents
+
+        # print(_edge1.reshape((-1, 3)).norm(), _edge2.reshape((-1, 3)).norm(), _edge3.reshape((-1, 3)).norm())
+
+        n = np.cross(_edge1/np.sqrt(np.sum(_edge1**2)), _edge2/np.sqrt(np.sum(_edge2**2)))
+        n /= np.sqrt(np.sum(n**2))
+        b = (v0 + v1 + v2) / 3.
+
+        if np.sum(n0**2) == 0 or np.sum(n1**2) == 0 or np.sum(n2**2) == 0:
+            raise Exception('Null normal detected (n0, n1, n2): %d %d %d' %
+                            (np.sum(n0**2) == 0, np.sum(n1**2) == 0, np.sum(n2**2) == 0))
+
+        # print(id_tri, n0.tolist(), n1.tolist(), n2.tolist())
+
+        if obj_centroid is not None:
+            n_outward = is_outward_normal(n, b, obj_centroid)
+            if not n_outward:
+                v0, v2, v1 = vertices
+                n0, n2, n1 = normals
+                _edge1, _edge2 = (v1 - v0, v2 - v0)
+
+        # Some vertices can lie on a degenerate or singular point where the normal vector cannot be uniquely defined.
+        # This might occur, for example, at sharp edges, corners, or vertices where multiple surfaces meet.
+        if np.sqrt(np.sum(n0 ** 2)) < 1e-6:
+            tests = Vector([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+            test_norm = tests.dot(t0.reshape((-1, 3)))
+            tests = tests[test_norm < .5]
+            candidates = tests.cross(t0.reshape((-1, 3))).normalised()
+            test_dot_n = candidates.dot(n.reshape((-1, 3)))
+            idx_candidate = np.abs(test_dot_n).argmax()
+            n0 = candidates[idx_candidate, :]
+        if np.sqrt(np.sum(n1 ** 2)) < 1e-6:
+            tests = Vector([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+            test_norm = tests.dot(t1.reshape((-1, 3)))
+            tests = tests[test_norm < .5]
+            candidates = tests.cross(t1.reshape((-1, 3))).normalised()
+            test_dot_n = candidates.dot(n.reshape((-1, 3)))
+            idx_candidate = np.abs(test_dot_n).argmax()
+            n1 = candidates[idx_candidate, :]
+        if np.sqrt(np.sum(n2 ** 2)) < 1e-6:
+            tests = Vector([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+            test_norm = tests.dot(t2.reshape((-1, 3)))
+            tests = tests[test_norm < .5]
+            candidates = tests.cross(t2.reshape((-1, 3))).normalised()
+            test_dot_n = candidates.dot(n.reshape((-1, 3)))
+            idx_candidate = np.abs(test_dot_n).argmax()
+            n2 = candidates[idx_candidate, :]
+
+        if np.dot(n0, n) < 0:
+            n0 = n0 * -1
+        if np.dot(n1, n) < 0:
+            n1 = n1 * -1
+        if np.dot(n2, n) < 0:
+            n2 = n2 * -1
+
+        self.normals = (n0, n1, n2)
+        self.curvatures = curvatures
+        self.vertices = (v0, v1, v2)
+
+        # store precomputed values
+        self.edges = (_edge1, _edge2)
+
+    @property
+    def cpp(self):
+        v0, v1, v2 = self.vertices
+        n0, n1, n2 = self.normals
+        c0, c1, c2 = self.curvatures
+        e0, e1 = self.edges
+        dioptre = self.dioptre
+        tri = qosm_core.Triangle(
+            vertices=(
+                qosm_core.Vec3(v0[0], v0[1], v0[2]),
+                qosm_core.Vec3(v1[0], v1[1], v1[2]),
+                qosm_core.Vec3(v2[0], v2[1], v2[2])),
+            normals=(
+                qosm_core.Vec3(n0[0], n0[1], n0[2]),
+                qosm_core.Vec3(n1[0], n1[1], n1[2]),
+                qosm_core.Vec3(n2[0], n2[1], n2[2])),
+            edges=(qosm_core.Vec3(e0[0], e0[1], e0[2]),
+                   qosm_core.Vec3(e1[0], e1[1], e1[2])),
+            curvatures=(c0, c1, c2),
+            dioptre=dioptre,
+            id_triangle=self.id_tri
+        )
+        tri.bounces = self.bounces
+        return tri