eic-opticks/CSG/CSGScan.cc

0001 /**
0002 CSGScan.cc
0003 ==========
0004
0005 Restructured this in preparaction for the CUDA/MOCK_CUDA impl
0006 by preparing all the rays ahead of time and
0007 then doing all the intersects together.
0008
0009 DONE : split off the CUDA/MOCK_CUDA part that can be parallelized
0010 DONE : add launch, pullback
0011
0012 TODO : try using CUDA __constant__ and cudaMemcpyToSymbol for the CSGParams
0013
0014 * HMM: it makes more sense to do that for the geometry pointers, not so much
0015   for the less constant rays and intersects
0016
0017
0018 **/
0019
0020 #include "NPX.h"
0021 #include "sstr.h"
0022
0023 #include "CSGFoundry.h"
0024 #include "CSGSolid.h"
0025 #include "CSGScan.h"
0026
0027 #include "CSGParams.h"
0028
0029 #include "CU.h"
0030
0031 /**
0032 CSGScan::CSGScan
0033 -----------------
0034
0035 h
0036    host pointer to CSGParams populated on host
0037 d
0038    host pointer to CSGParams populated on device
0039 d_d
0040    device pointer of the copied to device d instance
0041 c
0042    host pointer to the copied back CSGParams instance
0043
0044 **/
0045
0046
0047 CSGScan::CSGScan( const CSGFoundry* fd_, const CSGSolid* so_, const char* opts_  )
0048     :
0049     fd(fd_),
0050     prim0(fd->getPrim(0)),
0051     node0(fd->getNode(0)),
0052     so(so_),
0053     primIdx0(so->primOffset),
0054     primIdx1(so->primOffset+so->numPrim),
0055     primIdx(primIdx0),   //
0056     prim(prim0 + primIdx),
0057     nodeOffset(prim->nodeOffset()),
0058     node(node0 + nodeOffset),
0059     h(new CSGParams {}),
0060     d(new CSGParams {}),
0061     d_d(nullptr),
0062     c(new CSGParams {})
0063 {
0064     initGeom_h();
0065     initRays_h(opts_);
0066
0067     initGeom_d();
0068     initRays_d();
0069     initParams_d();
0070 }
0071
0072 void CSGScan::initGeom_h()
0073 {
0074     h->devp = false ;
0075     h->node = node ;
0076     h->plan = fd->getPlan(0) ;
0077     h->itra = fd->getItra(0) ;
0078 }
0079
0080
0081 void CSGScan::initRays_h(const char* opts_)
0082 {
0083     std::vector<std::string> opts ;
0084     sstr::Split(opts_, ',', opts );
0085
0086     std::vector<quad4> qq ;
0087     for(unsigned i=0 ; i < opts.size() ; i++)
0088     {
0089         const char* opt = opts[i].c_str();
0090         add_scan(qq, opt);
0091     }
0092
0093     h->num = qq.size() ;
0094     h->qq = new quad4[h->num];
0095     memcpy( (void*)h->qq, qq.data(), sizeof(quad4)*h->num ) ;
0096
0097     h->tt = new quad4[h->num];
0098 }
0099
0100
0101 void CSGScan::initGeom_d()
0102 {
0103     assert( fd->isUploaded() );
0104
0105     d->devp = true ;
0106     d->node = fd->d_node ;
0107     d->plan = fd->d_plan ;
0108     d->itra = fd->d_itra ;
0109 }
0110 void CSGScan::initRays_d()
0111 {
0112     d->num = h->num ;
0113     d->qq = CU::UploadArray<quad4>( h->qq, h->num ) ;
0114     d->tt = CU::AllocArray<quad4>( d->num ) ;
0115 }
0116 void CSGScan::initParams_d()
0117 {
0118     d_d = CU::UploadArray<CSGParams>( d, 1 ) ;
0119 }
0120
0121
0122
0123 void CSGScan::add_scan(std::vector<quad4>& qq, const char* opt)
0124 {
0125     if(strcmp(opt,"axis")==0)   add_axis_scan(qq);
0126     if(strcmp(opt,"circle")==0) add_circle_scan(qq);
0127     if(strcmp(opt,"rectangle")==0) add_rectangle_scan(qq);
0128 }
0129
0130 /**
0131 CSGScan::add_axis_scan
0132 ----------------------
0133
0134 Six rays from so->center_extent center along X,Y,Z,-X,-Y,-Z axes
0135
0136 **/
0137
0138 void CSGScan::add_axis_scan(std::vector<quad4>& qq)
0139 {
0140     float t_min = 0.f ;
0141     float4 ce = so->center_extent ;
0142     float3 origin = make_float3(ce);
0143
0144     std::vector<float3> dirs ;
0145     dirs.push_back( make_float3( 1.f, 0.f, 0.f));
0146     dirs.push_back( make_float3( 0.f, 1.f, 0.f));
0147     dirs.push_back( make_float3( 0.f, 0.f, 1.f));
0148
0149     dirs.push_back( make_float3(-1.f, 0.f, 0.f));
0150     dirs.push_back( make_float3( 0.f,-1.f, 0.f));
0151     dirs.push_back( make_float3( 0.f, 0.f,-1.f));
0152
0153     add_q(qq, t_min, origin, dirs );
0154 }
0155
0156 /**
0157 CSGScan::add_circle_scan
0158 -------------------------
0159
0160
0161
0162 **/
0163
0164 void CSGScan::add_circle_scan(std::vector<quad4>& qq)
0165 {
0166     float t_min = 0.f ;
0167     float4 ce = so->center_extent ;
0168     float3 center = make_float3( ce );
0169     float extent = ce.w ;
0170     float radius = 2.0f*extent ;
0171
0172     if(0) std::cout
0173         << "CSGScan::add_circle_scan"
0174         << " extent " << extent
0175         << " radius " << radius
0176         << std::endl
0177         ;
0178
0179     // M_PIf from sutil_vec_math.h
0180     for(float phi=0. ; phi <= M_PIf*2.0 ; phi+=M_PIf*2.0/1000.0 )
0181     {
0182         float3 origin = center + make_float3( radius*sin(phi), 0.f, radius*cos(phi) );
0183         float3 direction = make_float3( -sin(phi),  0.f, -cos(phi) );
0184         add_q(qq, t_min, origin, direction);
0185     }
0186 }
0187
0188 void CSGScan::add_rectangle_scan(std::vector<quad4>& qq)
0189 {
0190     float4 ce = so->center_extent ;
0191     float extent = ce.w ;
0192     float halfside = 2.0f*extent ;
0193     unsigned nxz = 100 ;
0194     unsigned ny = 10 ;
0195     float t_min = 0.f ;
0196
0197     if(0) std::cout
0198         << "CSGScan::add_rectangle_scan"
0199         << " extent " << extent
0200         << " halfside " << halfside
0201         << std::endl
0202         ;
0203
0204
0205     for(float y=-halfside ; y <= halfside ; y += halfside/float(ny) )
0206     {
0207         _add_rectangle_scan(qq, t_min, nxz, halfside,   y );
0208     }
0209 }
0210
0211 void CSGScan::_add_rectangle_scan(std::vector<quad4>& qq, float t_min, unsigned n, float halfside, float y )
0212 {
0213     // shooting up/down
0214
0215     float3 z_up   = make_float3( 0.f, 0.f,  1.f);
0216     float3 z_down = make_float3( 0.f, 0.f, -1.f);
0217
0218     float3 z_top = make_float3( 0.f, y,  halfside );
0219     float3 z_bot = make_float3( 0.f, y, -halfside );
0220
0221     // shooting left/right
0222
0223     float3 x_right = make_float3(  1.f, 0.f,  0.f);
0224     float3 x_left  = make_float3( -1.f, 0.f,  0.f);
0225
0226     float3 x_lhs = make_float3( -halfside, y,  0.f );
0227     float3 x_rhs = make_float3(  halfside, y,  0.f );
0228
0229     for(float v=-halfside ; v <= halfside ; v+= halfside/float(n) )
0230     {
0231         z_top.x = v ;
0232         z_bot.x = v ;
0233
0234         add_q(qq, t_min, z_top, z_down );
0235         add_q(qq, t_min, z_bot, z_up   );
0236
0237         x_lhs.z = v ;
0238         x_rhs.z = v ;
0239         add_q(qq, t_min, x_lhs, x_right );
0240         add_q(qq, t_min, x_rhs, x_left  );
0241     }
0242 }
0243
0244 void CSGScan::add_q(std::vector<quad4>& qq, const float t_min, const float3& ray_origin, const std::vector<float3>& dirs )
0245 {
0246     for(unsigned i=0 ; i < dirs.size() ; i++)
0247     {
0248         const float3& ray_direction = dirs[i] ;
0249         add_q(qq, t_min, ray_origin, ray_direction );
0250     }
0251 }
0252
0253 void CSGScan::add_q(std::vector<quad4>& qq, float t_min, const float3& ray_origin, const float3& ray_direction )
0254 {
0255     quad4 q = {} ;
0256
0257     q.q0.f = make_float4(ray_origin);
0258     q.q1.f = make_float4(ray_direction);
0259     q.q1.f.w = t_min ;
0260
0261     qq.push_back(q);
0262 }
0263
0264 void CSGScan::intersect_h()
0265 {
0266     for(int i=0 ; i < h->num ; i++)
0267     {
0268         h->intersect(i);
0269     }
0270 }
0271
0272
0273
0274 extern void CSGScan_intersect( dim3 numBlocks, dim3 threadsPerBlock, CSGParams* d );
0275
0276 void CSGScan::intersect_d()
0277 {
0278     dim3 numBlocks ;
0279     dim3 threadsPerBlock ;
0280     CU::ConfigureLaunch1D( numBlocks, threadsPerBlock, d->num, 512u );
0281
0282     CSGScan_intersect( numBlocks, threadsPerBlock, d_d ) ;
0283
0284     download();
0285 }
0286
0287
0288 void CSGScan::download()
0289 {
0290     c->num = d->num ;
0291     c->qq = CU::DownloadArray<quad4>( d->qq, d->num ) ;
0292     c->tt = CU::DownloadArray<quad4>( d->tt, d->num ) ;
0293     assert( d->devp == true ) ;
0294     c->devp = false ;
0295 }
0296
0297
0298
0299
0300 void CSGScan::dump( const quad4& t )  // stat
0301 {
0302     bool valid_isect = t.q0.i.w == 1 ;
0303
0304     const float4& isect = t.q3.f ;
0305     const float4& ray_origin  = t.q0.f ;
0306     const float4& ray_direction = t.q1.f ;
0307
0308     std::cout
0309         << std::setw(30) << so->label
0310         << " valid_isect " << valid_isect
0311         << " isect ( "
0312         << std::setw(10) << std::fixed << std::setprecision(3) << isect.x
0313         << std::setw(10) << std::fixed << std::setprecision(3) << isect.y
0314         << std::setw(10) << std::fixed << std::setprecision(3) << isect.z
0315         << std::setw(10) << std::fixed << std::setprecision(3) << isect.w
0316         << " ) "
0317         << " dir ( "
0318         << std::setw(10) << std::fixed << std::setprecision(3) << ray_direction.x
0319         << std::setw(10) << std::fixed << std::setprecision(3) << ray_direction.y
0320         << std::setw(10) << std::fixed << std::setprecision(3) << ray_direction.z
0321         << " ) "
0322         << " ori ( "
0323         << std::setw(10) << std::fixed << std::setprecision(3) << ray_origin.x
0324         << std::setw(10) << std::fixed << std::setprecision(3) << ray_origin.y
0325         << std::setw(10) << std::fixed << std::setprecision(3) << ray_origin.z
0326         << " ) "
0327         << std::endl
0328         ;
0329 }
0330
0331 std::string CSGScan::brief() const
0332 {
0333     std::stringstream ss ;
0334     ss << " h " << brief(h) << "\n" ;
0335     ss << " d " << brief(c) << "\n" ; // actually c copied to host from d
0336     std::string str = ss.str() ;
0337     return str ;
0338 }
0339
0340 std::string CSGScan::brief(CSGParams* s) const
0341 {
0342     int n_hit = s->num_valid_isect();
0343     std::stringstream ss ;
0344     ss
0345         << " num " << s->num
0346         << " n_hit " << n_hit
0347         << " (num-n_hit) " << (s->num-n_hit)
0348         ;
0349     std::string str = ss.str() ;
0350     return str ;
0351 }
0352
0353 NPFold* CSGScan::serialize_(CSGParams* s) const
0354 {
0355     NPFold* fold = new NPFold ;
0356     NP* _qq = NPX::ArrayFromData<float>( (float*)s->qq, s->num, 4, 4 ) ;
0357     NP* _tt = NPX::ArrayFromData<float>( (float*)s->tt, s->num, 4, 4 ) ;
0358     fold->add("qq", _qq );
0359     fold->add("tt", _tt );
0360     return fold ;
0361 }
0362
0363 NPFold* CSGScan::serialize() const
0364 {
0365     NPFold* fold = new NPFold ;
0366     fold->add_subfold("h", serialize_(h));
0367     fold->add_subfold("d", serialize_(c));   // c is d copied back to host
0368     return fold ;
0369 }
0370
0371 void CSGScan::save(const char* base) const
0372 {
0373    NPFold* fold = serialize();
0374    fold->save(base) ;
0375 }