eic-opticks/sysrap/sphotonlite.h

0001 #pragma once
0002
0003 /**
0004 sphotonlite.hh
0005 ===============
0006
0007 +----+-----------------------+----------------+-----------------------+----------------+------------------------------+
0008 | q  |      x                |      y         |     z                 |      w         |  notes                       |
0009 +====+=======================+================+=======================+================+==============================+
0010 |    | u:hitcount_identity   |  f:time        | u:lposcost_lposfphi   | u:flagmask     |                              |
0011 | q0 |                       |                |                       |                |                              |
0012 |    | off:0, 2              | off:4          | off:8,10              | off:12         |  off:byte offsets            |
0013 +----+-----------------------+----------------+-----------------------+----------------+------------------------------+
0014
0015
0016 **/
0017
0018
0019 #if defined(__CUDACC__) || defined(__CUDABE__)
0020 #    define SPHOTONLITE_METHOD __host__ __device__ __forceinline__
0021 #else
0022 #    define SPHOTONLITE_METHOD inline
0023 #endif
0024
0025
0026 #include <cstdint>
0027 #if defined(__CUDACC__) || defined(__CUDABE__)
0028 #else
0029    #include <iostream>
0030    #include <iomanip>
0031    #include <sstream>
0032    #include <cstring>
0033    #include <cassert>
0034    #include "OpticksPhoton.h"
0035    #include "OpticksPhoton.hh"
0036    #include "NP.hh"
0037    struct sphotonlite_selector ;
0038 #endif
0039
0040
0041 struct sphotonlite
0042 {
0043 #if defined(__CUDACC__) || defined(__CUDABE__)
0044 #else
0045     static constexpr const char* NAME = "sphotonlite" ;
0046 #endif
0047
0048     uint32_t  hitcount_identity ; // hi16:hitcount, lo16:identity
0049     float     time ;
0050     uint32_t  lposcost_lposfphi ; // hi16:lposcost, lo16:lposfphi
0051     uint32_t  flagmask ;
0052
0053     SPHOTONLITE_METHOD void init(unsigned _identity, float _time, unsigned _flagmask);
0054     SPHOTONLITE_METHOD void set_lpos(float lposcost, float lposfphi);
0055
0056     SPHOTONLITE_METHOD void set_hitcount_identity( unsigned hitcount, unsigned identity );
0057     SPHOTONLITE_METHOD void set_hitcount( unsigned hitcount );
0058     SPHOTONLITE_METHOD void set_identity( unsigned identity );
0059
0060     SPHOTONLITE_METHOD unsigned hitcount() const { return hitcount_identity >> 16 ; }
0061     SPHOTONLITE_METHOD unsigned identity() const { return hitcount_identity & 0xffffu ; }
0062     SPHOTONLITE_METHOD unsigned pmtid() const
0063     {
0064         unsigned id = identity() ;
0065         return id == 0 ? 0xffffffffu : id - 1 ;    // id:0 means not-a-sensor
0066     }
0067
0068     struct key_functor
0069     {
0070         float    tw;
0071         SPHOTONLITE_METHOD uint64_t operator()(const sphotonlite& p) const
0072         {
0073             unsigned id = p.identity() ;
0074             unsigned bucket = static_cast<unsigned>(p.time / tw);
0075             return (uint64_t(id) << 48) | uint64_t(bucket);
0076         }
0077     };
0078
0079     struct reduce_op
0080     {
0081         SPHOTONLITE_METHOD sphotonlite operator()(const sphotonlite& a, const sphotonlite& b) const
0082         {
0083             sphotonlite r = a;
0084             r.time = fminf(a.time, b.time);
0085             r.flagmask |= b.flagmask;
0086             unsigned hc = a.hitcount() + b.hitcount();
0087             r.set_hitcount_identity(hc, a.identity());
0088             return r;
0089         }
0090     };
0091
0092     // this requires any bit match, not all bit like the below sphotonlite_selector
0093     struct select_pred
0094     {
0095         unsigned mask;
0096         SPHOTONLITE_METHOD bool operator()(const sphotonlite& p) const { return (p.flagmask & mask) != 0; }
0097     };
0098
0099
0100
0101
0102
0103 #if defined(__CUDACC__) || defined(__CUDABE__)
0104 #else
0105
0106     bool operator==(const sphotonlite& other) const noexcept
0107     {
0108         return hitcount_identity  == other.hitcount_identity &&
0109                time               == other.time              &&
0110                lposcost_lposfphi  == other.lposcost_lposfphi &&
0111                flagmask           == other.flagmask ;
0112     }
0113     bool operator!=(const sphotonlite& other) const noexcept
0114     {
0115         return !(*this == other);
0116     }
0117
0118
0119     SPHOTONLITE_METHOD void get_lpos(float& lposcost, float& lposfphi) const ;
0120     SPHOTONLITE_METHOD void get_lpos_theta_phi(float& lpos_theta, float& lpos_phi) const ;
0121
0122     SPHOTONLITE_METHOD std::string lpos_() const ;
0123
0124
0125     SPHOTONLITE_METHOD std::string flagmask_() const ;
0126     SPHOTONLITE_METHOD std::string desc() const ;
0127
0128     SPHOTONLITE_METHOD static NP* zeros(    size_t num_photon);
0129     SPHOTONLITE_METHOD static bool expected( const NP* l );
0130
0131     SPHOTONLITE_METHOD static NP* demoarray(size_t num_photon);
0132     SPHOTONLITE_METHOD static NP* select( const NP* photonlite, const sphotonlite_selector* photonlite_selector );
0133     SPHOTONLITE_METHOD static void loadbin( std::vector<sphotonlite>& photonlite, const char* path );
0134
0135     SPHOTONLITE_METHOD static std::string Desc(const NP* a, size_t edge=20);
0136     SPHOTONLITE_METHOD static NP*         MockupForMergeTest(size_t ni);
0137
0138
0139     SPHOTONLITE_METHOD static std::string desc_diff( const sphotonlite* a, const sphotonlite* b, size_t ni );
0140
0141 #endif
0142
0143 };
0144
0145 /**
0146 sphotonlite_selector
0147 --------------------
0148
0149 All bits that are set within the *hitmask* are required to be set within the *flagmask* for a match
0150 This functor binds to the below predicate signature, where S=sphotonlite::
0151
0152      std::function<bool(const S*)>
0153
0154 Instances of this functor are used on CPU via the below sphotonlite::select
0155 and on GPU from QEvt::gatherHitLite_
0156
0157 **/
0158
0159 struct sphotonlite_selector
0160 {
0161     uint32_t hitmask ;
0162     sphotonlite_selector(uint32_t hitmask_) : hitmask(hitmask_) {};
0163     SPHOTONLITE_METHOD bool operator() (const sphotonlite& p) const { return ( p.flagmask  & hitmask ) == hitmask  ; }
0164     SPHOTONLITE_METHOD bool operator() (const sphotonlite* p) const { return ( p->flagmask & hitmask ) == hitmask  ; }
0165 };
0166
0167
0168
0169
0170
0171
0172 /**
0173 sphotonlite::init
0174 -------------------
0175
0176 Example, CSGOptiX/CSGOptiX7.cu::
0177
0178
0179     457     if( evt->photonlite )
0180     458     {
0181     459         sphotonlite l ;
0182     460         l.init( ctx.p.identity, ctx.p.time, ctx.p.flagmask );
0183     461         l.set_lpos(prd->lposcost(), prd->lposfphi() );
0184     462         evt->photonlite[idx] = l ;  // *idx* (not *photon_idx*) as needs to go from zero for photons from a slice of genstep array
0185     463     }
0186
0187
0188
0189 **/
0190
0191
0192 inline void sphotonlite::init(unsigned _identity, float _time, unsigned _flagmask)
0193 {
0194     hitcount_identity = ( 0x1u << 16 ) | ( _identity & 0xFFFFu ) ; // hitcount starts as one
0195     time = _time ;
0196     lposcost_lposfphi = 0u ;
0197     flagmask = _flagmask ;
0198 }
0199
0200
0201 inline void sphotonlite::set_hitcount_identity( unsigned hitcount, unsigned identity )
0202 {
0203     hitcount_identity = (( hitcount & 0xFFFFu ) << 16 ) | ( identity & 0xFFFFu ) ;
0204 }
0205 inline void sphotonlite::set_hitcount( unsigned hc )
0206 {
0207     hitcount_identity = ( hitcount_identity & 0x0000ffffu ) | (( 0x0000ffffu & hc ) << 16);
0208 }
0209 inline void sphotonlite::set_identity( unsigned id )
0210 {
0211     hitcount_identity = ( hitcount_identity & 0xffff0000u ) | (( 0x0000ffffu & id ) <<  0);
0212 }
0213
0214
0215
0216
0217
0218
0219 /**
0220 sphotonlite::set_lpos
0221 ----------------------
0222
0223 lposcost
0224      cosine of the intersect local position polar angle,
0225      within range 0.->1. [assumes only front hemi-sphere is sensitive]
0226
0227 lposfphi
0228      fraction of azimuthal angle in radians of local intersect position,
0229      scaled by 2pi to be within range 0.->1.
0230
0231 **/
0232
0233
0234 inline void sphotonlite::set_lpos(float lposcost, float lposfphi)
0235 {
0236 #if defined(__CUDACC__) || defined(__CUDABE__)
0237 #else
0238     assert( lposcost >= 0.f && lposcost <= 1.f );
0239     assert( lposfphi >= 0.f && lposfphi <= 1.f );
0240 #endif
0241     lposcost_lposfphi =
0242                         ((uint16_t)(lposcost * 0xffffu + 0.5f) << 16) |
0243                         ((uint16_t)(lposfphi * 0xffffu + 0.5f) <<  0)
0244                       ;
0245 }
0246
0247
0248 #if defined(__CUDACC__) || defined(__CUDABE__)
0249 #else
0250 inline void sphotonlite::get_lpos(float& lposcost, float& lposfphi) const
0251 {
0252     lposcost = float(lposcost_lposfphi >> 16    ) / 0xffffu ;
0253     lposfphi = float(lposcost_lposfphi & 0xffffu) / 0xffffu ;
0254
0255     assert( lposcost >= 0.f && lposcost <= 1.f );
0256     assert( lposfphi >= 0.f && lposfphi <= 1.f );
0257 }
0258
0259
0260 inline void sphotonlite::get_lpos_theta_phi(float& lpos_theta, float& lpos_phi) const
0261 {
0262     float lposcost = float(lposcost_lposfphi >> 16    ) / 0xffffu ;
0263     float lposfphi = float(lposcost_lposfphi & 0xffffu) / 0xffffu ;
0264
0265     lpos_theta = std::acos( lposcost );  // HMM: could store fthe instead of cost ?
0266     lpos_phi = phi_from_fphi( lposfphi );
0267
0268 }
0269
0270 inline std::string sphotonlite::lpos_() const
0271 {
0272     float lposcost, lposfphi ;
0273     get_lpos(lposcost, lposfphi);
0274
0275     float lpos_theta, lpos_phi ;
0276     get_lpos_theta_phi(lpos_theta, lpos_phi);
0277
0278
0279     std::stringstream ss ;
0280     ss
0281         << " cost " << std::setw(7) << std::fixed << std::setprecision(5) << lposcost
0282         << " theta " << std::setw(7) << std::fixed << std::setprecision(5) << lpos_theta
0283         << " fphi " << std::setw(7) << std::fixed << std::setprecision(5) << lposfphi
0284         << " phi " << std::setw(7) << std::fixed << std::setprecision(5) << lpos_phi
0285
0286
0287         ;
0288     std::string str = ss.str() ;
0289     return str ;
0290 }
0291
0292
0293
0294
0295
0296
0297 inline std::string sphotonlite::desc() const
0298 {
0299     std::stringstream ss ;
0300     ss << "[sphotonlite:"
0301        << " t " <<  std::setw(8) << time
0302        << " id " <<  std::setw(8) << identity()
0303        << " hc " <<  std::setw(8) << hitcount()
0304        << " lp " <<  lpos_()
0305        << " fm " <<  flagmask_()
0306        << "]"
0307        ;
0308
0309     std::string str = ss.str() ;
0310     return str ;
0311 }
0312
0313
0314
0315 /**
0316
0317
0318
0319     In [2]: t.demoarray.view(np.uint32)[:,0] >> 16   ## hitcount
0320     Out[2]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype=uint32)
0321
0322     In [3]: t.demoarray.view(np.uint32)[:,0] & 0xffff   ## identity
0323     Out[3]: array([   0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000], dtype=uint32)
0324
0325
0326     In [7]: (t.demoarray.view(np.uint32)[:,2] >> 16)/0xffff
0327     Out[7]: array([0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5])
0328
0329     n [6]: (t.demoarray.view(np.uint32)[:,2] & 0xffff)/0xffff
0330     Out[6]: array([0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6])
0331
0332
0333 **/
0334
0335
0336 inline NP* sphotonlite::zeros(size_t num_photon) // static
0337 {
0338     NP* l = NP::Make<uint32_t>(num_photon, 4);
0339     return l ;
0340 }
0341
0342 inline bool sphotonlite::expected( const NP* l ) // static
0343 {
0344     bool expected_type = l && l->uifc == 'u' && l->ebyte == 4 ;
0345     bool expected_shape =  l && l->has_shape(-1, 4) ;
0346     bool expected_arr = expected_type && expected_shape ;
0347
0348     if(!expected_arr) std::cerr
0349         << "sphotonlite::expected"
0350         << " expected_arr " << ( expected_arr ? "YES" : "NO " )
0351         << " expected_type " << ( expected_type ? "YES" : "NO " )
0352         << " expected_shape " << ( expected_shape ? "YES" : "NO " )
0353         << " l.sstr " << ( l ? l->sstr() : "-" )
0354         ;
0355     return expected_arr ;
0356 }
0357
0358
0359 inline NP* sphotonlite::demoarray(size_t num_photon) // static
0360 {
0361     NP* l = zeros(num_photon);
0362     sphotonlite* ll  = (sphotonlite*)l->bytes();
0363     for(size_t i=0 ; i < num_photon ; i++)
0364     {
0365         ll[i].set_hitcount_identity( i, i*1000 );
0366         ll[i].flagmask = EFFICIENCY_COLLECT | SURFACE_DETECT ;
0367         ll[i].time = i*0.1f ;
0368         ll[i].set_lpos( 0.5f, 0.6f );
0369     }
0370     return l ;
0371 }
0372
0373 inline std::string sphotonlite::flagmask_() const
0374 {
0375     return OpticksPhoton::FlagMaskLabel(flagmask) ;
0376 }
0377
0378 inline NP* sphotonlite::select( const NP* photonlite, const sphotonlite_selector* photonlite_selector )
0379 {
0380     NP* hitlite = photonlite ? photonlite->copy_if<uint32_t, sphotonlite>(*photonlite_selector) : nullptr ;
0381     return hitlite ;
0382 }
0383
0384
0385 inline void sphotonlite::loadbin( std::vector<sphotonlite>& photonlite, const char* path )
0386 {
0387     // 1. open path with pointer at-end "ate", use tellg to determine number of hits
0388
0389     std::ifstream f(path, std::ios::binary | std::ios::ate);
0390     assert(f && "sphotonlite::loadbin failed to open path");
0391     size_t bytes = f.tellg();
0392     size_t n = bytes / sizeof(sphotonlite);
0393     f.seekg(0);
0394
0395     // 2. read into photonlite vector
0396     photonlite.resize(n);
0397     f.read((char*)photonlite.data(), bytes);
0398 }
0399
0400
0401 inline std::string sphotonlite::Desc(const NP* a, size_t edge) // static
0402 {
0403     std::stringstream ss ;
0404     ss << "[sphotonlite::Desc a.sstr " << ( a ? a->sstr() : "-" ) << "\n" ;
0405     size_t ni = a->num_items();
0406     sphotonlite* ll = (sphotonlite*)a->bytes();
0407     for(size_t i=0 ; i < ni ; i++)
0408     {
0409         if(i < edge || i > (ni - edge)) ss << std::setw(6) << i << " : " << ll[i].desc() << "\n" ;
0410         else if( i == edge )  ss << std::setw(6) << "" << " : " << "..." << "\n" ;
0411     }
0412
0413     std::string str = ss.str() ;
0414     return str ;
0415 }
0416
0417 SPHOTONLITE_METHOD NP* sphotonlite::MockupForMergeTest(size_t ni) // static
0418 {
0419     NP* a = NP::Make<unsigned>(ni, 4);
0420     sphotonlite* aa = (sphotonlite*)a->bytes();
0421     for(size_t i=0 ; i < ni ; i++)
0422     {
0423         sphotonlite& p = aa[i];
0424         p.time = float( i % 100 )*0.1f ;
0425         p.set_hitcount_identity( 1, i % 100 );
0426         p.flagmask = i % 5 == 0 ? EFFICIENCY_COLLECT : EFFICIENCY_CULL ;
0427     }
0428     return a ;
0429 }
0430
0431
0432 inline std::string sphotonlite::desc_diff( const sphotonlite* a, const sphotonlite* b, size_t ni )
0433 {
0434     size_t tot_diff = 0 ;
0435     size_t tot_same = 0 ;
0436
0437     std::stringstream ss ;
0438     ss << "[sphotonlite::desc_diff\n"
0439        << " ni " << ni
0440        << " a " << ( a ? "YES" : "NO " )
0441        << " b " << ( b ? "YES" : "NO " )
0442        << "\n"
0443        ;
0444
0445     if( a != nullptr && b != nullptr )
0446     {
0447         for(size_t i = 0; i < ni; ++i)
0448         {
0449             if (a[i] != b[i])
0450             {
0451                 tot_diff +=1 ;
0452                 ss << "Index " << i << " differs:\n"
0453                           << "  a: " << a[i].desc() << "\n"
0454                           << "  b: " << b[i].desc() << "\n";
0455             }
0456             else
0457             {
0458                 tot_same +=1 ;
0459             }
0460         }
0461     }
0462
0463     ss << "]sphotonlite::desc_diff\n"
0464        << " tot_diff " << tot_diff
0465        << " tot_same " << tot_same
0466        << "\n"
0467        ;
0468
0469     std::string str = ss.str() ;
0470     return str ;
0471 }
0472
0473
0474
0475
0476
0477
0478
0479
0480 #endif
0481
0482
0483