eic-opticks/CSG/CSGGeometry.cc

0001 #include <csignal>
0002
0003 #include "SLOG.hh"
0004 #include "ssys.h"
0005 #include "spath.h"
0006 #include "sdirectory.h"
0007 #include "scuda.h"
0008 #include "squad.h"
0009 #include "sc4u.h"
0010
0011 #include "SCenterExtentGenstep.hh"
0012 #include "NP.hh"
0013
0014 #include "CSGQuery.h"
0015 #include "CSGDraw.h"
0016 #include "CSGFoundry.h"
0017 #include "CSGMaker.h"
0018 #include "CSGGeometry.h"
0019 #include "CSGRecord.h"
0020 #include "CSGGrid.h"
0021
0022 const plog::Severity CSGGeometry::LEVEL = SLOG::EnvLevel("CSGGeometry", "DEBUG") ;
0023
0024
0025 void CSGGeometry::operator()()
0026 {
0027     LOG(info) ;
0028
0029     bool dump_ = ssys::getenvbool("DUMP");
0030     if(dump_)
0031     {
0032         LOG(info) << " dump  " ;
0033         dump();
0034         LOG(info) << desc();
0035
0036     }
0037     else
0038     {
0039         LOG(info) << " centerExtentGenstepIntersect " ;
0040         centerExtentGenstepIntersect();
0041     }
0042 }
0043
0044
0045 const char* CSGGeometry::OutDir(const char* cfbase, const char* geom, const char* sopr)   // static
0046 {
0047     if(cfbase == nullptr )
0048     {
0049         LOG(LEVEL) << " require CFBASE to control output dir " ;
0050         return nullptr ;
0051     }
0052
0053     const char* reldir = geom ? geom : sopr ;
0054
0055     if(reldir == nullptr )
0056     {
0057         LOG(fatal) << " require geom or sopr to control reldir included in output dir " ;
0058         return nullptr ;
0059     }
0060
0061     const char* outdir = spath::Resolve(cfbase, "CSGIntersectSolidTest", reldir );
0062     sdirectory::MakeDirs(outdir);
0063
0064     return outdir ;
0065 }
0066
0067 /**
0068 CSGGeometry::CSGGeometry
0069 -------------------------
0070
0071 Can boot in three ways:
0072
0073 1. externally created CSGFoundry instance
0074 2. GEOM envvar identifying a test geometry, and resulting in creation of a CSGFoundry instance
0075 3. CFBASE envvar identifying directory containing a persisted CSGFoundry geometry that is loaded
0076
0077
0078 **/
0079
0080 CSGGeometry::CSGGeometry(const char* default_cfbase, const CSGFoundry* fd_)
0081     :
0082     default_geom(nullptr),
0083     default_sopr(nullptr),
0084     geom(ssys::getenvvar("CSGGeometry_GEOM", default_geom)),
0085     sopr(ssys::getenvvar("SOPR", default_sopr)),
0086     cfbase(ssys::getenvvar("CFBASE", default_cfbase)),
0087     outdir(OutDir(cfbase,geom,sopr)),
0088     name(nullptr),
0089     fd(fd_),
0090     q(nullptr),
0091     ce(nullptr),
0092     d(nullptr),
0093     sxyzw(ssys::getenvintvec("SXYZW",',')),
0094     sxyz(ssys::getenvintvec("SXYZ",',')),
0095     no_selection(sxyzw == nullptr && sxyz == nullptr),
0096     sx(0),
0097     sy(0),
0098     sz(0),
0099     sw(0),
0100     rc(0)
0101 {
0102     init();
0103 }
0104
0105 void CSGGeometry::init()
0106 {
0107     LOG(LEVEL) << " CSGGeometry_GEOM " << geom  ;
0108     init_fd();
0109
0110     if(fd == nullptr)
0111     {
0112         LOG(fatal) << " ABORT " ;
0113         rc = 101 ;
0114         return ;
0115     }
0116
0117     q = new CSGQuery(fd);
0118     q->dumpPrim("CSGGeometry::init");
0119     ce = new float4(q->select_prim->ce()) ;
0120     d = new CSGDraw(q, 'Z') ;
0121     init_selection();
0122 }
0123
0124 void CSGGeometry::init_fd()
0125 {
0126     if( fd == nullptr )
0127     {
0128         if( geom != nullptr && cfbase == nullptr)
0129         {
0130             name = strdup(geom);
0131             LOG(info) << "init from GEOM " << geom << " name " << name ;
0132             fd = CSGMaker::MakeGeom(geom) ;
0133         }
0134         else if(cfbase != nullptr)
0135         {
0136             name = spath::Basename(cfbase);
0137             LOG(info) << "init from CFBASE " << cfbase << " name " << name  ;
0138             fd = CSGFoundry::Load(cfbase, "CSGFoundry");
0139             LOG(info) << " fd.meta\n" << ( fd->hasMeta() ? fd->meta : " NO meta " ) ;
0140         }
0141         else
0142         {
0143             LOG(fatal) << " neither GEOM or CFBASE envvars are defined and fd pointer not provided : ABORT " ;
0144             //std::raise(SIGINT);
0145         }
0146     }
0147     else
0148     {
0149         LOG(LEVEL) << " booting from provided CSGFoundry pointer " ;
0150         cfbase = fd->getCFBase();
0151     }
0152 }
0153
0154 void CSGGeometry::init_selection()
0155 {
0156     if(sxyz)
0157     {
0158         bool sxyz_expect = sxyz->size() == 3 ;
0159         if(!sxyz_expect)
0160         {
0161             LOG(fatal) << "SXYZ envvar is provided but with size: " << sxyz->size() << " when 3 is expected" ;
0162             assert(0);
0163         }
0164         sx = (*sxyz)[0] ;
0165         sy = (*sxyz)[1] ;
0166         sz = (*sxyz)[2] ;
0167
0168         LOG(info) << "SXYZ (sx,sy,sz) " << "(" << sx << "," << sy << "," << sz << ")" ;
0169     }
0170     else if(sxyzw)
0171     {
0172         bool sxyzw_expect = sxyzw->size() == 4 ;
0173         if(!sxyzw_expect)
0174         {
0175             LOG(fatal) << "SXYZW envvar is provided but with size: " << sxyzw->size() << " when 4 is expected" ;
0176             assert(0);
0177         }
0178         sx = (*sxyzw)[0] ;
0179         sy = (*sxyzw)[1] ;
0180         sz = (*sxyzw)[2] ;
0181         sw = (*sxyzw)[3] ;
0182         LOG(info) << "SXYZW (sx,sy,sz,sw) " << "(" << sx << "," << sy << "," << sz << "," << sw << ")" ;
0183     }
0184     else
0185     {
0186         LOG(LEVEL) << " no SXYZ or SXYZW selection " ;
0187     }
0188 }
0189
0190
0191
0192 void CSGGeometry::saveSignedDistanceField() const
0193 {
0194     LOG_IF(error, q == nullptr ) << " q null : ABORT " ;
0195     if(!q) return ;
0196
0197     int resolution = ssys::getenvint("RESOLUTION", 25) ;
0198     LOG(info) << " name " << name << " RESOLUTION " << resolution ;
0199
0200     q->dumpPrim();
0201
0202     const CSGGrid* grid = q->scanPrim(resolution);
0203     assert( grid );
0204     grid->save(name);
0205 }
0206
0207
0208 /**
0209 CSGGeometry::centerExtentGenstepIntersect
0210 -------------------------------------------
0211
0212 When SELECT_ISECT envvar is set to the path of an isect.npy
0213 those intersects are rerun. This is convenient for debugging
0214 a selection of intersects.
0215
0216 **/
0217
0218 void CSGGeometry::centerExtentGenstepIntersect()
0219 {
0220     const char* path = ssys::getenvvar("SELECTED_ISECT");
0221     if(path == nullptr)
0222     {
0223         float t_min = ssys::getenvfloat("TMIN", 0.f );
0224         saveCenterExtentGenstepIntersect(t_min);
0225     }
0226     else
0227     {
0228         intersectSelected(path);
0229     }
0230
0231     LOG(info) << desc() ;
0232 }
0233
0234
0235 void CSGGeometry::saveCenterExtentGenstepIntersect(float t_min) const
0236 {
0237     //float4 ce = make_float4(0.f, 0.f, 0.f, 70.f );
0238     SCenterExtentGenstep* cegs = new SCenterExtentGenstep(ce) ;
0239     const std::vector<quad4>& pp = cegs->pp ;
0240     std::vector<quad4>& ii = cegs->ii ;
0241
0242     LOG(info) << "[ pp.size " << pp.size() << " t_min " << std::fixed << std::setw(10) << std::setprecision(4) << t_min  ;
0243
0244     C4U gsid ;
0245     quad4 isect ;
0246     unsigned num_ray(0);
0247
0248     for(unsigned i=0 ; i < pp.size() ; i++)
0249     {
0250         const quad4& p = pp[i];
0251         gsid.u = p.q3.u.w ;
0252
0253         int ix = gsid.c4.x ;
0254         int iy = gsid.c4.y ;
0255         int iz = gsid.c4.z ;
0256         int iw = gsid.c4.w ;
0257
0258         if( sxyz == nullptr && sxyzw == nullptr )   // no restriction
0259         {
0260             num_ray += 1 ;
0261             if(q->intersect(isect, t_min, p )) ii.push_back(isect);
0262         }
0263         else if( sxyz && ix == sx && iy == sy && iz == sz )    // restrict to single genstep
0264         {
0265             num_ray += 1 ;
0266             if(q->intersect(isect, t_min, p )) ii.push_back(isect);
0267         }
0268         else if( sxyzw && ix == sx && iy == sy && iz == sz && iw == sw )  // restrict to single photon
0269         {
0270             num_ray += 1 ;
0271             LOG(info) << "[ single photon selected" ;
0272             bool valid_isect = q->intersect(isect, t_min, p ) ;
0273             if(valid_isect) ii.push_back(isect);
0274             LOG(info) << std::endl << CSGQuery::Desc(isect, "single photon selected", &valid_isect ) ;
0275
0276 #ifdef DEBUG_RECORD
0277             std::string gsid_name = C4U_name( gsid.u , "gsid", '_' );
0278             const char* dir = spath::Resolve(outdir, "saveCenterExtentGenstepIntersect", gsid_name.c_str() );
0279             sdirectory::MakeDirs(dir,0);
0280
0281             LOG(info) << " DEBUG_RECORD saving to " << dir ;
0282             CSGRecord::Save(dir);
0283             CSGRecord::Dump(dir);
0284             CSGRecord::Clear();
0285 #else
0286             LOG(info) << " DEBUG_RECORD flag not enabled " ;
0287 #endif
0288             LOG(info) << "] single photon selected " ;
0289         }
0290     }
0291
0292     unsigned num_isect = ii.size() ;
0293     LOG(info)
0294         << " pp.size " << pp.size()
0295         << " num_ray " << num_ray
0296         << " ii.size " << num_isect
0297         << ( num_isect == 0 ? " WARNING : NO INTERSECTS " : " " )
0298         ;
0299
0300     if( no_selection )
0301     {
0302         LOG(info) << " saving to outdir " << outdir ;
0303         cegs->save(outdir);
0304     }
0305     else
0306     {
0307         LOG(info) << " NOT saving as an SXYZ or SXYZW selection is in use " ;
0308
0309     }
0310
0311     LOG(info) << "]" ;
0312 }
0313
0314 void CSGGeometry::intersectSelected(const char* path)
0315 {
0316     NP* a = NP::Load(path);
0317     LOG_IF(fatal, a == nullptr) << " FAILED to load from path " << path ;
0318     if( a == nullptr ) return ;
0319
0320     bool expected_shape = a->has_shape(-1,4,4) ;
0321     assert(expected_shape);
0322     if(!expected_shape) std::raise(SIGINT);
0323
0324     LOG(info) << " load SELECTED_ISECT from " << path << " a.sstr " << a->sstr() ;
0325
0326     unsigned num_isect = a->shape[0] ;
0327     std::vector<quad4> isects(num_isect) ;
0328     memcpy( isects.data(), a->bytes(),  sizeof(float)*16 );
0329
0330     for(unsigned i=0 ; i < num_isect ; i++)
0331     {
0332         const quad4& prev_isect = isects[i] ;
0333
0334         int4 gsid ;
0335         C4U_decode(gsid, prev_isect.q3.u.w );
0336
0337         quad4 isect ;
0338         bool valid_intersect = q->intersect_again(isect, prev_isect );
0339
0340         std::string gsid_name = C4U_name(gsid, "gsid", '_' );
0341
0342 #ifdef DEBUG_RECORD
0343         const char* dir = spath::Resolve(outdir, "intersectSelected", gsid_name.c_str() );
0344         sdirectory::MakeDirs(dir,0);
0345
0346         CSGRecord::Save(dir);
0347         if(i == 0) CSGRecord::Dump(dir);
0348         CSGRecord::Clear();
0349 #endif
0350
0351         std::cout
0352             << " gsid " << C4U_desc(gsid)
0353             << " gsid_name " << gsid_name
0354             << " valid_intersect " << valid_intersect
0355             << " t:prev_isect.q0.f.w "
0356             << std::setw(10) << std::fixed << std::setprecision(4) << prev_isect.q0.f.w
0357             << " t:isect.q0.f.w "
0358             << std::setw(10) << std::fixed << std::setprecision(4) << isect.q0.f.w
0359             << std::endl
0360             ;
0361     }
0362 }
0363
0364 void CSGGeometry::dump(const char* msg) const
0365 {
0366     LOG_IF(fatal, !fd ) << " fd null : ABORT " ;
0367     if(!fd) return ;
0368
0369     LOG(error) << "fd.dumpSolid" ;
0370     fd->dumpSolid();
0371     LOG(error) << "fd.dumpPrim" ;
0372     fd->dumpPrim();
0373     LOG(error) << "fd.dumpNode" ;
0374     fd->dumpNode();
0375     LOG(error) << "q.dump" ;
0376     q->dump(msg);
0377 }
0378
0379
0380 std::string CSGGeometry::desc()
0381 {
0382     return d ? d->desc() :  "ERROR-d-null" ;
0383 }
0384
0385 std::string CSGGeometry::Desc( const CSGFoundry* fd ) // static
0386 {
0387     CSGGeometry cg(nullptr, fd);
0388     return cg.desc();
0389 }
0390