DDCAD/src/ASSIMPWriter.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013
0014 /// Framework include files
0015 #include <DD4hep/Detector.h>
0016 #include <DD4hep/Printout.h>
0017 #include <DDCAD/ASSIMPWriter.h>
0018 #include <DDCAD/Utilities.h>
0019
0020 /// Open Asset Importer Library
0021 #include "assimp/postprocess.h"
0022 #include "assimp/Exporter.hpp"
0023 #include "assimp/scene.h"
0024
0025 /// ROOT include files
0026 #include <TBuffer3D.h>
0027 #include <TBuffer3DTypes.h>
0028 #include <TClass.h>
0029 #include <TGeoBoolNode.h>
0030 #include <TGeoMatrix.h>
0031 #include <CsgOps.h>
0032
0033 /// C/C++ include files
0034 #include <set>
0035
0036 using namespace dd4hep::cad;
0037
0038 using Vertex = Tessellated::Vertex_t;
0039
0040 namespace  {
0041
0042   void _collect(std::vector<std::pair<dd4hep::PlacedVolume,TGeoHMatrix*> >& cont,
0043                 bool recursive, const TGeoHMatrix& to_global, dd4hep::PlacedVolume pv)
0044   {
0045     dd4hep::Volume v = pv.volume();
0046     for(Int_t i=0; i<v->GetNdaughters(); ++i)  {
0047       dd4hep::PlacedVolume p   = v->GetNode(i);
0048       dd4hep::Solid        sol = p.volume().solid();
0049       bool                 use = sol->IsA() != TGeoShapeAssembly::Class();
0050       std::unique_ptr<TGeoHMatrix> mother(new TGeoHMatrix(to_global));
0051       mother->Multiply(p->GetMatrix());
0052
0053       if ( use )  {
0054         TGeoHMatrix* m = mother.release();
0055         cont.push_back(std::make_pair(p, m));
0056         if ( recursive )
0057           _collect(cont, recursive, *m, p);
0058       }
0059       else if ( recursive )  {
0060         _collect(cont, recursive, *mother, p);
0061       }
0062     }
0063   }
0064   bool equals(Vertex const &lhs, Vertex const &rhs)  {
0065     constexpr double kTolerance = 1.e-32;
0066     return TMath::Abs(lhs[0] - rhs[0]) < kTolerance &&
0067       TMath::Abs(lhs[1] - rhs[1]) < kTolerance &&
0068       TMath::Abs(lhs[2] - rhs[2]) < kTolerance;
0069   }
0070
0071   struct TessellateShape   {
0072   public:
0073     TessellateShape() = default;
0074     virtual ~TessellateShape() = default;
0075     RootCsg::TBaseMesh* make_mesh(TGeoShape* sh)  const;
0076     RootCsg::TBaseMesh* collect_composite(TGeoCompositeShape* sh)    const;
0077     std::unique_ptr<TGeoTessellated> build_mesh(int id, const std::string& name, TGeoShape* shape);
0078     std::unique_ptr<TGeoTessellated> tessellate_primitive(const std::string& name, dd4hep::Solid solid);
0079     std::unique_ptr<TGeoTessellated> close_tessellated(int id, TGeoShape* shape, int nskip, std::unique_ptr<TGeoTessellated>&& tes);
0080   };
0081
0082   std::unique_ptr<TGeoTessellated>
0083   TessellateShape::close_tessellated(int id, TGeoShape* shape, int nskip, std::unique_ptr<TGeoTessellated>&& tes)   {
0084     std::string nam = shape->GetName(), typ = "["+std::string(shape->IsA()->GetName())+"]";
0085     nam = nam.substr(0, nam.find("_0x"));
0086     tes->CloseShape(true, true, true);
0087     if ( nskip > 0 )   {
0088       dd4hep::printout(dd4hep::ALWAYS,
0089                        "ASSIMPWriter","+++ %3d %-48s %-24s Skipped %3ld/%-4d degenerate facets %4d vertices closed:%s defined:%s",
0090                        id, nam.c_str(), typ.c_str(), nskip, tes->GetNfacets(),  tes->GetNvertices(),
0091                        dd4hep::yes_no(tes->IsClosedBody()), dd4hep::yes_no(tes->IsDefined()));
0092     }
0093     else   {
0094       dd4hep::printout(dd4hep::ALWAYS,
0095                        "ASSIMPWriter","+++ %3d %-48s %-24s Tessellated %4d facets %4d vertices closed:%s defined:%s",
0096                        id, nam.c_str(), typ.c_str(),
0097                        tes->GetNfacets(), tes->GetNvertices(),
0098                        dd4hep::yes_no(tes->IsClosedBody()),
0099                        dd4hep::yes_no(tes->IsDefined()));
0100     }
0101     std::cout << std::flush;
0102     return std::move(tes);
0103   }
0104
0105   std::unique_ptr<TGeoTessellated> TessellateShape::build_mesh(int id, const std::string& name, TGeoShape* shape)      {
0106     auto mesh = std::unique_ptr<RootCsg::TBaseMesh>(this->make_mesh(shape));
0107     std::vector<Vertex> vertices;
0108     std::size_t nskip = 0;
0109
0110     vertices.reserve(mesh->NumberOfVertices());
0111     std::map<std::size_t,std::size_t> vtx_index_replacements;
0112     for( size_t ipoint = 0, npoints = mesh->NumberOfVertices(); ipoint < npoints; ++ipoint )   {
0113       long found = -1;
0114       const double* v = mesh->GetVertex(ipoint);
0115       Vertex vtx(v[0], v[1], v[2]);
0116       for(std::size_t i=0; i < vertices.size(); ++i)   {
0117         if ( equals(vertices[i],vtx) )  {
0118           vtx_index_replacements[ipoint] = found = i;
0119           break;
0120         }
0121       }
0122       if ( found < 0 )   {
0123         vtx_index_replacements[ipoint] = vertices.size();
0124         vertices.emplace_back(v[0], v[1], v[2]);
0125       }
0126     }
0127     std::size_t vtx_len = vertices.size();
0128     std::unique_ptr<TGeoTessellated> tes = std::make_unique<TGeoTessellated>(name.c_str(), vertices);
0129     for( std::size_t ipoly = 0, npols = mesh->NumberOfPolys(); ipoly < npols; ++ipoly)    {
0130       std::size_t npoints = mesh->SizeOfPoly(ipoly);
0131       if ( npoints >= 3 )   {
0132         printout(dd4hep::DEBUG,"ASSIMPWriter","+++ Got polygon with %ld points",npoints);
0133         ///
0134         /// 3-vertex polygons automatically translate to GL_TRIANGLES
0135         /// See Kronos documentation to glBegin / glEnd from the glu library:
0136         /// https://www.khronos.org/registry/OpenGL-Refpages/gl2.1/xhtml/glBegin.xml
0137         ///
0138         /// Otherwise:
0139 #if 1
0140         /// Apparently this is the correct choice:
0141         ///
0142         /// Interprete as FAN:  GL_TRIANGLE_FAN
0143         /// One triangle is defined for each vertex presented after the first two vertices.
0144         /// Vertices 1 , n + 1 , and n + 2 define triangle n.
0145         /// N - 2 triangles are drawn.
0146         std::size_t v0  = mesh->GetVertexIndex(ipoly, 0);
0147         std::size_t vv0 = vtx_index_replacements[v0];
0148         for( std::size_t ipoint = 0; ipoint < npoints-2; ++ipoint )   {
0149           std::size_t v1 = mesh->GetVertexIndex(ipoly, ipoint+1);
0150           std::size_t v2 = mesh->GetVertexIndex(ipoly, ipoint+2);
0151           std::size_t vv1 = vtx_index_replacements[v1];
0152           std::size_t vv2 = vtx_index_replacements[v2];
0153           if ( vv0 > vtx_len || vv1 > vtx_len || vv2 > vtx_len )  {
0154             ++nskip;
0155             continue;
0156           }
0157           if ( vv0 == vv1 || vv0 == vv2 || vv1 == vv2 )   {
0158             ++nskip;
0159             continue;
0160           }
0161           Vertex w[3] = {vertices[vv0],vertices[vv1],vertices[vv2]};
0162           if ( TGeoFacet::CompactFacet(w, 3) < 3 )   {
0163             ++nskip;
0164             continue;
0165           }
0166 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
0167           bool degenerated = dd4hep::cad::facetIsDegenerated({vertices[vv0],vertices[vv1],vertices[vv2]});
0168 #else
0169           bool degenerated = true;
0170           TGeoFacet f(&vertices, 3, vv0, vv1, vv2);
0171           f.ComputeNormal(degenerated);
0172 #endif
0173           if ( degenerated )    {
0174             ++nskip;
0175             continue;
0176           }
0177           tes->AddFacet(vv0, vv1, vv2);
0178         }
0179 #else
0180         /// Interprete as STRIP: GL_TRIANGLE_STRIP
0181         /// One triangle is defined for each vertex presented after the first two vertices.
0182         /// For odd n, vertices n, n + 1 , and n + 2 define triangle n.
0183         /// For even n, vertices n + 1 , n, and n + 2 define triangle n.
0184         /// N - 2 triangles are drawn.
0185         for( std::size_t ipoint = 0; ipoint < npoints-2; ++ipoint )   {
0186           vtx_t v0(mesh->GetVertex(mesh->GetVertexIndex(ipoly, ipoint)));
0187           vtx_t v1(mesh->GetVertex(mesh->GetVertexIndex(ipoly, ipoint+1)));
0188           vtx_t v2(mesh->GetVertex(mesh->GetVertexIndex(ipoly, ipoint+2)));
0189           ((ipoint%2) == 0) ? tes->AddFacet(v1, v0, v2) : tes->AddFacet(v0, v1, v2);
0190         }
0191 #endif
0192       }
0193     }
0194     return close_tessellated(id, shape, nskip, std::move(tes));
0195   }
0196
0197   RootCsg::TBaseMesh* TessellateShape::make_mesh(TGeoShape* sh)   const   {
0198     if (TGeoCompositeShape *shape = dynamic_cast<TGeoCompositeShape *>(sh))   {
0199       return collect_composite(shape);
0200     }
0201     UInt_t  flags = TBuffer3D::kCore|TBuffer3D::kBoundingBox|TBuffer3D::kRawSizes|TBuffer3D::kRaw|TBuffer3D::kShapeSpecific;
0202     const TBuffer3D& buffer = sh->GetBuffer3D(flags, kFALSE);
0203     return RootCsg::ConvertToMesh(buffer);
0204   }
0205
0206   RootCsg::TBaseMesh* TessellateShape::collect_composite(TGeoCompositeShape* sh)  const  {
0207     TGeoBoolNode* node  = sh->GetBoolNode();
0208     TGeoShape*    left  = node->GetLeftShape();
0209     TGeoShape*    right = node->GetRightShape();
0210     TGeoHMatrix*  glmat = (TGeoHMatrix*)TGeoShape::GetTransform();
0211     UInt_t        oper  = node->GetBooleanOperator();
0212     TGeoHMatrix   copy(*glmat); // keep a copy
0213
0214     // Do not wonder about this logic.
0215     // GetBuffer3D (->make_mesh) uses static variable fgTransform of TGeoShape!
0216     glmat->Multiply(node->GetLeftMatrix());
0217     auto left_mesh  = std::unique_ptr<RootCsg::TBaseMesh>(make_mesh(left));
0218     *glmat = &copy;
0219
0220     glmat->Multiply(node->GetRightMatrix());
0221     auto right_mesh = std::unique_ptr<RootCsg::TBaseMesh>(make_mesh(right));
0222     *glmat = &copy;
0223
0224     switch (oper) {
0225     case TGeoBoolNode::kGeoUnion:
0226       return RootCsg::BuildUnion(left_mesh.get(), right_mesh.get());
0227     case TGeoBoolNode::kGeoIntersection:
0228       return RootCsg::BuildIntersection(left_mesh.get(), right_mesh.get());
0229     case TGeoBoolNode::kGeoSubtraction:
0230       return RootCsg::BuildDifference(left_mesh.get(), right_mesh.get());
0231     default:
0232       Error("BuildComposite", "Wrong boolean operation code %d\n", oper);
0233       return 0;
0234     }
0235   }
0236
0237   std::unique_ptr<TGeoTessellated> TessellateShape::tessellate_primitive(const std::string& name, dd4hep::Solid solid)   {
0238     using  vtx_t = Vertex;
0239     const  TBuffer3D& buf3D = solid->GetBuffer3D(TBuffer3D::kAll, false);
0240     struct pol_t { int c, n; int segs[1]; } *pol = nullptr;
0241     struct seg_t { int c, _1, _2;         };
0242     const  seg_t* segs = (seg_t*)buf3D.fSegs;
0243     const  vtx_t* vtcs = (vtx_t*)buf3D.fPnts;
0244     std::size_t i, num_facet = 0;
0245     const  Int_t* q;
0246
0247     for( i=0, q=buf3D.fPols; i<buf3D.NbPols(); ++i, q += (2+pol->n))  {
0248       pol = (pol_t*)q;
0249       for( int j=0; j < pol->n-1; ++j ) num_facet += 2;
0250     }
0251
0252     std::unique_ptr<TGeoTessellated> tes = std::make_unique<TGeoTessellated>(name.c_str(), num_facet);
0253     q = buf3D.fPols;
0254     for( i=0, q=buf3D.fPols; i<buf3D.NbPols(); ++i)  {
0255       pol = (pol_t*)q;
0256       q += (2+pol->n);
0257       for( int j=0; j < pol->n; j += 2 )   {
0258         /* ------------------------------------------------------------
0259         //   Convert quadri-linear facet to 2 tri-linear facets
0260         //
0261         //    f1 +---------------+ v2/v3: f0
0262         //      /                /
0263         //     /                /
0264         //    /                /
0265         //    +---------------+
0266         //  v0             v1 v2/v3
0267         // --------------------------------------------------------- */
0268         const int    s1  = pol->segs[j], s2 = pol->segs[(j+1)%pol->n];
0269         const int    s[] = { segs[s1]._1, segs[s1]._2, segs[s2]._1, segs[s2]._2 };
0270         const vtx_t& v0  = vtcs[s[0]], &v1=vtcs[s[1]], &v2=vtcs[s[2]], &v3=vtcs[s[3]];
0271
0272         if ( s[0] == s[2] )   {       // Points are ( s[1], s[0], s[3] )
0273           tes->AddFacet(v1, v0, v3);
0274         }
0275         else if ( s[0] == s[3] )   {  // Points are ( s[1], s[0], s[2] )
0276           tes->AddFacet(v1, v0, v2);
0277         }
0278         else if ( s[1] == s[2] )   {  // Points are ( s[0], s[1], s[3] )
0279           tes->AddFacet(v0, v1, v3);
0280         }
0281         else   {                      // Points are ( s[0], s[1], s[2] )
0282           tes->AddFacet(v0, v1, v2);
0283         }
0284       }
0285     }
0286     return tes;
0287   }
0288 }
0289
0290 /// Write output file
0291 int ASSIMPWriter::write(const std::string& file_name,
0292                         const std::string& file_type,
0293                         const VolumePlacements& places,
0294                         bool   recursive,
0295                         double unit_scale)  const
0296 {
0297   std::vector<std::pair<PlacedVolume,TGeoHMatrix*> > placements;
0298   int  build_mode  = ((flags&0x1) != 0) ? 1 : 0;
0299   bool dump_facets = ((flags&0x2) != 0);
0300   std::vector<Material>  materials;
0301   TGeoHMatrix            toGlobal;
0302
0303   for( auto pv : places )
0304     _collect(placements, recursive, toGlobal, pv);
0305
0306   std::size_t num_mesh = placements.size();
0307
0308   aiScene scene;
0309   scene.mNumMaterials = 0;
0310   scene.mNumMeshes    = 0;
0311   scene.mMeshes       = new aiMesh* [num_mesh];
0312   scene.mMaterials    = new aiMaterial* [num_mesh];
0313
0314   aiNode *root        = new aiNode();
0315   scene.mRootNode     = root;
0316   root->mName.Set("<STL>");
0317   root->mNumMeshes    = 0;
0318   root->mNumChildren  = 0;
0319   root->mChildren     = new aiNode* [num_mesh];
0320   root->mMeshes       = 0;
0321   auto* geo_transform = TGeoShape::GetTransform();
0322
0323   TGeoHMatrix identity;
0324   for( std::size_t imesh=0; imesh < num_mesh; ++imesh )   {
0325     std::unique_ptr<TGeoHMatrix>     trafo(placements[imesh].second);
0326     std::unique_ptr<TGeoTessellated> shape_holder;
0327     PlacedVolume     pv  = placements[imesh].first;
0328     Volume           vol = pv.volume();
0329     Solid            sol = vol.solid();
0330     Material         mat = vol.material();
0331     TessellatedSolid tes = sol;
0332     aiString         node_name(vol.name());
0333
0334     identity.Clear();
0335     TGeoShape::SetTransform(&identity);
0336
0337     /// If the shape is already tessellated, no need to create another one!
0338     if ( !tes.isValid() )   {
0339       TessellateShape helper;
0340       auto* shape = dynamic_cast<TGeoCompositeShape*>(sol.ptr());
0341       if ( build_mode || shape )   {  // Always use this method!
0342         auto* paintVol = detector.manager().GetPaintVolume();
0343         detector.manager().SetPaintVolume(vol.ptr());
0344         shape_holder = helper.build_mesh(imesh, vol.name(), sol.ptr());
0345         detector.manager().SetPaintVolume(paintVol);
0346       }
0347       else   {
0348         shape_holder = helper.tessellate_primitive(vol.name(), sol);
0349       }
0350       tes = shape_holder.get();
0351     }
0352     if ( tes->GetNfacets() == 0 )   {
0353       continue;
0354     }
0355
0356     std::size_t num_vert = 0;
0357     for( long j=0, n=tes->GetNfacets(); j < n; ++j )
0358       num_vert += tes->GetFacet(j).GetNvert();
0359
0360     std::size_t index = std::numeric_limits<std::size_t>::max();
0361     for( std::size_t j=0; j<materials.size(); ++j )  {
0362       if( materials[j] == mat )   {
0363         index = j;
0364         break;
0365       }
0366     }
0367     if ( index > materials.size() )   {
0368       aiString name(mat.name());
0369       auto* ai_mat = new aiMaterial();
0370       index = materials.size();
0371       materials.push_back(mat);
0372       ai_mat->AddProperty(&name, AI_MATKEY_NAME);
0373       scene.mMaterials[scene.mNumMaterials] = ai_mat;
0374       ++scene.mNumMaterials;
0375     }
0376     aiMesh* mesh = new aiMesh;
0377     mesh->mName = node_name;
0378     mesh->mMaterialIndex = index;
0379     if ( vol.visAttributes().isValid() )   {
0380       float cr = 0e0, cg = 0e0, cb = 0e0, ca = 0e0;
0381       vol.visAttributes().argb(ca, cr, cg, cb);
0382       mesh->mColors[0] = new aiColor4D(cr, cg, cb, ca);
0383     }
0384     mesh->mFaces       = new aiFace[tes->GetNfacets()];
0385     mesh->mVertices    = new aiVector3D[num_vert];
0386     mesh->mNormals     = new aiVector3D[num_vert];
0387     mesh->mTangents    = nullptr;
0388     mesh->mBitangents  = nullptr;
0389     mesh->mNumFaces    = 0;
0390     mesh->mNumVertices = 0;
0391
0392     for( long j=0, n=tes->GetNfacets(); j < n; ++j )   {
0393       aiFace& face = mesh->mFaces[j];
0394       face.mNumIndices = 0;
0395       face.mIndices = nullptr;
0396     }
0397     Vertex vtx, tmp, norm;
0398     for( long j=0, nvx=0, n=tes->GetNfacets(); j < n; ++j )  {
0399       bool degenerated  = false;
0400       const auto& facet = tes->GetFacet(j);
0401 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
0402       tmp = tes->FacetComputeNormal(j, degenerated);
0403 #else
0404       tmp = facet.ComputeNormal(degenerated);
0405 #endif
0406       if ( !degenerated && facet.GetNvert() > 0 )   {
0407         aiFace& face  = mesh->mFaces[mesh->mNumFaces];
0408         double  u     = unit_scale;
0409
0410         face.mIndices = new unsigned int[facet.GetNvert()];
0411         trafo->LocalToMaster(tmp.fVec, norm.fVec);
0412         face.mNumIndices = 0;
0413         for( long k=0; k < facet.GetNvert(); ++k )  {
0414 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
0415           tmp = tes->GetVertex(facet[k]);
0416 #else
0417           tmp = facet.GetVertex(k);
0418 #endif
0419           trafo->LocalToMaster(tmp.fVec, vtx.fVec);
0420           face.mIndices[face.mNumIndices] = nvx;
0421           mesh->mNormals[nvx]  = aiVector3D(norm.x(), norm.y(), norm.z());
0422           mesh->mVertices[nvx] = aiVector3D(vtx.x()*u,vtx.y()*u,vtx.z()*u);
0423           ++mesh->mNumVertices;
0424           ++face.mNumIndices;
0425           ++nvx;
0426         }
0427         ++mesh->mNumFaces;
0428         if ( dump_facets )   {
0429           const auto* id = face.mIndices;
0430           const auto* vv = mesh->mVertices;
0431           ROOT::Geom::Vertex_t v1(vv[id[0]].x, vv[id[0]].y, vv[id[0]].z);
0432           ROOT::Geom::Vertex_t v2(vv[id[1]].x, vv[id[1]].y, vv[id[1]].z);
0433           ROOT::Geom::Vertex_t v3(vv[id[2]].x, vv[id[2]].y, vv[id[2]].z);
0434           std::string str = dd4hep::cad::streamVertices(v1, v2, v3);
0435           printout(ALWAYS,"ASSIMPWriter","++ Facet %4ld : %s", j, str.c_str());
0436         }
0437       }
0438       else   {
0439         printout(ALWAYS,"ASSIMPWriter",
0440                  "+++ Found degenerate facet while writing [Should not happen]");
0441       }
0442     }
0443
0444     /// Check if we have here a valid mesh
0445     if ( 0 == mesh->mNumFaces || 0 == mesh->mNumVertices )    {
0446       if ( mesh->mVertices ) delete [] mesh->mVertices;
0447       mesh->mVertices = nullptr;
0448       mesh->mNumVertices = 0;
0449       if ( mesh->mNormals ) delete [] mesh->mNormals;
0450       mesh->mNormals = nullptr;
0451       if ( mesh->mFaces ) delete [] mesh->mFaces;
0452       mesh->mFaces = nullptr;
0453       mesh->mNumFaces = 0;
0454       delete mesh;
0455       continue;
0456     }
0457
0458     scene.mMeshes[scene.mNumMeshes] = mesh;
0459
0460     aiNode *node      = new aiNode;
0461     node->mMeshes     = new unsigned int[node->mNumMeshes=1];
0462     node->mMeshes[0]  = scene.mNumMeshes;
0463     node->mParent     = root;
0464     node->mName.Set("<STL>");
0465
0466     root->mChildren[root->mNumChildren] = node;
0467     ++root->mNumChildren;
0468     ++scene.mNumMeshes;
0469   }
0470   TGeoShape::SetTransform(geo_transform);
0471   printout(ALWAYS,"ASSIMPWriter","+++ Analysed %ld out of %ld meshes.",
0472            scene.mNumMeshes, placements.size());
0473   if ( scene.mNumMeshes > 0 )   {
0474     Assimp::Exporter exporter;
0475     Assimp::ExportProperties *props = new Assimp::ExportProperties;
0476     props->SetPropertyBool(AI_CONFIG_EXPORT_POINT_CLOUDS, flags&EXPORT_POINT_CLOUDS ? true : false);
0477     exporter.Export(&scene, file_type.c_str(), file_name.c_str(), 0, props);
0478     return 1;
0479   }
0480   return 0;
0481 }