BRL-CAD
photonmap.c
Go to the documentation of this file.
1 /* P H O T O N M A P . C
2  * BRL-CAD
3  *
4  * Copyright (c) 2002-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 /** @file liboptical/photonmap.c
21  *
22  * Implementation of Photon Mapping
23  *
24  */
25 
26 #include "common.h"
27 
28 #include <limits.h>
29 #include <stdlib.h>
30 #include <time.h>
31 
32 #ifdef HAVE_UNISTD_H
33 # include <unistd.h>
34 #endif
35 #ifdef HAVE_ALARM
36 # include <signal.h>
37 #endif
38 
39 #include "bu/parallel.h"
40 #include "photonmap.h"
41 
42 
45 
46 struct PhotonMap *PMap[PM_MAPS];/* Photon Map (KD-TREE) */
47 struct Photon *Emit[PM_MAPS]; /* Emitted Photons */
48 struct Photon CurPh;
49 vect_t BBMin; /* Min Bounding Box */
50 vect_t BBMax; /* Max Bounding Box */
51 int Depth; /* Used to determine how many times the photon has propagated */
52 int PType; /* Used to determine the type of Photon: Direct, Indirect, Specular, Caustic */
53 int PInit;
54 int EPL; /* Emitted Photons For the Light */
55 int EPS[PM_MAPS]; /* Emitted Photons For the Light */
56 int ICSize;
57 double ScaleFactor;
58 struct IrradCache *IC; /* Irradiance Cache for Hypersampling */
59 char *Map; /* Used for Irradiance HyperSampling Cache */
60 int GPM_IH; /* Irradiance Hypersampling Toggle, 0=off, 1=on */
63 int GPM_RAYS; /* Number of Sample Rays for each Direction in Irradiance Hemi */
64 double GPM_ATOL; /* Angular Tolerance for Photon Gathering */
65 struct resource GPM_RTAB[MAX_PSW]; /* Resource Table for Multi-threading */
66 int HitG, HitB;
67 
68 
69 /* Split so that equal numbers are above and below the splitting plane */
70 int
71 FindMedian(struct Photon *List, int Num, int Axis)
72 {
73  int i;
74  fastf_t Min, Max, Mean;
75 
76  Min = Max = List[0].Pos[Axis];
77  for (i = 1; i < Num; i++) {
78  if (List[i].Pos[Axis] < Min)
79  Min = List[i].Pos[Axis];
80  if (List[i].Pos[Axis] > Max)
81  Max = List[i].Pos[Axis];
82  }
83  Mean = (Min+Max)/2.0;
84  i = 0;
85  while (List[i].Pos[Axis] < Mean && i < Num)
86  i++;
87 
88  return i;
89 }
90 
91 
92 /* Generate a KD-Tree from a Flat Array of Photons */
93 void
94 BuildTree(struct Photon *EList, int ESize, struct PNode *Root)
95 {
96  struct Photon *LList, *RList;
97  vect_t Min, Max;
98  int i, Axis, MedianIndex, LInd, RInd;
99 
100 
101  /* Allocate memory for left and right lists */
102  LList = (struct Photon*)bu_calloc(ESize, sizeof(struct Photon), "LList");
103  RList = (struct Photon*)bu_calloc(ESize, sizeof(struct Photon), "RList");
104 
105  /* Find the Bounding volume of the Current list of photons */
106  Min[0] = Max[0] = EList[0].Pos[0];
107  Min[1] = Max[1] = EList[0].Pos[1];
108  Min[2] = Max[2] = EList[0].Pos[2];
109  for (i = 1; i < ESize; i++) {
110  if (EList[i].Pos[0] < Min[0]) Min[0] = EList[i].Pos[0];
111  if (EList[i].Pos[1] < Min[1]) Min[1] = EList[i].Pos[1];
112  if (EList[i].Pos[2] < Min[2]) Min[2] = EList[i].Pos[2];
113  if (EList[i].Pos[0] > Max[0]) Max[0] = EList[i].Pos[0];
114  if (EList[i].Pos[1] > Max[1]) Max[1] = EList[i].Pos[1];
115  if (EList[i].Pos[2] > Max[2]) Max[2] = EList[i].Pos[2];
116  }
117 
118  /* Obtain splitting Axis, which is the largest dimension of the bounding volume */
119  Max[0] -= Min[0];
120  Max[1] -= Min[1];
121  Max[2] -= Min[2];
122  Axis = 0;
123  if (Max[1] > Max[0] && Max[1] > Max[2]) Axis = 1;
124  if (Max[2] > Max[0] && Max[2] > Max[1]) Axis = 2;
125 
126  /* Find Median Photon to split by. */
127  MedianIndex = FindMedian(EList, ESize, Axis);
128 
129  /* Build Left and Right Lists and make sure the Median Photon is not included in either list. */
130  LInd = RInd = 0;
131  for (i = 0; i < ESize; i++) {
132  if (i != MedianIndex) {
133  if (EList[i].Pos[Axis] < EList[MedianIndex].Pos[Axis]) {
134  LList[LInd++] = EList[i];
135  } else {
136  RList[RInd++] = EList[i];
137  }
138  /*
139  if (EList[i].Pos[Axis] < Median.Pos[Axis]) {
140  LList[LInd++] = EList[i];
141  } else {
142  RList[RInd++] = EList[i];
143  }
144  */
145  }
146  }
147 
148  /* Store the Median Photon into the KD-Tree. */
149  /* bu_log("insertKD: %.3f, %.3f, %.3f\n", EList[MedianIndex].Pos[0], EList[MedianIndex].Pos[1], EList[MedianIndex].Pos[2]);*/
150  Root->P = EList[MedianIndex];
151  Root->P.Axis = Axis;
152  Root->C = 0;
153 
154  /* With Left and Right if either contain any photons then repeat this process */
155  /* if (LInd) bu_log("Left Branch\n");*/
156  if (LInd) {
157  BU_ALLOC(Root->L, struct PNode);
158  Root->L->L = 0;
159  Root->L->R = 0;
160  BuildTree(LList, LInd, Root->L);
161  }
162  /* if (RInd) bu_log("Right Branch\n");*/
163  if (RInd) {
164  BU_ALLOC(Root->R, struct PNode);
165  Root->R->L = 0;
166  Root->R->R = 0;
167  BuildTree(RList, RInd, Root->R);
168  }
169 
170  bu_free(LList, "LList");
171  bu_free(RList, "RList");
172 }
173 
174 
175 void
176 LocatePhotons(struct PhotonSearch *Search, struct PNode *Root)
177 {
178  fastf_t Dist, TDist, angle, MDist;
179  int i, MaxInd, Axis;
180 
181  if (!Root)
182  return;
183 
184  Axis = Root->P.Axis;
185  Dist = Search->Pos[Axis] - Root->P.Pos[Axis];
186 
187  if (Dist < 0) {
188  /* Left of plane - search left subtree first */
189  LocatePhotons(Search, Root->L);
190  if (Dist*Dist < Search->RadSq)
191  LocatePhotons(Search, Root->R);
192  } else {
193  /* Right of plane - search right subtree first */
194  LocatePhotons(Search, Root->R);
195  if (Dist*Dist < Search->RadSq)
196  LocatePhotons(Search, Root->L);
197  }
198 
199  /* REPLACE, Find Distance between Root Photon and NP->Pos */
200  Dist = (Root->P.Pos[0] - Search->Pos[0])*(Root->P.Pos[0] - Search->Pos[0]) + (Root->P.Pos[1] - Search->Pos[1])*(Root->P.Pos[1] - Search->Pos[1]) + (Root->P.Pos[2] - Search->Pos[2])*(Root->P.Pos[2] - Search->Pos[2]);
201 
202  angle = VDOT(Search->Normal, Root->P.Normal);
203  if (Dist < Search->RadSq && angle > GPM_ATOL) {
204  /* Check that Result is within Radius and Angular Tolerance */
205  /* if (Dist < NP->RadSq) {*/
206  if (Search->Found < Search->Max) {
207  Search->List[Search->Found++].P = Root->P;
208  } else {
209  MDist = (Search->Pos[0] - Search->List[0].P.Pos[0])*(Search->Pos[0] - Search->List[0].P.Pos[0])+(Search->Pos[1] - Search->List[0].P.Pos[1])*(Search->Pos[1] - Search->List[0].P.Pos[1])+(Search->Pos[2] - Search->List[0].P.Pos[2])*(Search->Pos[2] - Search->List[0].P.Pos[2]);
210  MaxInd = 0;
211  for (i = 1; i < Search->Found; i++) {
212  TDist = (Search->Pos[0] - Search->List[i].P.Pos[0])*(Search->Pos[0] - Search->List[i].P.Pos[0])+(Search->Pos[1] - Search->List[i].P.Pos[1])*(Search->Pos[1] - Search->List[i].P.Pos[1])+(Search->Pos[2] - Search->List[i].P.Pos[2])*(Search->Pos[2] - Search->List[i].P.Pos[2]);
213  if (TDist > MDist) {
214  MDist = TDist;
215  MaxInd = i;
216  }
217  }
218 
219  if (Dist < MDist)
220  Search->List[MaxInd].P = Root->P;
221  }
222  }
223 }
224 
225 
226 /* Places photon into flat array that wwill form the final kd-tree. */
227 void
228 Store(point_t Pos, vect_t Dir, vect_t Normal, int map)
229 {
230  struct PhotonSearch Search;
231  int i;
232 
233  /* If Importance Mapping is enabled, Check to see if the Photon is in an area that is considered important, if not then disregard it */
234  if (map != PM_IMPORTANCE && PMap[PM_IMPORTANCE]->StoredPhotons) {
235  /* Do a KD-Tree lookup and if the photon is within a distance of sqrt(ScaleFactor) from the nearest importon then keep it, otherwise discard it */
236 
237  Search.RadSq = ScaleFactor;
238  Search.Found = 0;
239  Search.Max = 1;
240  Search.Pos[0] = Pos[0];
241  Search.Pos[1] = Pos[1];
242  Search.Pos[2] = Pos[2];
243  Search.Normal[0] = Normal[0];
244  Search.Normal[1] = Normal[1];
245  Search.Normal[2] = Normal[2];
246 
247  Search.List = (struct PSN*)bu_calloc(Search.Max, sizeof(struct PSN), "Search.List");
248  LocatePhotons(&Search, PMap[PM_IMPORTANCE]->Root);
249  bu_free(Search.List, "Search.List");
250 
251  if (!Search.Found) {
252  HitB++;
253  return;
254  }
255  }
256 
257 
258  if (PMap[map]->StoredPhotons < PMap[map]->MaxPhotons) {
259  HitG++;
260  for (i = 0; i < 3; i++) {
261  /* Store Position, Direction, and Power of Photon */
262  Emit[map][PMap[map]->StoredPhotons].Pos[i] = Pos[i];
263  Emit[map][PMap[map]->StoredPhotons].Dir[i] = Dir[i];
264  Emit[map][PMap[map]->StoredPhotons].Normal[i] = Normal[i];
265  Emit[map][PMap[map]->StoredPhotons].Power[i] = CurPh.Power[i];
266  }
267  PMap[map]->StoredPhotons++;
268  /*
269  if (map != PM_IMPORTANCE && PMap[PM_IMPORTANCE]->StoredPhotons < PMap[PM_IMPORTANCE]->MaxPhotons)
270  bu_log("Map2: %d, Size: %d\n", map, PMap[map]->StoredPhotons);
271  */
272  /*
273  if (map == PM_IMPORTANCE)
274  bu_log("Map: %d, Size: %d, [%.3f, %.3f, %.3f] [%.3f, %.3f, %.3f]\n", map, PMap[map]->StoredPhotons, Pos[0], Pos[1], Pos[2], CurPh.Power[0], CurPh.Power[1], CurPh.Power[2]);
275  */
276  }
277 
278  /*
279  bu_log("[%d][%d][%.3f, %.3f, %.3f]\n", map, PMap[map]->StoredPhotons, CurPh.Power[0], CurPh.Power[1], CurPh.Power[2]);
280  if (!(PMap[map]->StoredPhotons % 64))
281  bu_log("[%d][%d][%.3f, %.3f, %.3f]\n", map, PMap[map]->StoredPhotons, Pos[0], Pos[1], Pos[2]);
282  */
283 }
284 
285 
286 /* Compute a specular reflection */
287 void
288 SpecularReflect(vect_t normal, vect_t rdir)
289 {
290  vect_t d;
291  fastf_t dot;
292 
293  VSCALE(d, rdir, -1);
294  dot = VDOT(d, normal);
295 
296  if (dot < 0.0f) {
297  rdir[0] = rdir[1] = rdir[2] = 0;
298  } else {
299  VSCALE(rdir, normal, 2*dot);
300  VSUB2(rdir, rdir, d);
301  }
302 }
303 
304 
305 /* Compute a random reflected diffuse direction */
306 void
307 DiffuseReflect(vect_t normal, vect_t rdir)
308 {
309  /* Allow Photons to get a random direction at most 60 degrees to the normal */
310  do {
311 #ifndef HAVE_DRAND48
312  rdir[0] = 2.0*rand()/(double)RAND_MAX-1.0;
313  rdir[1] = 2.0*rand()/(double)RAND_MAX-1.0;
314  rdir[2] = 2.0*rand()/(double)RAND_MAX-1.0;
315 #else
316  rdir[0] = 2.0*drand48()-1.0;
317  rdir[1] = 2.0*drand48()-1.0;
318  rdir[2] = 2.0*drand48()-1.0;
319 #endif
320  VUNITIZE(rdir);
321  } while (VDOT(rdir, normal) < 0.5);
322 }
323 
324 
325 /* Compute refracted ray given Incident Ray, Normal, and 2 refraction indices */
326 int
327 Refract(vect_t I, vect_t N, fastf_t n1, fastf_t n2)
328 {
329  fastf_t n, c1, c2, radicand;
330  vect_t t, r;
331 
332  n = n1/n2;
333  c1 = -VDOT(I, N);
334  radicand = 1.0 - (n*n)*(1.0-c1*c1);
335  if (radicand < 0) {
336  /* Total Internal Reflection */
337  I[0] = I[1] = I[2] = 0;
338  return 0;
339  }
340  c2 = sqrt(radicand);
341 
342  VSCALE(r, I, n);
343  VSCALE(t, N, n*c1-c2);
344  VADD2(I, r, t);
345  return 1;
346 }
347 
348 
349 int
350 CheckMaterial(char *cmp, char *MS)
351 {
352  size_t i;
353 
354  if (MS) {
355  for (i = 0; i < strlen(cmp) && i < strlen(MS); i++)
356  if (MS[i] != cmp[i])
357  return 0;
358  return 1;
359  } else {
360  return 0;
361  }
362 }
363 
364 
365 /* This function parses the material string to obtain specular and refractive values */
366 void
367 GetMaterial(char *MS, vect_t spec, fastf_t *refi, fastf_t *transmit)
368 {
369  struct phong_specific *phong_sp;
370  struct bu_vls matparm = BU_VLS_INIT_ZERO;
371 
372  BU_GET(phong_sp, struct phong_specific);
373 
374  /* Initialize spec and refi */
375  spec[0] = spec[1] = spec[2] = *refi = *transmit = 0;
376  if (CheckMaterial("plastic", MS)) {
377  /* Checks that the first 7 chars match any of the characters found in plastic */
378  /* Plastic Shader */
379  phong_sp->magic = PL_MAGIC;
380  phong_sp->shine = 10;
381  phong_sp->wgt_specular = 0.7;
382  phong_sp->wgt_diffuse = 0.3;
383  phong_sp->transmit = 0.0;
384  phong_sp->reflect = 0.0;
385  phong_sp->refrac_index = 1.0;
386  phong_sp->extinction = 0.0;
387 
388  MS += 7;
389  bu_vls_printf(&matparm, "%s", MS);
390  if (bu_struct_parse(&matparm, phong_parse, (char *)phong_sp, NULL) < 0)
391  bu_log("Warning - bu_struct_parse failure (matparm, phone_parse, material = plastic) in GetMaterial!\n");
392  bu_vls_free(&matparm);
393 
394  /* bu_log("string: %s\n", MS);*/
395  /* bu_log("info[%d]: sh: %d, spec: %.3f, diff: %.3f, tran: %.3f, refl: %.3f, refr: %.3f\n", PMap->StoredPhotons, phong_sp->shine, phong_sp->wgt_specular, phong_sp->wgt_diffuse, phong_sp->transmit, phong_sp->reflect, phong_sp->refrac_index);*/
396  /* bu_log("spec[%d]: %.3f\n", PMap->StoredPhotons, phong_sp->wgt_specular);*/
397 
398  spec[0] = spec[1] = spec[2] = phong_sp->wgt_specular;
399  *refi = phong_sp->refrac_index;
400  spec[0] = spec[1] = spec[2] = 0.7;
401  /*
402  *refi = 1.0;
403  *transmit = 0.0;
404  */
405  /*
406  if (phong_sp->refrac_index != 1.0)
407  bu_log("refi: %.3f\n", phong_sp->refrac_index);
408  */
409  } else if (CheckMaterial("glass", MS)) {
410  /* Glass Shader */
411  phong_sp->magic = PL_MAGIC;
412  phong_sp->shine = 4;
413  phong_sp->wgt_specular = 0.7;
414  phong_sp->wgt_diffuse = 0.3;
415  phong_sp->transmit = 0.8;
416  phong_sp->reflect = 0.1;
417  phong_sp->refrac_index = 1.65;
418  phong_sp->extinction = 0.0;
419 
420  MS += 5; /* move pointer past "pm " (3 characters) */
421  bu_vls_printf(&matparm, "%s", MS);
422  if (bu_struct_parse(&matparm, phong_parse, (char *)phong_sp, NULL) < 0)
423  bu_log("Warning - bu_struct_parse failure (matparm, phone_parse, material = glass) in GetMaterial!\n");
424  bu_vls_free(&matparm);
425 
426  /* bu_log("string: %s\n", MS);*/
427  /* bu_log("info[%d]: sh: %d, spec: %.3f, diff: %.3f, tran: %.3f, refl: %.3f, refr: %.3f\n", PMap->StoredPhotons, phong_sp->shine, phong_sp->wgt_specular, phong_sp->wgt_diffuse, phong_sp->transmit, phong_sp->reflect, phong_sp->refrac_index);*/
428  /* bu_log("spec[%d]: %.3f\n", PMap->StoredPhotons, phong_sp->wgt_specular);*/
429 
430  spec[0] = spec[1] = spec[2] = phong_sp->wgt_specular;
431  *refi = phong_sp->refrac_index;
432  *transmit = phong_sp->transmit;
433 
434  /*
435  if (phong_sp->refrac_index != 1.0)
436  bu_log("refi: %.3f\n", phong_sp->refrac_index);
437  */
438  }
439 
440  BU_PUT(phong_sp, struct phong_specific);
441 }
442 
443 
444 static fastf_t
445 MaxFloat(fastf_t a, fastf_t b, fastf_t c)
446 {
447  return a > b ? a > c ? a : c : b > c ? b : c;
448 }
449 
450 
451 /* break PHit <-> HitRef cycle */
452 int PHit(struct application *ap, struct partition *PartHeadp, struct seg *UNUSED(finished_segs));
453 
454 int
455 HitRef(struct application *ap, struct partition *PartHeadp, struct seg *UNUSED(finished_segs))
456 {
457  struct partition *part;
458  vect_t pt, normal, spec;
459  fastf_t refi, transmit;
460 
461  ap->a_hit = PHit;
462  part = PartHeadp->pt_forw;
463 
464  /* Compute Intersection Point, Pt = o + td */
465  VJOIN1(pt, ap->a_ray.r_pt, part->pt_outhit->hit_dist, ap->a_ray.r_dir);
466 
467  /* Fetch Intersection Normal */
468  RT_HIT_NORMAL(normal, part->pt_inhit, part->pt_inseg->seg_stp, &(ap->a_ray), part->pt_inflip);
469  /* RT_HIT_NORMAL(normal, part->pt_outhit, part->pt_inseg->seg_stp, &(ap->a_ray), part->pt_inflip);*/
470  /* RT_HIT_NORMAL(normal, part->pt_outhit, part->pt_outseg->seg_stp, &(ap->a_ray), part->pt_outflip);*/
471 
472  /* Assign pt */
473  ap->a_ray.r_pt[0] = pt[0];
474  ap->a_ray.r_pt[1] = pt[1];
475  ap->a_ray.r_pt[2] = pt[2];
476 
477 
478  /* Fetch Material */
479  GetMaterial(part->pt_regionp->reg_mater.ma_shader, spec, &refi, &transmit);
480 
481 
482  if (Refract(ap->a_ray.r_dir, normal, refi, 1.0)) {
483  /*
484  bu_log("1D: %d, [%.3f, %.3f, %.3f], [%.3f, %.3f, %.3f], [%.3f, %.3f, %.3f]\n", Depth, pt[0], pt[1], pt[2], ap->a_ray.r_dir[0], ap->a_ray.r_dir[1], ap->a_ray.r_dir[2], normal[0], normal[1], normal[2]);
485  bu_log("p1: [%.3f, %.3f, %.3f]\n", part->pt_inhit->hit_point[0], part->pt_inhit->hit_point[1], part->pt_inhit->hit_point[2]);
486  bu_log("p2: [%.3f, %.3f, %.3f]\n", part->pt_outhit->hit_point[0], part->pt_outhit->hit_point[1], part->pt_outhit->hit_point[2]);
487  */
488  Depth++;
489  rt_shootray(ap);
490  } else {
491  bu_log("TIF\n");
492  }
493 
494  ap->a_onehit = 0;
495  return 1;
496 }
497 
498 
499 /* Callback for Photon Hit, The 'current' photon is Emit[PMap->StoredPhotons] */
500 int
501 PHit(struct application *ap, struct partition *PartHeadp, struct seg *UNUSED(finished_segs))
502 {
503  struct partition *part;
504  vect_t pt, normal, color, spec, power;
505  fastf_t refi, transmit, prob, prob_diff, prob_spec, prob_ref;
506  int hit;
507 
508 
509  /* Move ptr forward to next region and call Hit recursively until reaching a region this is either
510  not a light or none at all */
511  hit = 0;
512  for (BU_LIST_FOR(part, partition, (struct bu_list *)PartHeadp)) {
513  if (part != PartHeadp) {
514  hit++;
515  VJOIN1(pt, ap->a_ray.r_pt, part->pt_inhit->hit_dist, ap->a_ray.r_dir);
516  /* printf("pt[%d][%d]: --- [%.3f, %.3f, %.3f], %s\n", hit, CheckMaterial("light", part->pt_regionp->reg_mater.ma_shader), pt[0], pt[1], pt[2], part->pt_regionp->reg_mater.ma_shader);*/
517 
518  if (!CheckMaterial("light", part->pt_regionp->reg_mater.ma_shader)) {
519  /* bu_log(" Found object!\n");*/
520  break;
521  }
522  }
523  }
524 
525  if (part == PartHeadp)
526  return 0;
527 
528 
529  if (CheckMaterial("light", part->pt_regionp->reg_mater.ma_shader))
530  return 0;
531 
532 
533  /* Compute Intersection Point, Pt = o + td */
534  VJOIN1(pt, ap->a_ray.r_pt, part->pt_inhit->hit_dist, ap->a_ray.r_dir);
535 
536 
537  /* Generate Bounding Box for Scaling Phase */
538  if (PInit) {
539  BBMin[0] = BBMax[0] = pt[0];
540  BBMin[1] = BBMax[1] = pt[1];
541  BBMin[2] = BBMax[2] = pt[2];
542  PInit = 0;
543  } else {
544  if (pt[0] < BBMin[0])
545  BBMin[0] = pt[0];
546  if (pt[0] > BBMax[0])
547  BBMax[0] = pt[0];
548  if (pt[1] < BBMin[1])
549  BBMin[1] = pt[1];
550  if (pt[1] > BBMax[1])
551  BBMax[1] = pt[1];
552  if (pt[2] < BBMin[2])
553  BBMin[2] = pt[2];
554  if (pt[2] > BBMax[2])
555  BBMax[2] = pt[2];
556  }
557 
558  /* Fetch Intersection Normal */
559  RT_HIT_NORMAL(normal, part->pt_inhit, part->pt_inseg->seg_stp, &(ap->a_ray), part->pt_inflip);
560 
561  /* Fetch Material */
562  GetMaterial(part->pt_regionp->reg_mater.ma_shader, spec, &refi, &transmit);
563  /* bu_log("Spec: [%.3f, %.3f, %.3f], Refi: %.3f\n", spec[0], spec[1], spec[2], refi);*/
564 
565 
566  /* Compute Diffuse, Specular, and Caustics */
567  color[0] = part->pt_regionp->reg_mater.ma_color[0];
568  color[1] = part->pt_regionp->reg_mater.ma_color[1];
569  color[2] = part->pt_regionp->reg_mater.ma_color[2];
570 
571  prob_ref = MaxFloat(color[0]+spec[0], color[1]+spec[1], color[2]+spec[2]);
572  prob_diff = ((color[0]+color[1]+color[2])/(color[0]+color[1]+color[2]+spec[0]+spec[1]+spec[2]))*prob_ref;
573  prob_spec = prob_ref - prob_diff;
574 #ifndef HAVE_DRAND48
575  prob = rand()/(double)RAND_MAX;
576 #else
577  prob = drand48();
578 #endif
579 
580  /* bu_log("pr: %.3f, pd: %.3f, [%.3f, %.3f, %.3f] [%.3f, %.3f, %.3f]\n", prob_ref, prob_diff, color[0], color[1], color[2], spec[0], spec[1], spec[2]);*/
581  /* bu_log("prob: %.3f, prob_diff: %.3f, pd+ps: %.3f\n", prob, prob_diff, prob_diff+prob_spec);*/
582 
583  if (prob < 1.0 - transmit) {
584  if (prob < prob_diff) {
585  /* Store power of incident Photon */
586  power[0] = CurPh.Power[0];
587  power[1] = CurPh.Power[1];
588  power[2] = CurPh.Power[2];
589 
590 
591  /* Scale Power of reflected photon */
592  CurPh.Power[0] = power[0]*color[0]/prob_diff;
593  CurPh.Power[1] = power[1]*color[1]/prob_diff;
594  CurPh.Power[2] = power[2]*color[2]/prob_diff;
595 
596  /* Store Photon */
597  Store(pt, ap->a_ray.r_dir, normal, PType);
598 
599  /* Assign diffuse reflection direction */
600  DiffuseReflect(normal, ap->a_ray.r_dir);
601 
602  /* Assign pt */
603  ap->a_ray.r_pt[0] = pt[0];
604  ap->a_ray.r_pt[1] = pt[1];
605  ap->a_ray.r_pt[2] = pt[2];
606 
607  if (PType != PM_CAUSTIC) {
608  Depth++;
609  rt_shootray(ap);
610  }
611  } else if (prob >= prob_diff && prob < prob_diff + prob_spec) {
612  /* Store power of incident Photon */
613  power[0] = CurPh.Power[0];
614  power[1] = CurPh.Power[1];
615  power[2] = CurPh.Power[2];
616 
617  /* Scale power of reflected photon */
618  CurPh.Power[0] = power[0]*spec[0]/prob_spec;
619  CurPh.Power[1] = power[1]*spec[1]/prob_spec;
620  CurPh.Power[2] = power[2]*spec[2]/prob_spec;
621 
622  /* Reflective */
623  SpecularReflect(normal, ap->a_ray.r_dir);
624 
625  /* Assign pt */
626  ap->a_ray.r_pt[0] = pt[0];
627  ap->a_ray.r_pt[1] = pt[1];
628  ap->a_ray.r_pt[2] = pt[2];
629 
630  if (PType != PM_IMPORTANCE)
631  PType = PM_CAUSTIC;
632  Depth++;
633  rt_shootray(ap);
634  } else {
635  /* Store Photon */
636  Store(pt, ap->a_ray.r_dir, normal, PType);
637  }
638  } else {
639  if (refi > 1.0 && (PType == PM_CAUSTIC || Depth == 0)) {
640  if (PType != PM_IMPORTANCE)
641  PType = PM_CAUSTIC;
642 
643  /* Store power of incident Photon */
644  power[0] = CurPh.Power[0];
645  power[1] = CurPh.Power[1];
646  power[2] = CurPh.Power[2];
647 
648  /* Scale power of reflected photon */
649  CurPh.Power[0] = power[0]*spec[0]/prob_spec;
650  CurPh.Power[1] = power[1]*spec[1]/prob_spec;
651  CurPh.Power[2] = power[2]*spec[2]/prob_spec;
652 
653  /* Refractive or Reflective */
654  if (refi > 1.0 && prob < transmit) {
655  CurPh.Power[0] = power[0];
656  CurPh.Power[1] = power[1];
657  CurPh.Power[2] = power[2];
658 
659  if (!Refract(ap->a_ray.r_dir, normal, 1.0, refi))
660  printf("TIF0\n");
661 
662  ap->a_hit = HitRef;
663 
664  /*
665  bu_log("dir: [%.3f, %.3f, %.3f]\n", ap->a_ray.r_dir[0], ap->a_ray.r_dir[1], ap->a_ray.r_dir[2]);
666  bu_log("ref: [%.3f, %.3f, %.3f]\n", ap->a_ray.r_dir[0], ap->a_ray.r_dir[1], ap->a_ray.r_dir[2]);
667  bu_log("nor: [%.3f, %.3f, %.3f]\n", normal[0], normal[1], normal[2]);
668  bu_log("p0: [%.3f, %.3f, %.3f], [%.3f, %.3f, %.3f]\n", pt[0], pt[1], pt[2], ap->a_ray.r_dir[0], ap->a_ray.r_dir[1], ap->a_ray.r_dir[2]);
669  */
670  ap->a_onehit = 0;
671  } else {
672  SpecularReflect(normal, ap->a_ray.r_dir);
673  }
674 
675  /* Assign pt */
676  ap->a_ray.r_pt[0] = pt[0];
677  ap->a_ray.r_pt[1] = pt[1];
678  ap->a_ray.r_pt[2] = pt[2];
679 
680  /* bu_log("2D: %d, [%.3f, %.3f, %.3f], [%.3f, %.3f, %.3f], [%.3f, %.3f, %.3f]\n", Depth, pt[0], pt[1], pt[2], ap->a_ray.r_dir[0], ap->a_ray.r_dir[1], ap->a_ray.r_dir[2], normal[0], normal[1], normal[2]);*/
681  Depth++;
682  rt_shootray(ap);
683  }
684  }
685  return 1;
686 }
687 
688 
689 /* Callback for Photon Miss */
690 int
692 {
693  return 0;
694 }
695 
696 
697 /* ScalePhotonPower() is used to scale the power of all photons once they
698  * have been emitted from the light source. Scale = 1/(#emitted photons).
699  * Call this function after each light source is processed.
700  * This function also handles setting a default power for the photons based
701  * on the size of the scene, i.e. power of light source */
702 void
704 {
705  int i;
706 
707  for (i = 0; i < PMap[map]->StoredPhotons; i++) {
708  Emit[map][i].Power[0] *= ScaleFactor/(double)EPS[map];
709  Emit[map][i].Power[1] *= ScaleFactor/(double)EPS[map];
710  Emit[map][i].Power[2] *= ScaleFactor/(double)EPS[map];
711  }
712 }
713 
714 
715 /* Generate Importons and emit them into the scene from the eye position */
716 void
717 EmitImportonsRandom(struct application *ap, point_t eye_pos)
718 {
719  while (PMap[PM_IMPORTANCE]->StoredPhotons < PMap[PM_IMPORTANCE]->MaxPhotons) {
720  do {
721  /* Set Ray Direction to application ptr */
722 #ifndef HAVE_DRAND48
723  ap->a_ray.r_dir[0] = 2.0*rand()/(double)RAND_MAX-1.0;
724  ap->a_ray.r_dir[1] = 2.0*rand()/(double)RAND_MAX-1.0;
725  ap->a_ray.r_dir[2] = 2.0*rand()/(double)RAND_MAX-1.0;
726 #else
727  ap->a_ray.r_dir[0] = 2.0*drand48()-1.0;
728  ap->a_ray.r_dir[1] = 2.0*drand48()-1.0;
729  ap->a_ray.r_dir[2] = 2.0*drand48()-1.0;
730 #endif
731  } while (ap->a_ray.r_dir[0]*ap->a_ray.r_dir[0] + ap->a_ray.r_dir[1]*ap->a_ray.r_dir[1] + ap->a_ray.r_dir[2]*ap->a_ray.r_dir[2] > 1);
732 
733  /* Normalize Ray Direction */
734  VUNITIZE(ap->a_ray.r_dir);
735 
736  /* Set Ray Position to application ptr */
737  ap->a_ray.r_pt[0] = eye_pos[0];
738  ap->a_ray.r_pt[1] = eye_pos[1];
739  ap->a_ray.r_pt[2] = eye_pos[2];
740 
741  /* Shoot Importon into Scene */
742  CurPh.Power[0] = 0;
743  CurPh.Power[1] = 100000000;
744  CurPh.Power[2] = 0;
745 
746  Depth = 0;
747  PType = PM_IMPORTANCE;
748  rt_shootray(ap);
749  }
750 }
751 
752 
753 /* Emit a photons in a random direction based on a point light */
754 void
755 EmitPhotonsRandom(struct application *ap, double ScaleIndirect)
756 {
757  struct light_specific *lp;
758  int i;
759 
760  /*
761  for (i = 0; i < 8; i++)
762  bu_log("sample points: [%.3f, %.3f, %.3f]\n", lp->lt_sample_pts[i].lp_pt[0], lp->lt_sample_pts[i].lp_pt[1], lp->lt_sample_pts[i].lp_pt[2]);
763  */
764  while (1) {
765  for (BU_LIST_FOR(lp, light_specific, &(LightHead.l))) {
766  /* If the Global Photon Map Completes before the Caustics Map, then it probably means there are no caustic objects in the Scene */
767  if (PMap[PM_GLOBAL]->StoredPhotons == PMap[PM_GLOBAL]->MaxPhotons && (!PMap[PM_CAUSTIC]->StoredPhotons || PMap[PM_CAUSTIC]->StoredPhotons == PMap[PM_CAUSTIC]->MaxPhotons))
768  return;
769 
770  do {
771  /* do {*/
772  /* Set Ray Direction to application ptr */
773 #ifndef HAVE_DRAND48
774  ap->a_ray.r_dir[0] = 2.0*rand()/(double)RAND_MAX-1.0;
775  ap->a_ray.r_dir[1] = 2.0*rand()/(double)RAND_MAX-1.0;
776  ap->a_ray.r_dir[2] = 2.0*rand()/(double)RAND_MAX-1.0;
777 #else
778  ap->a_ray.r_dir[0] = 2.0*drand48()-1.0;
779  ap->a_ray.r_dir[1] = 2.0*drand48()-1.0;
780  ap->a_ray.r_dir[2] = 2.0*drand48()-1.0;
781 #endif
782  } while (ap->a_ray.r_dir[0]*ap->a_ray.r_dir[0] + ap->a_ray.r_dir[1]*ap->a_ray.r_dir[1] + ap->a_ray.r_dir[2]*ap->a_ray.r_dir[2] > 1);
783  /* Normalize Ray Direction */
784  VUNITIZE(ap->a_ray.r_dir);
785 
786  /* Set Ray Position to application ptr */
787  ap->a_ray.r_pt[0] = lp->lt_pos[0];
788  ap->a_ray.r_pt[1] = lp->lt_pos[1];
789  ap->a_ray.r_pt[2] = lp->lt_pos[2];
790 
791 
792  /* Shoot Photon into Scene, (4.0) is used to align phong's attenuation with photonic energies, it's a heuristic */
793  /*bu_log("Shooting Ray: [%.3f, %.3f, %.3f] [%.3f, %.3f, %.3f]\n", lp->lt_pos[0], lp->lt_pos[1], lp->lt_pos[2], x, y, z);*/
794  CurPh.Power[0] = 1000.0 * ScaleIndirect * lp->lt_intensity * lp->lt_color[0];
795  CurPh.Power[1] = 1000.0 * ScaleIndirect * lp->lt_intensity * lp->lt_color[1];
796  CurPh.Power[2] = 1000.0 * ScaleIndirect * lp->lt_intensity * lp->lt_color[2];
797 
798  Depth = 0;
799  PType = PM_GLOBAL;
800 
801  EPL++;
802  for (i = 0; i < PM_MAPS; i++)
803  if (PMap[i]->StoredPhotons < PMap[i]->MaxPhotons)
804  EPS[i]++;
805 
806  rt_shootray(ap);
807  /* bu_log("1: %d, 2: %d\n", PMap[PM_GLOBAL]->StoredPhotons, PMap[PM_CAUSTIC]->StoredPhotons);*/
808  }
809  }
810 }
811 
812 
813 void
814 SanityCheck(struct PNode *Root, int LR)
815 {
816  if (!Root)
817  return;
818 
819  bu_log("Pos[%d]: [%.3f, %.3f, %.3f]\n", LR, Root->P.Pos[0], Root->P.Pos[1], Root->P.Pos[2]);
820  SanityCheck(Root->L, 1);
821  SanityCheck(Root->R, 2);
822 }
823 
824 
825 fastf_t
827 {
828  return 0.918 * (1.0 - (1.0 - exp(-1.953*dist*dist/(2.0*rad*rad)))/(1.0 - exp(-1.953)));
829 }
830 
831 
832 fastf_t
834 {
835  return 1.0 - dist/rad;
836 }
837 
838 
839 void
840 GetEstimate(vect_t irrad, point_t pos, vect_t normal, fastf_t rad, int np, int map, double max_rad, int centog, int min_np)
841 {
842  struct PhotonSearch Search;
843  int i;
844  fastf_t tmp, dist, Filter, ScaleFilter;
845  vect_t t, Centroid;
846 
847 
848  irrad[0] = irrad[1] = irrad[2] = 0;
849  if (!PMap[map]->StoredPhotons) return;
850 
851  Search.Pos[0] = pos[0];
852  Search.Pos[1] = pos[1];
853  Search.Pos[2] = pos[2];
854 
855  Search.RadSq = rad*rad/4.0;
856  Search.Max = np < min_np ? min_np : np;
857  Search.Normal[0] = normal[0];
858  Search.Normal[1] = normal[1];
859  Search.Normal[2] = normal[2];
860 
861  Search.List = (struct PSN*)bu_calloc(Search.Max, sizeof(struct PSN), "PSN");
862  do {
863  Search.Found = 0;
864  Search.RadSq *= 4.0;
865  LocatePhotons(&Search, PMap[map]->Root);
866  if (!Search.Found && Search.RadSq > ScaleFactor*ScaleFactor/100.0)
867  break;
868  } while (Search.Found < Search.Max && Search.RadSq < max_rad*max_rad);
869 
870  /* bu_log("Found: %d\n", Search.Found);*/
871  if (Search.Found < min_np) {
872  bu_free(Search.List, "Search.List");
873  return;
874  }
875 
876  /* Calculate Max Distance */
877  Search.RadSq = 1;
878  Centroid[0] = Centroid[1] = Centroid[2] = 0;
879 
880  for (i = 0; i < Search.Found; i++) {
881  t[0] = Search.List[i].P.Pos[0] - pos[0];
882  t[1] = Search.List[i].P.Pos[1] - pos[1];
883  t[2] = Search.List[i].P.Pos[2] - pos[2];
884 
885  Centroid[0] += Search.List[i].P.Pos[0];
886  Centroid[1] += Search.List[i].P.Pos[1];
887  Centroid[2] += Search.List[i].P.Pos[2];
888 
889  dist = t[0]*t[0] + t[1]*t[1] + t[2]*t[2];
890  if (dist > Search.RadSq)
891  Search.RadSq = dist;
892  }
893 
894  if (Search.Found) {
895  Centroid[0] /= (double)Search.Found;
896  Centroid[1] /= (double)Search.Found;
897  Centroid[2] /= (double)Search.Found;
898  }
899 
900 
901  /* This needs a little debugging, splotches in moss cause tmp gets too small, will look at later, ||1 to turn it off */
902  if (!centog||1) {
903  Centroid[0] = pos[0];
904  Centroid[1] = pos[1];
905  Centroid[2] = pos[2];
906  ScaleFilter = 2.0;
907  } else {
908  ScaleFilter = 1.0;
909  }
910 
911  for (i = 0; i < Search.Found; i++) {
912  t[0] = Search.List[i].P.Pos[0] - pos[0];
913  t[1] = Search.List[i].P.Pos[1] - pos[1];
914  t[2] = Search.List[i].P.Pos[2] - pos[2];
915 
916  dist = t[0]*t[0] + t[1]*t[1] + t[2]*t[2];
917  /* Filter= 0.50;*/
918  /* Filter= ConeFilter(dist, NP.RadSq);*/
919  Filter = 0.5*GaussFilter(dist, Search.RadSq);
920 
921  irrad[0] += Search.List[i].P.Power[0]*Filter*ScaleFilter;
922  irrad[1] += Search.List[i].P.Power[1]*Filter*ScaleFilter;
923  irrad[2] += Search.List[i].P.Power[2]*Filter*ScaleFilter;
924  }
925 
926  tmp = M_PI*Search.RadSq;
927  t[0] = sqrt((Centroid[0] - pos[0])*(Centroid[0] - pos[0])+(Centroid[1] - pos[1])*(Centroid[1] - pos[1])+(Centroid[2] - pos[2])*(Centroid[2] - pos[2]));
928  tmp = M_PI*(sqrt(Search.RadSq)-t[0])*(sqrt(Search.RadSq)-t[0]);
929 
930 
931  irrad[0] /= tmp;
932  irrad[1] /= tmp;
933  irrad[2] /= tmp;
934 
935  /*
936  if (irrad[0] > 10 || irrad[1] > 10 || irrad[2] > 10) {
937  bu_log("found: %d, tmp: %.1f\n", Search.Found, tmp);
938  }
939  */
940  if (map == PM_CAUSTIC) {
941  tmp = (double)Search.Found/(double)Search.Max;
942  irrad[0] *= tmp;
943  irrad[1] *= tmp;
944  irrad[2] *= tmp;
945  }
946 
947  /*
948  irrad[0] *= M_1_PI / NP.RadSq;
949  irrad[1] *= M_1_PI / NP.RadSq;
950  irrad[2] *= M_1_PI / NP.RadSq;
951  */
952  bu_free(Search.List, "Search.List");
953  /* bu_log("Radius: %.3f, Max Phot: %d, Found: %d, Power: [%.4f, %.4f, %.4f], Pos: [%.3f, %.3f, %.3f]\n", sqrt(NP.RadSq), NP.Max, NP.Found, irrad[0], irrad[1], irrad[2], pos[0], pos[1], pos[2]);*/
954 }
955 
956 
957 int
958 ICHit(struct application *ap, struct partition *PartHeadp, struct seg *UNUSED(finished_segs))
959 {
960  struct partition *part;
961  point_t pt, normal;
962  vect_t C1, C2;
963 
964  part = PartHeadp->pt_forw;
965 
966  VJOIN1(pt, ap->a_ray.r_pt, part->pt_inhit->hit_dist, ap->a_ray.r_dir);
967 
968  RT_HIT_NORMAL(normal, part->pt_inhit, part->pt_inseg->seg_stp, &(ap->a_ray), part->pt_inflip);
969  /* GetEstimate(C1, pt, normal, ScaleFactor/10.0, PMap[PM_GLOBAL]->StoredPhotons / 100, PM_GLOBAL, 5, 1);*/
970  /* GetEstimate(C1, pt, normal, ScaleFactor/1000.0, 10.0*log(PMap[PM_GLOBAL]->StoredPhotons), PM_GLOBAL, ScaleFactor/5.0, 1);*/
971  GetEstimate(C1, pt, normal, ScaleFactor/1024.0, 128, PM_GLOBAL, ScaleFactor/8.0, 1, 15);
972  /* GetEstimate(C2, pt, normal, (int)(ScaleFactor/pow(2, (log(PMap[PM_CAUSTIC]->MaxPhotons/2)/log(4)))), PMap[PM_CAUSTIC]->MaxPhotons / 50, PM_CAUSTIC, 0, 0);*/
973  GetEstimate(C2, pt, normal, ScaleFactor/1024.0, PMap[PM_CAUSTIC]->MaxPhotons / 100, PM_CAUSTIC, ScaleFactor/256.0, 1, 15);
974  /* GetEstimate(IMColor2, pt, normal, (int)(ScaleFactor/100.0), PMap[PM_CAUSTIC]->MaxPhotons/50, PM_CAUSTIC, 1, 0);*/
975 
976  (*(vect_t*)ap->a_purpose)[0] += C1[0] + C2[0];
977  (*(vect_t*)ap->a_purpose)[1] += C1[1] + C2[1];
978  (*(vect_t*)ap->a_purpose)[2] += C1[2] + C2[2];
979 
980  return 1;
981 }
982 
983 
984 int
986 {
987  /* Set to Background/Ambient Color later */
988  return 0;
989 }
990 
991 
992 /* Convert a Polar vector into a euclidian vector:
993  * - Requires that an orthogonal basis, so generate one using the photon normal as the first vector.
994  * - The Normal passed to this function must be unitized.
995  * - This took me almost 3 hours to write, and it's tight.
996  */
997 void
998 Polar2Euclidian(vect_t Dir, vect_t Normal, double Theta, double Phi)
999 {
1000  int i;
1001  vect_t BasisX, BasisY;
1002 
1003 
1004  BasisX[0] = fabs(Normal[0]) < 0.001 ? 1.0 : fabs(Normal[0]) > 0.999 ? 0.0 : -(1.0/Normal[0]);
1005  BasisX[1] = fabs(Normal[1]) < 0.001 ? 1.0 : fabs(Normal[1]) > 0.999 ? 0.0 : (1.0/Normal[1]);
1006  BasisX[2] = 0.0;
1007  VUNITIZE(BasisX);
1008  VCROSS(BasisY, Normal, BasisX);
1009 
1010  for (i = 0; i < 3; i++)
1011  Dir[i] = sin(Theta)*cos(Phi)*BasisX[i] + sin(Theta)*sin(Phi)*BasisY[i] + cos(Theta)*Normal[i];
1012 }
1013 
1014 
1015 /*
1016  * Irradiance Calculation for a given position
1017  */
1018 void
1019 Irradiance(int pid, struct Photon *P, struct application *ap)
1020 {
1021  struct application *lap; /* local application instance */
1022  int i, j, M, N;
1023  double theta, phi, Coef;
1024 
1025  BU_ALLOC(lap, struct application);
1026  RT_APPLICATION_INIT(lap);
1027  lap->a_rt_i = ap->a_rt_i;
1028  lap->a_hit = ap->a_hit;
1029  lap->a_miss = ap->a_miss;
1030  lap->a_resource = &GPM_RTAB[pid];
1031  lap->a_logoverlap = ap->a_logoverlap;
1032 
1033  M = N = GPM_RAYS;
1034  P->Irrad[0] = P->Irrad[1] = P->Irrad[2] = 0.0;
1035  for (i = 1; i <= M; i++) {
1036  for (j = 1; j <= N; j++) {
1037 #ifndef HAVE_DRAND48
1038  theta = asin(sqrt((j-rand()/(double)RAND_MAX)/M));
1039  phi = (M_2PI)*((i-rand()/(double)RAND_MAX)/N);
1040 #else
1041  theta = asin(sqrt((j-drand48())/M));
1042  phi = (M_2PI)*((i-drand48())/N);
1043 #endif
1044 
1045  /* Assign pt */
1046  lap->a_ray.r_pt[0] = P->Pos[0];
1047  lap->a_ray.r_pt[1] = P->Pos[1];
1048  lap->a_ray.r_pt[2] = P->Pos[2];
1049 
1050  /* Assign Dir */
1051  Polar2Euclidian(lap->a_ray.r_dir, P->Normal, theta, phi);
1052 
1053  /* Utilize the purpose pointer as a pointer to the Irradiance Color */
1054  lap->a_purpose = (const char *)P->Irrad;
1055 
1056  /* bu_log("Vec: [%.3f, %.3f, %.3f]\n", ap->a_ray.r_dir[0], ap->a_ray.r_dir[1], ap->a_ray.r_dir[2]);*/
1057  rt_shootray(lap);
1058 
1059  /* bu_log("[%.3f, %.3f, %.3f] [%.3f, %.3f, %.3f] [%.3f, %.3f, %.3f]\n", P.Pos[0], P.Pos[1], P.Pos[2], P.Normal[0], P.Normal[1], P.Normal[2], IMColor[0], IMColor[1], IMColor[2]);*/
1060  }
1061  }
1062 
1063  Coef = 1.0/((double)(M*N));
1064  P->Irrad[0] *= Coef;
1065  P->Irrad[1] *= Coef;
1066  P->Irrad[2] *= Coef;
1067 
1068  bu_free(lap, "app");
1069 }
1070 
1071 
1072 /*
1073  * Irradiance Cache for Indirect Illumination
1074  * Go through each photon and use it for the position of the hemisphere
1075  * and then determine whether that should be included as a Cache Pt.
1076  */
1077 void
1078 BuildIrradianceCache(int pid, struct PNode *Node, struct application *ap)
1079 {
1080  if (!Node)
1081  return;
1082 
1083  /* Determine if this pt will be used by calculating a weight */
1084  bu_semaphore_acquire(PM_SEM);
1085  if (!Node->C) {
1086  ICSize++;
1087  Node->C++;
1088 #ifndef HAVE_ALARM
1089  if (!(ICSize%(PMap[PM_GLOBAL]->MaxPhotons/8)))
1090  bu_log(" Irradiance Cache Progress: %d%%\n", (int)(0.5+100.0*ICSize/PMap[PM_GLOBAL]->MaxPhotons));
1091 #endif
1092  bu_semaphore_release(PM_SEM);
1093  Irradiance(pid, &Node->P, ap);
1094  } else {
1095  bu_semaphore_release(PM_SEM);
1096  }
1097 
1098  BuildIrradianceCache(pid, Node->L, ap);
1099  BuildIrradianceCache(pid, Node->R, ap);
1100 }
1101 
1102 
1103 #ifdef HAVE_ALARM
1104 static int starttime = 0;
1105 void
1106 alarmhandler(int sig)
1107 {
1108  int t;
1109  float p, tl;
1110  if (sig != SIGALRM)
1111  bu_bomb("Funky signals\n");
1112  t = time(NULL) - starttime;
1113  p = (float)ICSize/PMap[PM_GLOBAL]->MaxPhotons + .015;
1114  tl = (float)t*1.0/p - t;
1115  bu_log(" Irradiance Cache Progress: %d%% Approximate time left: %f%f",
1116  (int)(100.0*p), (1.0/p-1.0)*(float)t, (float)t*1.0/p);
1117 #define BAH(s, w) if (tl > (s)) { float d = floor(tl / (((s) == 0)?1.0:(float)(s))); \
1118  tl -= d * (s); bu_log("%d "w, (int)d, d>1?"s":""); }
1119  BAH(60*60*24, "day%s, ");
1120  BAH(60*60, "hour%s, ");
1121  BAH(60, "minute%s, ");
1122  BAH(0, "second%s.");
1123 #undef BAH
1124  bu_log("\n");
1125  alarm(60);
1126  return;
1127 }
1128 #endif
1129 
1130 void
1131 IrradianceThread(int pid, void *arg)
1132 {
1133 #ifdef HAVE_ALARM
1134  starttime = time(NULL);
1135  signal(SIGALRM, alarmhandler);
1136  alarm(60);
1137  BuildIrradianceCache(pid, PMap[PM_GLOBAL]->Root, (struct application*)arg);
1138  alarm(0);
1139  starttime = 0;
1140 #else
1141  BuildIrradianceCache(pid, PMap[PM_GLOBAL]->Root, (struct application*)arg);
1142 #endif
1143 }
1144 
1145 
1146 void
1147 Initialize(int MAP, int MapSize)
1148 {
1149  BU_ALLOC(PMap[MAP], struct PhotonMap);
1150  PMap[MAP]->MaxPhotons = MapSize;
1151 
1152  BU_ALLOC(PMap[MAP]->Root, struct PNode);
1153  PMap[MAP]->StoredPhotons = 0;
1154  if (MapSize > 0)
1155  BU_ALLOC(Emit[MAP], struct Photon);
1156  else
1157  Emit[MAP] = NULL;
1158 }
1159 
1160 
1161 int
1162 LoadFile(char *pmfile)
1163 {
1164  size_t ret;
1165  FILE *FH;
1166  int I1 = 0, i;
1167  short S1;
1168  char C1;
1169 
1170  FH = fopen(pmfile, "rb");
1171  if (FH) {
1172  bu_log(" Reading Irradiance Cache File...\n");
1173  ret = fread(&S1, sizeof(short), 1, FH);
1174  if (ret != 1) {
1175  fclose(FH);
1176  bu_log("Error reading irradiance cache file (endian)\n");
1177  return 0;
1178  }
1179 
1180  bu_log("endian: %d\n", S1);
1181 
1182  ret = fread(&S1, sizeof(short), 1, FH);
1183  if (ret != 1) {
1184  fclose(FH);
1185  bu_log("Error reading irradiance cache file (revision)\n");
1186  return 0;
1187  }
1188  bu_log("revision: %d\n", S1);
1189 
1190  ret = fread(&ScaleFactor, sizeof(double), 1, FH);
1191  if (ret != 1) {
1192  fclose(FH);
1193  bu_log("Error reading irradiance cache file (scale factor)\n");
1194  return 0;
1195  }
1196  bu_log("Scale Factor: %.3f\n", ScaleFactor);
1197 
1198  /* Read in Map Type */
1199  ret = fread(&C1, sizeof(char), 1, FH);
1200  if (ret != 1) {
1201  fclose(FH);
1202  bu_log("Error reading irradiance cache file (map type)\n");
1203  return 0;
1204  }
1205  ret = fread(&I1, sizeof(int), 1, FH);
1206  if (ret != 1 || I1 < 0 || I1 > INT_MAX) {
1207  fclose(FH);
1208  bu_log("Error reading irradiance cache file (map type)\n");
1209  return 0;
1210  }
1211 
1212  Initialize(PM_GLOBAL, I1);
1213  bu_log("Reading Global: %d\n", I1);
1214  for (i = 0; i < I1; i++) {
1215  ret = fread(&Emit[PM_GLOBAL][i], sizeof(struct Photon), 1, FH);
1216  if (ret != 1) {
1217  fclose(FH);
1218  bu_log("Error reading irradiance cache file (global)\n");
1219  /* bu_log("Pos: [%.3f, %.3f, %.3f], Power: [%.3f, %.3f, %.3f]\n", Emit[PM_GLOBAL][i].Pos[0], Emit[PM_GLOBAL][i].Pos[1], Emit[PM_GLOBAL][i].Pos[2], Emit[PM_GLOBAL][i].Power[0], Emit[PM_GLOBAL][i].Power[1], Emit[PM_GLOBAL][i].Power[2]);*/
1220  /* bu_log("Pos: [%.3f, %.3f, %.3f], Irrad: [%.3f, %.3f, %.3f]\n", Emit[PM_GLOBAL][i].Pos[0], Emit[PM_GLOBAL][i].Pos[1], Emit[PM_GLOBAL][i].Pos[2], Emit[PM_GLOBAL][i].Irrad[0], Emit[PM_GLOBAL][i].Irrad[1], Emit[PM_GLOBAL][i].Irrad[2]);*/
1221  return 0;
1222  }
1223  }
1224 
1225  ret = fread(&C1, sizeof(char), 1, FH);
1226  if (ret != 1) {
1227  fclose(FH);
1228  bu_log("Error reading irradiance cache file (C1)\n");
1229  return 0;
1230  }
1231 
1232  ret = fread(&I1, sizeof(int), 1, FH);
1233  if (ret != 1 || I1 < 0 || I1 > INT_MAX) {
1234  fclose(FH);
1235  bu_log("Error reading irradiance cache file (l1)\n");
1236  return 0;
1237  }
1238 
1239  Initialize(PM_CAUSTIC, I1);
1240  bu_log("Reading Caustic: %d\n", I1);
1241  for (i = 0; i < I1; i++) {
1242  ret = fread(&Emit[PM_CAUSTIC][i], sizeof(struct Photon), 1, FH);
1243  if (ret != 1) {
1244  fclose(FH);
1245  bu_log("Error reading irradiance cache file (caustic)\n");
1246  return 0;
1247  }
1248  }
1249 
1250  PMap[PM_GLOBAL]->StoredPhotons = PMap[PM_GLOBAL]->MaxPhotons;
1251  BuildTree(Emit[PM_GLOBAL], PMap[PM_GLOBAL]->StoredPhotons, PMap[PM_GLOBAL]->Root);
1252 
1253  PMap[PM_CAUSTIC]->StoredPhotons = PMap[PM_CAUSTIC]->MaxPhotons;
1254  BuildTree(Emit[PM_CAUSTIC], PMap[PM_CAUSTIC]->StoredPhotons, PMap[PM_CAUSTIC]->Root);
1255  fclose(FH);
1256  return 1;
1257  }
1258 
1259  return 0;
1260 }
1261 
1262 
1263 void
1264 WritePhotons(struct PNode *Root, FILE *FH)
1265 {
1266  size_t ret;
1267  if (!Root)
1268  return;
1269 
1270  ret = fwrite(&Root->P, sizeof(struct Photon), 1, FH);
1271  if (ret != 1)
1272  bu_log("Unable to write photons\n");
1273  WritePhotons(Root->L, FH);
1274  WritePhotons(Root->R, FH);
1275 }
1276 
1277 
1278 void
1279 WritePhotonFile(char *pmfile)
1280 {
1281  FILE *FH;
1282  int I1;
1283  short S1;
1284  char C1;
1285  size_t ret;
1286 
1287  FH = fopen(pmfile, "wb");
1288  if (FH) {
1289  /* Write 2 Byte Endian Code and 2 Byte Revision Code */
1290  S1 = 1;
1291  ret = fwrite(&S1, sizeof(short), 1, FH);
1292  if (ret != 1)
1293  bu_log("Error writing irradiance cache file (endian)\n");
1294 
1295  S1 = 0;
1296  ret = fwrite(&S1, sizeof(short), 1, FH);
1297  if (ret != 1)
1298  bu_log("Error writing irradiance cache file (revision)\n");
1299 
1300  /* Write Scale Factor */
1301  bu_log("writing sf: %.3f\n", ScaleFactor);
1302  ret = fwrite(&ScaleFactor, sizeof(double), 1, FH);
1303  if (ret != 1)
1304  bu_log("Error writing irradiance cache file (scale factor)\n");
1305 
1306  /* === Write PM_GLOBAL Data === */
1307  C1 = PM_GLOBAL;
1308  ret = fwrite(&C1, sizeof(char), 1, FH);
1309  if (ret != 1)
1310  bu_log("Error writing irradiance cache file (PM_GLOBAL)\n");
1311 
1312  /* Write number of Photons */
1313  I1 = PMap[PM_GLOBAL]->StoredPhotons;
1314  ret = fwrite(&I1, sizeof(int), 1, FH);
1315  if (ret != 1)
1316  bu_log("Error writing irradiance cache file (photons)\n");
1317 
1318  /* Write each photon to file */
1319  if (PMap[PM_GLOBAL]->StoredPhotons)
1320  WritePhotons(PMap[PM_GLOBAL]->Root, FH);
1321 
1322  /* === Write PM_CAUSTIC Data === */
1323  C1 = PM_CAUSTIC;
1324  ret = fwrite(&C1, sizeof(char), 1, FH);
1325  if (ret != 1)
1326  bu_log("Error writing irradiance cache file (PM_CAUSTIC)\n");
1327 
1328  /* Write number of Photons */
1329  I1 = PMap[PM_CAUSTIC]->StoredPhotons;
1330  ret = fwrite(&I1, sizeof(int), 1, FH);
1331  if (ret != 1)
1332  bu_log("Error writing irradiance cache file (number of photons)\n");
1333 
1334  /* Write each photon to file */
1335  if (PMap[PM_CAUSTIC]->StoredPhotons)
1336  WritePhotons(PMap[PM_CAUSTIC]->Root, FH);
1337 
1338  fclose(FH);
1339  }
1340 }
1341 
1342 
1343 /*
1344  * Main Photon Mapping Function
1345  */
1346 void
1347 BuildPhotonMap(struct application *ap, point_t eye_pos, int cpus, int width, int height, int UNUSED(Hypersample), int GlobalPhotons, double CausticsPercent, int Rays, double AngularTolerance, int RandomSeed, int ImportanceMapping, int IrradianceHypersampling, int VisualizeIrradiance, double ScaleIndirect, char pmfile[255])
1348 {
1349  int i, MapSize[PM_MAPS];
1350  double ratio;
1351 
1352  PM_Visualize = VisualizeIrradiance;
1353  GPM_IH = IrradianceHypersampling;
1354  GPM_WIDTH = width;
1355  GPM_HEIGHT = height;
1356 
1357  /* If the user has specified a cache file then first check to see if there is any valid data within it,
1358  otherwise utilize the file to push the resulting irradiance cache data into for future use. */
1359  if (!LoadFile(pmfile)) {
1360  /*
1361  bu_log("pos: [%.3f, %.3f, %.3f]\n", eye_pos[0], eye_pos[1], eye_pos[2]);
1362  bu_log("I, V, Imp, H: %.3f, %d, %d, %d\n", LightIntensity, VisualizeIrradiance, ImportanceMapping, IrradianceHypersampling);
1363  */
1364  bu_log("Building Photon Map:\n");
1365 
1366 
1367  GPM_RAYS = Rays;
1368  GPM_ATOL = cos(AngularTolerance*DEG2RAD);
1369 
1370  PInit = 1;
1371 
1372 #ifndef HAVE_SRAND48
1373  srand(RandomSeed);
1374 #else
1375  srand48(RandomSeed);
1376 #endif
1377  /* bu_log("Photon Structure Size: %d\n", sizeof(struct PNode));*/
1378 
1379  /*
1380  bu_log("Checking application struct\n");
1381  RT_CK_APPLICATION(ap);
1382  */
1383 
1384  /* Initialize Emitted Photons for each map to 0 */
1385  EPL = 0;
1386  for (i = 0; i < PM_MAPS; i++)
1387  EPS[i] = 0;
1388 
1389  CausticsPercent /= 100.0;
1390  MapSize[PM_IMPORTANCE] = GlobalPhotons/8;
1391  MapSize[PM_GLOBAL] = (int)((1.0-CausticsPercent)*GlobalPhotons);
1392  MapSize[PM_CAUSTIC] = (int)(CausticsPercent*GlobalPhotons);
1393  MapSize[PM_SHADOW] = 0;
1394 
1395  /* bu_log("Caustic Photons: %d\n", MapSize[PM_CAUSTIC]);*/
1396  /* Allocate Memory for Photon Maps */
1397  Initialize(PM_GLOBAL, MapSize[PM_GLOBAL]);
1398  Initialize(PM_CAUSTIC, MapSize[PM_CAUSTIC]);
1399  Initialize(PM_SHADOW, MapSize[PM_SHADOW]);
1400  Initialize(PM_IMPORTANCE, MapSize[PM_IMPORTANCE]);
1401 
1402  /* Populate Application Structure */
1403  /* Set Recursion Level, Magic Number, Hit/Miss Callbacks, and Purpose */
1404  ap->a_level = 1;
1405  ap->a_onehit = 0;
1406  ap->a_ray.magic = RT_RAY_MAGIC;
1407  ap->a_hit = PHit;
1408  ap->a_miss = PMiss;
1410  ap->a_purpose = "Importance Mapping";
1411 
1412 
1413  if (ImportanceMapping) {
1414  bu_log(" Building Importance Map...\n");
1415  EmitImportonsRandom(ap, eye_pos);
1416  BuildTree(Emit[PM_IMPORTANCE], PMap[PM_IMPORTANCE]->StoredPhotons, PMap[PM_IMPORTANCE]->Root);
1417  ScaleFactor = MaxFloat(BBMax[0]-BBMin[0], BBMax[1]-BBMin[1], BBMax[2]-BBMin[2]);
1418  }
1419 
1420  HitG = HitB = 0;
1421  bu_log(" Emitting Photons...\n");
1422  EmitPhotonsRandom(ap, ScaleIndirect);
1423  /* EmitPhotonsRandom(ap, &(LightHead.l), LightIntensity);*/
1424 
1425  /* Generate Scale Factor */
1426  ScaleFactor = MaxFloat(BBMax[0]-BBMin[0], BBMax[1]-BBMin[1], BBMax[2]-BBMin[2]);
1427 
1428  bu_log("HitGB: %d, %d\n", HitG, HitB);
1429  bu_log("Scale Factor: %.3f\n", ScaleFactor);
1430  ratio = (double)HitG/((double)(HitG+HitB));
1431  bu_log("EPL: %d, Adjusted EPL: %d\n", (int)EPL, (int)(EPL*ratio));
1432  EPS[PM_GLOBAL] *= ratio;
1433  EPS[PM_CAUSTIC] *= ratio;
1434 
1435  /* Scale Photon Power */
1436  ScalePhotonPower(PM_GLOBAL);
1437  ScalePhotonPower(PM_CAUSTIC);
1438 
1439  /*
1440  for (i = 0; i < PMap->StoredPhotons; i++)
1441  bu_log("insertLS[%d]: %.3f, %.3f, %.3f\n", i, Emit[i].Pos[0], Emit[i].Pos[1], Emit[i].Pos[2]);
1442  */
1443 
1444 
1445  bu_log(" Building KD-Tree...\n");
1446  /* Balance KD-Tree */
1447  for (i = 0; i < 3; i++)
1448  if (PMap[i]->StoredPhotons)
1449  BuildTree(Emit[i], PMap[i]->StoredPhotons, PMap[i]->Root);
1450 
1451 
1452  bu_log(" Building Irradiance Cache...\n");
1453  ap->a_level = 1;
1454  ap->a_onehit = 0;
1455  ap->a_ray.magic = RT_RAY_MAGIC;
1456  ap->a_hit = ICHit;
1457  ap->a_miss = ICMiss;
1459  ICSize = 0;
1460 
1461  if (cpus > 1) {
1462  memset(GPM_RTAB, 0, sizeof(GPM_RTAB));
1463  for (i = 0; i < MAX_PSW; i++) {
1464  rt_init_resource(&GPM_RTAB[i], i, ap->a_rt_i);
1465  }
1466  bu_parallel(IrradianceThread, cpus, ap);
1467  } else {
1468  /* This will allow profiling for single threaded rendering */
1469  IrradianceThread(0, ap);
1470  }
1471 
1472  /* Allocate Memory for Irradiance Cache and Initialize Pixel Map */
1473  /* bu_log("Image Size: %d, %d\n", width, height);*/
1474  if (GPM_IH) {
1475  Map = (char*)bu_calloc(width*height, sizeof(char), "Map");
1476  IC = (struct IrradCache*)bu_malloc(sizeof(struct IrradCache)*width*height, "IrradCache");
1477  for (i = 0; i < width*height; i++) {
1478  BU_ALLOC(IC[i].List, struct IrradNode);
1479  IC[i].Num = 0;
1480  }
1481  }
1482 
1483 
1484  /*
1485  bu_log(" Sanity Check...\n");
1486  SanityCheck(PMap[PM_GLOBAL]->Root, 0);
1487  */
1488 
1489  WritePhotonFile(pmfile);
1490 
1491  for (i = 0; i < PM_MAPS; i++) {
1492  bu_free(Emit[i], "Photons");
1493  Emit[i] = NULL;
1494  }
1495 
1496  }
1497 }
1498 
1499 
1500 void
1501 Swap(struct PSN *a, struct PSN *b)
1502 {
1503  struct PSN c;
1504 
1505  /*
1506  c.P = a->P;
1507  c.Dist = a->Dist;
1508  a->P = b->P;
1509  a->Dist = b->Dist;
1510  b->P = c.P;
1511  b->Dist = c.Dist;
1512  */
1513  /* bu_log(" SWAP_IN: %.3f, %.3f\n", a->Dist, b->Dist);*/
1514  memcpy(&c, a, sizeof(struct PSN));
1515  memcpy(a, b, sizeof(struct PSN));
1516  memcpy(b, &c, sizeof(struct PSN));
1517  /* bu_log(" SWAP_OT: %.3f, %.3f\n", a->Dist, b->Dist);*/
1518 }
1519 
1520 
1521 /*
1522  After inserting a new node it must be brought upwards until both children
1523  are less than it.
1524 */
1525 void
1526 HeapUp(struct PhotonSearch *S, int ind)
1527 {
1528  int i;
1529 
1530  if (!ind)
1531  return;
1532 
1533  i = ((ind+1)-(ind+1)%2)/2-1;
1534  /* bu_log(" CHECK: %.3f > %.3f :: [%d] > [%d]\n", S->List[ind].Dist, S->List[i].Dist, ind, i);*/
1535  if (S->List[ind].Dist > S->List[i].Dist) {
1536  /* bu_log("SWAP_A: %.3f, %.3f\n", S->List[i].Dist, S->List[ind].Dist);*/
1537  Swap(&S->List[i], &S->List[ind]);
1538  /* bu_log("SWAP_B: %.3f, %.3f\n", S->List[i].Dist, S->List[ind].Dist);*/
1539  }
1540  HeapUp(S, i);
1541 }
1542 
1543 
1544 /*
1545  Sift the new Root node down, by choosing the child with the highest number
1546  since choosing a child with the highest number may reduce the number of
1547  recursions the number will have to propagate
1548 */
1549 void
1550 HeapDown(struct PhotonSearch *S, int ind)
1551 {
1552  int c;
1553 
1554  if (2*ind+1 > S->Found)
1555  return;
1556 
1557  c = 2*ind+1 < S->Found ? S->List[2*ind+2].Dist > S->List[2*ind+1].Dist ? 2*ind+2 : 2*ind+1 : 2*ind+1;
1558  /* bu_log(" c: %d\n", c);*/
1559 
1560  if (S->List[c].Dist > S->List[ind].Dist) {
1561  /* bu_log("SWAP_C: %.3f, %.3f :: %d, %d :: %d\n", S->List[c].Dist, S->List[ind].Dist, c, ind, S->Found);*/
1562  Swap(&S->List[c], &S->List[ind]);
1563  /* bu_log("SWAP_D: %.3f, %.3f :: %d, %d :: %d\n", S->List[c].Dist, S->List[ind].Dist, c, ind, S->Found);*/
1564  }
1565  HeapDown(S, c);
1566 }
1567 
1568 
1569 fastf_t
1570 Dist(point_t a, point_t b)
1571 {
1572  return (a[0]-b[0])*(a[0]-b[0]) + (a[1]-b[1])*(a[1]-b[1]) + (a[2]-b[2])*(a[2]-b[2]);
1573 }
1574 
1575 
1576 void
1577 IrradianceEstimate(struct application *ap, vect_t irrad, point_t pos, vect_t normal)
1578 {
1579  struct PhotonSearch Search;
1580  int i, idx;
1581  fastf_t dist, TotDist;
1582  vect_t t, cirrad;
1583 
1584 
1585  idx = 0;
1586  if (GPM_IH) {
1587  idx = ap->a_x + ap->a_y*GPM_WIDTH;
1588  /* See if there is a cached irradiance calculation for this point */
1589  for (i = 0; i < IC[idx].Num; i++) {
1590  dist = (pos[0]-IC[idx].List[i].Pos[0])*(pos[0]-IC[idx].List[i].Pos[0])+(pos[1]-IC[idx].List[i].Pos[1])*(pos[1]-IC[idx].List[i].Pos[1])+(pos[2]-IC[idx].List[i].Pos[2])*(pos[2]-IC[idx].List[i].Pos[2]);
1591  if (dist < (ScaleFactor/100.0)*(ScaleFactor*100.0)) {
1592  irrad[0] = IC[idx].List[i].RGB[0];
1593  irrad[1] = IC[idx].List[i].RGB[1];
1594  irrad[2] = IC[idx].List[i].RGB[2];
1595  return;
1596  }
1597  }
1598 
1599  /* There is no precomputed irradiance for this point, allocate space
1600  for a new one if necessary. */
1601  if (IC[idx].Num) {
1602  IC[idx].List = (struct IrradNode*)bu_realloc(IC[idx].List, sizeof(struct IrradNode)*(IC[idx].Num+1), "List");
1603  }
1604  }
1605 
1606  Search.Pos[0] = pos[0];
1607  Search.Pos[1] = pos[1];
1608  Search.Pos[2] = pos[2];
1609 
1610  /* NP.RadSq = (ScaleFactor/(10.0*pow(2, (log(PMap[PM_GLOBAL]->MaxPhotons/2)/log(4))))) * (ScaleFactor/(10.0*pow(2, (log(PMap[PM_GLOBAL]->MaxPhotons/2)/log(4)))));*/
1611  /* bu_log("SF: %.3f\n", (ScaleFactor/(10.0*pow(2, (log(PMap[PM_GLOBAL]->MaxPhotons/2)/log(4))))));*/
1612  /* bu_log("SF: %.3f\n", ScaleFactor/pow(PMap[PM_GLOBAL]->MaxPhotons, 0.5));*/
1613 
1614  /*
1615  Search.RadSq = 0.5*ScaleFactor/pow(PMap[PM_GLOBAL]->MaxPhotons, 0.5) * ScaleFactor/pow(PMap[PM_GLOBAL]->MaxPhotons, 0.5);
1616  Search.Max = pow(PMap[PM_GLOBAL]->StoredPhotons, 0.5);
1617  */
1618 
1619  Search.RadSq = (ScaleFactor/2048.0);
1620  Search.RadSq *= Search.RadSq;
1621  /* NP.RadSq = (4.0*ScaleFactor/PMap[PM_GLOBAL]->MaxPhotons) * (4.0*ScaleFactor/PMap[PM_GLOBAL]->MaxPhotons);*/
1622  /* NP.Max = 2.0*pow(PMap[PM_GLOBAL]->StoredPhotons, 0.5);*/
1623  /* Search.Max = PMap[PM_GLOBAL]->StoredPhotons / 50;*/
1624  Search.Max = 32;
1625 
1626  Search.Normal[0] = normal[0];
1627  Search.Normal[1] = normal[1];
1628  Search.Normal[2] = normal[2];
1629 
1630  Search.List = (struct PSN*)bu_calloc(Search.Max, sizeof(struct PSN), "Search.List");
1631  do {
1632  Search.Found = 0;
1633  Search.RadSq *= 4.0;
1634  LocatePhotons(&Search, PMap[PM_GLOBAL]->Root);
1635  } while (Search.Found < Search.Max && Search.RadSq < ScaleFactor * ScaleFactor / 64.0);
1636 
1637 
1638  irrad[0] = irrad[1] = irrad[2] = 0;
1639  TotDist = 0;
1640  for (i = 0; i < Search.Found; i++) {
1641  t[0] = Search.List[i].P.Pos[0] - pos[0];
1642  t[1] = Search.List[i].P.Pos[1] - pos[1];
1643  t[2] = Search.List[i].P.Pos[2] - pos[2];
1644  TotDist += t[0]*t[0] + t[1]*t[1] + t[2]*t[2];
1645  }
1646 
1647 
1648  for (i = 0; i < Search.Found; i++) {
1649  t[0] = Search.List[i].P.Pos[0] - pos[0];
1650  t[1] = Search.List[i].P.Pos[1] - pos[1];
1651  t[2] = Search.List[i].P.Pos[2] - pos[2];
1652 
1653  t[0] = (t[0]*t[0] + t[1]*t[1] + t[2]*t[2])/TotDist;
1654  /*
1655  irrad[0] += Search.List[i].P.Irrad[0] * t[0];
1656  irrad[1] += Search.List[i].P.Irrad[1] * t[0];
1657  irrad[2] += Search.List[i].P.Irrad[2] * t[0];
1658  */
1659  irrad[0] += Search.List[i].P.Irrad[0];
1660  irrad[1] += Search.List[i].P.Irrad[1];
1661  irrad[2] += Search.List[i].P.Irrad[2];
1662  }
1663  if (Search.Found) {
1664  irrad[0] /= (double)Search.Found;
1665  irrad[1] /= (double)Search.Found;
1666  irrad[2] /= (double)Search.Found;
1667  }
1668  bu_free(Search.List, "Search.List");
1669 
1670  /* GetEstimate(cirrad, pos, normal, (int)(ScaleFactor/100.0), PMap[PM_CAUSTIC]->MaxPhotons/50, PM_CAUSTIC, 1, 0);*/
1671  /* GetEstimate(cirrad, pos, normal, (int)(ScaleFactor/pow(2, (log(PMap[PM_CAUSTIC]->MaxPhotons/2)/log(4)))), PMap[PM_CAUSTIC]->MaxPhotons / 50, PM_CAUSTIC, 0, 0);*/
1672  GetEstimate(cirrad, pos, normal, ScaleFactor/1024.0, PMap[PM_CAUSTIC]->MaxPhotons / 100, PM_CAUSTIC, ScaleFactor/128.0, 1, 15);
1673 
1674  irrad[0] += cirrad[0];
1675  irrad[1] += cirrad[1];
1676  irrad[2] += cirrad[2];
1677 
1678 
1679  /* Visualize Green Importons */
1680  /*
1681  GetEstimate(cirrad, pos, normal, ScaleFactor/512.0, 1, PM_IMPORTANCE, ScaleFactor/256.0, 1, 1);
1682 
1683  irrad[0] += cirrad[0];
1684  irrad[1] += cirrad[1];
1685  irrad[2] += cirrad[2];
1686  */
1687 
1688  if (GPM_IH) {
1689  /* Store Irradiance */
1690  IC[idx].List[IC[idx].Num].RGB[0] = irrad[0];
1691  IC[idx].List[IC[idx].Num].RGB[1] = irrad[1];
1692  IC[idx].List[IC[idx].Num].RGB[2] = irrad[2];
1693 
1694  IC[idx].List[IC[idx].Num].Pos[0] = pos[0];
1695  IC[idx].List[IC[idx].Num].Pos[1] = pos[1];
1696  IC[idx].List[IC[idx].Num].Pos[2] = pos[2];
1697 
1698  IC[idx].Num++;
1699  }
1700 }
1701 
1702 
1703 /*
1704  * Local Variables:
1705  * mode: C
1706  * tab-width: 8
1707  * indent-tabs-mode: t
1708  * c-file-style: "stroustrup"
1709  * End:
1710  * ex: shiftwidth=4 tabstop=8
1711  */
struct xray a_ray
Actual ray to be shot.
Definition: raytrace.h:1583
#define BU_LIST_FOR(p, structure, hp)
Definition: list.h:365
fastf_t GaussFilter(fastf_t dist, fastf_t rad)
Definition: photonmap.c:826
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
struct region * pt_regionp
ptr to containing region
Definition: raytrace.h:580
struct hit * pt_outhit
OUT hit ptr.
Definition: raytrace.h:579
int HitB
Definition: photonmap.c:66
void WritePhotonFile(char *pmfile)
Definition: photonmap.c:1279
void BuildPhotonMap(struct application *ap, point_t eye_pos, int cpus, int width, int height, int Hypersample, int GlobalPhotons, double CausticsPercent, int Rays, double AngularTolerance, int RandomSeed, int ImportanceMapping, int IrradianceHypersampling, int VisualizeIrradiance, double ScaleIndirect, char pmfile[255])
Definition: photonmap.c:1347
Definition: list.h:118
#define RT_RAY_MAGIC
Definition: magic.h:165
void rt_silent_logoverlap(struct application *ap, const struct partition *pp, const struct bu_ptbl *regiontable, const struct partition *InputHdp)
Definition: bool.c:1122
char * Map
Definition: photonmap.c:59
struct soltab * seg_stp
pointer back to soltab
Definition: raytrace.h:372
int Refract(vect_t I, vect_t N, fastf_t n1, fastf_t n2)
Definition: photonmap.c:327
fastf_t Dist(point_t a, point_t b)
Definition: photonmap.c:1570
#define N
Definition: randmt.c:39
struct bu_structparse phong_parse[]
Definition: sh_plastic.c:54
void bu_semaphore_acquire(unsigned int i)
Definition: semaphore.c:180
#define M_PI
Definition: fft.h:35
Definition: raytrace.h:368
vect_t BBMax
Definition: photonmap.c:50
char pt_inflip
flip inhit->hit_normal
Definition: raytrace.h:581
void Initialize(int MAP, int MapSize)
Definition: photonmap.c:1147
Header file for the BRL-CAD common definitions.
long time(time_t *)
struct resource * a_resource
dynamic memory resources
Definition: raytrace.h:1591
int PType
Definition: photonmap.c:52
int a_onehit
flag to stop on first hit
Definition: raytrace.h:1586
void IrradianceThread(int pid, void *arg)
Definition: photonmap.c:1131
ustring width
int PHit(struct application *ap, struct partition *PartHeadp, struct seg *finished_segs)
Definition: photonmap.c:501
int(* a_hit)(struct application *, struct partition *, struct seg *)
called when shot hits model
Definition: raytrace.h:1584
struct hit * pt_inhit
IN hit pointer.
Definition: raytrace.h:577
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
#define Max(a, b)
Definition: multipoly.c:43
void rt_init_resource(struct resource *resp, int cpu_num, struct rt_i *rtip)
Definition: prep.c:585
void bu_vls_free(struct bu_vls *vp)
Definition: vls.c:248
fastf_t ConeFilter(fastf_t dist, fastf_t rad)
Definition: photonmap.c:833
void Irradiance(int pid, struct Photon *P, struct application *ap)
Definition: photonmap.c:1019
struct Photon CurPh
Definition: photonmap.c:48
struct rt_i * a_rt_i
this librt instance
Definition: raytrace.h:1588
void * memset(void *s, int c, size_t n)
int PM_Activated
Definition: photonmap.c:43
void Store(point_t Pos, vect_t Dir, vect_t Normal, int map)
Definition: photonmap.c:228
void EmitImportonsRandom(struct application *ap, point_t eye_pos)
Definition: photonmap.c:717
void BuildIrradianceCache(int pid, struct PNode *Node, struct application *ap)
Definition: photonmap.c:1078
void GetMaterial(char *MS, vect_t spec, fastf_t *refi, fastf_t *transmit)
Definition: photonmap.c:367
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
void * bu_calloc(size_t nelem, size_t elsize, const char *str)
Definition: malloc.c:321
int FindMedian(struct Photon *List, int Num, int Axis)
Definition: photonmap.c:71
int EPL
Definition: photonmap.c:54
int GPM_IH
Definition: photonmap.c:60
int a_x
Screen X of ray, if applicable.
Definition: raytrace.h:1596
#define BU_GET(_ptr, _type)
Definition: malloc.h:201
int ICMiss(struct application *ap)
Definition: photonmap.c:985
struct seg * pt_inseg
IN seg ptr (gives stp)
Definition: raytrace.h:576
int Depth
Definition: photonmap.c:51
int CheckMaterial(char *cmp, char *MS)
Definition: photonmap.c:350
char * ma_shader
shader name & parms
Definition: raytrace.h:527
int a_level
recursion level (for printing)
Definition: raytrace.h:1595
void * bu_realloc(void *ptr, size_t siz, const char *str)
const char * a_purpose
Debug string: purpose of ray.
Definition: raytrace.h:1598
#define Min(a, b)
Definition: multipoly.c:44
#define UNUSED(parameter)
Definition: common.h:239
#define BU_PUT(_ptr, _type)
Definition: malloc.h:215
void Polar2Euclidian(vect_t Dir, vect_t Normal, double Theta, double Phi)
Definition: photonmap.c:998
void HeapDown(struct PhotonSearch *S, int ind)
Definition: photonmap.c:1550
int GPM_HEIGHT
Definition: photonmap.c:62
double drand48(void)
Definition: heap.c:81
struct resource GPM_RTAB[MAX_PSW]
Definition: photonmap.c:65
void bu_semaphore_release(unsigned int i)
Definition: semaphore.c:218
vect_t r_dir
Direction of ray (UNIT Length)
Definition: raytrace.h:219
#define MAX_PSW
Definition: parallel.h:48
void WritePhotons(struct PNode *Root, FILE *FH)
Definition: photonmap.c:1264
void bu_parallel(void(*func)(int func_ncpu, void *func_data), int ncpu, void *data)
void GetEstimate(vect_t irrad, point_t pos, vect_t normal, fastf_t rad, int np, int map, double max_rad, int centog, int min_np)
Definition: photonmap.c:840
float ma_color[3]
explicit color: 0..1
Definition: raytrace.h:522
struct mater_info reg_mater
Real material information.
Definition: raytrace.h:546
int PInit
Definition: photonmap.c:53
void SpecularReflect(vect_t normal, vect_t rdir)
Definition: photonmap.c:288
void(* a_logoverlap)(struct application *, const struct partition *, const struct bu_ptbl *, const struct partition *)
called to log overlaps
Definition: raytrace.h:1594
point_t r_pt
Point at which ray starts.
Definition: raytrace.h:218
double GPM_ATOL
Definition: photonmap.c:64
int EPS[PM_MAPS]
Definition: photonmap.c:55
int hit(struct application *ap, struct partition *PartHeadp, struct seg *segs)
Definition: gqa.c:963
void SanityCheck(struct PNode *Root, int LR)
Definition: photonmap.c:814
int ICHit(struct application *ap, struct partition *PartHeadp, struct seg *finished_segs)
Definition: photonmap.c:958
void ScalePhotonPower(int map)
Definition: photonmap.c:703
int GPM_RAYS
Definition: photonmap.c:63
int bu_struct_parse(const struct bu_vls *in_vls, const struct bu_structparse *desc, const char *base, void *data)
Definition: parse.c:878
int a_y
Screen Y of ray, if applicable.
Definition: raytrace.h:1597
void bu_vls_printf(struct bu_vls *vls, const char *fmt,...) _BU_ATTR_PRINTF23
Definition: vls.c:694
void LocatePhotons(struct PhotonSearch *Search, struct PNode *Root)
Definition: photonmap.c:176
void Swap(struct PSN *a, struct PSN *b)
Definition: photonmap.c:1501
int LoadFile(char *pmfile)
Definition: photonmap.c:1162
#define PL_MAGIC
Definition: magic.h:207
int rt_shootray(struct application *ap)
void HeapUp(struct PhotonSearch *S, int ind)
Definition: photonmap.c:1526
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
int(* a_miss)(struct application *)
called when shot misses
Definition: raytrace.h:1585
double ScaleFactor
Definition: photonmap.c:57
#define BU_VLS_INIT_ZERO
Definition: vls.h:84
#define RT_HIT_NORMAL(_normal, _hitp, _stp, _unused, _flipflag)
Definition: raytrace.h:273
vect_t BBMin
Definition: photonmap.c:49
struct light_specific LightHead
Definition: sh_light.c:50
void DiffuseReflect(vect_t normal, vect_t rdir)
Definition: photonmap.c:307
struct IrradCache * IC
Definition: photonmap.c:58
int PMiss(struct application *ap)
Definition: photonmap.c:691
struct partition * pt_forw
forwards link
Definition: raytrace.h:574
int PM_Visualize
Definition: photonmap.c:44
#define RT_APPLICATION_INIT(_p)
Definition: raytrace.h:1676
fastf_t hit_dist
dist from r_pt to hit_point
Definition: raytrace.h:250
Definition: vls.h:56
int HitG
Definition: photonmap.c:66
void bu_bomb(const char *str) _BU_ATTR_NORETURN
Definition: bomb.c:91
int HitRef(struct application *ap, struct partition *PartHeadp, struct seg *finished_segs)
Definition: photonmap.c:455
double fastf_t
Definition: defines.h:300
void BuildTree(struct Photon *EList, int ESize, struct PNode *Root)
Definition: photonmap.c:94
int GPM_WIDTH
Definition: photonmap.c:61
void IrradianceEstimate(struct application *ap, vect_t irrad, point_t pos, vect_t normal)
Definition: photonmap.c:1577
int ICSize
Definition: photonmap.c:56
uint32_t magic
Definition: raytrace.h:216
struct PhotonMap * PMap[PM_MAPS]
Definition: photonmap.c:46
void EmitPhotonsRandom(struct application *ap, double ScaleIndirect)
Definition: photonmap.c:755
#define M
Definition: msr.c:52
struct Photon * Emit[PM_MAPS]
Definition: photonmap.c:47