eic-opticks/CSG/CSGQuery.cc

0001 #include "SLOG.hh"
0002 #include "SSys.hh"
0003 #include "SPath.hh"
0004 #include "sc4u.h"
0005
0006 #include "CSGFoundry.h"
0007 #include "SName.h"
0008 #include "CSGQuery.h"
0009 #include "CSGGrid.h"
0010
0011 #ifdef DEBUG_RECORD
0012 #include "CSGRecord.h"
0013 #endif
0014
0015 #ifdef DEBUG_CYLINDER
0016 #include "CSGDebug_Cylinder.hh"
0017 #endif
0018
0019
0020 #include "OpticksCSG.h"
0021
0022
0023
0024 #include "csg_intersect_leaf.h"
0025 #include "csg_intersect_node.h"
0026 #include "csg_intersect_tree.h"
0027
0028
0029 const plog::Severity CSGQuery::LEVEL = SLOG::EnvLevel("CSGQuery", "DEBUG") ;
0030
0031 const int CSGQuery::VERBOSE = SSys::getenvint("VERBOSE", 0);
0032
0033
0034 CSGQuery::CSGQuery( const CSGFoundry* fd_ )
0035     :
0036     fd(fd_),
0037     prim0(fd->getPrim(0)),
0038     node0(fd->getNode(0)),
0039     plan0(fd->getPlan(0)),
0040     itra0(fd->getItra(0)),
0041     select_prim(nullptr),
0042     select_nodeOffset(0),
0043     select_prim_numNode(0),
0044     select_root_node(nullptr),   // set by selectPrim
0045     select_root_typecode(CSG_ZERO),
0046     select_root_subNum(0),
0047     select_is_tree(true)
0048 {
0049     init();
0050 }
0051
0052 /**
0053 CSGQuery::init
0054 ----------------
0055
0056 MOI meshName:meshOrdinal:instanceIndex
0057    selects by meshname irrespective of which solid it appears in
0058
0059 SOPR
0060    selects more formally by solidIdx:primIdxRel
0061
0062 **/
0063
0064 void CSGQuery::init()
0065 {
0066     // select first prim of first solid unless SOPR overrides that
0067     int solidIdx = 0 ;
0068     int primIdxRel = 0 ;
0069
0070     const char* sopr = SSys::getenvvar("SOPR", "0:0" );
0071     SName::ParseSOPR(solidIdx, primIdxRel, sopr );
0072
0073     LOG(LEVEL) << " sopr " << sopr << " solidIdx " << solidIdx << " primIdxRel " << primIdxRel ;
0074     selectPrim(solidIdx, primIdxRel );
0075 }
0076
0077 void CSGQuery::selectPrim(unsigned solidIdx, unsigned primIdxRel )
0078 {
0079     const CSGPrim* pr = fd->getSolidPrim(solidIdx, primIdxRel);
0080     assert( pr );
0081     selectPrim(pr);
0082 }
0083
0084 /**
0085 CSGQuery::selectPrim
0086 ---------------------
0087
0088 Recall that each CSGSolid have multiple CSGPrim corresponding
0089 in Geant4 to G4VSolid root nodes.  Also each CSGPrim generally has
0090 multiple CSGNode corresponding to the CSG constituent G4VSolid.
0091
0092 **/
0093
0094 void CSGQuery::selectPrim( const CSGPrim* pr )
0095 {
0096     select_prim = pr ;
0097     select_prim_ce = pr->ce();
0098     select_nodeOffset = pr->nodeOffset() ;
0099     select_prim_numNode = pr->numNode() ;
0100     select_root_node = node0 + pr->nodeOffset() ;
0101     select_root_typecode = select_root_node->typecode();
0102     select_root_subNum = select_root_node->subNum();
0103
0104
0105
0106     select_is_tree = CSG::IsTree((OpticksCSG_t)select_root_typecode) ;
0107
0108     LOG_IF(info, VERBOSE > 0)
0109         << " select_prim " << select_prim
0110         << " select_nodeOffset " << select_nodeOffset
0111         << " select_prim_numNode " << select_prim_numNode
0112         << " select_root_node " << select_root_node
0113         << " select_root_typecode " << CSG::Name(select_root_typecode)
0114         << " select_root_subNum " << select_root_subNum
0115         << " getSelectedTreeHeight " << getSelectedTreeHeight()
0116         ;
0117
0118     LOG_IF(LEVEL, select_root_subNum == 0) << "select_root_subNum ZERO " ;  // hmm formerly fatal ? why ?
0119
0120 }
0121
0122 int CSGQuery::getSelectedType() const
0123 {
0124     return select_root_typecode ;
0125 }
0126 int CSGQuery::getSelectedTreeHeight() const
0127 {
0128     return select_is_tree ? TREE_HEIGHT(select_root_subNum) : -1 ;
0129 }
0130
0131 /**
0132 CSGQuery::getSelectedNode
0133 ---------------------------
0134
0135 **/
0136
0137 const CSGNode* CSGQuery::getSelectedNode( int nodeIdx ) const
0138 {
0139     const CSGNode* nd = nodeIdx < select_prim_numNode ? select_root_node + nodeIdx : nullptr  ;
0140     return nd ;
0141 }
0142
0143
0144 /**
0145 CSGQuery::getSelectedTreeNode
0146 ------------------------------
0147
0148 nodeIdx
0149    0-based tree index, root=0
0150
0151 **/
0152
0153 const CSGNode* CSGQuery::getSelectedTreeNode( int nodeIdx ) const
0154 {
0155     const CSGNode* nd = nodeIdx < select_root_subNum ? select_root_node + nodeIdx : nullptr  ;
0156     return nd ;
0157 }
0158
0159
0160 float CSGQuery::distance(const float3& position ) const
0161 {
0162     float sd = distance_prim( position, select_root_node, plan0, itra0 ) ;
0163     return sd ;
0164 }
0165
0166 float CSGQuery::operator()(const float3& position ) const
0167 {
0168     return distance(position);
0169 }
0170
0171 /**
0172 CSGQuery::intersect
0173 ----------------------
0174
0175 This is invoked by CSGGeometry::saveCenterExtentGenstepIntersect using photons
0176 generated by SCenterExtentGenstep
0177
0178 **/
0179
0180 bool CSGQuery::intersect( quad4& isect,  float t_min, const quad4& p ) const
0181 {
0182     // HMM: NOT USING SIMTRACE LAYOUT
0183     float3 ray_origin    = make_float3(  p.q0.f.x, p.q0.f.y, p.q0.f.z );    // HMM: this could use sphoton approach avoiding temporaries
0184     float3 ray_direction = make_float3(  p.q1.f.x, p.q1.f.y, p.q1.f.z );
0185     unsigned gsid = p.q3.u.w ;
0186     return intersect( isect, t_min, ray_origin, ray_direction, gsid );
0187 }
0188
0189
0190
0191
0192
0193
0194 /**
0195 CSGQuery::intersect
0196 --------------------
0197
0198 CPU exercise of CUDA intersect code
0199
0200 isect.q0.f.xyx
0201     surface normal at the intersect
0202
0203 isect.q0.f.w
0204     t, ray_direction multiple at the intersect
0205
0206 isect.q1.f.xyz
0207     intersect position
0208
0209 isect.q1.f.w
0210     surface distance at intersect position.
0211
0212     This is expected to be close to zero,
0213     deviations from zero can be used to identify some
0214     forms of spurious intersects.
0215
0216     Note however that the most common cause of spurious intersects,
0217     coincident faces, does not show up this way as the surface distance is
0218     zero for both constituents across the common face.
0219
0220 isect.q2.f.xyz
0221     ray_origin
0222
0223 isect.q3.f.xyz
0224     ray_direction
0225
0226
0227 TODO: compare the isect with simtrace, it would be convenient to have same layout
0228
0229 **/
0230
0231 bool CSGQuery::intersect( quad4& isect,  float t_min, const float3& ray_origin, const float3& ray_direction, unsigned gsid ) const
0232 {
0233     isect.zero();
0234
0235     isect.q2.f.x = ray_origin.x ;
0236     isect.q2.f.y = ray_origin.y ;
0237     isect.q2.f.z = ray_origin.z ;
0238     isect.q2.f.w = t_min ;
0239
0240     isect.q3.f.x = ray_direction.x ;
0241     isect.q3.f.y = ray_direction.y ;
0242     isect.q3.f.z = ray_direction.z ;
0243     isect.q3.u.w = gsid ;
0244
0245 #ifdef DEBUG
0246     std::cout << "CSGQuery::intersect  ray_origin " << ray_origin << " ray_direction " << ray_direction << std::endl ;
0247 #endif
0248
0249     bool dump = false ;
0250     bool valid_intersect = intersect_prim( isect.q0.f, select_root_node, plan0, itra0, t_min, ray_origin, ray_direction, dump ) ;
0251     if( valid_intersect )
0252     {
0253         float t = isect.q0.f.w ;
0254         float3 ipos = ray_origin + t*ray_direction ;
0255
0256         float sd = distance(ipos) ;
0257
0258         isect.q1.f.x = ipos.x ;
0259         isect.q1.f.y = ipos.y ;
0260         isect.q1.f.z = ipos.z ;
0261         isect.q1.f.w = sd ;     // HMM: overwrite of tmin ?
0262     }
0263     return valid_intersect ;
0264 }
0265
0266
0267
0268
0269
0270
0271 /**
0272 CSGQuery::distance
0273 -------------------
0274
0275 Used for debugging distance issues
0276
0277 **/
0278
0279 void CSGQuery::distance( quad4& isect,  const float3& ray_origin ) const
0280 {
0281     isect.zero();
0282     isect.q1.f.w = distance(ray_origin) ;  // HMM thats an overwrite
0283
0284     isect.q2.f.x = ray_origin.x ;
0285     isect.q2.f.y = ray_origin.y ;
0286     isect.q2.f.z = ray_origin.z ;
0287     isect.q2.f.w = 0.f ;
0288 }
0289
0290
0291 const float CSGQuery::SD_CUT = -1e-3 ;
0292
0293 bool CSGQuery::IsSpurious( const quad4& isect ) // static
0294 {
0295     float sd = isect.q1.f.w ;
0296     bool spurious = sd < SD_CUT ;
0297     return spurious ;
0298 }
0299
0300
0301 std::string CSGQuery::Label() // static
0302 {
0303     std::stringstream ss ;
0304     ss << "CSGQuery::Label " ;
0305
0306 #ifdef DEBUG
0307     ss << " DEBUG" ;
0308 #else
0309     ss << " not-DEBUG" ;
0310 #endif
0311
0312 #ifdef DEBUG_RECORD
0313     ss << " DEBUG_RECORD" ;
0314 #else
0315     ss << " not-DEBUG_RECORD" ;
0316 #endif
0317
0318 #ifdef DEBUG_CYLINDER
0319     ss << " DEBUG_CYLINDER" ;
0320 #else
0321     ss << " not-DEBUG_CYLINDER" ;
0322 #endif
0323
0324
0325
0326     std::string s = ss.str() ;
0327     return s ;
0328 }
0329
0330
0331 /**
0332
0333
0334
0335 **/
0336
0337
0338 std::string CSGQuery::Desc( const quad4& isect, const char* label, bool* valid_intersect  )  // static
0339 {
0340     float sd = isect.q1.f.w ;
0341     bool spurious = IsSpurious(isect);
0342     std::stringstream ss ;
0343     ss
0344          << std::setw(30) << label
0345          << " "
0346          << ( valid_intersect ? ( *valid_intersect ? "HIT" : "MISS" )  : " " )
0347          << std::endl
0348          << std::setw(30) << " q0 norm t "
0349          << "("
0350          << std::setw(10) << std::fixed << std::setprecision(4) << isect.q0.f.x
0351          << std::setw(10) << std::fixed << std::setprecision(4) << isect.q0.f.y
0352          << std::setw(10) << std::fixed << std::setprecision(4) << isect.q0.f.z
0353          << std::setw(10) << std::fixed << std::setprecision(4) << isect.q0.f.w
0354          << ")"
0355          << std::endl
0356          << std::setw(30) << " q1 ipos sd "
0357          << "("
0358          << std::setw(10) << std::fixed << std::setprecision(4) << isect.q1.f.x
0359          << std::setw(10) << std::fixed << std::setprecision(4) << isect.q1.f.y
0360          << std::setw(10) << std::fixed << std::setprecision(4) << isect.q1.f.z
0361          << std::setw(10) << std::fixed << std::setprecision(4) << sd
0362          << ")"
0363          << ( spurious ? " SPURIOUS INTERSECT " : "-" )
0364          << " sd < SD_CUT : "
0365          << std::setw(10) << std::fixed << std::setprecision(4) <<  SD_CUT
0366          << std::endl
0367          << std::setw(30) << " q2 ray_ori t_min "
0368          << "("
0369          << std::setw(10) << std::fixed << std::setprecision(4) << isect.q2.f.x
0370          << std::setw(10) << std::fixed << std::setprecision(4) << isect.q2.f.y
0371          << std::setw(10) << std::fixed << std::setprecision(4) << isect.q2.f.z
0372          //<< std::setw(10) << std::fixed << std::setprecision(4) << isect.q2.f.w        // problem as simtrace has boundary here
0373          << ")"
0374          << std::endl
0375          << std::setw(30) << " q3 ray_dir gsid "
0376          << "("
0377          << std::setw(10) << std::fixed << std::setprecision(4) << isect.q3.f.x
0378          << std::setw(10) << std::fixed << std::setprecision(4) << isect.q3.f.y
0379          << std::setw(10) << std::fixed << std::setprecision(4) << isect.q3.f.z
0380          << std::setw(20) << C4U_desc(isect.q3.u.w)
0381          << ")"
0382          << std::endl
0383          ;
0384
0385     std::string s = ss.str() ;
0386     return s ;
0387 }
0388
0389 /**
0390 CSGQuery::simtrace
0391 --------------------
0392
0393 **/
0394
0395 bool CSGQuery::simtrace( quad4& p ) const
0396 {
0397     float3* ray_origin = p.v2() ;
0398     float3* ray_direction = p.v3() ;
0399     float t_min = p.q1.f.w ;
0400
0401     // the 1st float4 argumnent gives surface normal at intersect and distance
0402     bool dump = false ;
0403     bool valid_intersect = intersect_prim( p.q0.f, select_root_node, plan0, itra0, t_min, *ray_origin, *ray_direction, dump ) ;
0404     if( valid_intersect )
0405     {
0406         float t = p.q0.f.w ;
0407         float3 ipos = (*ray_origin) + t*(*ray_direction) ;
0408         p.q1.f.x = ipos.x ;
0409         p.q1.f.y = ipos.y ;
0410         p.q1.f.z = ipos.z ;
0411         //p.q1.f.w = distance(ipos) ;     // HMM: overwrite of tmin is problematic
0412     }
0413     return valid_intersect ;
0414 }
0415
0416
0417 void CSGQuery::post(const char* outdir)
0418 {
0419 #ifdef DEBUG_CYLINDER
0420     if(outdir)
0421     {
0422         const std::vector<CSGDebug_Cylinder>& record = CSGDebug_Cylinder::record ;
0423         LOG(info) << " CSGDebug_Cylinder::record.size " << record.size() << " outdir " << outdir ;
0424         CSGDebug_Cylinder::Save(outdir);
0425     }
0426 #endif
0427 }
0428
0429
0430
0431 /**
0432 CSGQuery::intersect_again
0433 --------------------------
0434
0435
0436 **/
0437 bool CSGQuery::intersect_again( quad4& isect, const quad4& prev ) const
0438 {
0439     float3* ray_origin = prev.v2() ;
0440     float3* ray_direction = prev.v3() ;
0441     float t_min = prev.q1.f.w ;    // t_min was formerly  prev.q2.f.w
0442
0443     /*
0444     std::cout
0445         << "CSGQuery::intersect_again"
0446         << " ray_origin " << *ray_origin
0447         << " ray_direction " << *ray_direction
0448         << " t_min " << t_min
0449         << std::endl
0450         ;
0451
0452     */
0453
0454     //if(t_min != 0.f) LOG(error) << " t_min " << t_min ;  // not zero eg from propagate epsilon 0.05
0455     //assert( t_min == 0.f );
0456     unsigned gsid = prev.q3.u.w ;
0457
0458 #ifdef DEBUG_RECORD
0459     CSGRecord::SetEnabled(true);
0460 #endif
0461
0462     bool valid_intersect = intersect( isect, t_min, *ray_origin, *ray_direction, gsid );
0463
0464 #ifdef DEBUG_RECORD
0465     CSGRecord::SetEnabled(false);
0466 #endif
0467
0468    /*
0469     LOG(info)
0470         << std::endl
0471         << Desc(prev,  "prev_isect" )
0472         << Desc(isect, "isect" )
0473         ;
0474
0475     */
0476
0477     return valid_intersect ;
0478 }
0479
0480
0481 /**
0482 CSGQuery::scanPrim
0483 --------------------
0484
0485 Used by sdf_geochain.sh CSGGeometry::saveSignedDistanceField
0486
0487 **/
0488
0489 CSGGrid* CSGQuery::scanPrim(int resolution) const
0490 {
0491     const CSGPrim* pr = select_prim ;
0492     if( pr == nullptr )
0493     {
0494         LOG(fatal) << " no prim is selected " ;
0495         return nullptr ;
0496     }
0497
0498     const float4 ce =  pr->ce() ;
0499     LOG(info) << " ce " << ce << " resolution " << resolution ;
0500
0501     CSGGrid* grid = new CSGGrid( ce, resolution, resolution, resolution );
0502     grid->scan(*this) ;
0503     return grid ;
0504 }
0505
0506
0507
0508 std::string CSGQuery::descPrim() const
0509 {
0510     std::stringstream ss ;
0511     for(int nodeIdx=select_nodeOffset ; nodeIdx < select_nodeOffset+select_prim_numNode ; nodeIdx++)
0512     {
0513         const CSGNode* nd = node0 + nodeIdx ;
0514         ss << nd->desc() << std::endl ;
0515     }
0516     std::string s = ss.str();
0517     return s ;
0518 }
0519
0520
0521 void CSGQuery::dumpPrim(const char* msg) const
0522 {
0523     if( select_prim_numNode == 0 ) return ;
0524     LOG(LEVEL)
0525           << msg
0526           << " select_prim_numNode " << select_prim_numNode
0527           << " select_nodeOffset " << select_nodeOffset
0528           << std::endl
0529           << descPrim()
0530           ;
0531
0532 }
0533
0534 void CSGQuery::dump(const char* msg) const
0535 {
0536     LOG(info) << msg ;
0537     dumpPrim();
0538 }
0539
0540