O2BVHSurfaceSolid.h

Go to the documentation of this file.

// Copyright 2019-2026 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
#ifndef ALICEO2_CADSUPPORT_O2BVHSURFACESOLID_
#define ALICEO2_CADSUPPORT_O2BVHSURFACESOLID_
 
#include "TGeoBBox.h"
 
#include <array>
#include <iosfwd>
#include <vector>
 
class TBuffer3D;
 
namespace o2
{
namespace cad
{
 
struct BVHSurfaceCurveRecord {
  int kind = 0; 
  double lineStart[2] = {0., 0.};
  double lineEnd[2] = {0., 0.};
  double center[2] = {0., 0.};
  double radius = 0.;
  double startAngle = 0.;
  double endAngle = 0.;
  int degree = 0;              
  std::vector<double> poles;   
  std::vector<double> weights; 
  std::vector<double> knots;   
};
 
struct BVHSurfaceRecord {
  enum Kind { PlanarPolygon = 0,
              CurvedPlanar = 1,
              Cylindrical = 2,
              Spherical = 3,
              Conical = 4,
              Toroidal = 5 };
 
  int kind = PlanarPolygon;
  double origin[3] = {0., 0., 0.}; 
  double axisA[3] = {0., 0., 0.};  
  double axisB[3] = {0., 0., 0.};  
 
  std::vector<double> scalars;
 
  bool innerWall = false;
  bool trimmed = false; 
 
  std::vector<double> polygonPoints;
  std::vector<BVHSurfaceCurveRecord> curves;
  std::vector<int> wireSizes;
 
  std::vector<unsigned int> boundaryEdgeIds;
  std::vector<unsigned char> boundaryEdgeFlags;
 
  static int expectedScalarCount(int recordKind);
};
 
class O2BVHSurfaceSolid : public TGeoBBox
{
 public:
  using Point2D = std::array<double, 2>;
  using Point3D = std::array<double, 3>;
 
  O2BVHSurfaceSolid();
  explicit O2BVHSurfaceSolid(const char* name);
  ~O2BVHSurfaceSolid() override;
 
  O2BVHSurfaceSolid(const O2BVHSurfaceSolid&) = delete;
  O2BVHSurfaceSolid& operator=(const O2BVHSurfaceSolid&) = delete;
 
  bool AddPlanarSurface(const Point3D& origin, const Point3D& axisU, const Point3D& axisV,
                        const std::vector<Point2D>& outerWire,
                        const std::vector<std::vector<Point2D>>& innerWires = {});
 
  struct PlanarBoundaryCurve {
    enum Kind { Line,
                Arc,
                BSpline };
    Kind kind = Line;
    Point2D lineStart{0., 0.};
    Point2D lineEnd{0., 0.};
    Point2D center{0., 0.};
    double radius = 0.;
    double startAngle = 0.;
    double endAngle = 0.;
    int degree = 0;              
    std::vector<Point2D> poles;  
    std::vector<double> weights; 
    std::vector<double> knots;   
 
    static PlanarBoundaryCurve makeLine(const Point2D& start, const Point2D& end)
    {
      PlanarBoundaryCurve curve;
      curve.kind = Line;
      curve.lineStart = start;
      curve.lineEnd = end;
      return curve;
    }
    static PlanarBoundaryCurve makeArc(const Point2D& c, double r, double start, double end)
    {
      PlanarBoundaryCurve curve;
      curve.kind = Arc;
      curve.center = c;
      curve.radius = r;
      curve.startAngle = start;
      curve.endAngle = end;
      return curve;
    }
    static PlanarBoundaryCurve makeBSpline(int splineDegree, std::vector<Point2D> splinePoles,
                                           std::vector<double> splineWeights, std::vector<double> splineKnots)
    {
      PlanarBoundaryCurve curve;
      curve.kind = BSpline;
      curve.degree = splineDegree;
      curve.poles = std::move(splinePoles);
      curve.weights = std::move(splineWeights);
      curve.knots = std::move(splineKnots);
      return curve;
    }
  };
 
  bool AddCurvedPlanarSurface(const Point3D& origin, const Point3D& axisU, const Point3D& axisV,
                              const std::vector<PlanarBoundaryCurve>& outerWire,
                              const std::vector<std::vector<PlanarBoundaryCurve>>& innerWires = {});
 
  bool AddCylindricalSurface(const Point3D& centerPoint, const Point3D& axis, const Point3D& referenceAxisU,
                             double radius, double heightMin, double heightMax, double phiStart = 0.,
                             double phiSweep = 6.283185307179586, bool innerWall = false);
 
  bool AddCylindricalSurface(const Point3D& centerPoint, const Point3D& axis, const Point3D& referenceAxisU,
                             double radius, double heightMin, double heightMax, double phiStart, double phiSweep,
                             bool innerWall, const std::vector<PlanarBoundaryCurve>& outerTrim,
                             const std::vector<std::vector<PlanarBoundaryCurve>>& innerTrims = {});
 
  bool AddSphericalSurface(const Point3D& center, const Point3D& polarAxis, const Point3D& referenceAxisU,
                           double radius, double thetaMin = 0., double thetaMax = 3.141592653589793,
                           double phiStart = 0., double phiSweep = 6.283185307179586, bool innerWall = false);
 
  bool AddSphericalSurface(const Point3D& center, const Point3D& polarAxis, const Point3D& referenceAxisU,
                           double radius, double thetaMin, double thetaMax, double phiStart, double phiSweep,
                           bool innerWall, const std::vector<PlanarBoundaryCurve>& outerTrim,
                           const std::vector<std::vector<PlanarBoundaryCurve>>& innerTrims = {});
 
  bool AddConicalSurface(const Point3D& centerPoint, const Point3D& axis, const Point3D& referenceAxisU,
                         double radiusAtMin, double radiusAtMax, double heightMin, double heightMax,
                         double phiStart = 0., double phiSweep = 6.283185307179586, bool innerWall = false);
 
  bool AddConicalSurface(const Point3D& centerPoint, const Point3D& axis, const Point3D& referenceAxisU,
                         double radiusAtMin, double radiusAtMax, double heightMin, double heightMax, double phiStart,
                         double phiSweep, bool innerWall, const std::vector<PlanarBoundaryCurve>& outerTrim,
                         const std::vector<std::vector<PlanarBoundaryCurve>>& innerTrims = {});
 
  bool AddToroidalSurface(const Point3D& centerPoint, const Point3D& axis, const Point3D& referenceAxisU,
                          double majorRadius, double minorRadius, double phiStart = 0.,
                          double phiSweep = 6.283185307179586, double tubeStart = 0.,
                          double tubeSweep = 6.283185307179586, bool innerWall = false);
 
  bool AddToroidalSurface(const Point3D& centerPoint, const Point3D& axis, const Point3D& referenceAxisU,
                          double majorRadius, double minorRadius, double phiStart, double phiSweep, double tubeStart,
                          double tubeSweep, bool innerWall, const std::vector<PlanarBoundaryCurve>& outerTrim,
                          const std::vector<std::vector<PlanarBoundaryCurve>>& innerTrims = {});
 
  enum BoundaryEdgeFlag : unsigned char {
    kEdgeReversed = 1u << 0,   
    kEdgeDegenerate = 1u << 1, 
    kEdgeAnchored = 1u << 2    
  };
 
  bool SetSurfaceBoundaryEdges(int surfaceIndex, const std::vector<unsigned int>& edgeIds,
                               const std::vector<unsigned char>& edgeFlags);
 
  void CloseShape(bool check = true);
 
  int GetNsurfaces() const;
  bool IsDefined() const;
 
  void SetModelTolerance(double toleranceCm);
  double GetModelTolerance() const { return fModelTolerance; }
 
  bool HasBVH() const;
  bool GetBVHRootBounds(Point3D& lower, Point3D& upper) const;
  int CountBVHRayCandidates(const Point3D& point, const Point3D& direction) const;
 
  static void SetRayTMaxPruning(bool enable);
  static bool GetRayTMaxPruning();
 
  static void ResetRayCandidateCounter();
  static long long GetRayCandidateCount();
 
  static void ResetSafetyCandidateCounter();
  static long long GetSafetyCandidateCount();
 
  static void SetSafetyBoundUnsoundForTest(bool enable);
  static bool GetSafetyBoundUnsoundForTest();
 
  struct ContainsCrossing {
    double distance = 0.;        
    double normalAlignment = 0.; 
    bool onTrimBoundary = false;
  };
 
  void DescribeContainsCrossings(const Point3D& point, std::vector<ContainsCrossing>& bvhCrossings,
                                 std::vector<ContainsCrossing>& loopCrossings) const;
 
  void DescribeContainsCrossings(const Point3D& point, const Point3D& direction,
                                 std::vector<ContainsCrossing>& bvhCrossings,
                                 std::vector<ContainsCrossing>& loopCrossings) const;
 
  bool IsClosed() const;
  bool IsOrientationConsistent() const;
 
  enum class NavigationReliability {
    Undetermined = 0, 
    Reliable,         
    ReversedFaces,    
    OpenSurfaceSet,   
    NonManifold       
  };
 
  NavigationReliability GetNavigationReliability() const;
  bool IsNavigable() const;
  static const char* GetNavigationReliabilityName(NavigationReliability reliability);
 
  int GetBoundaryEdgeCount() const;
  int GetNonManifoldEdgeCount() const;
  int GetReversedEdgeCount() const;
 
  double GetMaxRimIsolation() const;
  bool HasEdgeIdentity() const;
  int GetSourceEdgeCount() const;
  int GetSharedSourceEdgeCount() const;
  int GetBoundarySourceEdgeCount() const;
  int GetNonManifoldSourceEdgeCount() const;
  int GetReversedSourceEdgeCount() const;
  int GetDegenerateSourceEdgeCount() const;
  double GetMaxSharedEdgeDeviation() const;
  int GetMeasuredSharedEdgeCount() const;
  int GetUnmeasuredSharedEdgeCount() const;
  double GetRimChordResolution() const;
  double GetRimMatchTolerance() const;
  double GetTotalRimLength() const;
  double GetUnmatchedRimLength() const;
  int GetRimCount() const;
  int GetMatchedRimCount() const;
  int GetBoundaryRimCount() const;
  int GetNonManifoldRimCount() const;
  int GetReversedRimCount() const;
 
  struct RimReport {
    int surface = -1;      
    int rimOnSurface = -1; 
    bool closed = false;   
    int chords = 0;
    int unmatchedChords = 0;     
    double length = 0.;          
    double unmatchedLength = 0.; 
    double maxIsolation = 0.;
    std::array<double, 3> maxIsolationPoint{{0., 0., 0.}};
    int maxIsolationFace = -1;
    NavigationReliability state = NavigationReliability::Undetermined;
  };
  const std::vector<RimReport>& GetRimReports() const;
 
  void GetSurfaceCapacityContributions(std::vector<double>& contributions) const;
 
  void ComputeBBox() override;
 
  int DistancetoPrimitive(int, int) override { return 99999; }
  const TBuffer3D& GetBuffer3D(int reqSections, Bool_t localFrame) const override;
  void GetMeshNumbers(int& nvert, int& nsegs, int& npols) const override;
  int GetNmeshVertices() const override;
 
  Bool_t GetPointsOnSegments(Int_t npoints, Double_t* array) const override;
 
  static constexpr double kSurfacePointTolerance = 1.e-11;
 
  void InspectShape() const override {}
  TBuffer3D* MakeBuffer3D() const override;
  void Print(Option_t* option = "") const override;
  void SavePrimitive(std::ostream&, Option_t*) override {}
  void SetPoints(double* points) const override;
  void SetPoints(Float_t* points) const override;
  void SetSegsAndPols(TBuffer3D& buff) const override;
  void Sizeof3D() const override {}
 
  Double_t DistFromOutside(const Double_t* point, const Double_t* dir, Int_t iact = 1,
                           Double_t step = TGeoShape::Big(), Double_t* safe = nullptr) const override;
  Double_t DistFromInside(const Double_t* point, const Double_t* dir, Int_t iact = 1,
                          Double_t step = TGeoShape::Big(), Double_t* safe = nullptr) const override;
  bool Contains(const Double_t* point) const override;
  bool Contains_Loop(const Double_t* point) const;
  bool ContainsAlongDirection(const Double_t* point, const Double_t* direction) const;
  Double_t DistFromOutside_Loop(const Double_t* point, const Double_t* dir,
                                Double_t stepmax = TGeoShape::Big()) const;
  Double_t DistFromInside_Loop(const Double_t* point, const Double_t* dir,
                               Double_t stepmax = TGeoShape::Big()) const;
  Double_t Safety(const Double_t* point, Bool_t in = kTRUE) const override;
  void ComputeNormal(const Double_t* point, const Double_t* dir, Double_t* norm) const override;
  Double_t Safety_Loop(const Double_t* point, Bool_t in = kTRUE) const;
  void ComputeNormal_Loop(const Double_t* point, const Double_t* dir, Double_t* norm) const;
  Double_t Capacity() const override;
 
  const std::vector<BVHSurfaceRecord>& GetSurfaceRecords() const { return fRecords; }
 
 private:
  bool containsByParity(const Double_t* point, bool useBVH) const;
 
  void computeNormalFrom(const Double_t* point, const Double_t* dir, Double_t* norm, bool useLoop) const;
 
  bool RebuildFromRecords();
 
  bool ProjectOntoPatch(int surfaceIndex, double* point) const;
 
  struct Impl;
  Impl* fImpl = nullptr; 
 
  std::vector<BVHSurfaceRecord> fRecords;
 
  double fModelTolerance = 0.;
 
  ClassDefOverride(O2BVHSurfaceSolid, 3) // BVH surface-bounded shape class
};
 
} // namespace cad
} // namespace o2
 
#endif