1 /*
2  * Copyright (C) 2020 Linux Studio Plugins Project <https://lsp-plug.in/>
3  *           (C) 2020 Vladimir Sadovnikov <sadko4u@gmail.com>
4  *
5  * This file is part of lsp-plugins
6  * Created on: 9 апр. 2017 г.
7  *
8  * lsp-plugins is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU Lesser General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * any later version.
12  *
13  * lsp-plugins is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  * GNU Lesser General Public License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public License
19  * along with lsp-plugins. If not, see <https://www.gnu.org/licenses/>.
20  */
21 
22 #include <core/alloc.h>
23 #include <core/3d/common.h>
24 #include <core/3d/RayTrace3D.h>
25 #include <stdlib.h>
26 #include <sys/time.h>
27 
28 #define SAMPLE_QUANTITY     512
29 #define TASK_LO_THRESH      0x2000
30 #define TASK_HI_THRESH      0x4000
31 
32 namespace lsp
33 {
34     static const size_t bbox_map[] =
35     {
36         0, 1, 2,
37         0, 2, 3,
38         6, 5, 4,
39         6, 4, 7,
40         1, 0, 4,
41         1, 4, 5,
42         3, 2, 6,
43         3, 6, 7,
44         1, 5, 2,
45         2, 5, 6,
46         0, 3, 4,
47         3, 7, 4
48     };
49 
50 
TaskThread(RayTrace3D * trace)51     RayTrace3D::TaskThread::TaskThread(RayTrace3D *trace)
52     {
53         this->trace     = trace;
54         heavy_state     = S_SCAN_OBJECTS;
55     }
56 
~TaskThread()57     RayTrace3D::TaskThread::~TaskThread()
58     {
59         // Cleanup capture state
60         for (size_t i=0; i<bindings.size(); ++i)
61         {
62             rt_binding_t *b     = bindings.at(i);
63             if (b == NULL)
64                 continue;
65             for (size_t j=0; j<b->bindings.size(); ++j)
66             {
67                 // Cleanup sample
68                 sample_t *samp  = b->bindings.at(j);
69                 if (samp->sample != NULL)
70                 {
71                     samp->sample->destroy();
72                     delete samp->sample;
73                     samp->sample = NULL;
74                 }
75             }
76             delete b;
77         }
78 
79         destroy_objects(&objects);
80         bindings.flush();
81     }
82 
run()83     status_t RayTrace3D::TaskThread::run()
84     {
85         // Initialize DSP context
86         dsp::context_t ctx;
87         dsp::start(&ctx);
88 
89         // Enter the main loop
90         status_t res = main_loop();
91         destroy_tasks(&tasks);
92         destroy_objects(&objects);
93 
94         // Finalize DSP context and return result
95         dsp::finish(&ctx);
96         return res;
97     }
98 
main_loop()99     status_t RayTrace3D::TaskThread::main_loop()
100     {
101         rt_context_t *ctx   = NULL;
102         bool report         = false;
103         status_t res        = STATUS_OK;
104 
105         // Perform main loop of raytracing
106         while (true)
107         {
108             // Check cancellation flag
109             if ((trace->bCancelled) || (trace->bFailed))
110             {
111                 res = STATUS_CANCELLED;
112                 break;
113             }
114 
115             // Try to fetch new task from internal queue
116             if (!tasks.pop(&ctx))
117             {
118                 // Are there any tasks left in global space?
119                 trace->lkTasks.lock();
120                 if (!trace->vTasks.pop(&ctx))
121                 {
122                     trace->lkTasks.unlock();
123                     break;
124                 }
125 
126                 // Update statistics
127                 if (trace->nQueueSize > trace->vTasks.size())
128                 {
129                     report      = true;
130                     trace->nQueueSize  = trace->vTasks.size();
131                 }
132                 ++stats.root_tasks;
133                 trace->lkTasks.unlock();
134             }
135             else
136                 ++stats.local_tasks;
137 
138             if (ctx == NULL)
139                 break;
140 
141             // Process context state
142             res     = process_context(ctx);
143 
144             // Report status if required
145             if ((res == STATUS_OK) && (report))
146             {
147                 report      = false;
148 
149                 trace->lkTasks.lock();
150                 float prg   = float(trace->nProgressPoints) / float(trace->nProgressMax);
151                 lsp_trace("Reporting progress %d/%d = %.2f%%", int(trace->nProgressPoints), int(trace->nProgressMax), prg * 100.0f);
152                 ++trace->nProgressPoints;
153                 res         = trace->report_progress(prg);
154                 trace->lkTasks.unlock();
155             }
156 
157             if (res != STATUS_OK)
158             {
159                 trace->bFailed  = true; // Report fail status if at least one thread has failed
160                 break;
161             }
162         }
163 
164         return res;
165     }
166 
submit_task(rt_context_t * ctx)167     status_t RayTrace3D::TaskThread::submit_task(rt_context_t *ctx)
168     {
169         // 'Heavy' state and pretty high number of pending tasks - submit task to global queue
170         if (ctx->state == heavy_state)
171         {
172             if (trace->vTasks.size() < TASK_LO_THRESH)
173             {
174                 trace->lkTasks.lock();
175                 status_t res = (trace->vTasks.push(ctx)) ? STATUS_OK : STATUS_NO_MEM;
176                 trace->lkTasks.unlock();
177                 return res;
178             }
179         }
180 
181         // Otherwise, submit to local task queue
182         return (tasks.push(ctx)) ? STATUS_OK : STATUS_NO_MEM;
183     }
184 
process_context(rt_context_t * ctx)185     status_t RayTrace3D::TaskThread::process_context(rt_context_t *ctx)
186     {
187         status_t res;
188 
189         switch (ctx->state)
190         {
191             case S_SCAN_OBJECTS:
192                 ++stats.calls_scan;
193                 res     = scan_objects(ctx);
194                 break;
195             case S_SPLIT:
196                 ++stats.calls_split;
197                 res     = split_view(ctx);
198                 break;
199             case S_CULL_BACK:
200                 ++stats.calls_cullback;
201                 res     = cullback_view(ctx);
202                 break;
203             case S_REFLECT:
204                 ++stats.calls_reflect;
205                 res     = reflect_view(ctx);
206                 break;
207             default:
208                 res = STATUS_BAD_STATE;
209                 break;
210         }
211 
212         // Analyze status
213         RT_TRACE(trace->pDebug,
214             if (res == STATUS_BREAK_POINT)
215             {
216                 trace->pDebug->ignored.swap(&ctx->ignored);
217                 trace->pDebug->trace.swap(&ctx->trace);
218             }
219         )
220 
221         // Force context to be deleted
222         if (res != STATUS_OK)
223             delete ctx;
224 
225         return res;
226     }
227 
generate_tasks(cvector<rt_context_t> * tasks,float initial)228     status_t RayTrace3D::TaskThread::generate_tasks(cvector<rt_context_t> *tasks, float initial)
229     {
230         status_t res;
231 
232         for (size_t i=0,n=trace->vSources.size(); i<n; ++i)
233         {
234             rt_source_settings_t *src = trace->vSources.get(i);
235             if (src == NULL)
236                 return STATUS_CORRUPTED;
237 
238             // Generate source mesh
239             cstorage<rt_group_t> groups;
240             res = rt_gen_source_mesh(groups, src);
241             if (res != STATUS_OK)
242                 return res;
243 
244             // Generate sources
245             matrix3d_t tm = src->pos;
246 
247             for (size_t ti=0, n=groups.size(); ti<n; ++ti)
248             {
249                 rt_group_t *grp = groups.at(ti);
250                 if (grp == NULL)
251                     continue;
252 
253                 rt_context_t *ctx   = new rt_context_t();
254                 if (ctx == NULL)
255                     return STATUS_NO_MEM;
256 
257                 RT_TRACE(trace->pDebug, ctx->set_debug_context(trace->pDebug); )
258 
259                 dsp::apply_matrix3d_mp2(&ctx->view.s, &grp->s, &tm);
260                 dsp::apply_matrix3d_mp2(&ctx->view.p[0], &grp->p[0], &tm);
261                 dsp::apply_matrix3d_mp2(&ctx->view.p[1], &grp->p[1], &tm);
262                 dsp::apply_matrix3d_mp2(&ctx->view.p[2], &grp->p[2], &tm);
263 
264                 ctx->state          = S_SCAN_OBJECTS;
265                 ctx->view.location  = 1.0f;
266                 ctx->view.oid       = -1;
267                 ctx->view.face      = -1;
268                 ctx->view.speed     = SOUND_SPEED_M_S;
269 
270                 ctx->view.amplitude = src->amplitude;
271                 ctx->view.time[0]   = 0.0f;
272                 ctx->view.time[1]   = 0.0f;
273                 ctx->view.time[2]   = 0.0f;
274 
275                 if (!tasks->add(ctx))
276                 {
277                     delete ctx;
278                     return STATUS_NO_MEM;
279                 }
280 
281                 RT_TRACE_BREAK(trace->pDebug,
282                     lsp_trace("Generated %d raytrace contexts for source %d", int(groups.size()), int(i));
283                     for (size_t i=0, n=tasks->size(); i<n; ++i)
284                     {
285                         rt_context_t *ctx = tasks->at(i);
286                         trace->pDebug->trace.add_view_1c(&ctx->view, &C3D_MAGENTA);
287                     }
288                 );
289             }
290         }
291 
292         return STATUS_OK;
293     }
294 
check_object(rt_context_t * ctx,Object3D * obj,const matrix3d_t * m)295     status_t RayTrace3D::TaskThread::check_object(rt_context_t *ctx, Object3D *obj, const matrix3d_t *m)
296     {
297         // Ensure that we need to perform additional checks
298         if (obj->num_triangles() < 16)
299             return STATUS_OK;
300 
301         // Prepare bounding-box check
302         bound_box3d_t box = *(obj->bound_box());
303         for (size_t j=0; j<8; ++j)
304             dsp::apply_matrix3d_mp1(&box.p[j], m);
305 
306         RT_TRACE_BREAK(trace->pDebug,
307             lsp_trace("Testing bound box");
308 
309             ctx->trace.add_view_1c(&ctx->view, &C3D_MAGENTA);
310             v_vertex3d_t v[3];
311             for (size_t j=0, m = sizeof(bbox_map)/sizeof(size_t); j < m; )
312             {
313                 v[0].p      = box.p[bbox_map[j++]];
314                 v[0].c      = C3D_YELLOW;
315                 v[1].p      = box.p[bbox_map[j++]];
316                 v[1].c      = C3D_YELLOW;
317                 v[2].p      = box.p[bbox_map[j++]];
318                 v[2].c      = C3D_YELLOW;
319 
320                 dsp::calc_normal3d_p3(&v[0].n, &v[0].p, &v[1].p, &v[2].p);
321                 v[1].n      = v[0].n;
322                 v[2].n      = v[0].n;
323 
324                 ctx->trace.add_triangle(v);
325             }
326         )
327 
328         // Perform simple bounding-box check
329         bool res    = check_bound_box(&box, &ctx->view);
330 
331 #ifdef LSP_RT_TRACE
332         if (!res)
333         {
334             RT_TRACE(trace->pDebug,
335                 for (size_t j=0,n=obj->num_triangles(); j<n; ++j)
336                 {
337                     obj_triangle_t *st = obj->triangle(j);
338 
339                     v_triangle3d_t t;
340                     dsp::apply_matrix3d_mp2(&t.p[0], st->v[0], m);
341                     dsp::apply_matrix3d_mp2(&t.p[1], st->v[1], m);
342                     dsp::apply_matrix3d_mp2(&t.p[2], st->v[2], m);
343 
344                     dsp::apply_matrix3d_mv2(&t.n[0], st->n[0], m);
345                     dsp::apply_matrix3d_mv2(&t.n[1], st->n[1], m);
346                     dsp::apply_matrix3d_mv2(&t.n[2], st->n[2], m);
347 
348                     ctx->ignored.add(&t);
349                 }
350             );
351         }
352 #endif /* LSP_RT_TRACE */
353 
354         return (res) ? STATUS_OK : STATUS_SKIP;
355     }
356 
generate_root_mesh()357     status_t RayTrace3D::TaskThread::generate_root_mesh()
358     {
359         status_t res;
360         size_t obj_id = 0;
361 
362         // Clear contents of the root context
363         rt_mesh_t root;
364         RT_TRACE(trace->pDebug, root.set_debug_context(trace->pDebug, &trace->pDebug->trace); );
365 
366         // Add capture objects as fake icosphere objects
367         for (size_t i=0, n=trace->vCaptures.size(); i<n; ++i, ++obj_id)
368         {
369             capture_t *cap      = trace->vCaptures.get(i);
370             if (cap == NULL)
371                 return STATUS_BAD_STATE;
372             if ((res = generate_capture_mesh(obj_id, cap)) != STATUS_OK)
373                 return res;
374         }
375 
376         // Add scene objects
377         for (size_t i=0, oid=obj_id, n=trace->pScene->num_objects(); i<n; ++i, ++oid)
378         {
379             // Get object
380             Object3D *obj       = trace->pScene->object(i);
381             if (obj == NULL)
382                 return STATUS_BAD_STATE;
383             else if (!obj->is_visible()) // Skip invisible objects
384                 continue;
385 
386             // Get material
387             rt_material_t *m    = trace->vMaterials.get(i);
388             if (m == NULL)
389                 return STATUS_BAD_STATE;
390 
391             // Add object to context
392             res         = root.add_object(obj, oid, obj->matrix(), m);
393             if (res != STATUS_OK)
394                 return res;
395         }
396 
397         RT_TRACE_BREAK(trace->pDebug,
398             lsp_trace("Prepared scene (%d triangles)", int(root.triangle.size()));
399             for (size_t i=0,n=root.triangle.size(); i<n; ++i)
400                 trace->pDebug->trace.add_triangle_3c(root.triangle.get(i), &C3D_RED, &C3D_GREEN, &C3D_BLUE);
401         );
402 
403         // Solve conflicts between all objects
404         res = root.solve_conflicts();
405         if (res != STATUS_OK)
406             return res;
407 
408         lsp_trace("Overall mesh statistics: %d vertexes, %d edges, %d triangles",
409                 int(root.vertex.size()), int(root.edge.size()), int(root.triangle.size()));
410 
411         // Generate object meshes
412         destroy_objects(&objects);
413         for (size_t i=0, n=trace->pScene->num_objects(); i<n; ++i, ++obj_id)
414         {
415             // Get object
416             Object3D *obj       = trace->pScene->object(i);
417             if (obj == NULL)
418                 return STATUS_BAD_STATE;
419             else if (!obj->is_visible()) // Skip invisible objects
420                 continue;
421 
422             // Create object
423             rt_object_t *rt = new rt_object_t();
424             if (rt == NULL)
425                 return STATUS_NO_MEM;
426             else if (!objects.add(rt)) {
427                 delete rt;
428                 return STATUS_NO_MEM;
429             }
430 
431             // Compute object's bounding box
432             obj->calc_bound_box();
433             if ((res = generate_object_mesh(obj_id, rt, &root, obj, obj->matrix())) != STATUS_OK)
434                 return res;
435         }
436 
437         RT_TRACE(trace->pDebug,
438             if (!trace->pScene->validate())
439                 return STATUS_CORRUPTED;
440         )
441 
442         RT_TRACE_BREAK(trace->pDebug,
443             lsp_trace("Added capture objects (%d triangles)", int(root.triangle.size()));
444             for (size_t i=0,n=root.triangle.size(); i<n; ++i)
445                 trace->pDebug->trace.add_triangle_3c(root.triangle.get(i), &C3D_RED, &C3D_GREEN, &C3D_BLUE);
446         );
447 
448         return res;
449     }
450 
generate_capture_mesh(size_t id,capture_t * c)451     status_t RayTrace3D::TaskThread::generate_capture_mesh(size_t id, capture_t *c)
452     {
453         status_t res;
454         cstorage<raw_triangle_t> mesh;
455 
456         // Generate capture mesh
457         if ((res = rt_gen_capture_mesh(mesh, c)) != STATUS_OK)
458             return res;
459 
460         // Generate bound box for capture and it's mesh
461         bound_box3d_t *b    = &c->bbox;
462         float r             = c->radius;
463         dsp::init_point_xyz(&b->p[0], -r, r, r);
464         dsp::init_point_xyz(&b->p[1], -r, -r, r);
465         dsp::init_point_xyz(&b->p[2], r, -r, r);
466         dsp::init_point_xyz(&b->p[3], r, r, r);
467         dsp::init_point_xyz(&b->p[4], -r, r, -r);
468         dsp::init_point_xyz(&b->p[5], -r, -r, -r);
469         dsp::init_point_xyz(&b->p[6], r, -r, -r);
470         dsp::init_point_xyz(&b->p[7], r, r, -r);
471 
472         // Apply changes to bounding box
473         for (size_t j=0; j<8; ++j)
474             dsp::apply_matrix3d_mp1(&b->p[j], &c->pos);
475 
476         // Convert data
477         size_t count = mesh.size();
478         raw_triangle_t *src = mesh.get_array();
479         rt_triangle_t *dst = c->mesh.append_n(count);
480         if (dst == NULL)
481             return STATUS_NO_MEM;
482 
483         for (size_t j=0; j < count; ++j, ++src, ++dst)
484         {
485             // Generate points and fill additional fields
486             dsp::apply_matrix3d_mp2(&dst->v[0], &src->v[0], &c->pos);
487             dsp::apply_matrix3d_mp2(&dst->v[1], &src->v[1], &c->pos);
488             dsp::apply_matrix3d_mp2(&dst->v[2], &src->v[2], &c->pos);
489             dsp::calc_plane_pv(&dst->n, src->v);
490 
491             dst->oid        = id;
492             dst->face       = j;
493             dst->m          = NULL;
494         }
495 
496         return res;
497     }
498 
generate_object_mesh(ssize_t id,rt_object_t * o,rt_mesh_t * src,Object3D * obj,const matrix3d_t * m)499     status_t RayTrace3D::TaskThread::generate_object_mesh(ssize_t id, rt_object_t *o, rt_mesh_t *src, Object3D *obj, const matrix3d_t *m)
500     {
501         rtx_edge_t *e;
502 
503         // Initialize ptag
504         for (size_t i=0, n=src->edge.size(); i<n; ++i)
505             src->edge.get(i)->itag  = -1;
506 
507         // Build set of triangles and edges
508         ssize_t itag = 0;
509         for (size_t i=0, n=src->triangle.size(); i<n; ++i)
510         {
511             // Fetch triangle
512             const rtm_triangle_t *t = src->triangle.get(i);
513             if (t->oid != id)
514                 continue;
515 
516             // Add triangle
517             rtx_triangle_t *rt = o->mesh.add();
518             if (rt == NULL)
519                 return STATUS_NO_MEM;
520 
521             rt->v[0]    = *(t->v[0]);
522             rt->v[1]    = *(t->v[1]);
523             rt->v[2]    = *(t->v[2]);
524             rt->n       = t->n;
525             rt->oid     = t->oid;
526             rt->face    = t->face;
527             rt->m       = t->m;
528 
529             // Add edges to if they are still not present
530             for (size_t j=0; j<3; ++j)
531             {
532                 rtm_edge_t *se  = t->e[j];
533                 rt->e[j]        = reinterpret_cast<rtx_edge_t *>(se);
534 
535                 if (se->itag < 0)
536                 {
537                     if ((e = o->plan.add()) == NULL)
538                         return STATUS_NO_MEM;
539                     e->v[0]         = *(se->v[0]);
540                     e->v[1]         = *(se->v[1]);
541                     se->itag        = itag++;
542                 }
543             }
544         }
545 
546         // Patch edge pointers according to stored indexes
547         for (size_t i=0, n=o->mesh.size(); i<n; ++i)
548         {
549             rtx_triangle_t *rt = o->mesh.at(i);
550             for (size_t j=0; j<3; ++j)
551             {
552                 rtm_edge_t *se  = reinterpret_cast<rtm_edge_t *>(rt->e[j]);
553                 rt->e[j]        = o->plan.at(se->itag);
554             }
555         }
556 
557         // Apply changes to bound box
558         const obj_boundbox_t *bbox = obj->bound_box();
559         for (size_t i=0; i<8; ++i)
560             dsp::apply_matrix3d_mp2(&o->bbox.p[i], &bbox->p[i], m);
561 
562         return STATUS_OK;
563     }
564 
scan_objects(rt_context_t * ctx)565     status_t RayTrace3D::TaskThread::scan_objects(rt_context_t *ctx)
566     {
567         status_t res = STATUS_OK;
568 
569         RT_TRACE_BREAK(trace->pDebug,
570             lsp_trace("Scanning objects...");
571             ctx->trace.add_view_1c(&ctx->view, &C3D_MAGENTA);
572         )
573 
574         // Initialize context's view
575         ctx->init_view();
576 
577         // Initialize scan variables
578         size_t max_objs     = trace->pScene->num_objects() + trace->vCaptures.size();
579         size_t segs         = (max_objs + (sizeof(size_t) * 8) - 1) / (sizeof(size_t) * 8);
580         size_t *objs        = reinterpret_cast<size_t *>(alloca(sizeof(size_t) * segs));
581         for (size_t i=0; i<segs; ++i)
582             objs[i]             = 0;
583 
584         size_t n_objs       = 0;
585         size_t obj_id       = 0;
586 
587         // Add captures as opaque objects
588         obj_id = 0;
589         for (size_t i=0, n=trace->vCaptures.size(); i<n; ++i, ++obj_id)
590         {
591             capture_t *cap      = trace->vCaptures.at(i);
592             if (cap == NULL)
593                 return STATUS_BAD_STATE;
594 
595             RT_TRACE_BREAK(trace->pDebug,
596                 lsp_trace("Testing capture bound box");
597 
598                 ctx->trace.add_view_1c(&ctx->view, &C3D_MAGENTA);
599                 v_vertex3d_t v[3];
600                 for (size_t j=0, m = sizeof(bbox_map)/sizeof(size_t); j < m; )
601                 {
602                     v[0].p      = cap->bbox.p[bbox_map[j++]];
603                     v[0].c      = C3D_YELLOW;
604                     v[1].p      = cap->bbox.p[bbox_map[j++]];
605                     v[1].c      = C3D_YELLOW;
606                     v[2].p      = cap->bbox.p[bbox_map[j++]];
607                     v[2].c      = C3D_YELLOW;
608 
609                     dsp::calc_normal3d_p3(&v[0].n, &v[0].p, &v[1].p, &v[2].p);
610                     v[1].n      = v[0].n;
611                     v[2].n      = v[0].n;
612 
613                     ctx->trace.add_triangle(v);
614                 }
615             )
616 
617             if (!check_bound_box(&cap->bbox, &ctx->view))
618             {
619                 RT_TRACE_BREAK(trace->pDebug,
620                     ctx->trace.add_view_1c(&ctx->view, &C3D_MAGENTA);
621                     for (size_t j=0,m=cap->mesh.size(); j<m; ++j)
622                     {
623                         rt_triangle_t *st = cap->mesh.at(j);
624 
625                         v_triangle3d_t t;
626                         t.p[0]  = st->v[0];
627                         t.p[1]  = st->v[1];
628                         t.p[2]  = st->v[2];
629                         t.n[0]  = st->n;
630                         t.n[1]  = st->n;
631                         t.n[2]  = st->n;
632 
633                         ctx->ignored.add(&t);
634                     }
635                 );
636                 continue;
637             }
638 
639             // Add capture as opaque object
640             if ((res = ctx->add_opaque_object(cap->mesh.get_array(), cap->mesh.size())) != STATUS_OK)
641                 return res;
642         }
643 
644         // Iterate all object and add to the context if the object is potentially participating the ray tracing algorithm
645         for (size_t i=0, n=objects.size(); i<n; ++i)
646         {
647             rt_object_t *rt = objects.at(i);
648             if (rt == NULL)
649                 return STATUS_BAD_STATE;
650 
651             // Check bound box
652             bool check = rt->mesh.size() > 16;
653             if (check)
654             {
655                 RT_TRACE_BREAK(trace->pDebug,
656                     lsp_trace("Testing bound box");
657 
658                     ctx->trace.add_view_1c(&ctx->view, &C3D_MAGENTA);
659                     v_vertex3d_t v[3];
660                     for (size_t j=0, m = sizeof(bbox_map)/sizeof(size_t); j < m; )
661                     {
662                         v[0].p      = rt->bbox.p[bbox_map[j++]];
663                         v[0].c      = C3D_YELLOW;
664                         v[1].p      = rt->bbox.p[bbox_map[j++]];
665                         v[1].c      = C3D_YELLOW;
666                         v[2].p      = rt->bbox.p[bbox_map[j++]];
667                         v[2].c      = C3D_YELLOW;
668 
669                         dsp::calc_normal3d_p3(&v[0].n, &v[0].p, &v[1].p, &v[2].p);
670                         v[1].n      = v[0].n;
671                         v[2].n      = v[0].n;
672 
673                         ctx->trace.add_triangle(v);
674                     }
675                 )
676 
677                 // Perform simple bounding-box check
678                 check = !check_bound_box(&rt->bbox, &ctx->view);
679             }
680             if (check)
681             {
682                 RT_TRACE(trace->pDebug,
683                     for (size_t j=0,m=rt->mesh.size(); j<m; ++j)
684                     {
685                         rtx_triangle_t *st = rt->mesh.at(j);
686 
687                         v_triangle3d_t t;
688                         t.p[0]  = st->v[0];
689                         t.p[1]  = st->v[1];
690                         t.p[2]  = st->v[2];
691                         t.n[0]  = st->n;
692                         t.n[1]  = st->n;
693                         t.n[2]  = st->n;
694 
695                         ctx->ignored.add(&t);
696                     }
697                 );
698                 continue;
699             }
700 
701             // Add object to context
702             res = ctx->add_object(rt->mesh.get_array(), rt->plan.get_array(), rt->mesh.size(), rt->plan.size());
703             if (res != STATUS_OK)
704                 return res;
705             ++n_objs;
706         }
707 
708         RT_TRACE_BREAK(trace->pDebug,
709             lsp_trace("Fetched %d objects", int(n_objs));
710             ctx->trace.add_view_1c(&ctx->view, &C3D_MAGENTA);
711             for (size_t i=0,n=ctx->triangle.size(); i<n; ++i)
712                 ctx->trace.add_triangle_3c(ctx->triangle.get(i), &C3D_RED, &C3D_GREEN, &C3D_BLUE);
713 
714             for (size_t i=0, n=ctx->plan.items.size(); i<n; ++i)
715                 ctx->trace.add_segment(ctx->plan.items.get(i), &C3D_YELLOW);
716         )
717 
718         // Update state
719         if (!ctx->plan.is_empty())
720             ctx->state  = S_SPLIT;
721         else if (ctx->triangle.size() <= 0)
722         {
723             delete ctx;
724             return STATUS_OK;
725         }
726         else
727             ctx->state  = S_REFLECT;
728 
729         return submit_task(ctx);
730     }
731 
cull_view(rt_context_t * ctx)732     status_t RayTrace3D::TaskThread::cull_view(rt_context_t *ctx)
733     {
734         status_t res = ctx->cull_view();
735         if (res != STATUS_OK)
736             return res;
737 
738         if (!ctx->plan.is_empty())
739             ctx->state  = S_SPLIT;
740         else if (ctx->triangle.size() <= 0)
741         {
742             delete ctx;
743             return STATUS_OK;
744         }
745         else
746             ctx->state  = S_REFLECT;
747 
748         return submit_task(ctx);
749     }
750 
split_view(rt_context_t * ctx)751     status_t RayTrace3D::TaskThread::split_view(rt_context_t *ctx)
752     {
753         rt_context_t out;
754         RT_TRACE(trace->pDebug, out.set_debug_context(trace->pDebug); );
755 
756         // Perform binary split
757         status_t res = ctx->edge_split(&out);
758         if (res == STATUS_NOT_FOUND)
759         {
760             ctx->state      = S_CULL_BACK;
761             return submit_task(ctx);
762         }
763         else if (res != STATUS_OK)
764             return res;
765 
766         // Analyze state of current and out context
767         if (ctx->triangle.size() > 0)
768         {
769             // Analyze state of 'out' context
770             if (out.triangle.size() > 0)
771             {
772                 // Allocate additional context and add to task list
773                 rt_context_t *nctx = new rt_context_t(&ctx->view, (out.triangle.size() > 1) ? S_SPLIT : S_REFLECT);
774                 if (nctx == NULL)
775                     return STATUS_NO_MEM;
776 
777                 RT_TRACE(trace->pDebug,
778                     nctx->set_debug_context(trace->pDebug);
779 
780                     nctx->ignored.add_all(&ctx->ignored);
781                     nctx->trace.add_all(&ctx->trace);
782 
783                     for (size_t i=0,n=ctx->triangle.size(); i<n; ++i)
784                         nctx->ignore(ctx->triangle.get(i));
785                     for (size_t i=0,n=out.triangle.size(); i<n; ++i)
786                         ctx->ignore(out.triangle.get(i));
787                 );
788 
789                 nctx->swap(&out);
790 
791                 // Submit task
792                 res = submit_task(nctx);
793                 if (res != STATUS_OK)
794                 {
795                     delete nctx;
796                     return res;
797                 }
798             }
799 
800             // Update context state
801             ctx->state  = (ctx->plan.is_empty()) ? S_REFLECT : S_SPLIT;
802             return submit_task(ctx);
803         }
804         else if (out.triangle.size() > 0)
805         {
806             ctx->swap(&out);
807             ctx->state  = (ctx->plan.is_empty()) ? S_REFLECT : S_SPLIT;
808 
809             return submit_task(ctx);
810         }
811 
812         delete ctx;
813         return STATUS_OK;
814     }
815 
cullback_view(rt_context_t * ctx)816     status_t RayTrace3D::TaskThread::cullback_view(rt_context_t *ctx)
817     {
818         RT_TRACE_BREAK(trace->pDebug,
819             lsp_trace("Performing depth test...");
820             ctx->trace.add_view_1c(&ctx->view, &C3D_MAGENTA);
821             for (size_t i=0,n=ctx->triangle.size(); i<n; ++i)
822                 ctx->trace.add_triangle_1c(ctx->triangle.get(i), &C3D_RED);
823         )
824 
825         // Perform cullback
826         status_t res = ctx->depth_test();
827         if (res != STATUS_OK)
828             return res;
829 
830         RT_TRACE_BREAK(trace->pDebug,
831             lsp_trace("After depth test %d triangles", int(ctx->triangle.size()));
832         ctx->trace.add_view_1c(&ctx->view, &C3D_MAGENTA);
833             for (size_t i=0,n=ctx->triangle.size(); i<n; ++i)
834                 ctx->trace.add_triangle_1c(ctx->triangle.get(i), &C3D_GREEN);
835         );
836 
837         if (ctx->triangle.size() <= 0)
838         {
839             delete ctx;
840             return STATUS_OK;
841         }
842 
843         ctx->state  = S_REFLECT;
844         return submit_task(ctx);
845     }
846 
847 #ifdef LSP_TRACE
invalid_state_hook()848     void invalid_state_hook()
849     {
850         lsp_error("Invalid state");
851     }
852 #endif
853 
reflect_view(rt_context_t * ctx)854     status_t RayTrace3D::TaskThread::reflect_view(rt_context_t *ctx)
855     {
856         rt_context_t *rc;
857         rt_view_t sv, v, rv, tv;        // source view, view, captured view, reflected view, transparent trace
858         vector3d_t vpl;                 // trace plane, split plane
859         point3d_t p[3];                 // Projection points
860         float d[3], t[3];               // distance, time
861         float a[3], A, kd;              // particular area, area, dispersion coefficient
862 
863         sv      = ctx->view;
864         A       = dsp::calc_area_pv(sv.p);
865 
866         if (A <= trace->fTolerance)
867         {
868             delete ctx;
869             return STATUS_OK;
870         }
871         float revA = 1.0f / A;
872 
873         dsp::calc_plane_pv(&vpl, ctx->view.p);
874         status_t res    = STATUS_OK;
875 
876         for (size_t i=0,n=ctx->triangle.size(); i<n; ++i)
877         {
878             // Get current triangle to perform reflection
879             rt_triangle_t *ct = ctx->triangle.get(i);
880 
881             // Estimate co-location of source point and reflecting triangle
882             float distance  = sv.s.x * ct->n.dx + sv.s.y * ct->n.dy + sv.s.z * ct->n.dz + ct->n.dw;
883 
884             // We need to skip the interaction with triangle due two reasons:
885             //   1. The ray group reaches the face with non-matching normal
886             //   2. The ray group reaches the face with matching normal but invalid object identifier
887             if (distance > 0.0f) // we're entering from outside the object inside the object
888             {
889                 // The source point should be above the triangle
890                 if (sv.location <= 0.0f)
891                     continue;
892             }
893             else if (distance < 0.0f)
894             {
895                 // The source point should be below triangle and match the object identifier
896                 if ((sv.location >= 0.0f) || (sv.oid != ct->oid))
897                     continue;
898             }
899             else
900                 continue;
901 
902             // Estimate the start time for each trace point using barycentric coordinates
903             // Additionally filter invalid triangles
904             bool valid = true;
905             for (size_t j=0; j<3; ++j)
906             {
907                 dsp::calc_split_point_p2v1(&p[j], &sv.s, &ct->v[j], &vpl);    // Project triangle point to trace point
908                 d[j]        = dsp::calc_distance_p2(&p[j], &ct->v[j]);        // Compute distance between projected point and triangle point
909 
910                 a[0]        = dsp::calc_area_p3(&p[j], &sv.p[1], &sv.p[2]);   // Compute area 0
911                 a[1]        = dsp::calc_area_p3(&p[j], &sv.p[0], &sv.p[2]);   // Compute area 1
912                 a[2]        = dsp::calc_area_p3(&p[j], &sv.p[0], &sv.p[1]);   // Compute area 2
913 
914                 float dA    = A - (a[0] + a[1] + a[2]); // Point should lay within the view
915                 if ((dA <= -trace->fTolerance) || (dA >= trace->fTolerance))
916                 {
917                     valid = false;
918                     break;
919                 }
920 
921                 t[j]        = (sv.time[0] * a[0] + sv.time[1] * a[1] + sv.time[2] * a[2]) * revA; // Compute projected point's time
922                 v.time[j]   = t[j] + (d[j] / sv.speed);
923                 #ifdef LSP_TRACE
924                     if (v.time[j] > 30.0f)
925                         invalid_state_hook();
926                 #endif
927             }
928 
929             if (!valid)
930                 continue;
931 
932             // Compute area of projected triangle
933             float area  = dsp::calc_area_pv(p);
934             if (area <= trace->fDetalization)
935                 continue;
936 
937             // Determine the direction from which came the wave front
938             v.oid       = ct->oid;
939             v.face      = ct->face;
940             v.s         = sv.s;
941             v.amplitude = sv.amplitude * sqrtf(area * revA);
942             v.location  = sv.location;
943             v.speed     = sv.speed;
944             v.rnum      = sv.rnum;
945             v.p[0]      = ct->v[0];
946             v.p[1]      = ct->v[1];
947             v.p[2]      = ct->v[2];
948 
949             RT_TRACE_BREAK(trace->pDebug,
950                 lsp_trace("Projection points distance: {%f, %f, %f}", d[0], d[1], d[2]);
951                 lsp_trace("View points time: {%f, %f, %f}", sv.time[0], sv.time[1], sv.time[2]);
952                 lsp_trace("Projection points time: {%f, %f, %f}", t[0], t[1], t[2]);
953                 lsp_trace("Target points time: {%f, %f, %f}", v.time[0], v.time[1], v.time[2]);
954                 lsp_trace("Amplitude: %e -> %e", sv.amplitude, v.amplitude);
955                 lsp_trace("Distance between source point and triangle: %f", distance);
956                 lsp_trace("oid = %d, caprues.size=%d", int(ct->oid), int(trace->vCaptures.size()));
957 
958                 ctx->trace.add_triangle_1c(ct, &C3D_YELLOW);
959                 ctx->trace.add_plane_3pn1c(&ct->v[0], &ct->v[1], &ct->v[2], &ct->n, &C3D_GREEN);
960                 ctx->trace.add_view_1c(&sv, &C3D_MAGENTA);
961                 ctx->trace.add_segment(&p[0], &ct->v[0], &C3D_RED);
962                 ctx->trace.add_segment(&p[1], &ct->v[1], &C3D_GREEN);
963                 ctx->trace.add_segment(&p[2], &ct->v[2], &C3D_BLUE);
964             )
965 
966             // Perform capture/reflect
967             capture_t *cap  = trace->vCaptures.get(ct->oid);
968             if (cap != NULL)
969             {
970                 rt_binding_t *b = bindings.get(ct->oid);
971                 if (b != NULL)
972                 {
973                     // Perform synchronized capturing
974                     ++stats.calls_capture;
975                     #ifdef LSP_RT_TRACE
976                         res = capture(cap, &b->bindings, &v, &ctx->trace);
977                     #else
978                         res = capture(cap, &b->bindings, &v);
979                     #endif /* LSP_RT_TRACE */
980                 }
981                 else
982                     res = STATUS_CORRUPTED;
983             }
984             else
985             {
986                 // Get material
987                 rt_material_t *m    = ct->m;
988 
989                 // Compute reflected and refracted views
990                 rv          = v;
991                 tv          = v;
992 
993                 if (distance > 0.0f)
994                 {
995                     v.amplitude    *= (1.0f - m->absorption[0]);
996 
997                     kd              = (1.0f + 1.0f / m->diffusion[0]) * distance;
998                     rv.amplitude    = v.amplitude * (m->transparency[0] - 1.0f); // Sign negated
999                     rv.s.x         -= kd * ct->n.dx;
1000                     rv.s.y         -= kd * ct->n.dy;
1001                     rv.s.z         -= kd * ct->n.dz;
1002                     rv.rnum         = v.rnum + 1;       // Increment reflection number
1003 
1004                     kd              = (m->permeability/m->dispersion[0] - 1.0f) * distance;
1005                     tv.amplitude    = v.amplitude * m->transparency[0];
1006                     tv.speed       *= m->permeability;
1007                     tv.s.x         += kd * ct->n.dx;
1008                     tv.s.y         += kd * ct->n.dy;
1009                     tv.s.z         += kd * ct->n.dz;
1010                     tv.location     = - v.location;     // Invert location of transparent trace
1011 
1012                     #ifdef LSP_TRACE
1013                         if (tv.speed > 10000.0f)
1014                             invalid_state_hook();
1015                         else if (tv.speed < 10.0f)
1016                             invalid_state_hook();
1017                     #endif
1018 
1019                     RT_TRACE_BREAK(trace->pDebug,
1020                         lsp_trace("Outside->inside reflect_view");
1021                         lsp_trace("Amplitude: reflected=%e, refracted=%e", rv.amplitude, tv.amplitude);
1022                         lsp_trace("Material: absorption=%f, transparency=%f, permeability=%f, diffusion=%f, dispersion=%f",
1023                                 m->absorption[0], m->transparency[0], m->permeability, m->diffusion[0], m->dispersion[0]);
1024 
1025                         ctx->trace.add_view_1c(&sv, &C3D_RED);
1026                         ctx->trace.add_view_1c(&rv, &C3D_GREEN);
1027                         ctx->trace.add_view_1c(&tv, &C3D_CYAN);
1028                     )
1029                 }
1030                 else
1031                 {
1032                     v.amplitude    *= (1.0f - m->absorption[1]);
1033 
1034                     kd              = (1.0f + 1.0f / m->diffusion[1]) * distance;
1035                     rv.amplitude    = v.amplitude * (m->transparency[1] - 1.0f); // Sign negated
1036                     rv.s.x         -= kd * ct->n.dx;
1037                     rv.s.y         -= kd * ct->n.dy;
1038                     rv.s.z         -= kd * ct->n.dz;
1039                     rv.rnum         = v.rnum + 1;       // Increment reflection number
1040 
1041                     kd              = (1.0f/(m->dispersion[1]*m->permeability) - 1.0f) * distance;
1042                     tv.amplitude    = v.amplitude * m->transparency[1];
1043                     tv.speed       /= m->permeability;
1044                     tv.s.x         += kd * ct->n.dx;
1045                     tv.s.y         += kd * ct->n.dy;
1046                     tv.s.z         += kd * ct->n.dz;
1047                     tv.location     = - v.location;     // Invert location of transparent trace
1048 
1049                     #ifdef LSP_TRACE
1050                         if (tv.speed > 10000.0f)
1051                             invalid_state_hook();
1052                         else if (tv.speed < 10.0f)
1053                             invalid_state_hook();
1054                     #endif
1055 
1056                     RT_TRACE_BREAK(trace->pDebug,
1057                         lsp_trace("Inside->outside reflect_view");
1058                         lsp_trace("Amplitude: reflected=%e, refracted=%e", rv.amplitude, tv.amplitude);
1059                         lsp_trace("Material: absorption=%f, transparency=%f, permeability=%f, diffusion=%f, dispersion=%f",
1060                                 m->absorption[1], m->transparency[1], m->permeability, m->diffusion[1], m->dispersion[1]);
1061 
1062                         ctx->trace.add_view_1c(&sv, &C3D_RED);
1063                         ctx->trace.add_view_1c(&rv, &C3D_GREEN);
1064                         ctx->trace.add_view_1c(&tv, &C3D_CYAN);
1065                     )
1066                 }
1067 
1068                 // Create reflection context
1069                 if ((rv.amplitude <= -trace->fEnergyThresh) || (rv.amplitude >= trace->fEnergyThresh))
1070                 {
1071                     // Revert the order of triangle because direction has changed to opposite
1072                     rv.p[1]     = v.p[2];
1073                     rv.p[2]     = v.p[1];
1074 
1075                     if ((rc = new rt_context_t(&rv, S_SCAN_OBJECTS)) != NULL)
1076                     {
1077                         RT_TRACE(trace->pDebug,
1078                             rc->set_debug_context(trace->pDebug);
1079                         );
1080 
1081                         if ((res = submit_task(rc)) != STATUS_OK)
1082                             delete rc;
1083                     }
1084                     else
1085                         res = STATUS_NO_MEM;
1086                 }
1087 
1088                 // Create refraction context
1089                 if ((tv.amplitude <= -trace->fEnergyThresh) || (tv.amplitude >= trace->fEnergyThresh))
1090                 {
1091                     if ((rc = new rt_context_t(&tv, S_SCAN_OBJECTS)) != NULL)
1092                     {
1093                         RT_TRACE(trace->pDebug,
1094                             rc->set_debug_context(trace->pDebug);
1095                         );
1096 
1097                         if ((res = submit_task(rc)) != STATUS_OK)
1098                             delete rc;
1099                     }
1100                     else
1101                         res = STATUS_NO_MEM;
1102                 }
1103             }
1104 
1105             // Analyze result
1106             if (res != STATUS_OK)
1107                 break;
1108         }
1109 
1110         // Drop current context on OK status
1111         if (res == STATUS_OK)
1112             delete ctx;
1113         return res;
1114     }
1115 
1116 #ifdef LSP_RT_TRACE
capture(capture_t * capture,cstorage<sample_t> * bindings,const rt_view_t * v,View3D * view)1117     status_t RayTrace3D::TaskThread::capture(capture_t *capture, cstorage<sample_t> *bindings, const rt_view_t *v, View3D *view)
1118 #else
1119     status_t RayTrace3D::TaskThread::capture(capture_t *capture, cstorage<sample_t> *bindings, const rt_view_t *v)
1120 #endif /* LSP_RT_TRACE */
1121     {
1122 //        lsp_trace("Capture:\n"
1123 //                "  coord  = {{%f, %f, %f}, {%f, %f, %f}, {%f, %f, %f}}\n"
1124 //                "  time   = {%f, %f, %f}\n"
1125 //                "  energy = %f",
1126 //                v->p[0].x, v->p[0].y, v->p[0].z,
1127 //                v->p[1].x, v->p[1].y, v->p[1].z,
1128 //                v->p[2].x, v->p[2].y, v->p[2].z,
1129 //                v->time[0], v->time[1], v->time[2],
1130 //                v->energy
1131 //            );
1132 
1133         // Compute the area of triangle
1134         float v_area = dsp::calc_area_pv(v->p);
1135         if (v_area <= trace->fDetalization) // Prevent from becoming NaNs
1136             return STATUS_OK;
1137 
1138         vector3d_t cv, pv;
1139         float afactor       = v->amplitude / sqrtf(v_area); // The norming energy factor
1140         dsp::unit_vector_p1pv(&cv, &v->s, v->p);
1141         pv                  = capture->direction;
1142         float kcos          = cv.dx*pv.dx + cv.dy*pv.dy + cv.dz * pv.dz; // cos(a)
1143 
1144         // Analyze capture type
1145         switch (capture->type)
1146         {
1147             case RT_AC_CARDIO:
1148                 afactor    *= 0.5f * (1.0f - kcos); // 0.5 * (1 + cos(a))
1149                 break;
1150 
1151             case RT_AC_SCARDIO:
1152                 afactor    *= 2*fabs(0.5 - kcos)/3.0f; // 2*(0.5 + cos(a)) / 3
1153                 break;
1154 
1155             case RT_AC_HCARDIO:
1156                 afactor    *= 0.8*fabs(0.25 - kcos); // 4*(0.25 + cos(a)) / 5
1157                 break;
1158 
1159             case RT_AC_BIDIR:
1160                 afactor    *= kcos;    // factor = factor * cos(a)
1161                 break;
1162 
1163             case RT_AC_EIGHT:
1164                 afactor    *= kcos*kcos;  // factor = factor * cos(a)^2
1165                 break;
1166 
1167             case RT_AC_OMNI: // factor = factor * 1
1168             default:
1169                 break;
1170         }
1171 
1172         // Estimate distance and time parameters for source point
1173         vector3d_t ds[3];
1174         raw_triangle_t src;
1175         float ts[3], tsn[3];
1176 
1177         for (size_t i=0; i<3; ++i)
1178         {
1179             src.v[i]    = v->p[i];
1180             dsp::init_vector_p2(&ds[i], &v->s, &src.v[i]);  // Delta vector
1181             float dist  = dsp::calc_distance_v1(&ds[i]);
1182             ts[i]       = v->time[i] - dist / v->speed;     // Time at the source point
1183             tsn[i]      = v->time[i] * trace->nSampleRate;  // Sample number
1184         }
1185 
1186         // Estimate the culling sample number
1187         ssize_t csn;
1188         if ((tsn[0] < tsn[1]) && (tsn[0] < tsn[2]))
1189             csn     = tsn[0];
1190         else
1191             csn     = (tsn[1] < tsn[2]) ? tsn[1] : tsn[2];
1192 #ifdef LSP_TRACE
1193 //        ssize_t ssn = csn;                                      // Start culling sample number
1194 #endif
1195         ++csn;
1196 
1197         // Perform integration
1198         vector3d_t spl;
1199         raw_triangle_t in[2], out[2];
1200         size_t n_in, n_out;
1201         float prev_area     = 0.0f;                             // The area of polygon at previous step
1202         point3d_t p[3];
1203 
1204         do {
1205             // Estimate the culling plane points
1206             float ctime     = float(csn) / trace->nSampleRate;
1207             for (size_t i=0; i<3; ++i)
1208             {
1209                 float factor    = (ctime  - ts[i]) / (v->time[i] - ts[i]);
1210                 p[i].x  = v->s.x + ds[i].dx * factor;
1211                 p[i].y  = v->s.y + ds[i].dy * factor;
1212                 p[i].z  = v->s.z + ds[i].dz * factor;
1213                 p[i].w  = 1.0f;
1214             }
1215 
1216             // Compute culling plane
1217             dsp::calc_oriented_plane_pv(&spl, &v->s, p);
1218 
1219             RT_TRACE_BREAK(trace->pDebug,
1220                 lsp_trace("Integrating triangle...");
1221                 view->add_view_1c(v, &C3D_MAGENTA);
1222                 view->add_triangle_pv1c(src.v, &C3D_ORANGE);
1223                 view->add_plane_pvn1c(p, &spl, &C3D_YELLOW);
1224             )
1225 
1226             // Perform split
1227             n_out = 0, n_in = 0;
1228             dsp::split_triangle_raw(out, &n_out, in, &n_in, &spl, &src);
1229 
1230             // Compute the area of triangle and ensure that it is greater than previous value
1231             float in_area       = 0.0f;
1232             for (size_t i=0; i<n_in; ++i)
1233                 in_area          += dsp::calc_area_pv(in[i].v);
1234             if (in_area > prev_area)
1235             {
1236                 // Compute the amplitude
1237                 float amplitude = sqrtf(in_area - prev_area) * afactor;
1238                 prev_area       = in_area;
1239 
1240                 RT_TRACE_BREAK(trace->pDebug,
1241                     lsp_trace("After split in=%d (GREEN), out=%d (RED) sample=%d, amplitude=%e",
1242                             int(n_in), int(n_out),
1243                             int(csn-1), amplitude
1244                         );
1245                     view->add_view_1c(v, &C3D_MAGENTA);
1246                     view->add_plane_pvn1c(p, &spl, &C3D_YELLOW);
1247                     for (size_t i=0; i<n_in; ++i)
1248                         view->add_triangle_pv1c(in[i].v, &C3D_GREEN);
1249                     for (size_t i=0; i<n_out; ++i)
1250                         view->add_triangle_pv1c(out[i].v, &C3D_RED);
1251                 )
1252 
1253                 // Deploy energy value to the sample
1254                 if (csn > 0)
1255                 {
1256                     // Lock capture data
1257 //                    trace->lkCapture.lock();
1258 
1259                     // Append sample to each matching capture
1260                     for (size_t ci=0, cn=bindings->size(); ci<cn; ++ci)
1261                     {
1262                         sample_t *s = bindings->at(ci);
1263 
1264                         // Skip reflection not in range
1265                         if ((s->r_min >= 0) && (v->rnum < s->r_min))
1266                             continue;
1267                         else if ((s->r_max >= 0) && (v->rnum > s->r_max))
1268                             continue;
1269 
1270                         // Ensure that we need to resize sample
1271                         size_t len = s->sample->length();
1272                         if (len <= size_t(csn))
1273                         {
1274                             // Need to resize sample?
1275                             if (s->sample->max_length() <= size_t(csn))
1276                             {
1277                                 len     = (csn + 1 + SAMPLE_QUANTITY) / SAMPLE_QUANTITY;
1278                                 len    *= SAMPLE_QUANTITY;
1279 
1280                                 lsp_trace("v->time = {%e, %e, %e}", v->time[0], v->time[1], v->time[2]);
1281                                 lsp_trace("ctime = %e, tsn = {%e, %e, %e}", ctime, tsn[0], tsn[1], tsn[2]);
1282                                 lsp_trace("spl = {%e, %e, %e, %e}",
1283                                     spl.dx, spl.dy, spl.dz, spl.dw);
1284                                 lsp_trace("Requesting sample resize: csn=0x%llx, len=0x%llx, channels=%d",
1285                                     (long long)csn, (long long)len, int(s->sample->channels())
1286                                     );
1287                                 #ifdef LSP_TRACE
1288                                     if (len > 0x100000) // TODO: This is currently impossible, added for debugging, remove in future
1289                                         invalid_state_hook();
1290                                 #endif
1291                                 if (!s->sample->resize(s->sample->channels(), len, len))
1292                                 {
1293                                     #ifdef LSP_TRACE
1294                                         invalid_state_hook();
1295                                     #endif
1296                                     return STATUS_NO_MEM;
1297                                 }
1298                             }
1299 
1300                             // Update sample length
1301                             s->sample->setLength(csn + 1);
1302                         }
1303 
1304                         // Deploy sample to curent channel
1305                         float *buf  = s->sample->getBuffer(s->channel);
1306                         buf[csn-1] += amplitude;
1307                     }
1308 
1309                     // Unlock capture data
1310 //                    trace->lkCapture.unlock();
1311                 }
1312             }
1313 
1314             ++csn; // Increment culling sample number for next iteration
1315         } while (n_out > 0);
1316 
1317 //        lsp_trace("Samples %d-%d -> area=%e amplitude=%e kcos=%f rnum=%d",
1318 //                int(ssn), int(csn-1), v_area, v->amplitude, kcos, int(v->rnum));
1319 
1320         return STATUS_OK;
1321     }
1322 
prepare_main_loop(float initial)1323     status_t RayTrace3D::TaskThread::prepare_main_loop(float initial)
1324     {
1325         // Cleanup stats
1326         clear_stats(&stats);
1327 
1328         // Report progress as 0%
1329         status_t res    = trace->report_progress(0.0f);
1330         if (res != STATUS_OK)
1331             return res;
1332         else if (trace->bCancelled)
1333             return STATUS_CANCELLED;
1334 
1335         RT_TRACE_BREAK(trace->pDebug,
1336             for (size_t i=0,n=trace->pScene->num_objects(); i<n; ++i)
1337             {
1338                 Object3D *obj = trace->pScene->get_object(i);
1339                 if (!obj->is_visible())
1340                     continue;
1341                 for (size_t j=0,m=obj->num_triangles(); j<m; ++j)
1342                     trace->pDebug->trace.add_triangle(obj->triangle(j), &C3D_RED, &C3D_GREEN, &C3D_BLUE);
1343             }
1344         );
1345 
1346         // Prepare root context and captures
1347         res         = generate_root_mesh();
1348         if (res == STATUS_OK)
1349             res         = prepare_captures();
1350 
1351         if (res != STATUS_OK)
1352             return res;
1353         else if (trace->bCancelled)
1354             return STATUS_CANCELLED;
1355 
1356         // Generate ray-tracing tasks
1357         rt_context_t *ctx = NULL;
1358         cvector<rt_context_t> estimate;
1359         res         = generate_tasks(&estimate, initial);
1360         if (res != STATUS_OK)
1361         {
1362             destroy_tasks(&estimate);
1363             return res;
1364         }
1365         else if (trace->bCancelled)
1366         {
1367             destroy_tasks(&estimate);
1368             return STATUS_CANCELLED;
1369         }
1370 
1371         // Estimate the progress by doing set of steps
1372         heavy_state = -1; // This guarantees that all tasks will be submitted to local task queue
1373         do
1374         {
1375             while (estimate.size() > 0)
1376             {
1377                 // Check cancellation flag
1378                 if (trace->bCancelled)
1379                 {
1380                     destroy_tasks(&tasks);
1381                     destroy_tasks(&estimate);
1382                     return STATUS_CANCELLED;
1383                 }
1384 
1385                 // Get new task
1386                 if (!estimate.pop(&ctx))
1387                 {
1388                     destroy_tasks(&tasks);
1389                     destroy_tasks(&estimate);
1390                     return STATUS_CORRUPTED;
1391                 }
1392 
1393                 // Process new task but store into
1394                 ++stats.root_tasks;
1395                 res     = process_context(ctx);
1396                 if (res != STATUS_OK)
1397                 {
1398                     destroy_tasks(&tasks);
1399                     destroy_tasks(&estimate);
1400                     return res;
1401                 }
1402             }
1403 
1404             // Perform swap: empty local tasks, fill 'estimate' with them
1405             estimate.swap_data(&tasks);
1406         } while ((estimate.size() > 0) && (estimate.size() < TASK_LO_THRESH));
1407 
1408         heavy_state         = S_SCAN_OBJECTS; // Enable global task queue for this thread
1409         trace->vTasks.swap_data(&estimate); // Now all generated tasks are global
1410 
1411         // Values to report progress
1412         trace->nProgressPoints  = 1;
1413         trace->nQueueSize       = trace->vTasks.size();
1414         trace->nProgressMax     = trace->nQueueSize + 2;
1415 
1416         // Report progress
1417         res         = trace->report_progress(float(trace->nProgressPoints++) / float(trace->nProgressMax));
1418         if (res != STATUS_OK)
1419         {
1420             destroy_tasks(&trace->vTasks);
1421             return res;
1422         }
1423         else if (trace->bCancelled)
1424         {
1425             destroy_tasks(&trace->vTasks);
1426             return STATUS_CANCELLED;
1427         }
1428 
1429         return STATUS_OK;
1430     }
1431 
prepare_captures()1432     status_t RayTrace3D::TaskThread::prepare_captures()
1433     {
1434         // Copy capture state
1435         for (size_t i=0; i<trace->vCaptures.size(); ++i)
1436         {
1437             // Allocate capture
1438             capture_t *scap = trace->vCaptures.get(i);
1439             rt_binding_t *b = new rt_binding_t();
1440             if (b == NULL)
1441                 return STATUS_NO_MEM;
1442             else if (!bindings.add(b))
1443             {
1444                 delete b;
1445                 return STATUS_NO_MEM;
1446             }
1447 
1448             // Copy bindings
1449             for (size_t j=0; j<scap->bindings.size(); ++j)
1450             {
1451                 // Allocate binding
1452                 sample_t *ssamp = scap->bindings.get(j);
1453                 sample_t *dsamp = b->bindings.add();
1454                 if (dsamp == NULL)
1455                     return STATUS_NO_MEM;
1456 
1457                 dsamp->sample   = NULL;
1458                 dsamp->channel  = ssamp->channel;
1459                 dsamp->r_min    = ssamp->r_min;
1460                 dsamp->r_max    = ssamp->r_max;
1461 
1462                 // Allocate sample and link
1463                 Sample *xsamp   = ssamp->sample;
1464                 Sample *tsamp   = new Sample();
1465                 if (tsamp == NULL)
1466                     return STATUS_NO_MEM;
1467                 else if (!tsamp->init(xsamp->channels(), xsamp->max_length(), xsamp->length()))
1468                 {
1469                     tsamp->destroy();
1470                     delete tsamp;
1471                     return STATUS_NO_MEM;
1472                 }
1473                 dsamp->sample   = tsamp;
1474             }
1475         }
1476 
1477         return STATUS_OK;
1478     }
1479 
prepare_supplementary_loop(TaskThread * t)1480     status_t RayTrace3D::TaskThread::prepare_supplementary_loop(TaskThread *t)
1481     {
1482         // Cleanup statistics
1483         clear_stats(&stats);
1484 
1485         // Prepare captures and data context
1486         status_t res = prepare_captures();
1487         if (res == STATUS_OK)
1488             res = copy_objects(&t->objects);
1489 
1490         return res;
1491     }
1492 
copy_objects(cvector<rt_object_t> * src)1493     status_t RayTrace3D::TaskThread::copy_objects(cvector<rt_object_t> *src)
1494     {
1495         rtx_edge_t *se, *de;
1496         rtx_triangle_t *dt;
1497 
1498         for (size_t i=0, n=src->size(); i<n; ++i)
1499         {
1500             // Get source object
1501             rt_object_t *s = src->at(i);
1502             if (s == NULL)
1503                 return STATUS_CORRUPTED;
1504 
1505             // Create destination object
1506             rt_object_t *d = new rt_object_t;
1507             if (d == NULL)
1508                 return STATUS_NO_MEM;
1509             else if (!objects.add(d))
1510             {
1511                 delete d;
1512                 return STATUS_NO_MEM;
1513             }
1514 
1515             // Copy edges and triangles
1516             if (!d->plan.add_all(&s->plan))
1517                 return STATUS_NO_MEM;
1518             if (!d->mesh.add_all(&s->mesh))
1519                 return STATUS_NO_MEM;
1520 
1521             // Patch pointers
1522             se = s->plan.get_array();
1523             de = d->plan.get_array();
1524             dt = d->mesh.get_array();
1525             for (size_t j=0, m=d->mesh.size(); j<m; ++j, ++dt)
1526             {
1527                 dt->e[0] = &de[dt->e[0] - se];
1528                 dt->e[1] = &de[dt->e[1] - se];
1529                 dt->e[2] = &de[dt->e[2] - se];
1530             }
1531 
1532             // Copy bound box
1533             d->bbox     = s->bbox;
1534         }
1535 
1536         return STATUS_OK;
1537     }
1538 
merge_result()1539     status_t RayTrace3D::TaskThread::merge_result()
1540     {
1541         cvector<capture_t> &dst = trace->vCaptures;
1542         if (dst.size() != bindings.size())
1543             return STATUS_CORRUPTED;
1544 
1545         for (size_t i=0; i<dst.size(); ++i)
1546         {
1547             rt_binding_t   *csrc    = bindings.at(i);
1548             capture_t      *cdst    = dst.at(i);
1549 
1550             if (csrc->bindings.size() != cdst->bindings.size())
1551                 return STATUS_CORRUPTED;
1552 
1553             for (size_t j=0; j<csrc->bindings.size(); ++j)
1554             {
1555                 sample_t *ssrc  = csrc->bindings.at(j);
1556                 sample_t *sdst  = cdst->bindings.at(j);
1557 
1558                 if ((ssrc->sample == NULL) || (sdst->sample == NULL))
1559                     return STATUS_CORRUPTED;
1560 
1561                 // Compute new sample size
1562                 size_t nc = ssrc->sample->channels();
1563                 if (nc != sdst->sample->channels())
1564                     return STATUS_CORRUPTED;
1565 
1566                 bool resize = false;
1567                 size_t maxlen = sdst->sample->max_length();
1568                 if (maxlen < ssrc->sample->max_length())
1569                 {
1570                     maxlen  = ssrc->sample->max_length();
1571                     resize  = true;
1572                 }
1573 
1574                 size_t len = sdst->sample->length();
1575                 if (len < ssrc->sample->length())
1576                 {
1577                     len     = ssrc->sample->length();
1578                     resize  = true;
1579                 }
1580 
1581                 if (maxlen < len)
1582                     maxlen  = len;
1583 
1584                 // Check that we need resize
1585                 if (resize)
1586                 {
1587                     if (!sdst->sample->resize(nc, maxlen, len))
1588                         return STATUS_NO_MEM;
1589                 }
1590 
1591                 // Apply changes to the target sample
1592                 for (size_t k=0; k<nc; ++k)
1593                 {
1594 //                    lsp_trace("Merge %p (%lld samples) -> %p (%lld samples), channel=%d/%d",
1595 //                            ssrc->sample->getBuffer(k),
1596 //                            (long long)(ssrc->sample->length()),
1597 //                            sdst->sample->getBuffer(k),
1598 //                            (long long)(sdst->sample->length()),
1599 //                            int(k), int(nc)
1600 //                            );
1601                     dsp::add2(
1602                             sdst->sample->getBuffer(k),
1603                             ssrc->sample->getBuffer(k),
1604                             ssrc->sample->length()
1605                         );
1606                 }
1607             }
1608         }
1609 
1610         return STATUS_OK;
1611     }
1612 
RayTrace3D()1613     RayTrace3D::RayTrace3D()
1614     {
1615         pScene          = NULL;
1616         pProgress       = NULL;
1617         pProgressData   = NULL;
1618         nSampleRate     = DEFAULT_SAMPLE_RATE;
1619         pDebug          = NULL;
1620         fEnergyThresh   = 1e-6f;
1621         fTolerance      = 1e-5f;
1622         fDetalization   = 1e-10f;
1623         bNormalize      = true;
1624         bCancelled      = false;
1625         nQueueSize      = 0;
1626         nProgressPoints = 0;
1627         nProgressMax    = 0;
1628     }
1629 
~RayTrace3D()1630     RayTrace3D::~RayTrace3D()
1631     {
1632         destroy(true);
1633     }
1634 
clear_stats(stats_t * stats)1635     void RayTrace3D::clear_stats(stats_t *stats)
1636     {
1637         stats->root_tasks       = 0;
1638         stats->local_tasks      = 0;
1639         stats->calls_scan       = 0;
1640         stats->calls_cull       = 0;
1641         stats->calls_split      = 0;
1642         stats->calls_cullback   = 0;
1643         stats->calls_reflect    = 0;
1644         stats->calls_capture    = 0;
1645     }
1646 
dump_stats(const char * label,const stats_t * stats)1647     void RayTrace3D::dump_stats(const char *label, const stats_t *stats)
1648     {
1649         lsp_trace("%s:\n"
1650                 "  root tasks processed     : %lld\n"
1651                 "  local tasks processed    : %lld\n"
1652                 "  scan_objects             : %lld\n"
1653                 "  cull_view                : %lld\n"
1654                 "  split_view               : %lld\n"
1655                 "  cullback_view            : %lld\n"
1656                 "  reflect_view             : %lld\n"
1657                 "  capture                  : %lld\n",
1658             label,
1659             (long long)stats->root_tasks,
1660             (long long)stats->local_tasks,
1661             (long long)stats->calls_scan,
1662             (long long)stats->calls_cull,
1663             (long long)stats->calls_split,
1664             (long long)stats->calls_cullback,
1665             (long long)stats->calls_reflect,
1666             (long long)stats->calls_capture
1667         );
1668     }
1669 
merge_stats(stats_t * dst,const stats_t * src)1670     void RayTrace3D::merge_stats(stats_t *dst, const stats_t *src)
1671     {
1672         dst->root_tasks        += src->root_tasks;
1673         dst->local_tasks       += src->local_tasks;
1674         dst->calls_scan        += src->calls_scan;
1675         dst->calls_cull        += src->calls_cull;
1676         dst->calls_split       += src->calls_split;
1677         dst->calls_cullback    += src->calls_cullback;
1678         dst->calls_reflect     += src->calls_reflect;
1679         dst->calls_capture     += src->calls_capture;
1680     }
1681 
destroy_tasks(cvector<rt_context_t> * tasks)1682     void RayTrace3D::destroy_tasks(cvector<rt_context_t> *tasks)
1683     {
1684         for (size_t i=0, n=tasks->size(); i<n; ++i)
1685         {
1686             rt_context_t *ctx   = tasks->get(i);
1687             if (ctx != NULL)
1688                 delete ctx;
1689         }
1690 
1691         tasks->flush();
1692     }
1693 
destroy_objects(cvector<rt_object_t> * objects)1694     void RayTrace3D::destroy_objects(cvector<rt_object_t> *objects)
1695     {
1696         for (size_t i=0, n=objects->size(); i<n; ++i)
1697         {
1698             rt_object_t *obj = objects->get(i);
1699             if (obj != NULL)
1700             {
1701                 obj->mesh.flush();
1702                 obj->plan.flush();
1703                 delete obj;
1704             }
1705         }
1706 
1707         objects->flush();
1708     }
1709 
remove_scene(bool destroy)1710     void RayTrace3D::remove_scene(bool destroy)
1711     {
1712         if (pScene != NULL)
1713         {
1714             if (destroy)
1715             {
1716                 pScene->destroy();
1717                 delete pScene;
1718             }
1719             pScene = NULL;
1720         }
1721     }
1722 
resize_materials(size_t objects)1723     status_t RayTrace3D::resize_materials(size_t objects)
1724     {
1725         size_t size = vMaterials.size();
1726 
1727         if (objects < size)
1728         {
1729             if (!vMaterials.remove_n(objects, size - objects))
1730                 return STATUS_UNKNOWN_ERR;
1731         }
1732         else if (objects > size)
1733         {
1734             if (!vMaterials.append_n(objects - size))
1735                 return STATUS_NO_MEM;
1736 
1737             while (size < objects)
1738             {
1739                 rt_material_t *m    = vMaterials.get(size++);
1740                 if (m == NULL)
1741                     return STATUS_UNKNOWN_ERR;
1742 
1743                 // By default, we set the material to 'Concrete'
1744                 m->absorption[0]    = 0.02f;
1745                 m->diffusion[0]     = 1.0f;
1746                 m->dispersion[0]    = 1.0f;
1747                 m->transparency[0]  = 0.48f;
1748 
1749                 m->absorption[1]    = 0.0f;
1750                 m->diffusion[1]     = 1.0f;
1751                 m->dispersion[1]    = 1.0f;
1752                 m->transparency[1]  = 0.52f;
1753 
1754                 m->permeability     = 12.88f;
1755             }
1756         }
1757 
1758         return STATUS_OK;
1759     }
1760 
init()1761     status_t RayTrace3D::init()
1762     {
1763         return STATUS_OK;
1764     }
1765 
destroy(bool recursive)1766     void RayTrace3D::destroy(bool recursive)
1767     {
1768         destroy_tasks(&vTasks);
1769         clear_progress_callback();
1770         remove_scene(recursive);
1771 
1772         for (size_t i=0, n=vCaptures.size(); i<n; ++i)
1773         {
1774             capture_t *cap = vCaptures.get(i);
1775             if (cap != NULL)
1776             {
1777                 cap->bindings.flush();
1778                 delete cap;
1779             }
1780         }
1781         vCaptures.flush();
1782 
1783         vMaterials.flush();
1784         vSources.flush();
1785         vCaptures.flush();
1786     }
1787 
add_source(const rt_source_settings_t * settings)1788     status_t RayTrace3D::add_source(const rt_source_settings_t *settings)
1789     {
1790         if (settings == NULL)
1791             return STATUS_BAD_ARGUMENTS;
1792 
1793         rt_source_settings_t *src = vSources.add();
1794         if (src == NULL)
1795             return STATUS_NO_MEM;
1796 
1797         *src        = *settings;
1798 
1799         return STATUS_OK;
1800     }
1801 
add_capture(const rt_capture_settings_t * settings)1802     ssize_t RayTrace3D::add_capture(const rt_capture_settings_t *settings)
1803     {
1804         if (settings == NULL)
1805             return STATUS_BAD_ARGUMENTS;
1806 
1807         capture_t *cap      = new capture_t();
1808         if (cap == NULL)
1809             return -STATUS_NO_MEM;
1810 
1811         size_t idx          = vCaptures.size();
1812         if (!vCaptures.add(cap))
1813         {
1814             delete cap;
1815             return -STATUS_NO_MEM;
1816         }
1817 
1818         cap->pos            = settings->pos;
1819         dsp::init_vector_dxyz(&cap->direction, 1.0f, 0.0f, 0.0f);
1820         cap->radius         = settings->radius;
1821         cap->type           = settings->type;
1822 
1823         dsp::apply_matrix3d_mv1(&cap->direction, &cap->pos);
1824         dsp::normalize_vector(&cap->direction);
1825 
1826         return idx;
1827     }
1828 
bind_capture(size_t id,Sample * sample,size_t channel,ssize_t r_min,ssize_t r_max)1829     status_t RayTrace3D::bind_capture(size_t id, Sample *sample, size_t channel, ssize_t r_min, ssize_t r_max)
1830     {
1831         capture_t *cap = vCaptures.get(id);
1832         if (cap == NULL)
1833             return STATUS_INVALID_VALUE;
1834 
1835         sample_t *s     = cap->bindings.add();
1836         if (s == NULL)
1837             return STATUS_NO_MEM;
1838 
1839         s->sample       = sample;
1840         s->channel      = channel;
1841         s->r_min        = r_min;
1842         s->r_max        = r_max;
1843 
1844         return STATUS_OK;
1845     }
1846 
set_scene(Scene3D * scene,bool destroy)1847     status_t RayTrace3D::set_scene(Scene3D *scene, bool destroy)
1848     {
1849         status_t res = resize_materials(scene->num_objects());
1850         if (res != STATUS_OK)
1851             return res;
1852 
1853         // Destroy scene
1854         remove_scene(destroy);
1855         pScene      = scene;
1856         return STATUS_OK;
1857     }
1858 
set_material(size_t idx,const rt_material_t * material)1859     status_t RayTrace3D::set_material(size_t idx, const rt_material_t *material)
1860     {
1861         rt_material_t *m = vMaterials.get(idx);
1862         if (m == NULL)
1863             return STATUS_INVALID_VALUE;
1864         *m = *material;
1865         return STATUS_OK;
1866     }
1867 
get_material(rt_material_t * material,size_t idx)1868     status_t RayTrace3D::get_material(rt_material_t *material, size_t idx)
1869     {
1870         if (material == NULL)
1871             return STATUS_BAD_ARGUMENTS;
1872         rt_material_t *m = vMaterials.get(idx);
1873         if (m == NULL)
1874             return STATUS_INVALID_VALUE;
1875 
1876         *material = *m;
1877         return STATUS_OK;
1878     }
1879 
set_progress_callback(rt_progress_t callback,void * data)1880     status_t RayTrace3D::set_progress_callback(rt_progress_t callback, void *data)
1881     {
1882         if (callback == NULL)
1883             return clear_progress_callback();
1884 
1885         pProgress       = callback;
1886         pProgressData   = data;
1887         return STATUS_OK;
1888     }
1889 
clear_progress_callback()1890     status_t RayTrace3D::clear_progress_callback()
1891     {
1892         pProgress       = NULL;
1893         pProgressData   = NULL;
1894         return STATUS_OK;
1895     }
1896 
report_progress(float progress)1897     status_t RayTrace3D::report_progress(float progress)
1898     {
1899         if (pProgress == NULL)
1900             return STATUS_OK;
1901         return pProgress(progress, pProgressData);
1902     }
1903 
do_process(size_t threads,float initial)1904     status_t RayTrace3D::do_process(size_t threads, float initial)
1905     {
1906         status_t res = STATUS_OK;
1907         bCancelled   = false;
1908         bFailed      = false;
1909 
1910         // Get time of execution start
1911 #ifdef LSP_TRACE
1912         struct timespec tstart;
1913         clock_gettime(CLOCK_REALTIME, &tstart);
1914 #endif
1915 
1916         // Create main thread
1917         TaskThread *root = new TaskThread(this);
1918         if (root == NULL)
1919             return STATUS_NO_MEM;
1920 
1921         // Launch prepare_main_loop in root thread's context
1922         res    = root->prepare_main_loop(initial);
1923         if (res != STATUS_OK)
1924         {
1925             delete root;
1926             return res;
1927         }
1928 
1929         // Launch supplementary threads
1930         cvector<TaskThread> workers;
1931         if (vTasks.size() > 0)
1932         {
1933             for (size_t i=1; i<threads; ++i)
1934             {
1935                 // Create thread object
1936                 TaskThread *t   = new TaskThread(this);
1937                 if ((t == NULL) || (!workers.add(t)))
1938                 {
1939                     if (t != NULL)
1940                         delete t;
1941                     res = STATUS_NO_MEM;
1942                     break;
1943                 }
1944 
1945                 // Sync thread data
1946                 res = t->prepare_supplementary_loop(root);
1947                 if (res != STATUS_OK)
1948                     break;
1949 
1950                 // Launch thread
1951                 res = t->start();
1952                 if (res != STATUS_OK)
1953                     break;
1954             }
1955         }
1956 
1957         // If successful status, perform main loop
1958         if (res == STATUS_OK)
1959             res     = root->run();
1960         else
1961             bFailed = true;
1962 
1963         // Wait for supplementary threads
1964         for (size_t i=0,n=workers.size(); i<n; ++i)
1965         {
1966             // Wait for thread completion
1967             TaskThread *t = workers.get(i);
1968             t->join();
1969             if (res == STATUS_OK)
1970                 res     = t->get_result(); // Update execution status
1971         }
1972 
1973         // Get root thread statistics
1974         stats_t overall;
1975         clear_stats(&overall);
1976         merge_stats(&overall, root->get_stats());
1977         root->merge_result();
1978         if (res != STATUS_BREAK_POINT)
1979             dump_stats("Main thread statistics", root->get_stats());
1980 
1981         // Output thread stats and destroy threads
1982         for (size_t i=0,n=workers.size(); i<n; ++i)
1983         {
1984             // Post-process each thread
1985             TaskThread *t = workers.get(i);
1986             t->merge_result();
1987 
1988             // Merge and output statistics
1989             LSPString s;
1990             s.fmt_utf8("Supplementary thread %d statistics", int(i));
1991             merge_stats(&overall, t->get_stats());
1992             if (res != STATUS_BREAK_POINT)
1993                 dump_stats(s.get_utf8(), t->get_stats());
1994 
1995             // Detroy thread object
1996             delete t;
1997         }
1998         delete root;
1999         workers.flush();
2000 
2001         // Dump overall statistics
2002         if (res != STATUS_BREAK_POINT)
2003         {
2004             // Get time of execution end
2005 #ifdef LSP_TRACE
2006             struct timespec tend;
2007             clock_gettime(CLOCK_REALTIME, &tend);
2008             double etime = double(tend.tv_sec - tstart.tv_sec) + double(tend.tv_nsec - tend.tv_nsec) * 1e-6;
2009 #endif
2010 
2011             dump_stats("Overall statistics", &overall);
2012             lsp_trace("Overall execution time:      %f s", etime);
2013         }
2014 
2015         // Destroy all tasks
2016         destroy_tasks(&vTasks);
2017         if (res != STATUS_OK)
2018             return res;
2019 
2020         // Normalize output
2021         if (bNormalize)
2022             normalize_output();
2023 
2024         float prg   = float(nProgressPoints) / float(nProgressMax);
2025         lsp_trace("Reporting progress %d/%d = %.2f%%", int(nProgressPoints), int(nProgressMax), prg * 100.0f);
2026         ++nProgressPoints;
2027 
2028         return report_progress(prg);
2029     }
2030 
process(size_t threads,float initial)2031     status_t RayTrace3D::process(size_t threads, float initial)
2032     {
2033         // We need to initialize DSP context for processing
2034         dsp::context_t ctx;
2035         dsp::start(&ctx);
2036 
2037         // Call processing routine
2038         status_t res = do_process(threads, initial);
2039 
2040         // Finalize DSP context and exit
2041         dsp::finish(&ctx);
2042         return res;
2043     }
2044 
is_already_passed(const sample_t * bind)2045     bool RayTrace3D::is_already_passed(const sample_t *bind)
2046     {
2047         for (size_t i=0; i<vCaptures.size(); ++i)
2048         {
2049             capture_t *cap = vCaptures.at(i);
2050             for (size_t j=0; j<cap->bindings.size(); ++j)
2051             {
2052                 sample_t *s = cap->bindings.at(j);
2053                 // We reached current position?
2054                 if (s == bind)
2055                     return false;
2056                 // Check that the same sample and channel have been processed earlier than 'bind'
2057                 if ((s->sample == bind->sample) && (s->channel == bind->channel))
2058                     return true;
2059             }
2060         }
2061         return false;
2062     }
2063 
normalize_output()2064     void RayTrace3D::normalize_output()
2065     {
2066         float max_gain = 0.0f;
2067         cstorage<sample_t> plan;
2068 
2069         // Estimate the maximum output gain for each capture's binding
2070         for (size_t i=0; i<vCaptures.size(); ++i)
2071         {
2072             capture_t *cap = vCaptures.at(i);
2073             for (size_t j=0; j<cap->bindings.size(); ++j)
2074             {
2075                 sample_t *s = cap->bindings.at(j);
2076                 if (is_already_passed(s)) // Protect from measuring gain more than once
2077                     continue;
2078 
2079                 // Estimate the maximum gain
2080                 float c_gain = dsp::abs_max(s->sample->getBuffer(s->channel), s->sample->length());
2081                 if (max_gain < c_gain)
2082                     max_gain = c_gain;
2083             }
2084         }
2085 
2086         // Now we know the maximum gain
2087         if (max_gain == 0.0f)
2088             return;
2089         max_gain = 1.0f / max_gain; // Now it's a norming factor
2090 
2091         // Perform the gain adjustment
2092         for (size_t i=0; i<vCaptures.size(); ++i)
2093         {
2094             capture_t *cap = vCaptures.at(i);
2095             for (size_t j=0; j<cap->bindings.size(); ++j)
2096             {
2097                 sample_t *s = cap->bindings.at(j);
2098                 if (is_already_passed(s)) // Protect from adjusting gain more than once
2099                     continue;
2100 
2101                 // Apply the norming factor
2102                 dsp::mul_k2(s->sample->getBuffer(s->channel), max_gain, s->sample->length());
2103             }
2104         }
2105     }
2106 
check_bound_box(const bound_box3d_t * bbox,const rt_view_t * view)2107     bool RayTrace3D::check_bound_box(const bound_box3d_t *bbox, const rt_view_t *view)
2108     {
2109         const vector3d_t *pl;
2110         raw_triangle_t buf1[16], buf2[16], *in, *out;
2111         size_t nin, nout;
2112 
2113         // Cull each triangle of bounding box with four scissor planes
2114         for (size_t i=0, m = sizeof(bbox_map)/sizeof(size_t); i < m; )
2115         {
2116             // Initialize input
2117             in          = buf1;
2118             out         = buf2;
2119             nin         = 1;
2120 
2121             in->v[0]    = bbox->p[bbox_map[i++]];
2122             in->v[1]    = bbox->p[bbox_map[i++]];
2123             in->v[2]    = bbox->p[bbox_map[i++]];
2124             pl          = view->pl;
2125 
2126             // Cull triangle with planes
2127             for (size_t j=0; j<4; ++j, ++pl)
2128             {
2129                 // Reset counters
2130                 nout    = 0;
2131                 for (size_t k=0; k < nin; ++k, ++in)
2132                     dsp::cull_triangle_raw(out, &nout, pl, in);
2133 
2134                 // Interrupt cycle if there is no data to process
2135                 if (!nout)
2136                     break;
2137 
2138                 // Update state
2139                 nin     = nout;
2140                 if (j & 1)
2141                     in = buf1, out = buf2;
2142                 else
2143                     in = buf2, out = buf1;
2144             }
2145 
2146             if (nout)
2147                 break;
2148         }
2149 
2150         return nout;
2151     }
2152 
2153 
2154 } /* namespace lsp */
2155