File indexing completed on 2026-04-09 07:49:37
0001 #include <limits>
0002
0003 #include "scuda.h"
0004 #include "squad.h"
0005 #include "sqat4.h"
0006 #include "stran.h"
0007
0008 #ifdef WITH_OLD_FRAME
0009 #include "sframe.h"
0010 #else
0011 #include "sfr.h"
0012 #endif
0013
0014 #include "sc4u.h"
0015 #include "ssincos.h"
0016 #include "ssys.h"
0017
0018 #include "SLOG.hh"
0019 #include "SRng.hh"
0020 #include "SGenstep.h"
0021 #include "OpticksGenstep.h"
0022 #include "SFrameGenstep.hh"
0023 #include "NP.hh"
0024
0025
0026 const plog::Severity SFrameGenstep::LEVEL = SLOG::EnvLevel("SFrameGenstep", "DEBUG" );
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041 void SFrameGenstep::CE_OFFSET(std::vector<float3>& ce_offset, const float4& ce )
0042 {
0043 const char* ekey = "CE_OFFSET" ;
0044 const char* val = ssys::getenvvar(ekey);
0045
0046 bool is_CE_string = val == nullptr ? false : ( strcmp(val, "CE")== 0 || strcmp(val, "ce")== 0 ) ;
0047 float3 offset = make_float3(0.f, 0.f, 0.f );
0048
0049
0050 if(is_CE_string)
0051 {
0052 offset.x = ce.x ;
0053 offset.y = ce.y ;
0054 offset.z = ce.z ;
0055 ce_offset.push_back(offset);
0056 }
0057 else
0058 {
0059 std::vector<float>* fvec = ssys::getenvfloatvec(ekey, "0,0,0");
0060 unsigned num_values = fvec->size() ;
0061 bool multiple_of_3 = num_values % 3 == 0 ;
0062 unsigned num_offset = num_values/3 ;
0063 assert(fvec);
0064
0065 LOG(LEVEL)
0066 << " ekey " << ekey
0067 << " num_values " << num_values
0068 << " multiple_of_3 " << multiple_of_3
0069 << " num_offset " << num_offset
0070 ;
0071
0072 LOG_IF(fatal, !multiple_of_3) << " not multiple_of_3 num_values " << num_values << " ekey " << ekey ;
0073 assert(multiple_of_3);
0074
0075 for(unsigned i=0 ; i < num_offset ; i++)
0076 {
0077 offset.x = (*fvec)[i*3+0] ;
0078 offset.y = (*fvec)[i*3+1] ;
0079 offset.z = (*fvec)[i*3+2] ;
0080 ce_offset.push_back(offset);
0081 }
0082 }
0083
0084 bool with_offset = ce_offset.size() > 0 ;
0085 LOG(LEVEL)
0086 << "ekey " << ekey
0087 << " val " << val
0088 << " is_CE_string " << is_CE_string
0089 << " ce_offset.size " << ce_offset.size()
0090 << " ce " << ce
0091 << " with_offset " << ( with_offset ? "YES" : "NO " )
0092 ;
0093
0094 LOG(LEVEL) << Desc(ce_offset) ;
0095
0096 LOG_IF(fatal, !with_offset) << " at least one ce_offset is required " ;
0097 assert(with_offset);
0098
0099
0100 }
0101
0102 std::string SFrameGenstep::Desc(const std::vector<float3>& ce_offset )
0103 {
0104 std::stringstream ss ;
0105 ss << "SFrameGenstep::Desc ce_offset.size " << ce_offset.size() << std::endl ;
0106 for(unsigned i=0 ; i < ce_offset.size() ; i++)
0107 {
0108 const float3& offset = ce_offset[i] ;
0109 ss << std::setw(4) << i << " : " << offset << std::endl ;
0110 }
0111
0112 std::string s = ss.str();
0113 return s ;
0114 }
0115
0116 std::string SFrameGenstep::Desc(const std::vector<int>& cegs )
0117 {
0118 std::stringstream ss ;
0119 ss << " size " << cegs.size() << "[" ;
0120 for(unsigned i=0 ; i < cegs.size() ; i++) ss << cegs[i] << " " ;
0121 ss << "]" ;
0122 std::string s = ss.str();
0123 return s ;
0124 }
0125
0126
0127
0128 void SFrameGenstep::GetGridConfig_(std::vector<int>& cegs, const char* ekey, char delim, const char* fallback )
0129 {
0130 ssys::getenvintvec(ekey, cegs, delim, fallback );
0131 LOG(LEVEL)
0132 << " ekey " << ( ekey ? ekey : "-" )
0133 << " fallback " << ( fallback ? fallback : "-" )
0134 << " delim " << delim
0135 << " cegs.size " << cegs.size()
0136 ;
0137
0138
0139 StandardizeCEGS(cegs);
0140 assert( cegs.size() == 0 || cegs.size() == 8 );
0141 LOG(LEVEL) << " ekey " << ekey << " Desc " << Desc(cegs) ;
0142 }
0143
0144 std::string SFrameGenstep::GetGridConfig(std::vector<int>& cegs)
0145 {
0146 const char* CEGS_fallback = CEGS_XZ ;
0147 GetGridConfig_(cegs, CEGS, CEGS_delim, CEGS_fallback );
0148
0149 std::stringstream ss ;
0150 assert( cegs.size() == 8 );
0151 for(int i=0 ; i < 8 ; i++ ) ss << cegs[i] << ( i < 8 - 1 ? "," : "" ) ;
0152 std::string str = ss.str();
0153 return str ;
0154 }
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179 const char* SFrameGenstep::CEGS_XY = "16:9:0:1000" ;
0180 const char* SFrameGenstep::CEGS_XZ = "16:0:9:1000" ;
0181
0182
0183 bool SFrameGenstep::HasConfigEnv()
0184 {
0185 return ssys::hasenv_("CEGS");
0186 }
0187
0188
0189
0190
0191
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
0208
0209
0210
0211
0212
0213
0214
0215
0216
0217
0218
0219
0220
0221
0222 #ifdef WITH_OLD_FRAME
0223 NP* SFrameGenstep::MakeCenterExtentGenstep_FromFrame(sframe& fr)
0224 {
0225 const float4& ce = fr.ce ;
0226 Tran<double>* geotran = fr.getTransform();
0227 char* _GRIDSCALE = getenv("GRIDSCALE") ;
0228 float gridscale = ssys::getenvfloat("GRIDSCALE", 0.1f ) ;
0229 int prim = -1 ;
0230
0231 LOG(LEVEL)
0232 << " _GRIDSCALE [" << ( _GRIDSCALE ? _GRIDSCALE : "-" ) << "]"
0233 << " GRIDSCALE " << gridscale
0234 ;
0235
0236 std::vector<int> cegs ;
0237 std::string str_cegs = GetGridConfig(cegs);
0238 fr.set_grid(cegs, gridscale);
0239
0240 NP* gs = MakeCenterExtentGenstep_From_CE_geotran( ce, cegs, gridscale, geotran, prim );
0241
0242 gs->set_meta<int>("midx", fr.midx() );
0243 gs->set_meta<int>("mord", fr.mord() );
0244 gs->set_meta<int>("gord", fr.gord() );
0245 gs->set_meta<float>("gridscale", fr.gridscale() );
0246 gs->set_meta<std::string>("cegs", str_cegs );
0247
0248 return gs ;
0249 }
0250 #else
0251 NP* SFrameGenstep::MakeCenterExtentGenstep_FromFrame(sfr& fr)
0252 {
0253 float4 ce = make_float4( fr.ce.x, fr.ce.y, fr.ce.z, fr.ce.w );
0254 Tran<double>* geotran = fr.getTransform();
0255 char* _GRIDSCALE = getenv("GRIDSCALE") ;
0256 float gridscale = ssys::getenvfloat("GRIDSCALE", 0.1f ) ;
0257 int prim = fr.get_prim() ;
0258
0259 LOG(LEVEL)
0260 << " _GRIDSCALE [" << ( _GRIDSCALE ? _GRIDSCALE : "-" ) << "]"
0261 << " GRIDSCALE " << gridscale
0262 ;
0263
0264
0265 std::vector<int> cegs ;
0266 std::string str_cegs = GetGridConfig(cegs);
0267 fr.set_gridscale(gridscale);
0268
0269 LOG(LEVEL)
0270 << " cegs.size " << cegs.size()
0271 ;
0272
0273
0274 NP* gs = MakeCenterExtentGenstep_From_CE_geotran( ce, cegs, gridscale, geotran, prim );
0275
0276 gs->set_meta<float>("gridscale", fr.get_gridscale() );
0277 gs->set_meta<std::string>("cegs", str_cegs );
0278
0279 return gs ;
0280 }
0281 #endif
0282
0283 NP* SFrameGenstep::MakeCenterExtentGenstep_From_CE_geotran(const float4& ce, const std::vector<int>& cegs, float gridscale, const Tran<double>* geotran, int prim)
0284 {
0285 std::vector<float>* cegs_radial_range = ssys::getenv_vec<float>(CEGS_RADIAL_RANGE, nullptr, CEGS_delim );
0286
0287 std::vector<float3> ce_offset ;
0288 CE_OFFSET(ce_offset, ce);
0289
0290 LOG(LEVEL)
0291 << " ce " << ce
0292 << " ce_offset.size " << ce_offset.size()
0293 ;
0294
0295 bool with_offset = ce_offset.size() > 0 ;
0296 LOG_IF(fatal, !with_offset) << "ce_offset vector of float3 needs at least one entry " ;
0297 assert( with_offset );
0298
0299 assert(!getenv("CE_SCALE") && "CE_SCALE ENVVAR IS NO LONGER USED :CHANGE YOUR SCRIPT" );
0300
0301 bool ce_scale_off = ssys::getenvbool("CE_SCALE_OFF") ;
0302
0303
0304
0305 std::vector<NP*> gsl ;
0306 NP* gs_base = MakeCenterExtentGenstep(ce, cegs, gridscale, geotran, ce_offset, !ce_scale_off, cegs_radial_range );
0307 gsl.push_back(gs_base) ;
0308
0309 std::vector<std::string> keys = {{
0310 "CEHIGH_0",
0311 "CEHIGH_1",
0312 "CEHIGH_2",
0313 "CEHIGH_3",
0314 "CEHIGH_4",
0315 "CEHIGH_5",
0316 "CEHIGH_6",
0317 "CEHIGH_7",
0318 "CEHIGH_8",
0319 "CEHIGH_9"
0320 }} ;
0321
0322 const char* CEHIGH_fallback = "" ;
0323
0324
0325 for(unsigned i=0 ; i < keys.size() ; i++)
0326 {
0327 const char* key = keys[i].c_str() ;
0328 std::vector<int> cehigh ;
0329 GetGridConfig_(cehigh, key, CEGS_delim, CEHIGH_fallback );
0330 LOG(LEVEL) << " key " << key << " cehigh.size " << cehigh.size() ;
0331 if(cehigh.size() == 8)
0332 {
0333 NP* gs_cehigh = MakeCenterExtentGenstep(ce, cehigh, gridscale, geotran, ce_offset, !ce_scale_off, cegs_radial_range );
0334 gsl.push_back(gs_cehigh) ;
0335 }
0336 }
0337
0338
0339 NP* CEGS_NPY = ssys::hasenv_("CEGS_NPY") ? NP::Load("$CEGS_NPY") : nullptr ;
0340 NP* gs_CEGS_NPY = CEGS_NPY ? Make_CEGS_NPY_Genstep(CEGS_NPY, geotran) : nullptr ;
0341 if(gs_CEGS_NPY) gsl.push_back(gs_CEGS_NPY);
0342
0343
0344 Maybe_Add_PRIOR_SIMTRACE_Genstep( gsl, prim );
0345
0346 LOG(LEVEL) << " gsl.size " << gsl.size() ;
0347
0348 LOG(LEVEL) << "[ NP::Concatenate " ;
0349 NP* gs = NP::Concatenate(gsl) ;
0350 LOG(LEVEL) << "] NP::Concatenate " ;
0351
0352 gs->set_meta<int>("ce_scale", int(!ce_scale_off) );
0353
0354 return gs ;
0355 }
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365
0366
0367
0368
0369
0370
0371
0372
0373
0374
0375
0376
0377
0378
0379
0380
0381
0382
0383
0384
0385
0386
0387
0388
0389
0390
0391
0392 void SFrameGenstep::Maybe_Add_PRIOR_SIMTRACE_Genstep( std::vector<NP*>& gsl, int prim )
0393 {
0394 bool with_PRIOR_SIMTRACE = ssys::hasenv_("PRIOR_SIMTRACE") ;
0395 NP* PRIOR_SIMTRACE = with_PRIOR_SIMTRACE ? NP::Load("$PRIOR_SIMTRACE") : nullptr ;
0396
0397 int PRIOR_SIMTRACE_PRIM = ssys::getenvint("PRIOR_SIMTRACE_PRIM", -2 );
0398 int uprim = prim ;
0399 if( PRIOR_SIMTRACE_PRIM == -2 )
0400 {
0401 uprim = prim ;
0402 }
0403 else if( PRIOR_SIMTRACE_PRIM == -1 )
0404 {
0405 uprim = -1 ;
0406 }
0407 else if( PRIOR_SIMTRACE_PRIM > -1 )
0408 {
0409 uprim = PRIOR_SIMTRACE_PRIM ;
0410 }
0411
0412 NP* gs_PRIOR_SIMTRACE = PRIOR_SIMTRACE ? Make_PRIOR_SIMTRACE_Genstep( PRIOR_SIMTRACE, uprim ) : nullptr ;
0413 if(gs_PRIOR_SIMTRACE) gsl.push_back(gs_PRIOR_SIMTRACE);
0414 LOG(info)
0415 << "with_PRIOR_SIMTRACE " << ( with_PRIOR_SIMTRACE ? "YES" : "NO " )
0416 << " prim " << prim
0417 << " uprim " << uprim
0418 << " PRIOR_SIMTRACE_PRIM " << PRIOR_SIMTRACE_PRIM
0419 << " PRIOR_SIMTRACE " << ( PRIOR_SIMTRACE ? PRIOR_SIMTRACE->sstr() : "-" )
0420 << " gs_PRIOR_SIMTRACE " << ( gs_PRIOR_SIMTRACE ? gs_PRIOR_SIMTRACE->sstr() : "-" )
0421 ;
0422 }
0423
0424
0425
0426
0427
0428
0429
0430
0431 NP* SFrameGenstep::Make_PRIOR_SIMTRACE_Genstep( const NP* _simtrace, int prim )
0432 {
0433 const NP* _ontoprim = prim > -1 ? quad4::select_prim( _simtrace, prim ) : _simtrace ;
0434 int ni = _ontoprim->shape[0] ;
0435 const quad4* ontoprim = (const quad4*)(_ontoprim->bytes());
0436
0437 quad6 gs ;
0438 gs.zero();
0439
0440 glm::tmat4x4<double> identity(1.);
0441 qat4* qc = Tran<double>::ConvertFrom( identity ) ;
0442
0443 std::vector<quad6> gensteps ;
0444
0445 float dx = 0.f;
0446 float dy = 0.f;
0447 float dz = 0.f;
0448
0449 float s = 1.f ;
0450
0451 for(int i=0 ; i < ni ; i++)
0452 {
0453 const quad4& q = ontoprim[i];
0454 float x = q.q1.f.x ;
0455 float y = q.q1.f.y ;
0456 float z = q.q1.f.z ;
0457
0458 for(int j=0 ; j < 6 ; j++)
0459 {
0460 switch(j)
0461 {
0462 case 0: dx = -s ; dy = 0.f ; dz = 0.f ; break ;
0463 case 1: dx = +s ; dy = 0.f ; dz = 0.f ; break ;
0464 case 2: dx = 0.f ; dy = -s ; dz = 0.f ; break ;
0465 case 3: dx = 0.f ; dy = +s ; dz = 0.f ; break ;
0466 case 4: dx = 0.f ; dy = 0.f ; dz = -s ; break ;
0467 case 5: dx = 0.f ; dy = 0.f ; dz = s ; break ;
0468 }
0469
0470 gs.q1.f.x = x + dx ;
0471 gs.q1.f.y = y + dy ;
0472 gs.q1.f.z = z + dz ;
0473 gs.q1.f.w = 1.f ;
0474
0475 qc->write(gs);
0476
0477 int gridaxes = 0 ;
0478 int gsid = 0 ;
0479 int photons_per_genstep = 100 ;
0480 SGenstep::ConfigureGenstep(gs, OpticksGenstep_FRAME, gridaxes, gsid, photons_per_genstep );
0481
0482 gensteps.push_back(gs);
0483 }
0484 }
0485
0486 return SGenstep::MakeArray(gensteps);
0487 }
0488
0489
0490
0491
0492
0493
0494
0495
0496
0497
0498
0499 NP* SFrameGenstep::Make_CEGS_NPY_Genstep( const NP* t, const Tran<double>* geotran )
0500 {
0501 assert( t && t->has_shape(-1,3) && t->uifc == 'f' && t->ebyte == 4 );
0502 int ni = t->shape[0] ;
0503 int nj = t->shape[1] ;
0504 assert( nj == 3 );
0505
0506 const float* tt = t->cvalues<float>();
0507
0508 qat4* qc = Tran<double>::ConvertFrom( geotran->t ) ;
0509
0510 quad6 gs ;
0511 gs.zero();
0512
0513 std::vector<quad6> gensteps ;
0514
0515 for(int i=0 ; i < ni ; i++)
0516 {
0517 gs.q1.f.x = tt[nj*i+0] ;
0518 gs.q1.f.y = tt[nj*i+1] ;
0519 gs.q1.f.z = tt[nj*i+2] ;
0520 gs.q1.f.w = 1.f ;
0521
0522 qc->write(gs);
0523
0524 int gridaxes = 0 ;
0525 int gsid = 0 ;
0526 int photons_per_genstep = 100 ;
0527 SGenstep::ConfigureGenstep(gs, OpticksGenstep_FRAME, gridaxes, gsid, photons_per_genstep );
0528
0529 gensteps.push_back(gs);
0530 }
0531 return SGenstep::MakeArray(gensteps);
0532 }
0533
0534
0535
0536
0537
0538
0539
0540
0541
0542
0543
0544
0545
0546
0547
0548
0549
0550
0551
0552
0553
0554
0555
0556
0557
0558
0559
0560
0561
0562
0563
0564
0565
0566
0567
0568
0569
0570
0571
0572
0573
0574
0575
0576
0577
0578
0579
0580
0581
0582
0583
0584
0585
0586
0587
0588
0589
0590
0591
0592
0593
0594
0595
0596
0597
0598
0599
0600
0601
0602
0603
0604 NP* SFrameGenstep::MakeCenterExtentGenstep(
0605 const float4& ce,
0606 const std::vector<int>& cegs,
0607 float gridscale,
0608 const Tran<double>* geotran,
0609 const std::vector<float3>& ce_offset,
0610 bool ce_scale,
0611 std::vector<float>* cegs_radial_range)
0612 {
0613
0614 quad6 gs ;
0615 gs.zero();
0616
0617 assert( cegs.size() == 8 );
0618
0619 int high = cegs[7] ;
0620 assert( high >= 1 && high <= 8);
0621 double scale = double(gridscale)/double(high) ;
0622
0623 int ix0 = cegs[0]*high ;
0624 int ix1 = cegs[1]*high ;
0625 int iy0 = cegs[2]*high ;
0626 int iy1 = cegs[3]*high ;
0627 int iz0 = cegs[4]*high ;
0628 int iz1 = cegs[5]*high ;
0629
0630 int photons_per_genstep = cegs[6] ;
0631
0632 int nx = (ix1 - ix0)/2 ;
0633 int ny = (iy1 - iy0)/2 ;
0634 int nz = (iz1 - iz0)/2 ;
0635
0636 int gridaxes = SGenstep::GridAxes(nx, ny, nz);
0637 int num_offset = int(ce_offset.size()) ;
0638
0639 LOG(LEVEL)
0640 << " num_offset " << num_offset
0641 << " ce_scale " << ce_scale
0642 << " nx " << nx
0643 << " ny " << ny
0644 << " nz " << nz
0645 << " GridAxes " << gridaxes
0646 << " GridAxesName " << SGenstep::GridAxesName(gridaxes)
0647 << " high " << high
0648 << " gridscale " << gridscale
0649 << " scale " << scale
0650 << " photons_per_genstep " << photons_per_genstep
0651 ;
0652
0653 double local_scale = ce_scale ? scale*ce.w : scale ;
0654
0655
0656
0657
0658 [[maybe_unused]] unsigned photon_offset = 0 ;
0659 std::vector<quad6> gensteps ;
0660
0661 float rmin = 0. ;
0662 float rmax = std::numeric_limits<float>::max();
0663 if(cegs_radial_range)
0664 {
0665 size_t crr = cegs_radial_range->size();
0666 if(crr > 0) rmin = (*cegs_radial_range)[0] ;
0667 if(crr > 1) rmax = (*cegs_radial_range)[1] ;
0668
0669 LOG(info)
0670 << " [" << CEGS_RADIAL_RANGE << "] "
0671 << " crr " << crr
0672 << " rmin " << std::fixed << std::setw(10) << std::setprecision(2) << rmin
0673 << " rmax " << std::fixed << std::setw(10) << std::setprecision(2) << rmax
0674 << " ce.w " << std::fixed << std::setw(10) << std::setprecision(2) << ce.w
0675 << " scale " << std::fixed << std::setw(10) << std::setprecision(2) << scale
0676 << " local_scale " << std::fixed << std::setw(10) << std::setprecision(2) << local_scale
0677 << "\n"
0678 ;
0679 }
0680
0681
0682
0683 for(int ip=0 ; ip < num_offset ; ip++)
0684 {
0685 const float3& offset = ce_offset[ip] ;
0686
0687 gs.q1.f.x = offset.x ;
0688 gs.q1.f.y = offset.y ;
0689 gs.q1.f.z = offset.z ;
0690 gs.q1.f.w = 1.f ;
0691
0692 for(int ix=ix0 ; ix < ix1+1 ; ix++ )
0693 for(int iy=iy0 ; iy < iy1+1 ; iy++ )
0694 for(int iz=iz0 ; iz < iz1+1 ; iz++ )
0695 {
0696 double tx = double(ix)*local_scale ;
0697 double ty = double(iy)*local_scale ;
0698 double tz = double(iz)*local_scale ;
0699
0700 double radius = std::sqrt(tx*tx + ty*ty + tz*tz);
0701 bool in_radial_range = cegs_radial_range ? ( radius >= rmin && radius <= rmax ) : true ;
0702 if(!in_radial_range) continue ;
0703
0704
0705
0706 const Tran<double>* grid_shift = Tran<double>::make_translate( tx, ty, tz );
0707
0708 bool reverse = false ;
0709 const Tran<double>* transform = Tran<double>::product( geotran, grid_shift, reverse );
0710
0711 qat4* qc = Tran<double>::ConvertFrom( transform->t ) ;
0712
0713 unsigned gsid = SGenstep::GenstepID(ix,iy,iz,ip) ;
0714
0715 SGenstep::ConfigureGenstep(gs, OpticksGenstep_FRAME, gridaxes, gsid, photons_per_genstep );
0716
0717 qc->write(gs);
0718
0719 gensteps.push_back(gs);
0720 photon_offset += std::abs(photons_per_genstep) ;
0721 }
0722 }
0723
0724 LOG(LEVEL)
0725 << " num_offset " << num_offset
0726 << " gensteps.size " << gensteps.size()
0727 ;
0728 return SGenstep::MakeArray(gensteps);
0729 }
0730
0731
0732
0733
0734
0735
0736
0737
0738
0739
0740
0741
0742
0743
0744
0745
0746
0747
0748
0749
0750
0751
0752
0753
0754
0755
0756
0757
0758
0759
0760
0761
0762
0763
0764
0765
0766
0767
0768
0769
0770
0771
0772
0773
0774
0775
0776
0777
0778
0779
0780
0781
0782
0783
0784
0785
0786
0787
0788
0789
0790
0791
0792
0793
0794
0795
0796
0797
0798
0799
0800
0801
0802
0803
0804
0805
0806
0807
0808
0809
0810
0811
0812
0813
0814
0815
0816
0817
0818
0819
0820 void SFrameGenstep::StandardizeCEGS( std::vector<int>& cegs )
0821 {
0822 if( cegs.size() == 0 ) return ;
0823 if( cegs.size() == 8 ) return ;
0824
0825 bool expect = cegs.size() == 4 || cegs.size() == 7 ;
0826 LOG_IF(error, !expect) << " unexpected cegs.size " << cegs.size() ;
0827 assert( expect );
0828
0829 int ix0 = 0 ;
0830 int ix1 = 0 ;
0831 int iy0 = 0 ;
0832 int iy1 = 0 ;
0833 int iz0 = 0 ;
0834 int iz1 = 0 ;
0835 int photons_per_genstep = 0 ;
0836
0837 if( cegs.size() == 4 )
0838 {
0839 ix0 = -cegs[0] ; ix1 = cegs[0] ;
0840 iy0 = -cegs[1] ; iy1 = cegs[1] ;
0841 iz0 = -cegs[2] ; iz1 = cegs[2] ;
0842 photons_per_genstep = cegs[3] ;
0843 }
0844 else if( cegs.size() == 7 )
0845 {
0846 int nx = std::abs(cegs[0]) ;
0847 int ny = std::abs(cegs[1]) ;
0848 int nz = std::abs(cegs[2]) ;
0849 int dx = cegs[3] ;
0850 int dy = cegs[4] ;
0851 int dz = cegs[5] ;
0852 photons_per_genstep = cegs[6] ;
0853
0854 ix0 = -nx + dx ;
0855 iy0 = -ny + dy ;
0856 iz0 = -nz + dz ;
0857 ix1 = nx + dx ;
0858 iy1 = ny + dy ;
0859 iz1 = nz + dz ;
0860 }
0861
0862 cegs.resize(8) ;
0863 cegs[0] = ix0 ;
0864 cegs[1] = ix1 ;
0865 cegs[2] = iy0 ;
0866 cegs[3] = iy1 ;
0867 cegs[4] = iz0 ;
0868 cegs[5] = iz1 ;
0869 cegs[6] = photons_per_genstep ;
0870 cegs[7] = 1 ;
0871
0872
0873
0874
0875 unsigned grid_points = (ix1-ix0+1)*(iy1-iy0+1)*(iz1-iz0+1) ;
0876 unsigned tot_photons = grid_points*std::abs(photons_per_genstep) ;
0877
0878 LOG(LEVEL)
0879 << " CEGS "
0880 << " ix0 ix1 " << ix0 << " " << ix1
0881 << " iy0 iy1 " << iy0 << " " << iy1
0882 << " iz0 iz1 " << iz0 << " " << iz1
0883 << " photons_per_genstep " << photons_per_genstep
0884 << " grid_points (ix1-ix0+1)*(iy1-iy0+1)*(iz1-iz0+1) " << grid_points
0885 << " tot_photons (grid_points*photons_per_genstep) " << tot_photons
0886 ;
0887 }
0888
0889
0890
0891
0892
0893
0894
0895
0896
0897
0898
0899 void SFrameGenstep::GetBoundingBox( float3& mn, float3& mx, const float4& ce, const std::vector<int>& standardized_cegs, float gridscale, const float3& ce_offset )
0900 {
0901 assert( standardized_cegs.size() == 7 ) ;
0902
0903 int ix0 = standardized_cegs[0] ;
0904 int ix1 = standardized_cegs[1] ;
0905 int iy0 = standardized_cegs[2] ;
0906 int iy1 = standardized_cegs[3] ;
0907 int iz0 = standardized_cegs[4] ;
0908 int iz1 = standardized_cegs[5] ;
0909 int photons_per_genstep = standardized_cegs[6] ;
0910
0911
0912 float x0 = float(ix0)*gridscale*ce.w + ce_offset.x ;
0913 float x1 = float(ix1)*gridscale*ce.w + ce_offset.x ;
0914
0915 float y0 = float(iy0)*gridscale*ce.w + ce_offset.y ;
0916 float y1 = float(iy1)*gridscale*ce.w + ce_offset.y ;
0917
0918 float z0 = float(iz0)*gridscale*ce.w + ce_offset.z ;
0919 float z1 = float(iz1)*gridscale*ce.w + ce_offset.z ;
0920
0921 mn.x = x0 ;
0922 mx.x = x1 ;
0923
0924 mn.y = y0 ;
0925 mx.y = y1 ;
0926
0927 mn.z = z0 ;
0928 mx.z = z1 ;
0929
0930 LOG(LEVEL)
0931 << " ce_offset " << ce_offset
0932 << " x0 " << std::setw(10) << std::fixed << std::setprecision(3) << x0
0933 << " x1 " << std::setw(10) << std::fixed << std::setprecision(3) << x1
0934 << " y0 " << std::setw(10) << std::fixed << std::setprecision(3) << y0
0935 << " y1 " << std::setw(10) << std::fixed << std::setprecision(3) << y1
0936 << " z0 " << std::setw(10) << std::fixed << std::setprecision(3) << z0
0937 << " z1 " << std::setw(10) << std::fixed << std::setprecision(3) << z1
0938 << " photons_per_genstep " << photons_per_genstep
0939 << " gridscale " << std::setw(10) << std::fixed << std::setprecision(3) << gridscale
0940 << " ce.w(extent) " << std::setw(10) << std::fixed << std::setprecision(3) << ce.w
0941 ;
0942 }
0943
0944
0945
0946
0947
0948
0949
0950
0951
0952
0953
0954
0955
0956
0957
0958
0959
0960
0961
0962 NP* SFrameGenstep::GenerateCenterExtentGenstepPhotons_( const NP* gsa, float gridscale )
0963 {
0964 std::vector<quad4> pp ;
0965 GenerateCenterExtentGenstepPhotons( pp, gsa, gridscale );
0966
0967 NP* ppa = NP::Make<float>( pp.size(), 4, 4 );
0968 memcpy( ppa->bytes(), (float*)pp.data(), ppa->arr_bytes() );
0969 return ppa ;
0970 }
0971
0972
0973
0974
0975
0976
0977
0978
0979
0980
0981
0982
0983
0984
0985
0986
0987
0988
0989
0990
0991
0992
0993 void SFrameGenstep::GenerateCenterExtentGenstepPhotons( std::vector<quad4>& pp, const NP* gsa, float gridscale )
0994 {
0995 LOG(LEVEL) << " gsa " << gsa->sstr() ;
0996
0997 assert( gsa->shape.size() == 3 && gsa->shape[1] == 6 && gsa->shape[2] == 4 );
0998 assert( gsa->has_shape(-1,6,4) );
0999
1000 std::vector<quad6> gsv(gsa->shape[0]) ;
1001 memcpy( gsv.data(), gsa->bytes(), gsa->arr_bytes() );
1002
1003 quad4 p ;
1004 p.zero();
1005
1006 unsigned seed = 0 ;
1007 SRng<float> rng(seed) ;
1008
1009
1010 float3 paradir ;
1011 qvals(paradir, "PARADIR", "0,0,0" );
1012 bool with_paradir = dot(paradir,paradir) > 0.f ;
1013 if(with_paradir)
1014 {
1015 paradir = normalize(paradir);
1016 LOG(LEVEL) << " PARADIR enabled " << paradir ;
1017 }
1018 else
1019 {
1020 LOG(LEVEL) << " PARADIR NOT-enabled " ;
1021 }
1022
1023
1024 for(unsigned i=0 ; i < gsv.size() ; i++)
1025 {
1026 const quad6& gs = gsv[i];
1027 qat4 qt(gs) ;
1028
1029 C4U gsid ;
1030
1031 int gencode = gs.q0.i.x ;
1032 int gridaxes = gs.q0.i.y ;
1033 gsid.u = gs.q0.u.z ;
1034 int num_photons_ = gs.q0.i.w ;
1035
1036 unsigned num_photons = std::abs(num_photons_);
1037
1038 bool expect = gencode == OpticksGenstep_TORCH || gencode == OpticksGenstep_FRAME ;
1039
1040 LOG_IF(error, !expect)
1041 << "unexpected gencod " << gencode
1042 << " OpticksGenstep_::Name " << OpticksGenstep_::Name(gencode) ;
1043
1044 assert(expect);
1045
1046
1047
1048 double u0, u1 ;
1049 double phi, sinPhi, cosPhi ;
1050 double sinTheta, cosTheta ;
1051
1052
1053
1054 for(unsigned j=0 ; j < num_photons ; j++)
1055 {
1056 u0 = num_photons_ < 0 ? double(j)/double(num_photons-1) : rng() ;
1057
1058 phi = 2.*M_PIf*u0 ;
1059 ssincos(phi,sinPhi,cosPhi);
1060
1061
1062
1063 u1 = rng();
1064 cosTheta = u1 ;
1065 sinTheta = sqrtf(1.0-u1*u1);
1066
1067 if( with_paradir )
1068 {
1069
1070 p.q0.f.x = 0.f ;
1071 p.q0.f.y = u0*float(num_photons-1)*gridscale ;
1072 p.q0.f.z = 0.f ;
1073 p.q0.f.w = 1.f ;
1074
1075 p.q1.f.x = paradir.x ;
1076 p.q1.f.y = paradir.y ;
1077 p.q1.f.z = paradir.z ;
1078 p.q1.f.w = 0.f ;
1079 }
1080 else
1081 {
1082
1083 p.q0.f.x = gs.q1.f.x ;
1084 p.q0.f.y = gs.q1.f.y ;
1085 p.q0.f.z = gs.q1.f.z ;
1086 p.q0.f.w = 1.f ;
1087
1088 SetGridPlaneDirection( p.q1.f, gridaxes, cosPhi, sinPhi, cosTheta, sinTheta );
1089 }
1090
1091
1092 qt.right_multiply_inplace( p.q0.f, 1.f );
1093 qt.right_multiply_inplace( p.q1.f, 0.f );
1094
1095 unsigned char ucj = (j < 255 ? j : 255 ) ;
1096 gsid.c4.w = ucj ;
1097 p.q3.u.w = gsid.u ;
1098
1099 pp.push_back(p) ;
1100 }
1101 }
1102
1103 LOG(LEVEL) << " pp.size " << pp.size() ;
1104 }
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116 void SFrameGenstep::GenerateSimtracePhotons( std::vector<quad4>& simtrace, const std::vector<quad6>& genstep )
1117 {
1118 LOG(LEVEL)
1119 << " genstep.size " << genstep.size()
1120 << " simtrace.size " << simtrace.size()
1121 ;
1122
1123 unsigned seed = 0 ;
1124 SRng<float> rng(seed) ;
1125
1126 bool simtrace_layout = true ;
1127 unsigned count = 0 ;
1128
1129 for(unsigned i=0 ; i < genstep.size() ; i++)
1130 {
1131 const quad6& gs = genstep[i];
1132 C4U gsid ;
1133 int gencode = gs.q0.i.x ;
1134 int gridaxes = gs.q0.i.y ;
1135 gsid.u = gs.q0.u.z ;
1136 int num_photons_ = gs.q0.i.w ;
1137
1138 qat4 qt(gs) ;
1139
1140 unsigned num_photons = std::abs(num_photons_);
1141 bool expect = gencode == OpticksGenstep_TORCH || gencode == OpticksGenstep_FRAME ;
1142
1143 LOG_IF(error, !expect)
1144 << "unexpected gencode " << gencode
1145 << " OpticksGenstep_::Name " << OpticksGenstep_::Name(gencode) ;
1146
1147 assert(expect);
1148
1149
1150 double u0, u1 ;
1151 double phi, sinPhi, cosPhi ;
1152 double sinTheta, cosTheta ;
1153
1154 for(unsigned j=0 ; j < num_photons ; j++)
1155 {
1156 u0 = num_photons_ < 0 ? double(j)/double(num_photons-1) : rng() ;
1157
1158 phi = 2.*M_PIf*u0 ;
1159 ssincos(phi,sinPhi,cosPhi);
1160
1161
1162 u1 = rng();
1163 cosTheta = u1 ;
1164 sinTheta = sqrtf(1.0-u1*u1);
1165
1166 float4 ori ;
1167 ori.x = gs.q1.f.x ;
1168 ori.y = gs.q1.f.y ;
1169 ori.z = gs.q1.f.z ;
1170 ori.w = 1.f ;
1171
1172 float4 dir ;
1173 SetGridPlaneDirection(dir, gridaxes, cosPhi, sinPhi, cosTheta, sinTheta );
1174
1175
1176 qt.right_multiply_inplace( ori, 1.f );
1177 qt.right_multiply_inplace( dir, 0.f );
1178
1179 unsigned char ucj = (j < 255 ? j : 255 ) ;
1180 gsid.c4.w = ucj ;
1181
1182 unsigned identity = gsid.u ;
1183
1184 quad4& p = simtrace[count] ;
1185 p.zero();
1186
1187 if(simtrace_layout)
1188 {
1189 p.q0.f = {0.f,0.f,0.f,0.f} ;
1190 p.q1.f = {0.f,0.f,0.f,0.f} ;
1191 p.q2.f = ori ;
1192 p.q3.f = dir ;
1193 p.q3.u.w = identity ;
1194 }
1195 else
1196 {
1197 p.q0.f = ori ;
1198 p.q1.f = dir ;
1199 p.q2.f = {0.f,0.f,0.f,0.f} ;
1200 p.q3.f = {0.f,0.f,0.f,0.f} ;
1201 p.q3.u.w = identity ;
1202 }
1203 count += 1 ;
1204 }
1205 }
1206 LOG(LEVEL) << " simtrace.size " << simtrace.size() ;
1207 }
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219 void SFrameGenstep::SetGridPlaneDirection( float4& dir, int gridaxes, double cosPhi, double sinPhi, double cosTheta, double sinTheta )
1220 {
1221 if( gridaxes == YZ )
1222 {
1223 dir.x = 0.f ;
1224 dir.y = float(cosPhi) ;
1225 dir.z = float(sinPhi) ;
1226 dir.w = 0.f ;
1227 }
1228 else if( gridaxes == XZ )
1229 {
1230 dir.x = float(cosPhi) ;
1231 dir.y = 0.f ;
1232 dir.z = float(sinPhi) ;
1233 dir.w = 0.f ;
1234 }
1235 else if( gridaxes == XY )
1236 {
1237 dir.x = float(cosPhi) ;
1238 dir.y = float(sinPhi) ;
1239 dir.z = 0.f ;
1240 dir.w = 0.f ;
1241 }
1242 else if( gridaxes == XYZ )
1243 {
1244 dir.x = float(sinTheta * cosPhi) ;
1245 dir.y = float(sinTheta * sinPhi) ;
1246 dir.z = float(cosTheta) ;
1247 dir.w = 0.f ;
1248 }
1249 else
1250 {
1251 LOG(fatal) << " invalid gridaxes value " << gridaxes ;
1252 assert(0);
1253 }
1254 }
1255