Analytic Geometry =================== Principal ---------- Triangle meshes are a convenient initial step to the GPU as all geometry can be treated with the same code. Special treatment of important geometry (PMTs) however is expected to have large performance gains. Ray intersection with CSG solids boils down to analytic solving quadratic/cubic polynomials. There is a technique to handle union intersections by applying boolean operations to intersection segments of the sub volumes. Aligning three geometries ----------------------------- Actually not 2 but 3 geometries: OpenGL GMergedMesh triangulated geometry, to first order only impacts visualization not propagation/tracing, but the solid identity infomation including boundary indices comes from GMergedMesh. Geometry is modified by *--box* option according to *--boxconfig*. The orignal instance-1 GMergedMesh is combined with the GSolid fabricated by GTestBox OptiX Triangulated translation of GMergedMesh triangulated into OptiX OptiX Analytic somewhat manual addon entering via GPmt and created in python by *pmt-parts* Containment Box Handling -------------------------- Actually the parameters of the containment bbox should be obtained from GMergedMesh.1, it should be enlarged according to input parameter and the resulting gbbox should then be fed down via GGeo::modifyGeometry to: GTestBox to fabricate the GSolid GPmt to create a *Part* struct to tack on to the parts buffer * Avoids mal-duplication. * allows pslice to apply to content and not the containment box (but pslicing is a low level Analytic Part debug thing, as it cannot be applied at triangulated level) GPmt Operation --------------- * reads *pmt-parts* written float parts buffer, deriving a solid buffer from it using solid or node indices which are embedded into the parts buffer * these five indices 0:4 hail from the detdesc tree parse of *hemi-pmt.xml* and these must match the indices from the GMergedMesh identity buffer * for debug want to skip some solids, eg to see behavior with a solid lump of Pyrex/Vacuum/Bialkali but at the moment the identity would then be misaligned ? Photocathode ------------- Geant4/detdesc model is with seperate very thin spherical parts, including a shared boundary with the inside of the pyrex. For surface based geometry coincident boundaries are unhealthy, so instead model it similar to how Pyrex/Vacuum modelling of the pyrex envelope is done. Again this might entail adding solids, which will mess up identity. Seems Test Box Debugging -------------------- :: ggv-box(){ ggv --box \ --animtimemax 7 \ --boxconfig "size=2;boundary=MineralOil/Rock//;" \ --torchconfig "source=0,0,400;target=0,0,0;radius=150;zenith_azimuth=1,0,1,0" \ $* } * need to sort out analytic boundary identity labelling, the missers think they are in Pyrex, when should be MO [FIXED] * GGeo/GMergedMesh should be orchestrating the analytic PMT for commality, currently OGeo/GPmt just grabbing from GCache [FIXED] * where do i set the boundary of the analytic test box ? again need to bring control of triangulated and analytic together to avoid confusion [FIXED] From OptiX point of view there 5 (or 6 with the container) primitives. These need to line up with the triangulated solids for identity to work. Each primitive has a small numbers of parts (up to 4). Total of 16 parts. When put the container at the end in pmt-parts the material mapping works better as that aligns with GMergedMesh::combine order. This is brittle, will fail when slicing. :: [2015-Nov-05 19:50:49.364081]:info: OGeo::makeAnalyticGeometry identity buffer BufferId -1 BufferTarget 0 NumBytes 96 ItemSize 16 NumElements 4 NumItems 6 NumElementsTotal 24 ( 0) 3199 47 27 0 ( 1) 3200 46 28 0 ( 2) 3201 43 29 3 ( 3) 3202 44 30 0 ( 4) 3203 45 30 0 ( 5) 5 1000 123 0 :: [2015-11-05 19:53:00.317306] [0x000007fff7448031] [info] GBndLib::dump ( 0) im: Vacuum om: Vacuum is: os: ... ( 27) im: Pyrex om: MineralOil is: os: ( 28) im: Vacuum om: Pyrex is: os: ( 29) im: Bialkali om: Vacuum is: os:lvPmtHemiCathodeSensorSurface ( 30) im: OpaqueVacuum om: Vacuum is: os: ... (122) im: RadRock om: Rock is: os: Implementing container making C++ side ? :: simon:pmt blyth$ ggv --pmt 0:15 [2015-Nov-05 20:44:09.782958]:info: 0:/usr/local/env/optix/ggeo/bin/GPmtTest [2015-Nov-05 20:44:09.783831]:info: 1:0:15 [2015-Nov-05 20:44:09.784031]:info: NPY::make_slice from 16 -> 15 slice NSlice 0 : 15 : 1 [2015-Nov-05 20:44:09.784205]:info: GPmt::loadFromCache slicing partBuf origBuf 16,4,4 partBuf 15,4,4 GPmt::make_container pbb min -101.168 -101.168 -23.838 max 101.168 101.168 56.000 ... GPmt::make_container pbb min -27.500 -27.500 -164.500 max 27.500 27.500 1.500 GPmt::make_container bb min -101.168 -101.168 -169.000 max 101.168 101.168 131.000 GPmt::make_container bb factor 3.0 min -551.168 -551.168 -619.000 max 551.168 551.168 581.000 [2015-Nov-05 20:44:09.784475]:info: parts shape: 15,4,4 0.0000 0.0000 69.0000 102.0000 :: simon:pmt blyth$ ggv --pmt 0:16 [2015-Nov-05 20:44:54.266290]:info: 0:/usr/local/env/optix/ggeo/bin/GPmtTest [2015-Nov-05 20:44:54.266963]:info: 1:0:16 [2015-Nov-05 20:44:54.267173]:info: NPY::make_slice from 16 -> 16 slice NSlice 0 : 16 : 1 [2015-Nov-05 20:44:54.267336]:info: GPmt::loadFromCache slicing partBuf origBuf 16,4,4 partBuf 16,4,4 GPmt::make_container pbb min -101.168 -101.168 -23.838 max 101.168 101.168 56.000 GPmt::make_container pbb min -101.168 -101.168 56.000 max 101.168 101.168 100.070 GPmt::make_container pbb min -84.540 -84.540 100.070 max 84.540 84.540 131.000 ... GPmt::make_container pbb min -98.143 -98.143 -30.000 max 98.143 98.143 56.000 GPmt::make_container pbb min -97.151 -97.151 -29.000 max 97.151 97.151 56.131 GPmt::make_container pbb min -27.500 -27.500 -164.500 max 27.500 27.500 1.500 GPmt::make_container pbb min -551.168 -551.168 -619.000 max 551.168 551.168 581.000 GPmt::make_container bb min -551.168 -551.168 -619.000 max 551.168 551.168 581.000 GPmt::make_container bb factor 3.0 min -2351.168 -2351.168 -2419.000 max 2351.168 2351.168 2381.000 [2015-Nov-05 20:44:54.267608]:info: parts shape: 16,4,4 Fixing box normals ------------------- After fixing ray box normals, get very pretty Lambertian render of PMT in box with *ggv-pmt* ie:: ggv-pmt () { ggv.sh --tracer --restrictmesh 1 --analyticmesh 1 --islice 0 --target 3199 $* } But the OptiX mode of *ggv-box* is far less pretty with nasty black faces, thats with:: ggv-box () { ggv --box --animtimemax 7 --boxconfig "size=2;boundary=MineralOil/Rock//;" --torchconfig "source=0,0,400;target=0,0,0;radius=102;zenith_azimuth=1,0,1,0" $* } Also photon reflections show non-symmetric behaviour, discriminating againt two of the box faces. How is that possible ? * different code in propagator and tracer ? * different geometry ? * normal issue or iimpinging other geometry ? :: ggv.sh --tracer --analyticmesh 1 --islice 0 --target 3199 $* # not restricting to instanced, see pretty render of analytic PMT with no extra box ? ggv.sh --tracer --islice 0 --target 3199 $* # triangulated PMT After fixing *ggv-box* mismatch, changing to *size=3* get the pretty render in OptiX mode and symmetric reflections:: ggv-box () { ggv --box --animtimemax 7 --boxconfig "size=3;boundary=MineralOil/Rock//;" --torchconfig "source=0,0,400;target=0,0,0;radius=102;zenith_azimuth=1,0,1,0" $* } * Explain that ? * Also, still material colors seem wrong. Face Slicing ------------- :: ggv-pmt --fslice 0:720 ggv-pmt --fslice 720:1392 ggv-pmt --fslice 1392:2352 ggv-pmt --fslice 2352:2832 ggv-pmt --fslice 2832:2928 # selecting faces of single solids, nodeinfo.npy provides the face index ranges :: In [1]: ni = np.load("GMergedMesh/1/nodeinfo.npy") In [2]: ni Out[2]: array([[ 720, 362, 3199, 3155], [ 672, 338, 3200, 3199], [ 960, 482, 3201, 3200], [ 480, 242, 3202, 3200], [ 96, 50, 3203, 3200]], dtype=uint32) In [3]: np.cumsum(ni[:,0]) Out[3]: array([ 720, 1392, 2352, 2832, 2928], dtype=uint64) Photocathode ------------- :: pmt-parts # move to writing full partition file, and pslicing as needed ggv-pmt --fslice 1392:2352 --pslice 8:10 ggv-pmt --fslice 1392:2352 --pslice 8:12 # after add inner spheres First and Second Solids, Pyrex and contained vacuum ------------------------------------------------------- OptiX render is as would expect, with pyrex and vacuum very thinly separated, to make the inner volume visible adjust near to control the ray trace epsilon OpenGL render not as would expect, much fatter to the back. As if pushed out by the dynode ? :: pmt-parts 0:8 ggv-pmt --fslice 0:1392 Tubs Issue FIXED, was caused by cylinder poking outside its bbox ------------------------------------------------------------------- * enable ENDCAP_P only in pmt-/dd.py and regen with :: pmt-parts 3:4 * setup coloring in cu/pinhole_camera.cu :: 100 // BGRA 101 uchar4 color = prd.flag == HP_PCAP_I ? RED : make_color( prd.result ); * get expected behavior for outer and inner HP_PCAP_O and HP_PCAP_I * PCAP endcap is to the right(in default initial ggv-pmt viewpoint) * doing the same for QCAP see view dependent shape mis-behaviour, but disabling the partition_union resetting of bbox avoids it * the problem was the bbox was clipped in at the 3spehere interseciton plane but ZSize was not changed * from point of view of cylinder rendering the relevant PQ vector is not (0,0,sizeZ) but rather (0,0,clipped_sizeZ) :: 194 static __device__ 195 void intersect_ztubs(quad& q0, quad& q1, quad& q2, quad& q3, const uint4& identity ) 196 { 197 /* 198 Position shift below is to match between different cylinder Z origin conventions 199 200 * Ericson calc implemented below has cylinder origin at endcap P 201 * detdesc/G4 Tubs has cylinder origin in the center 202 203 */ 204 float sizeZ = q1.f.x ; 205 float z0 = q0.f.z - sizeZ/2.f ; 206 float3 position = make_float3( q0.f.x, q0.f.y, z0 ); // 0,0,-169. 207 float clipped_sizeZ = q3.f.z - q2.f.z ; 208 209 float radius = q0.f.w ; 210 int flags = q1.i.w ; 211 212 bool PCAP = flags & ENDCAP_P ; 213 bool QCAP = flags & ENDCAP_Q ; 214 215 //rtPrintf("intersect_ztubs position %10.4f %10.4f %10.4f \n", position.x, position.y, position.z ); 216 //rtPrintf("intersect_ztubs flags %d PCAP %d QCAP %d \n", flags, PCAP, QCAP); 217 218 float3 m = ray.origin - position ; 219 float3 n = ray.direction ; 220 float3 d = make_float3(0.f, 0.f, clipped_sizeZ ); 221 222 float rr = radius*radius ; 223 float3 dnorm = normalize(d); 224 Just Tubs ---------- Some funny straight lines as rotate around:: pmt-parts 3:4 # just tubs ggv-pmt Either a bug or maybe optical illusion due to: * perspective projection * no depth/inside/outside queues Perhaps Z cut happening in wrong frame ? TODO: * get orthographic projection working for OptiX raygen * matplotlib projection plot of points of the mesh :: In [4]: v = np.load("GMergedMesh/1/vertices.npy") In [5]: v Out[5]: array([[ 0. , 0. , 131. ], [ 33.905, 0. , 126.536], [ 32.75 , 8.775, 126.536], ..., [ 26.563, -7.118, 1.5 ], [ 0. , 0. , 1.5 ], [ 0. , 0. , -164.5 ]], dtype=float32) In [6]: v.shape Out[6]: (1474, 3) In [7]: ni[:,1].sum() ## sum of vertices, it matches as these are fixed meshes with no dupes Out[7]: 1474 In [10]: i = np.load("GMergedMesh/1/indices.npy").reshape(-1,3) In [11]: i.shape Out[11]: (2928, 3) In [15]: np.unique(i[:720]).min() Out[15]: 0 In [16]: np.unique(i[:720]).max() Out[16]: 361 n [12]: ni[:,0].sum() Out[12]: 2928 In [19]: np.unique(i[:720]).size # hmm no need for doing indices look up into the vertices, its all contiguous Out[19]: 362 Just Tracing a single instance -------------------------------- Using OTracerTest with the below is much faster than with full context (including all those propagate buffers) and full geometry:: pmt-parts 0:4 # 3sphere + tubs ggv --tracer --restrictmesh 1 --analyticmesh 1 --islice 0 --target 3199 ggv-pmt # abbreviation for above ggv-allpmt --stack $((1024 + 512)) # stack can be reduced a bit with just the tracer ggv --tracer --restrictmesh 1 --analyticmesh 1 ggv-allpmt Plumbing check ---------------- :: ggv --restrictmesh 1 --analyticmesh 1 --torchconfig "radius=300;frame=3199;source=0,0,1000;target=0,0,0" How to OptiX intersect with CSG solid ? ----------------------------------------- :: simon:OptiX_380_sdk blyth$ find . -name '*.cu' -exec grep -l intersect {} \; ./ambocc/parallelogram.cu ./ambocc/sphere.cu ./buffersOfBuffers/parallelogram.cu ./buffersOfBuffers/sphere_texcoord.cu ./cook/clearcoat.cu ./cook/dof_camera.cu ./cook/parallelogram.cu ./cook/sphere.cu ./cook/sphere_texcoord.cu ./cuda/triangle_mesh.cu ./cuda/triangle_mesh_small.cu ./device_exceptions/device_exceptions.cu ./displacement/geometry_programs.cu ./glass/glass.cu ./glass/triangle_mesh_iterative.cu ./heightfield/heightfield.cu ./hybridShadows/triangle_mesh_fat.cu ./isgReflections/parallelogram.cu ./isgReflections/triangle_mesh_fat.cu ./isgShadows/triangle_mesh_fat.cu ./julia/block_floor.cu ./julia/julia.cu ... simon:OptiX_380_sdk blyth$ find . -type f -exec grep -l union {} \; ./julia/block_floor.cu ./julia/distance_field.h Julia sample has lots of non-trivial intersection examples julia/block_floor.cu:: 538 RT_PROGRAM void intersect(int primIdx) 539 { 540 object_factory::Object obj; 541 object_factory::make_object(obj, ray.direction); 542 543 // first check for intersection between the ray and aabb 544 optix::Ray tmp_ray = ray; 545 if(intersect_aabb(tmp_ray, obj)) { 546 float epsilon = 1.25e-3f; 547 float max_epsilon = 2.5e-2f; 548 549 float3 hit_point; 550 float t = adaptive_sphere_trace<1000>(tmp_ray, make_distance_to_primitive(obj), hit_point, epsilon, max_epsilon); 551 if(t < tmp_ray.tmax) 552 { 553 if(rtPotentialIntersection(t)) julia/distance_field.h:: 216 // The union of two primitives 217 template 218 class PrimitiveUnion 219 { 220 public: 221 // null constructor creates an undefined DistanceUnion 222 HD_DECL 223 PrimitiveUnion(void){} 224 225 HD_DECL 226 PrimitiveUnion(Primitive1 p1, Primitive2 p2):m_prim1(p1),m_prim2(p2){} 227 228 HD_DECL 229 float distance(const float3 &x) const 230 { 231 return fminf(m_prim1.distance(x), m_prim2.distance(x)); 232 } ... shadeTree/parallelogram.cu:: 37 RT_PROGRAM void intersect(int primIdx) 38 { 39 float3 n = make_float3( plane ); 40 float dt = dot(ray.direction, n ); 41 float t = (plane.w - dot(n, ray.origin))/dt; 42 if( t > ray.tmin && t < ray.tmax ) { 43 float3 p = ray.origin + ray.direction * t; 44 float3 vi = p - anchor; 45 float a1 = dot(v1, vi); 46 if(a1 >= 0 && a1 <= 1){ 47 float a2 = dot(v2, vi); 48 if(a2 >= 0 && a2 <= 1){ 49 if( rtPotentialIntersection( t ) ) { 50 geometric_normal = n; 51 shading_normal = n; 52 uv = make_float2(a1, a2); 53 rtReportIntersection( 0 ); 54 } 55 } 56 } 57 } 58 } tutorial.cpp:: 238 float4 make_plane( float3 n, float3 p ) 239 { 240 n = normalize(n); 241 float d = -dot(n, p); 242 return make_float4( n, d ); 243 } tutorial10.cu:: 313 // 314 // Intersection program for programmable convex hull primitive /// /// https://en.wikipedia.org/wiki/Line–plane_intersection /// http://geomalgorithms.com/index.html /// 315 // 316 rtBuffer planes; 317 RT_PROGRAM void chull_intersect(int primIdx) 318 { 319 int n = planes.size(); 320 float t0 = -FLT_MAX; 321 float t1 = FLT_MAX; 322 float3 t0_normal = make_float3(0); 323 float3 t1_normal = make_float3(0); 324 for(int i = 0; i < n && t0 < t1; ++i ) { 325 float4 plane = planes[i]; 326 float3 n = make_float3(plane); 327 float d = plane.w; 328 329 float denom = dot(n, ray.direction); 330 float t = -(d + dot(n, ray.origin))/denom; /// /// Plane eqn, p0 is point in plane, n is normal /// (p - p0).n = 0 /// /// Line /// p = ray.origin + t * ray.direction /// /// Intersect /// /// (ray.origin + t * ray.direction - p0 ).n = 0 /// /// dot(n, ray.origin) + t * dot(n, ray.direction) - dot(p0, n) = 0 /// /// dot(p0,n) - dot(n, ray.origin) /// t = -------------------------------- /// dot(n, ray.direction) /// /// 331 if( denom < 0){ 332 // enter 333 if(t > t0){ 334 t0 = t; 335 t0_normal = n; 336 } 337 } else { 338 //exit 339 if(t < t1){ 340 t1 = t; 341 t1_normal = n; 342 } 343 } 344 } 345 346 if(t0 > t1) 347 return; 348 349 if(rtPotentialIntersection( t0 )){ 350 shading_normal = geometric_normal = t0_normal; 351 rtReportIntersection(0); 352 } else if(rtPotentialIntersection( t1 )){ 353 shading_normal = geometric_normal = t1_normal; 354 rtReportIntersection(0); 355 } 356 } How to proceed ? ------------------ * on revisiting G4DAE include GDML G4 CSG model description together with the triangulated COLLADA detdesc PMT is involved ------------------------ Complicated assemblies of CSG solids. Implementing analytic is non-trivial. G5:/home/blyth/local/env/dyb/NuWa-trunk/dybgaudi/Detector/XmlDetDesc/DDDB/PMT/geometry.xml:: 08 09 10 11 12 13 14 15 16 .. dybgaudi/Detector/XmlDetDesc/DDDB/PMT/hemi-pmt.xml:: 37 38 39 40 41 43 44 46 47 48 50 51 52 53 56 57 58 59 61 62 :: 118 119 120 121 122 126 131 132 133