| 1 | 
greg | 
2.1 | 
/*  | 
| 2 | 
  | 
  | 
   ================================================================== | 
| 3 | 
  | 
  | 
   Photon map main module | 
| 4 | 
  | 
  | 
 | 
| 5 | 
  | 
  | 
   Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) | 
| 6 | 
  | 
  | 
   (c) Fraunhofer Institute for Solar Energy Systems, | 
| 7 | 
rschregle | 
2.4 | 
   (c) Lucerne University of Applied Sciences and Arts, | 
| 8 | 
  | 
  | 
   supported by the Swiss National Science Foundation (SNSF, #147053) | 
| 9 | 
greg | 
2.1 | 
   ================================================================== | 
| 10 | 
  | 
  | 
    | 
| 11 | 
rschregle | 
2.4 | 
   $Id: pmap.c,v 2.3 2015/04/22 20:28:16 rschregle Exp $ | 
| 12 | 
greg | 
2.1 | 
*/ | 
| 13 | 
  | 
  | 
 | 
| 14 | 
  | 
  | 
 | 
| 15 | 
  | 
  | 
 | 
| 16 | 
  | 
  | 
#include "pmap.h" | 
| 17 | 
  | 
  | 
#include "pmapmat.h" | 
| 18 | 
  | 
  | 
#include "pmapsrc.h" | 
| 19 | 
  | 
  | 
#include "pmaprand.h" | 
| 20 | 
  | 
  | 
#include "pmapio.h" | 
| 21 | 
  | 
  | 
#include "pmapbias.h" | 
| 22 | 
  | 
  | 
#include "pmapdiag.h" | 
| 23 | 
  | 
  | 
#include "otypes.h" | 
| 24 | 
  | 
  | 
#include <time.h> | 
| 25 | 
  | 
  | 
#include <sys/stat.h> | 
| 26 | 
  | 
  | 
 | 
| 27 | 
  | 
  | 
 | 
| 28 | 
  | 
  | 
 | 
| 29 | 
  | 
  | 
extern char *octname; | 
| 30 | 
  | 
  | 
 | 
| 31 | 
rschregle | 
2.4 | 
static char PmapRevision [] = "$Revision: 2.3 $"; | 
| 32 | 
greg | 
2.1 | 
 | 
| 33 | 
  | 
  | 
 | 
| 34 | 
  | 
  | 
 | 
| 35 | 
  | 
  | 
/* Photon map lookup functions per type */ | 
| 36 | 
  | 
  | 
void (*pmapLookup [NUM_PMAP_TYPES])(PhotonMap*, RAY*, COLOR) = { | 
| 37 | 
  | 
  | 
   photonDensity, photonPreCompDensity, photonDensity, volumePhotonDensity, | 
| 38 | 
  | 
  | 
   photonDensity, NULL | 
| 39 | 
  | 
  | 
}; | 
| 40 | 
  | 
  | 
 | 
| 41 | 
  | 
  | 
 | 
| 42 | 
  | 
  | 
 | 
| 43 | 
  | 
  | 
void colorNorm (COLOR c) | 
| 44 | 
  | 
  | 
/* Normalise colour channels to average of 1 */ | 
| 45 | 
  | 
  | 
{ | 
| 46 | 
  | 
  | 
   const float avg = colorAvg(c); | 
| 47 | 
  | 
  | 
    | 
| 48 | 
  | 
  | 
   if (!avg)  | 
| 49 | 
  | 
  | 
      return; | 
| 50 | 
  | 
  | 
       | 
| 51 | 
  | 
  | 
   c [0] /= avg; | 
| 52 | 
  | 
  | 
   c [1] /= avg; | 
| 53 | 
  | 
  | 
   c [2] /= avg; | 
| 54 | 
  | 
  | 
} | 
| 55 | 
  | 
  | 
 | 
| 56 | 
  | 
  | 
 | 
| 57 | 
  | 
  | 
 | 
| 58 | 
  | 
  | 
void loadPmaps (PhotonMap **pmaps, const PhotonMapParams *parm) | 
| 59 | 
  | 
  | 
{ | 
| 60 | 
  | 
  | 
   unsigned t; | 
| 61 | 
  | 
  | 
   struct stat octstat, pmstat; | 
| 62 | 
  | 
  | 
   PhotonMap *pm; | 
| 63 | 
  | 
  | 
   PhotonMapType type; | 
| 64 | 
  | 
  | 
    | 
| 65 | 
  | 
  | 
   for (t = 0; t < NUM_PMAP_TYPES; t++) | 
| 66 | 
  | 
  | 
      if (setPmapParam(&pm, parm + t)) {          | 
| 67 | 
  | 
  | 
         /* Check if photon map newer than octree */ | 
| 68 | 
rschregle | 
2.4 | 
         if (pm -> fileName && octname && | 
| 69 | 
  | 
  | 
             !stat(pm -> fileName, &pmstat) && !stat(octname, &octstat) &&  | 
| 70 | 
greg | 
2.1 | 
             octstat.st_mtime > pmstat.st_mtime) { | 
| 71 | 
  | 
  | 
            sprintf(errmsg, "photon map in file %s may be stale",  | 
| 72 | 
  | 
  | 
                    pm -> fileName); | 
| 73 | 
  | 
  | 
            error(USER, errmsg); | 
| 74 | 
  | 
  | 
         } | 
| 75 | 
  | 
  | 
          | 
| 76 | 
  | 
  | 
         /* Load photon map from file and get its type */ | 
| 77 | 
  | 
  | 
         if ((type = loadPhotonMap(pm, pm -> fileName)) == PMAP_TYPE_NONE) | 
| 78 | 
  | 
  | 
            error(USER, "failed loading photon map"); | 
| 79 | 
  | 
  | 
             | 
| 80 | 
  | 
  | 
         /* Assign to appropriate photon map type (deleting previously | 
| 81 | 
  | 
  | 
          * loaded photon map of same type if necessary) */ | 
| 82 | 
  | 
  | 
         if (pmaps [type]) { | 
| 83 | 
  | 
  | 
            deletePhotons(pmaps [type]); | 
| 84 | 
  | 
  | 
            free(pmaps [type]); | 
| 85 | 
  | 
  | 
         } | 
| 86 | 
  | 
  | 
         pmaps [type] = pm; | 
| 87 | 
  | 
  | 
          | 
| 88 | 
  | 
  | 
         /* Check for invalid density estimate bandwidth */                             | 
| 89 | 
  | 
  | 
         if (pm -> maxGather > pm -> heapSize) { | 
| 90 | 
  | 
  | 
            error(WARNING, "adjusting density estimate bandwidth"); | 
| 91 | 
  | 
  | 
            pm -> minGather = pm -> maxGather = pm -> heapSize; | 
| 92 | 
  | 
  | 
         } | 
| 93 | 
  | 
  | 
      } | 
| 94 | 
  | 
  | 
} | 
| 95 | 
  | 
  | 
 | 
| 96 | 
  | 
  | 
 | 
| 97 | 
  | 
  | 
 | 
| 98 | 
  | 
  | 
void savePmaps (const PhotonMap **pmaps, int argc, char **argv) | 
| 99 | 
  | 
  | 
{ | 
| 100 | 
  | 
  | 
   unsigned t; | 
| 101 | 
  | 
  | 
    | 
| 102 | 
  | 
  | 
   for (t = 0; t < NUM_PMAP_TYPES; t++) { | 
| 103 | 
  | 
  | 
      if (pmaps [t]) | 
| 104 | 
  | 
  | 
         savePhotonMap(pmaps [t], pmaps [t] -> fileName, t, argc, argv); | 
| 105 | 
  | 
  | 
   } | 
| 106 | 
  | 
  | 
}                    | 
| 107 | 
  | 
  | 
 | 
| 108 | 
  | 
  | 
 | 
| 109 | 
  | 
  | 
    | 
| 110 | 
  | 
  | 
void cleanUpPmaps (PhotonMap **pmaps) | 
| 111 | 
  | 
  | 
{ | 
| 112 | 
  | 
  | 
   unsigned t; | 
| 113 | 
  | 
  | 
    | 
| 114 | 
  | 
  | 
   for (t = 0; t < NUM_PMAP_TYPES; t++) { | 
| 115 | 
  | 
  | 
      if (pmaps [t]) { | 
| 116 | 
  | 
  | 
         deletePhotons(pmaps [t]); | 
| 117 | 
  | 
  | 
         free(pmaps [t]); | 
| 118 | 
  | 
  | 
      } | 
| 119 | 
  | 
  | 
   } | 
| 120 | 
  | 
  | 
} | 
| 121 | 
  | 
  | 
 | 
| 122 | 
  | 
  | 
 | 
| 123 | 
  | 
  | 
      | 
| 124 | 
  | 
  | 
static int photonParticipate (RAY *ray) | 
| 125 | 
  | 
  | 
/* Trace photon through participating medium. Returns 1 if passed through, | 
| 126 | 
  | 
  | 
   or 0 if absorbed and $*%&ed. Analogon to rayparticipate(). */ | 
| 127 | 
  | 
  | 
{ | 
| 128 | 
  | 
  | 
   int i; | 
| 129 | 
  | 
  | 
   RREAL cosTheta, cosPhi, du, dv; | 
| 130 | 
  | 
  | 
   const float cext = colorAvg(ray -> cext), | 
| 131 | 
  | 
  | 
               albedo = colorAvg(ray -> albedo); | 
| 132 | 
  | 
  | 
   FVECT u, v; | 
| 133 | 
  | 
  | 
   COLOR cvext; | 
| 134 | 
  | 
  | 
 | 
| 135 | 
  | 
  | 
   /* Mean free distance until interaction with medium */ | 
| 136 | 
  | 
  | 
   ray -> rmax = -log(pmapRandom(mediumState)) / cext; | 
| 137 | 
  | 
  | 
    | 
| 138 | 
  | 
  | 
   while (!localhit(ray, &thescene)) { | 
| 139 | 
  | 
  | 
      setcolor(cvext, exp(-ray -> rmax * ray -> cext [0]), | 
| 140 | 
  | 
  | 
                      exp(-ray -> rmax * ray -> cext [1]), | 
| 141 | 
  | 
  | 
                      exp(-ray -> rmax * ray -> cext [2])); | 
| 142 | 
  | 
  | 
                       | 
| 143 | 
  | 
  | 
      /* Modify ray color and normalise */ | 
| 144 | 
  | 
  | 
      multcolor(ray -> rcol, cvext); | 
| 145 | 
  | 
  | 
      colorNorm(ray -> rcol); | 
| 146 | 
  | 
  | 
      VCOPY(ray -> rorg, ray -> rop); | 
| 147 | 
  | 
  | 
       | 
| 148 | 
  | 
  | 
      if (albedo > FTINY) | 
| 149 | 
  | 
  | 
         /* Add to volume photon map */ | 
| 150 | 
  | 
  | 
         if (ray -> rlvl > 0) addPhoton(volumePmap, ray); | 
| 151 | 
  | 
  | 
          | 
| 152 | 
  | 
  | 
      /* Absorbed? */ | 
| 153 | 
  | 
  | 
      if (pmapRandom(rouletteState) > albedo) return 0; | 
| 154 | 
  | 
  | 
       | 
| 155 | 
  | 
  | 
      /* Colour bleeding without attenuation (?) */ | 
| 156 | 
  | 
  | 
      multcolor(ray -> rcol, ray -> albedo); | 
| 157 | 
  | 
  | 
      scalecolor(ray -> rcol, 1 / albedo);     | 
| 158 | 
  | 
  | 
       | 
| 159 | 
  | 
  | 
      /* Scatter photon */ | 
| 160 | 
  | 
  | 
      cosTheta = ray -> gecc <= FTINY ? 2 * pmapRandom(scatterState) - 1 | 
| 161 | 
  | 
  | 
                                      : 1 / (2 * ray -> gecc) *  | 
| 162 | 
  | 
  | 
                                            (1 + ray -> gecc * ray -> gecc -  | 
| 163 | 
  | 
  | 
                                               (1 - ray -> gecc * ray -> gecc) /  | 
| 164 | 
  | 
  | 
                                               (1 - ray -> gecc + 2 * ray -> gecc * | 
| 165 | 
  | 
  | 
                                                  pmapRandom(scatterState))); | 
| 166 | 
  | 
  | 
                                                   | 
| 167 | 
  | 
  | 
      cosPhi = cos(2 * PI * pmapRandom(scatterState)); | 
| 168 | 
  | 
  | 
      du = dv = sqrt(1 - cosTheta * cosTheta);   /* sin(theta) */ | 
| 169 | 
  | 
  | 
      du *= cosPhi; | 
| 170 | 
  | 
  | 
      dv *= sqrt(1 - cosPhi * cosPhi);           /* sin(phi) */ | 
| 171 | 
  | 
  | 
       | 
| 172 | 
  | 
  | 
      /* Get axes u & v perpendicular to photon direction */ | 
| 173 | 
  | 
  | 
      i = 0; | 
| 174 | 
  | 
  | 
      do { | 
| 175 | 
  | 
  | 
         v [0] = v [1] = v [2] = 0; | 
| 176 | 
  | 
  | 
         v [i++] = 1; | 
| 177 | 
  | 
  | 
         fcross(u, v, ray -> rdir); | 
| 178 | 
  | 
  | 
      } while (normalize(u) < FTINY); | 
| 179 | 
  | 
  | 
      fcross(v, ray -> rdir, u); | 
| 180 | 
  | 
  | 
       | 
| 181 | 
  | 
  | 
      for (i = 0; i < 3; i++) | 
| 182 | 
  | 
  | 
         ray -> rdir [i] = du * u [i] + dv * v [i] +  | 
| 183 | 
  | 
  | 
                           cosTheta * ray -> rdir [i]; | 
| 184 | 
  | 
  | 
      ray -> rlvl++; | 
| 185 | 
  | 
  | 
      ray -> rmax = -log(pmapRandom(mediumState)) / cext; | 
| 186 | 
  | 
  | 
   }   | 
| 187 | 
  | 
  | 
     | 
| 188 | 
  | 
  | 
   setcolor(cvext, exp(-ray -> rot * ray -> cext [0]), | 
| 189 | 
  | 
  | 
                   exp(-ray -> rot * ray -> cext [1]), | 
| 190 | 
  | 
  | 
                   exp(-ray -> rot * ray -> cext [2])); | 
| 191 | 
  | 
  | 
                    | 
| 192 | 
  | 
  | 
   /* Modify ray color and normalise */ | 
| 193 | 
  | 
  | 
   multcolor(ray -> rcol, cvext); | 
| 194 | 
  | 
  | 
   colorNorm(ray -> rcol); | 
| 195 | 
  | 
  | 
    | 
| 196 | 
  | 
  | 
   /* Passed through medium */   | 
| 197 | 
  | 
  | 
   return 1; | 
| 198 | 
  | 
  | 
} | 
| 199 | 
  | 
  | 
 | 
| 200 | 
  | 
  | 
 | 
| 201 | 
  | 
  | 
 | 
| 202 | 
  | 
  | 
void tracePhoton (RAY *ray) | 
| 203 | 
  | 
  | 
/* Follow photon as it bounces around... */ | 
| 204 | 
  | 
  | 
{ | 
| 205 | 
  | 
  | 
   long mod; | 
| 206 | 
  | 
  | 
   OBJREC* mat; | 
| 207 | 
  | 
  | 
  | 
| 208 | 
  | 
  | 
   if (ray -> rlvl > photonMaxBounce) { | 
| 209 | 
  | 
  | 
      error(WARNING, "runaway photon!"); | 
| 210 | 
  | 
  | 
      return; | 
| 211 | 
  | 
  | 
   } | 
| 212 | 
  | 
  | 
    | 
| 213 | 
  | 
  | 
   if (colorAvg(ray -> cext) > FTINY && !photonParticipate(ray))  | 
| 214 | 
  | 
  | 
      return; | 
| 215 | 
  | 
  | 
       | 
| 216 | 
  | 
  | 
   if (localhit(ray, &thescene)) { | 
| 217 | 
  | 
  | 
      mod = ray -> ro -> omod; | 
| 218 | 
  | 
  | 
       | 
| 219 | 
  | 
  | 
      if ((ray -> clipset && inset(ray -> clipset, mod)) || mod == OVOID) { | 
| 220 | 
  | 
  | 
         /* Transfer ray if modifier is VOID or clipped within antimatta */ | 
| 221 | 
  | 
  | 
         RAY tray; | 
| 222 | 
  | 
  | 
         photonRay(ray, &tray, PMAP_XFER, NULL); | 
| 223 | 
  | 
  | 
         tracePhoton(&tray); | 
| 224 | 
  | 
  | 
      } | 
| 225 | 
  | 
  | 
      else { | 
| 226 | 
  | 
  | 
         /* Scatter for modifier material */ | 
| 227 | 
  | 
  | 
         mat = objptr(mod); | 
| 228 | 
  | 
  | 
         photonScatter [mat -> otype] (mat, ray); | 
| 229 | 
  | 
  | 
      } | 
| 230 | 
  | 
  | 
   } | 
| 231 | 
  | 
  | 
} | 
| 232 | 
  | 
  | 
 | 
| 233 | 
  | 
  | 
 | 
| 234 | 
  | 
  | 
 | 
| 235 | 
  | 
  | 
static void preComputeGlobal (PhotonMap *pmap) | 
| 236 | 
  | 
  | 
/* Precompute irradiance from global photons for final gathering using  | 
| 237 | 
  | 
  | 
   the first finalGather * pmap -> heapSize photons in the heap. Returns  | 
| 238 | 
  | 
  | 
   new heap with precomputed photons. */ | 
| 239 | 
  | 
  | 
{ | 
| 240 | 
  | 
  | 
   unsigned long i, nuHeapSize; | 
| 241 | 
  | 
  | 
   unsigned j; | 
| 242 | 
  | 
  | 
   Photon *nuHeap, *p; | 
| 243 | 
  | 
  | 
   COLOR irrad; | 
| 244 | 
  | 
  | 
   RAY ray; | 
| 245 | 
  | 
  | 
   float nuMinPos [3], nuMaxPos [3]; | 
| 246 | 
  | 
  | 
 | 
| 247 | 
  | 
  | 
   repComplete = nuHeapSize = finalGather * pmap -> heapSize; | 
| 248 | 
  | 
  | 
    | 
| 249 | 
  | 
  | 
   if (photonRepTime) { | 
| 250 | 
  | 
  | 
      sprintf(errmsg,  | 
| 251 | 
  | 
  | 
              "Precomputing irradiance for %ld global photons...\n",  | 
| 252 | 
  | 
  | 
              nuHeapSize); | 
| 253 | 
  | 
  | 
      eputs(errmsg); | 
| 254 | 
  | 
  | 
      fflush(stderr); | 
| 255 | 
  | 
  | 
   } | 
| 256 | 
  | 
  | 
    | 
| 257 | 
  | 
  | 
   p = nuHeap = (Photon*)malloc(nuHeapSize * sizeof(Photon)); | 
| 258 | 
  | 
  | 
   if (!nuHeap)  | 
| 259 | 
  | 
  | 
      error(USER, "can't allocate photon heap"); | 
| 260 | 
  | 
  | 
       | 
| 261 | 
  | 
  | 
   for (j = 0; j <= 2; j++) { | 
| 262 | 
  | 
  | 
      nuMinPos [j] = FHUGE; | 
| 263 | 
  | 
  | 
      nuMaxPos [j] = -FHUGE; | 
| 264 | 
  | 
  | 
   } | 
| 265 | 
  | 
  | 
    | 
| 266 | 
  | 
  | 
   /* Record start time, baby */ | 
| 267 | 
  | 
  | 
   repStartTime = time(NULL); | 
| 268 | 
rschregle | 
2.3 | 
   #ifdef SIGCONT | 
| 269 | 
  | 
  | 
      signal(SIGCONT, pmapPreCompReport); | 
| 270 | 
  | 
  | 
   #endif | 
| 271 | 
greg | 
2.1 | 
   repProgress = 0; | 
| 272 | 
greg | 
2.2 | 
   memcpy(nuHeap, pmap -> heap, nuHeapSize * sizeof(Photon)); | 
| 273 | 
greg | 
2.1 | 
    | 
| 274 | 
  | 
  | 
   for (i = 0, p = nuHeap; i < nuHeapSize; i++, p++) { | 
| 275 | 
  | 
  | 
      ray.ro = NULL; | 
| 276 | 
  | 
  | 
      VCOPY(ray.rop, p -> pos); | 
| 277 | 
  | 
  | 
       | 
| 278 | 
  | 
  | 
      /* Update min and max positions & set ray normal */ | 
| 279 | 
  | 
  | 
      for (j = 0; j < 3; j++) { | 
| 280 | 
  | 
  | 
         if (p -> pos [j] < nuMinPos [j]) nuMinPos [j] = p -> pos [j]; | 
| 281 | 
  | 
  | 
         if (p -> pos [j] > nuMaxPos [j]) nuMaxPos [j] = p -> pos [j]; | 
| 282 | 
  | 
  | 
         ray.ron [j] = p -> norm [j] / 127.0; | 
| 283 | 
  | 
  | 
      } | 
| 284 | 
  | 
  | 
       | 
| 285 | 
  | 
  | 
      photonDensity(pmap, &ray, irrad); | 
| 286 | 
  | 
  | 
      setPhotonFlux(p, irrad); | 
| 287 | 
  | 
  | 
      repProgress++; | 
| 288 | 
  | 
  | 
       | 
| 289 | 
  | 
  | 
      if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) | 
| 290 | 
  | 
  | 
         pmapPreCompReport(); | 
| 291 | 
rschregle | 
2.3 | 
      #ifdef SIGCONT | 
| 292 | 
  | 
  | 
         else signal(SIGCONT, pmapPreCompReport); | 
| 293 | 
  | 
  | 
      #endif | 
| 294 | 
greg | 
2.1 | 
   } | 
| 295 | 
  | 
  | 
    | 
| 296 | 
rschregle | 
2.3 | 
   #ifdef SIGCONT    | 
| 297 | 
  | 
  | 
      signal(SIGCONT, SIG_DFL); | 
| 298 | 
  | 
  | 
   #endif | 
| 299 | 
greg | 
2.1 | 
    | 
| 300 | 
  | 
  | 
   /* Replace & rebuild heap */ | 
| 301 | 
  | 
  | 
   free(pmap -> heap); | 
| 302 | 
  | 
  | 
   pmap -> heap = nuHeap; | 
| 303 | 
  | 
  | 
   pmap -> heapSize = pmap -> heapEnd = nuHeapSize; | 
| 304 | 
  | 
  | 
   VCOPY(pmap -> minPos, nuMinPos); | 
| 305 | 
  | 
  | 
   VCOPY(pmap -> maxPos, nuMaxPos); | 
| 306 | 
  | 
  | 
    | 
| 307 | 
  | 
  | 
   if (photonRepTime) { | 
| 308 | 
  | 
  | 
      eputs("Rebuilding global photon heap...\n"); | 
| 309 | 
  | 
  | 
      fflush(stderr); | 
| 310 | 
  | 
  | 
   } | 
| 311 | 
  | 
  | 
    | 
| 312 | 
  | 
  | 
   balancePhotons(pmap, NULL); | 
| 313 | 
  | 
  | 
} | 
| 314 | 
  | 
  | 
 | 
| 315 | 
  | 
  | 
 | 
| 316 | 
  | 
  | 
 | 
| 317 | 
  | 
  | 
void distribPhotons (PhotonMap **pmaps) | 
| 318 | 
  | 
  | 
{ | 
| 319 | 
  | 
  | 
   EmissionMap emap; | 
| 320 | 
  | 
  | 
   char errmsg2 [128]; | 
| 321 | 
  | 
  | 
   unsigned t, srcIdx, passCnt = 0, prePassCnt = 0; | 
| 322 | 
  | 
  | 
   double totalFlux = 0; | 
| 323 | 
  | 
  | 
   PhotonMap *pm; | 
| 324 | 
  | 
  | 
    | 
| 325 | 
  | 
  | 
   for (t = 0; t < NUM_PMAP_TYPES && !photonMaps [t]; t++); | 
| 326 | 
  | 
  | 
   if (t >= NUM_PMAP_TYPES) | 
| 327 | 
  | 
  | 
      error(USER, "no photon maps defined"); | 
| 328 | 
  | 
  | 
       | 
| 329 | 
  | 
  | 
   if (!nsources) | 
| 330 | 
  | 
  | 
      error(USER, "no light sources"); | 
| 331 | 
  | 
  | 
 | 
| 332 | 
  | 
  | 
   /* =================================================================== | 
| 333 | 
  | 
  | 
    * INITIALISATION - Set up emission and scattering funcs | 
| 334 | 
  | 
  | 
    * =================================================================== */ | 
| 335 | 
  | 
  | 
   emap.samples = NULL; | 
| 336 | 
  | 
  | 
   emap.maxPartitions = MAXSPART; | 
| 337 | 
  | 
  | 
   emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1); | 
| 338 | 
  | 
  | 
   if (!emap.partitions) | 
| 339 | 
  | 
  | 
      error(INTERNAL, "can't allocate source partitions"); | 
| 340 | 
  | 
  | 
       | 
| 341 | 
  | 
  | 
   /* Initialise all defined photon maps */ | 
| 342 | 
  | 
  | 
   for (t = 0; t < NUM_PMAP_TYPES; t++) | 
| 343 | 
  | 
  | 
      initPhotonMap(photonMaps [t], t); | 
| 344 | 
  | 
  | 
 | 
| 345 | 
  | 
  | 
   initPhotonEmissionFuncs(); | 
| 346 | 
  | 
  | 
   initPhotonScatterFuncs(); | 
| 347 | 
  | 
  | 
    | 
| 348 | 
  | 
  | 
   /* Get photon ports if specified */ | 
| 349 | 
  | 
  | 
   if (ambincl == 1)  | 
| 350 | 
  | 
  | 
      getPhotonPorts(); | 
| 351 | 
  | 
  | 
 | 
| 352 | 
  | 
  | 
   /* Get photon sensor modifiers */ | 
| 353 | 
  | 
  | 
   getPhotonSensors(photonSensorList); | 
| 354 | 
  | 
  | 
    | 
| 355 | 
  | 
  | 
   /* Seed RNGs for photon distribution */ | 
| 356 | 
  | 
  | 
   pmapSeed(randSeed, partState); | 
| 357 | 
  | 
  | 
   pmapSeed(randSeed, emitState); | 
| 358 | 
  | 
  | 
   pmapSeed(randSeed, cntState); | 
| 359 | 
  | 
  | 
   pmapSeed(randSeed, mediumState); | 
| 360 | 
  | 
  | 
   pmapSeed(randSeed, scatterState); | 
| 361 | 
  | 
  | 
   pmapSeed(randSeed, rouletteState); | 
| 362 | 
  | 
  | 
    | 
| 363 | 
  | 
  | 
   if (photonRepTime) | 
| 364 | 
  | 
  | 
      eputs("\n"); | 
| 365 | 
  | 
  | 
    | 
| 366 | 
  | 
  | 
   /* =================================================================== | 
| 367 | 
  | 
  | 
    * FLUX INTEGRATION - Get total photon flux from light sources | 
| 368 | 
  | 
  | 
    * =================================================================== */ | 
| 369 | 
  | 
  | 
   for (srcIdx = 0; srcIdx < nsources; srcIdx++) {          | 
| 370 | 
  | 
  | 
      unsigned portCnt = 0; | 
| 371 | 
  | 
  | 
      emap.src = source + srcIdx;  | 
| 372 | 
  | 
  | 
       | 
| 373 | 
  | 
  | 
      do { | 
| 374 | 
  | 
  | 
         emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt  | 
| 375 | 
  | 
  | 
                                                   : NULL; | 
| 376 | 
  | 
  | 
         photonPartition [emap.src -> so -> otype] (&emap); | 
| 377 | 
  | 
  | 
          | 
| 378 | 
  | 
  | 
         if (photonRepTime) { | 
| 379 | 
  | 
  | 
            sprintf(errmsg, "Integrating flux from source %s ",  | 
| 380 | 
  | 
  | 
                    source [srcIdx].so -> oname); | 
| 381 | 
  | 
  | 
                     | 
| 382 | 
  | 
  | 
            if (emap.port) { | 
| 383 | 
  | 
  | 
               sprintf(errmsg2, "via port %s ",  | 
| 384 | 
  | 
  | 
                       photonPorts [portCnt].so -> oname); | 
| 385 | 
  | 
  | 
               strcat(errmsg, errmsg2); | 
| 386 | 
  | 
  | 
            } | 
| 387 | 
  | 
  | 
             | 
| 388 | 
  | 
  | 
            sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions); | 
| 389 | 
  | 
  | 
            strcat(errmsg, errmsg2); | 
| 390 | 
  | 
  | 
            eputs(errmsg); | 
| 391 | 
  | 
  | 
            fflush(stderr); | 
| 392 | 
  | 
  | 
         } | 
| 393 | 
  | 
  | 
          | 
| 394 | 
  | 
  | 
         for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; | 
| 395 | 
  | 
  | 
              emap.partitionCnt++) { | 
| 396 | 
  | 
  | 
            initPhotonEmission(&emap, pdfSamples); | 
| 397 | 
  | 
  | 
            totalFlux += colorAvg(emap.partFlux); | 
| 398 | 
  | 
  | 
         } | 
| 399 | 
  | 
  | 
          | 
| 400 | 
  | 
  | 
         portCnt++; | 
| 401 | 
  | 
  | 
      } while (portCnt < numPhotonPorts); | 
| 402 | 
  | 
  | 
   } | 
| 403 | 
  | 
  | 
 | 
| 404 | 
  | 
  | 
   if (totalFlux < FTINY) | 
| 405 | 
  | 
  | 
      error(USER, "zero flux from light sources"); | 
| 406 | 
  | 
  | 
 | 
| 407 | 
  | 
  | 
   /* Record start time and enable progress report signal handler */ | 
| 408 | 
  | 
  | 
   repStartTime = time(NULL); | 
| 409 | 
rschregle | 
2.3 | 
   #ifdef SIGCONT | 
| 410 | 
  | 
  | 
      signal(SIGCONT, pmapDistribReport); | 
| 411 | 
  | 
  | 
   #endif | 
| 412 | 
greg | 
2.1 | 
   repProgress = prePassCnt = 0; | 
| 413 | 
  | 
  | 
    | 
| 414 | 
  | 
  | 
   if (photonRepTime) | 
| 415 | 
  | 
  | 
      eputs("\n"); | 
| 416 | 
  | 
  | 
    | 
| 417 | 
  | 
  | 
   /* =================================================================== | 
| 418 | 
  | 
  | 
    * 2-PASS PHOTON DISTRIBUTION | 
| 419 | 
  | 
  | 
    * Pass 1 (pre):  emit fraction of target photon count | 
| 420 | 
  | 
  | 
    * Pass 2 (main): based on outcome of pass 1, estimate remaining number | 
| 421 | 
  | 
  | 
    *                of photons to emit to approximate target count | 
| 422 | 
  | 
  | 
    * =================================================================== */ | 
| 423 | 
  | 
  | 
   do { | 
| 424 | 
  | 
  | 
      double numEmit; | 
| 425 | 
  | 
  | 
       | 
| 426 | 
  | 
  | 
      if (!passCnt) {    | 
| 427 | 
  | 
  | 
         /* INIT PASS 1 */ | 
| 428 | 
  | 
  | 
         /* Skip if no photons contributed after sufficient iterations; make | 
| 429 | 
  | 
  | 
          * it clear to user which photon maps are missing so (s)he can | 
| 430 | 
  | 
  | 
          * check the scene geometry and materials */ | 
| 431 | 
  | 
  | 
         if (++prePassCnt > maxPreDistrib) { | 
| 432 | 
  | 
  | 
            sprintf(errmsg, "too many prepasses"); | 
| 433 | 
  | 
  | 
 | 
| 434 | 
  | 
  | 
            for (t = 0; t < NUM_PMAP_TYPES; t++) | 
| 435 | 
  | 
  | 
               if (photonMaps [t] && !photonMaps [t] -> heapEnd) { | 
| 436 | 
  | 
  | 
                  sprintf(errmsg2, ", no %s photons stored", pmapName [t]); | 
| 437 | 
  | 
  | 
                  strcat(errmsg, errmsg2); | 
| 438 | 
  | 
  | 
               } | 
| 439 | 
  | 
  | 
             | 
| 440 | 
  | 
  | 
            error(USER, errmsg); | 
| 441 | 
  | 
  | 
            break; | 
| 442 | 
  | 
  | 
         } | 
| 443 | 
  | 
  | 
 | 
| 444 | 
  | 
  | 
         /* Num to emit is fraction of minimum target count */ | 
| 445 | 
  | 
  | 
         numEmit = FHUGE; | 
| 446 | 
  | 
  | 
          | 
| 447 | 
  | 
  | 
         for (t = 0; t < NUM_PMAP_TYPES; t++) | 
| 448 | 
  | 
  | 
            if (photonMaps [t]) | 
| 449 | 
  | 
  | 
               numEmit = min(photonMaps [t] -> distribTarget, numEmit); | 
| 450 | 
  | 
  | 
                | 
| 451 | 
  | 
  | 
         numEmit *= preDistrib; | 
| 452 | 
  | 
  | 
      } | 
| 453 | 
  | 
  | 
 | 
| 454 | 
  | 
  | 
      else {             | 
| 455 | 
  | 
  | 
         /* INIT PASS 2 */ | 
| 456 | 
  | 
  | 
         /* Based on the outcome of the predistribution we can now estimate | 
| 457 | 
  | 
  | 
          * how many more photons we have to emit for each photon map to | 
| 458 | 
  | 
  | 
          * meet its respective target count. This value is clamped to 0 in | 
| 459 | 
  | 
  | 
          * case the target has already been exceeded in the pass 1. Note | 
| 460 | 
  | 
  | 
          * repProgress is the number of photons emitted thus far, while | 
| 461 | 
  | 
  | 
          * heapEnd is the number of photons stored in each photon map. */ | 
| 462 | 
  | 
  | 
         double maxDistribRatio = 0; | 
| 463 | 
  | 
  | 
 | 
| 464 | 
  | 
  | 
         /* Set the distribution ratio for each map; this indicates how many | 
| 465 | 
  | 
  | 
          * photons of each respective type are stored per emitted photon, | 
| 466 | 
  | 
  | 
          * and is used as probability for storing a photon by addPhoton(). | 
| 467 | 
  | 
  | 
          * Since this biases the photon density, addPhoton() promotes the | 
| 468 | 
  | 
  | 
          * flux of stored photons to compensate. */ | 
| 469 | 
  | 
  | 
         for (t = 0; t < NUM_PMAP_TYPES; t++) | 
| 470 | 
  | 
  | 
            if ((pm = photonMaps [t])) { | 
| 471 | 
  | 
  | 
               pm -> distribRatio = (double)pm -> distribTarget /  | 
| 472 | 
  | 
  | 
                                    pm -> heapEnd - 1; | 
| 473 | 
  | 
  | 
 | 
| 474 | 
  | 
  | 
               /* Check if photon map "overflowed", i.e. exceeded its target | 
| 475 | 
  | 
  | 
                * count in the prepass; correcting the photon flux via the | 
| 476 | 
  | 
  | 
                * distribution ratio is no longer possible, as no more | 
| 477 | 
  | 
  | 
                * photons of this type will be stored, so notify the user | 
| 478 | 
  | 
  | 
                * rather than deliver incorrect results.  | 
| 479 | 
  | 
  | 
                * In future we should handle this more intelligently by | 
| 480 | 
  | 
  | 
                * using the photonFlux in each photon map to individually | 
| 481 | 
  | 
  | 
                * correct the flux after distribution. */ | 
| 482 | 
  | 
  | 
               if (pm -> distribRatio <= FTINY) { | 
| 483 | 
  | 
  | 
                  sprintf(errmsg,  | 
| 484 | 
  | 
  | 
                          "%s photon map overflow in prepass, reduce -apD", | 
| 485 | 
  | 
  | 
                          pmapName [t]); | 
| 486 | 
  | 
  | 
                  error(INTERNAL, errmsg); | 
| 487 | 
  | 
  | 
               } | 
| 488 | 
  | 
  | 
                   | 
| 489 | 
  | 
  | 
               maxDistribRatio = max(pm -> distribRatio, maxDistribRatio); | 
| 490 | 
  | 
  | 
            } | 
| 491 | 
  | 
  | 
          | 
| 492 | 
  | 
  | 
         /* Normalise distribution ratios and calculate number of photons to | 
| 493 | 
  | 
  | 
          * emit in main pass */ | 
| 494 | 
  | 
  | 
         for (t = 0; t < NUM_PMAP_TYPES; t++) | 
| 495 | 
  | 
  | 
            if ((pm = photonMaps [t])) | 
| 496 | 
  | 
  | 
               pm -> distribRatio /= maxDistribRatio; | 
| 497 | 
  | 
  | 
                | 
| 498 | 
  | 
  | 
         if ((numEmit = repProgress * maxDistribRatio) < FTINY) | 
| 499 | 
  | 
  | 
            /* No photons left to distribute in main pass */ | 
| 500 | 
  | 
  | 
            break; | 
| 501 | 
  | 
  | 
      } | 
| 502 | 
  | 
  | 
       | 
| 503 | 
  | 
  | 
      /* Set completion count for progress report */ | 
| 504 | 
  | 
  | 
      repComplete = numEmit + repProgress;                              | 
| 505 | 
  | 
  | 
       | 
| 506 | 
  | 
  | 
      /* PHOTON DISTRIBUTION LOOP */ | 
| 507 | 
  | 
  | 
      for (srcIdx = 0; srcIdx < nsources; srcIdx++) { | 
| 508 | 
  | 
  | 
         unsigned portCnt = 0; | 
| 509 | 
  | 
  | 
         emap.src = source + srcIdx; | 
| 510 | 
  | 
  | 
                   | 
| 511 | 
  | 
  | 
         do { | 
| 512 | 
  | 
  | 
            emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt  | 
| 513 | 
  | 
  | 
                                                      : NULL; | 
| 514 | 
  | 
  | 
            photonPartition [emap.src -> so -> otype] (&emap); | 
| 515 | 
  | 
  | 
             | 
| 516 | 
  | 
  | 
            if (photonRepTime) { | 
| 517 | 
  | 
  | 
               if (!passCnt) | 
| 518 | 
  | 
  | 
                  sprintf(errmsg, "PREPASS %d on source %s ",  | 
| 519 | 
  | 
  | 
                          prePassCnt, source [srcIdx].so -> oname); | 
| 520 | 
  | 
  | 
               else  | 
| 521 | 
  | 
  | 
                  sprintf(errmsg, "MAIN PASS on source %s ", | 
| 522 | 
  | 
  | 
                          source [srcIdx].so -> oname); | 
| 523 | 
  | 
  | 
                        | 
| 524 | 
  | 
  | 
               if (emap.port) { | 
| 525 | 
  | 
  | 
                  sprintf(errmsg2, "via port %s ",  | 
| 526 | 
  | 
  | 
                          photonPorts [portCnt].so -> oname); | 
| 527 | 
  | 
  | 
                  strcat(errmsg, errmsg2); | 
| 528 | 
  | 
  | 
               } | 
| 529 | 
  | 
  | 
                | 
| 530 | 
  | 
  | 
               sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions); | 
| 531 | 
  | 
  | 
               strcat(errmsg, errmsg2); | 
| 532 | 
  | 
  | 
               eputs(errmsg); | 
| 533 | 
  | 
  | 
               fflush(stderr); | 
| 534 | 
  | 
  | 
            } | 
| 535 | 
  | 
  | 
             | 
| 536 | 
  | 
  | 
            for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; | 
| 537 | 
  | 
  | 
                 emap.partitionCnt++) {              | 
| 538 | 
  | 
  | 
               double partNumEmit; | 
| 539 | 
  | 
  | 
               unsigned long partEmitCnt; | 
| 540 | 
  | 
  | 
                | 
| 541 | 
  | 
  | 
               /* Get photon origin within current source partishunn and | 
| 542 | 
  | 
  | 
                * build emission map */ | 
| 543 | 
  | 
  | 
               photonOrigin [emap.src -> so -> otype] (&emap); | 
| 544 | 
  | 
  | 
               initPhotonEmission(&emap, pdfSamples); | 
| 545 | 
  | 
  | 
                | 
| 546 | 
  | 
  | 
               /* Number of photons to emit from ziss partishunn -- | 
| 547 | 
  | 
  | 
                * proportional to flux; photon ray weight and scalar flux | 
| 548 | 
  | 
  | 
                * are uniform (the latter only varying in RGB). */ | 
| 549 | 
  | 
  | 
               partNumEmit = numEmit * colorAvg(emap.partFlux) / totalFlux; | 
| 550 | 
  | 
  | 
               partEmitCnt = (unsigned long)partNumEmit; | 
| 551 | 
  | 
  | 
                | 
| 552 | 
  | 
  | 
               /* Probabilistically account for fractional photons */ | 
| 553 | 
  | 
  | 
               if (pmapRandom(cntState) < partNumEmit - partEmitCnt) | 
| 554 | 
  | 
  | 
                  partEmitCnt++; | 
| 555 | 
  | 
  | 
 | 
| 556 | 
  | 
  | 
               /* Integer counter avoids FP rounding errors */ | 
| 557 | 
  | 
  | 
               while (partEmitCnt--) { | 
| 558 | 
  | 
  | 
                  RAY photonRay; | 
| 559 | 
  | 
  | 
                   | 
| 560 | 
  | 
  | 
                  /* Emit photon based on PDF and trace through scene until | 
| 561 | 
  | 
  | 
                   * absorbed/leaked */ | 
| 562 | 
  | 
  | 
                  emitPhoton(&emap, &photonRay); | 
| 563 | 
  | 
  | 
                  tracePhoton(&photonRay); | 
| 564 | 
  | 
  | 
                   | 
| 565 | 
  | 
  | 
                  /* Record progress */ | 
| 566 | 
  | 
  | 
                  repProgress++; | 
| 567 | 
  | 
  | 
                   | 
| 568 | 
  | 
  | 
                  if (photonRepTime > 0 &&  | 
| 569 | 
  | 
  | 
                      time(NULL) >= repLastTime + photonRepTime) | 
| 570 | 
  | 
  | 
                     pmapDistribReport(); | 
| 571 | 
rschregle | 
2.3 | 
                  #ifdef SIGCONT | 
| 572 | 
greg | 
2.1 | 
                     else signal(SIGCONT, pmapDistribReport); | 
| 573 | 
  | 
  | 
                  #endif | 
| 574 | 
  | 
  | 
               } | 
| 575 | 
  | 
  | 
            } | 
| 576 | 
  | 
  | 
                         | 
| 577 | 
  | 
  | 
            portCnt++; | 
| 578 | 
  | 
  | 
         } while (portCnt < numPhotonPorts); | 
| 579 | 
  | 
  | 
      } | 
| 580 | 
  | 
  | 
       | 
| 581 | 
  | 
  | 
      for (t = 0; t < NUM_PMAP_TYPES; t++) | 
| 582 | 
  | 
  | 
         if (photonMaps [t] && !photonMaps [t] -> heapEnd) { | 
| 583 | 
  | 
  | 
            /* Double preDistrib in case a photon map is empty and redo | 
| 584 | 
  | 
  | 
             * pass 1 --> possibility of infinite loop for pathological | 
| 585 | 
  | 
  | 
             * scenes (e.g. absorbing materials) */ | 
| 586 | 
  | 
  | 
            preDistrib *= 2; | 
| 587 | 
  | 
  | 
            break; | 
| 588 | 
  | 
  | 
         } | 
| 589 | 
  | 
  | 
       | 
| 590 | 
  | 
  | 
      if (t >= NUM_PMAP_TYPES) { | 
| 591 | 
  | 
  | 
         /* No empty photon maps found; now do pass 2 */ | 
| 592 | 
  | 
  | 
         passCnt++; | 
| 593 | 
  | 
  | 
         if (photonRepTime) | 
| 594 | 
  | 
  | 
            eputs("\n"); | 
| 595 | 
  | 
  | 
      } | 
| 596 | 
  | 
  | 
   } while (passCnt < 2); | 
| 597 | 
  | 
  | 
 | 
| 598 | 
  | 
  | 
   /* =================================================================== | 
| 599 | 
  | 
  | 
    * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc. | 
| 600 | 
  | 
  | 
    * =================================================================== */ | 
| 601 | 
rschregle | 
2.3 | 
   #ifdef SIGCONT     | 
| 602 | 
  | 
  | 
      signal(SIGCONT, SIG_DFL); | 
| 603 | 
  | 
  | 
   #endif | 
| 604 | 
greg | 
2.1 | 
   free(emap.samples); | 
| 605 | 
  | 
  | 
    | 
| 606 | 
  | 
  | 
   /* Set photon flux (repProgress is total num emitted) */ | 
| 607 | 
  | 
  | 
   totalFlux /= repProgress; | 
| 608 | 
  | 
  | 
    | 
| 609 | 
  | 
  | 
   for (t = 0; t < NUM_PMAP_TYPES; t++) | 
| 610 | 
  | 
  | 
      if (photonMaps [t]) { | 
| 611 | 
  | 
  | 
         if (photonRepTime) { | 
| 612 | 
  | 
  | 
            sprintf(errmsg, "\nBuilding %s photon map...\n", pmapName [t]); | 
| 613 | 
  | 
  | 
            eputs(errmsg); | 
| 614 | 
  | 
  | 
            fflush(stderr); | 
| 615 | 
  | 
  | 
         } | 
| 616 | 
  | 
  | 
       | 
| 617 | 
  | 
  | 
         balancePhotons(photonMaps [t], &totalFlux); | 
| 618 | 
  | 
  | 
      } | 
| 619 | 
  | 
  | 
       | 
| 620 | 
  | 
  | 
   /* Precompute photon irradiance if necessary */ | 
| 621 | 
  | 
  | 
   if (preCompPmap) | 
| 622 | 
  | 
  | 
      preComputeGlobal(preCompPmap); | 
| 623 | 
  | 
  | 
} | 
| 624 | 
  | 
  | 
 | 
| 625 | 
  | 
  | 
 | 
| 626 | 
  | 
  | 
 | 
| 627 | 
  | 
  | 
void photonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad) | 
| 628 | 
  | 
  | 
/* Photon density estimate. Returns irradiance at ray -> rop. */ | 
| 629 | 
  | 
  | 
{ | 
| 630 | 
  | 
  | 
   unsigned i; | 
| 631 | 
  | 
  | 
   PhotonSQNode *sq; | 
| 632 | 
  | 
  | 
   float r; | 
| 633 | 
  | 
  | 
   COLOR flux; | 
| 634 | 
  | 
  | 
  | 
| 635 | 
  | 
  | 
   setcolor(irrad, 0, 0, 0); | 
| 636 | 
  | 
  | 
 | 
| 637 | 
  | 
  | 
   if (!pmap -> maxGather)  | 
| 638 | 
  | 
  | 
      return; | 
| 639 | 
  | 
  | 
       | 
| 640 | 
  | 
  | 
   /* Ignore sources */ | 
| 641 | 
  | 
  | 
   if (ray -> ro)  | 
| 642 | 
  | 
  | 
      if (islight(objptr(ray -> ro -> omod) -> otype))  | 
| 643 | 
  | 
  | 
         return; | 
| 644 | 
  | 
  | 
          | 
| 645 | 
  | 
  | 
   pmap -> squeueEnd = 0; | 
| 646 | 
  | 
  | 
   findPhotons(pmap, ray); | 
| 647 | 
  | 
  | 
    | 
| 648 | 
  | 
  | 
   /* Need at least 2 photons */ | 
| 649 | 
  | 
  | 
   if (pmap -> squeueEnd < 2) { | 
| 650 | 
  | 
  | 
      #ifdef PMAP_NONEFOUND    | 
| 651 | 
  | 
  | 
         sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)",  | 
| 652 | 
  | 
  | 
                 ray -> ro ? ray -> ro -> oname : "<null>", | 
| 653 | 
  | 
  | 
                 ray -> rop [0], ray -> rop [1], ray -> rop [2]); | 
| 654 | 
  | 
  | 
         error(WARNING, errmsg); | 
| 655 | 
  | 
  | 
      #endif       | 
| 656 | 
  | 
  | 
 | 
| 657 | 
  | 
  | 
      return; | 
| 658 | 
  | 
  | 
   } | 
| 659 | 
  | 
  | 
       | 
| 660 | 
  | 
  | 
   if (pmap -> minGather == pmap -> maxGather) { | 
| 661 | 
  | 
  | 
      /* No bias compensation. Just do a plain vanilla estimate */ | 
| 662 | 
  | 
  | 
      sq = pmap -> squeue + 1; | 
| 663 | 
  | 
  | 
       | 
| 664 | 
  | 
  | 
      /* Average radius between furthest two photons to improve accuracy */       | 
| 665 | 
  | 
  | 
      r = max(sq -> dist, (sq + 1) -> dist); | 
| 666 | 
  | 
  | 
      r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r));    | 
| 667 | 
  | 
  | 
       | 
| 668 | 
  | 
  | 
      /* Skip the extra photon */ | 
| 669 | 
  | 
  | 
      for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) { | 
| 670 | 
  | 
  | 
         getPhotonFlux(sq -> photon, flux);          | 
| 671 | 
  | 
  | 
#ifdef PMAP_EPANECHNIKOV | 
| 672 | 
  | 
  | 
         /* Apply Epanechnikov kernel to photon flux (dists are squared) */ | 
| 673 | 
  | 
  | 
         scalecolor(flux, 2 * (1 - sq -> dist / r)); | 
| 674 | 
  | 
  | 
#endif          | 
| 675 | 
  | 
  | 
         addcolor(irrad, flux); | 
| 676 | 
  | 
  | 
      } | 
| 677 | 
  | 
  | 
       | 
| 678 | 
  | 
  | 
      /* Divide by search area PI * r^2, 1 / PI required as ambient  | 
| 679 | 
  | 
  | 
         normalisation factor */          | 
| 680 | 
  | 
  | 
      scalecolor(irrad, 1 / (PI * PI * r));  | 
| 681 | 
  | 
  | 
       | 
| 682 | 
  | 
  | 
      return; | 
| 683 | 
  | 
  | 
   } | 
| 684 | 
  | 
  | 
   else  | 
| 685 | 
  | 
  | 
      /* Apply bias compensation to density estimate */ | 
| 686 | 
  | 
  | 
      biasComp(pmap, irrad); | 
| 687 | 
  | 
  | 
} | 
| 688 | 
  | 
  | 
 | 
| 689 | 
  | 
  | 
 | 
| 690 | 
  | 
  | 
 | 
| 691 | 
  | 
  | 
void photonPreCompDensity (PhotonMap *pmap, RAY *r, COLOR irrad) | 
| 692 | 
  | 
  | 
/* Returns precomputed photon density estimate at ray -> rop. */ | 
| 693 | 
  | 
  | 
{ | 
| 694 | 
  | 
  | 
   Photon *p; | 
| 695 | 
  | 
  | 
    | 
| 696 | 
  | 
  | 
   setcolor(irrad, 0, 0, 0); | 
| 697 | 
  | 
  | 
 | 
| 698 | 
  | 
  | 
   /* Ignore sources */ | 
| 699 | 
  | 
  | 
   if (r -> ro && islight(objptr(r -> ro -> omod) -> otype))  | 
| 700 | 
  | 
  | 
      return; | 
| 701 | 
  | 
  | 
       | 
| 702 | 
  | 
  | 
   if ((p = find1Photon(preCompPmap, r)))  | 
| 703 | 
  | 
  | 
      getPhotonFlux(p, irrad); | 
| 704 | 
  | 
  | 
} | 
| 705 | 
  | 
  | 
 | 
| 706 | 
  | 
  | 
 | 
| 707 | 
  | 
  | 
 | 
| 708 | 
  | 
  | 
void volumePhotonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad) | 
| 709 | 
  | 
  | 
/* Photon volume density estimate. Returns irradiance at ray -> rop. */ | 
| 710 | 
  | 
  | 
{ | 
| 711 | 
  | 
  | 
   unsigned i; | 
| 712 | 
  | 
  | 
   PhotonSQNode *sq; | 
| 713 | 
  | 
  | 
   float gecc2, r, ph; | 
| 714 | 
  | 
  | 
   COLOR flux; | 
| 715 | 
  | 
  | 
 | 
| 716 | 
  | 
  | 
   setcolor(irrad, 0, 0, 0); | 
| 717 | 
  | 
  | 
    | 
| 718 | 
  | 
  | 
   if (!pmap -> maxGather)  | 
| 719 | 
  | 
  | 
      return; | 
| 720 | 
  | 
  | 
       | 
| 721 | 
  | 
  | 
   pmap -> squeueEnd = 0; | 
| 722 | 
  | 
  | 
   findPhotons(pmap, ray); | 
| 723 | 
  | 
  | 
    | 
| 724 | 
  | 
  | 
   /* Need at least 2 photons */ | 
| 725 | 
  | 
  | 
   if (pmap -> squeueEnd < 2)  | 
| 726 | 
  | 
  | 
      return; | 
| 727 | 
  | 
  | 
       | 
| 728 | 
  | 
  | 
   if (pmap -> minGather == pmap -> maxGather) { | 
| 729 | 
  | 
  | 
      /* No bias compensation. Just do a plain vanilla estimate */ | 
| 730 | 
  | 
  | 
      gecc2 = ray -> gecc * ray -> gecc; | 
| 731 | 
  | 
  | 
      sq = pmap -> squeue + 1; | 
| 732 | 
  | 
  | 
       | 
| 733 | 
  | 
  | 
      /* Average radius between furthest two photons to improve accuracy */       | 
| 734 | 
  | 
  | 
      r = max(sq -> dist, (sq + 1) -> dist); | 
| 735 | 
  | 
  | 
      r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r));    | 
| 736 | 
  | 
  | 
       | 
| 737 | 
  | 
  | 
      /* Skip the extra photon */ | 
| 738 | 
  | 
  | 
      for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) { | 
| 739 | 
  | 
  | 
         /* Compute phase function for inscattering from photon */ | 
| 740 | 
  | 
  | 
         if (gecc2 <= FTINY)  | 
| 741 | 
  | 
  | 
            ph = 1; | 
| 742 | 
  | 
  | 
         else { | 
| 743 | 
  | 
  | 
            ph = DOT(ray -> rdir, sq -> photon -> norm) / 127; | 
| 744 | 
  | 
  | 
            ph = 1 + gecc2 - 2 * ray -> gecc * ph; | 
| 745 | 
  | 
  | 
            ph = (1 - gecc2) / (ph * sqrt(ph)); | 
| 746 | 
  | 
  | 
         } | 
| 747 | 
  | 
  | 
          | 
| 748 | 
  | 
  | 
         getPhotonFlux(sq -> photon, flux); | 
| 749 | 
  | 
  | 
         scalecolor(flux, ph); | 
| 750 | 
  | 
  | 
         addcolor(irrad, flux); | 
| 751 | 
  | 
  | 
      } | 
| 752 | 
  | 
  | 
       | 
| 753 | 
  | 
  | 
      /* Divide by search volume 4 / 3 * PI * r^3 and phase function | 
| 754 | 
  | 
  | 
         normalization factor 1 / (4 * PI) */ | 
| 755 | 
  | 
  | 
      scalecolor(irrad, 3 / (16 * PI * PI * r * sqrt(r))); | 
| 756 | 
  | 
  | 
       | 
| 757 | 
  | 
  | 
      return; | 
| 758 | 
  | 
  | 
   } | 
| 759 | 
  | 
  | 
    | 
| 760 | 
  | 
  | 
   else  | 
| 761 | 
  | 
  | 
      /* Apply bias compensation to density estimate */ | 
| 762 | 
  | 
  | 
      volumeBiasComp(pmap, ray, irrad); | 
| 763 | 
  | 
  | 
} |