| 1 | 
/*  | 
| 2 | 
   ====================================================================== | 
| 3 | 
   In-core kd-tree for photon map | 
| 4 | 
 | 
| 5 | 
   Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) | 
| 6 | 
   (c) Fraunhofer Institute for Solar Energy Systems, | 
| 7 | 
   (c) Lucerne University of Applied Sciences and Arts, | 
| 8 | 
       supported by the Swiss National Science Foundation (SNSF, #147053) | 
| 9 | 
   ====================================================================== | 
| 10 | 
    | 
| 11 | 
   $Id: pmapkdt.c,v 1.1 2016/05/18 08:22:45 rschregle Exp $ | 
| 12 | 
*/ | 
| 13 | 
 | 
| 14 | 
 | 
| 15 | 
 | 
| 16 | 
#include "pmapdata.h"   /* Includes pmapkdt.h */ | 
| 17 | 
#include "source.h" | 
| 18 | 
 | 
| 19 | 
 | 
| 20 | 
 | 
| 21 | 
 | 
| 22 | 
void kdT_Null (PhotonKdTree *kdt) | 
| 23 | 
{ | 
| 24 | 
   kdt -> nodes = NULL; | 
| 25 | 
} | 
| 26 | 
 | 
| 27 | 
 | 
| 28 | 
 | 
| 29 | 
static unsigned long kdT_MedianPartition (const Photon *heap,  | 
| 30 | 
                                          unsigned long *heapIdx, | 
| 31 | 
                                          unsigned long *heapXdi,  | 
| 32 | 
                                          unsigned long left,  | 
| 33 | 
                                          unsigned long right, unsigned dim) | 
| 34 | 
/* Returns index to median in heap from indices left to right  | 
| 35 | 
   (inclusive) in dimension dim. The heap is partitioned relative to  | 
| 36 | 
   median using a quicksort algorithm. The heap indices in heapIdx are | 
| 37 | 
   sorted rather than the heap itself. */ | 
| 38 | 
{ | 
| 39 | 
   const float    *p; | 
| 40 | 
   unsigned long  l, r, lg2, n2, m, n = right - left + 1; | 
| 41 | 
   unsigned       d; | 
| 42 | 
    | 
| 43 | 
   /* Round down n to nearest power of 2 */ | 
| 44 | 
   for (lg2 = 0, n2 = n; n2 > 1; n2 >>= 1, ++lg2); | 
| 45 | 
   n2 = 1 << lg2; | 
| 46 | 
    | 
| 47 | 
   /* Determine median position; this takes into account the fact that  | 
| 48 | 
      only the last level in the heap can be partially empty, and that  | 
| 49 | 
      it fills from left to right */ | 
| 50 | 
   m = left + ((n - n2) > (n2 >> 1) - 1 ? n2 - 1 : n - (n2 >> 1)); | 
| 51 | 
    | 
| 52 | 
   while (right > left) { | 
| 53 | 
      /* Pivot node */ | 
| 54 | 
      p = heap [heapIdx [right]].pos; | 
| 55 | 
      l = left; | 
| 56 | 
      r = right - 1; | 
| 57 | 
       | 
| 58 | 
      /* l & r converge, swapping elements out of order with respect to  | 
| 59 | 
         pivot node. Identical keys are resolved by cycling through | 
| 60 | 
         dim. The convergence point is then the pivot's position. */ | 
| 61 | 
      do { | 
| 62 | 
         while (l <= r) { | 
| 63 | 
            d = dim; | 
| 64 | 
             | 
| 65 | 
            while (heap [heapIdx [l]].pos [d] == p [d]) { | 
| 66 | 
               d = (d + 1) % 3; | 
| 67 | 
                | 
| 68 | 
               if (d == dim) { | 
| 69 | 
                  /* Ignore dupes? */ | 
| 70 | 
                  error(WARNING, "duplicate keys in photon heap"); | 
| 71 | 
                  l++; | 
| 72 | 
                  break; | 
| 73 | 
               } | 
| 74 | 
            } | 
| 75 | 
             | 
| 76 | 
            if (heap [heapIdx [l]].pos [d] < p [d])  | 
| 77 | 
               l++; | 
| 78 | 
            else break; | 
| 79 | 
         } | 
| 80 | 
          | 
| 81 | 
         while (r > l) { | 
| 82 | 
            d = dim; | 
| 83 | 
             | 
| 84 | 
            while (heap [heapIdx [r]].pos [d] == p [d]) { | 
| 85 | 
               d = (d + 1) % 3; | 
| 86 | 
                | 
| 87 | 
               if (d == dim) { | 
| 88 | 
                  /* Ignore dupes? */ | 
| 89 | 
                  error(WARNING, "duplicate keys in photon heap"); | 
| 90 | 
                  r--; | 
| 91 | 
                  break; | 
| 92 | 
               } | 
| 93 | 
            } | 
| 94 | 
             | 
| 95 | 
            if (heap [heapIdx [r]].pos [d] > p [d])  | 
| 96 | 
               r--; | 
| 97 | 
            else break; | 
| 98 | 
         } | 
| 99 | 
          | 
| 100 | 
         /* Swap indices (not the nodes they point to) */ | 
| 101 | 
         n2 = heapIdx [l]; | 
| 102 | 
         heapIdx [l] = heapIdx [r]; | 
| 103 | 
         heapIdx [r] = n2; | 
| 104 | 
         /* Update reverse indices */ | 
| 105 | 
         heapXdi [heapIdx [l]] = l; | 
| 106 | 
         heapXdi [n2] = r; | 
| 107 | 
      } while (l < r); | 
| 108 | 
       | 
| 109 | 
      /* Swap indices of convergence and pivot nodes */ | 
| 110 | 
      heapIdx [r] = heapIdx [l]; | 
| 111 | 
      heapIdx [l] = heapIdx [right]; | 
| 112 | 
      heapIdx [right] = n2; | 
| 113 | 
      /* Update reverse indices */ | 
| 114 | 
      heapXdi [heapIdx [r]] = r; | 
| 115 | 
      heapXdi [heapIdx [l]] = l; | 
| 116 | 
      heapXdi [n2] = right; | 
| 117 | 
       | 
| 118 | 
      if (l >= m)  | 
| 119 | 
         right = l - 1; | 
| 120 | 
      if (l <= m)  | 
| 121 | 
         left = l + 1; | 
| 122 | 
   } | 
| 123 | 
    | 
| 124 | 
   /* Once left & right have converged at m, we have found the median */ | 
| 125 | 
   return m; | 
| 126 | 
} | 
| 127 | 
 | 
| 128 | 
 | 
| 129 | 
 | 
| 130 | 
static void kdT_Build (Photon *heap, unsigned long *heapIdx, | 
| 131 | 
                       unsigned long *heapXdi, const float min [3],  | 
| 132 | 
                       const float max [3], unsigned long left,  | 
| 133 | 
                       unsigned long right, unsigned long root) | 
| 134 | 
/* Recursive part of balancePhotons(..). Builds heap from subarray | 
| 135 | 
   defined by indices left and right. min and max are the minimum resp.  | 
| 136 | 
   maximum photon positions in the array. root is the index of the | 
| 137 | 
   current subtree's root, which corresponds to the median's 1-based  | 
| 138 | 
   index in the heap. heapIdx are the balanced heap indices. The heap | 
| 139 | 
   is accessed indirectly through these. heapXdi are the reverse indices  | 
| 140 | 
   from the heap to heapIdx so that heapXdi [heapIdx [i]] = i. */ | 
| 141 | 
{ | 
| 142 | 
   float                maxLeft [3], minRight [3]; | 
| 143 | 
   Photon               rootNode; | 
| 144 | 
   unsigned             d; | 
| 145 | 
    | 
| 146 | 
   /* Choose median for dimension with largest spread and partition  | 
| 147 | 
      accordingly */ | 
| 148 | 
   const float          d0 = max [0] - min [0],  | 
| 149 | 
                        d1 = max [1] - min [1],  | 
| 150 | 
                        d2 = max [2] - min [2]; | 
| 151 | 
   const unsigned char  dim = d0 > d1 ? d0 > d2 ? 0 : 2  | 
| 152 | 
                                      : d1 > d2 ? 1 : 2; | 
| 153 | 
   const unsigned long  median = left == right  | 
| 154 | 
                                 ? left  | 
| 155 | 
                                 : kdT_MedianPartition(heap, heapIdx, heapXdi, | 
| 156 | 
                                                       left, right, dim); | 
| 157 | 
    | 
| 158 | 
   /* Place median at root of current subtree. This consists of swapping  | 
| 159 | 
      the median and the root nodes and updating the heap indices */ | 
| 160 | 
   memcpy(&rootNode, heap + heapIdx [median], sizeof(Photon)); | 
| 161 | 
   memcpy(heap + heapIdx [median], heap + root - 1, sizeof(Photon)); | 
| 162 | 
   rootNode.discr = dim; | 
| 163 | 
   memcpy(heap + root - 1, &rootNode, sizeof(Photon)); | 
| 164 | 
   heapIdx [heapXdi [root - 1]] = heapIdx [median]; | 
| 165 | 
   heapXdi [heapIdx [median]] = heapXdi [root - 1]; | 
| 166 | 
   heapIdx [median] = root - 1; | 
| 167 | 
   heapXdi [root - 1] = median; | 
| 168 | 
    | 
| 169 | 
   /* Update bounds for left and right subtrees and recurse on them */ | 
| 170 | 
   for (d = 0; d <= 2; d++) | 
| 171 | 
      if (d == dim)  | 
| 172 | 
         maxLeft [d] = minRight [d] = rootNode.pos [d]; | 
| 173 | 
      else { | 
| 174 | 
         maxLeft [d] = max [d]; | 
| 175 | 
         minRight [d] = min [d]; | 
| 176 | 
      } | 
| 177 | 
       | 
| 178 | 
   if (left < median)  | 
| 179 | 
      kdT_Build(heap, heapIdx, heapXdi, min, maxLeft, left, median - 1,  | 
| 180 | 
                root << 1); | 
| 181 | 
                 | 
| 182 | 
   if (right > median)  | 
| 183 | 
      kdT_Build(heap, heapIdx, heapXdi, minRight, max, median + 1, right,  | 
| 184 | 
                (root << 1) + 1); | 
| 185 | 
} | 
| 186 | 
 | 
| 187 | 
 | 
| 188 | 
 | 
| 189 | 
void kdT_BuildPhotonMap (struct PhotonMap *pmap) | 
| 190 | 
{ | 
| 191 | 
   Photon         *nodes; | 
| 192 | 
   unsigned long  i; | 
| 193 | 
   unsigned long  *heapIdx,        /* Photon index array */ | 
| 194 | 
                  *heapXdi;        /* Reverse index to heapIdx */ | 
| 195 | 
    | 
| 196 | 
   /* Allocate kd-tree nodes and load photons from heap file */ | 
| 197 | 
   if (!(nodes = calloc(pmap -> numPhotons, sizeof(Photon)))) | 
| 198 | 
      error(SYSTEM, "failed in-core heap allocation in kdT_BuildPhotonMap"); | 
| 199 | 
      | 
| 200 | 
   rewind(pmap -> heap); | 
| 201 | 
   i = fread(nodes, sizeof(Photon), pmap -> numPhotons, pmap -> heap); | 
| 202 | 
   if (i != | 
| 203 | 
       pmap -> numPhotons) | 
| 204 | 
      error(SYSTEM, "failed loading photon heap in kdT_BuildPhotonMap"); | 
| 205 | 
       | 
| 206 | 
   pmap -> store.nodes = nodes; | 
| 207 | 
   heapIdx = calloc(pmap -> numPhotons, sizeof(unsigned long)); | 
| 208 | 
   heapXdi = calloc(pmap -> numPhotons, sizeof(unsigned long)); | 
| 209 | 
   if (!heapIdx || !heapXdi) | 
| 210 | 
      error(SYSTEM, "failed heap index allocation in kdT_BuildPhotonMap"); | 
| 211 | 
          | 
| 212 | 
   /* Initialize index arrays */ | 
| 213 | 
   for (i = 0; i < pmap -> numPhotons; i++) | 
| 214 | 
      heapXdi [i] = heapIdx [i] = i; | 
| 215 | 
       | 
| 216 | 
   /* Build kd-tree */ | 
| 217 | 
   kdT_Build(nodes, heapIdx, heapXdi, pmap -> minPos, pmap -> maxPos, 0,  | 
| 218 | 
             pmap -> numPhotons - 1, 1); | 
| 219 | 
                 | 
| 220 | 
   /* Cleanup */ | 
| 221 | 
   free(heapIdx); | 
| 222 | 
   free(heapXdi); | 
| 223 | 
} | 
| 224 | 
 | 
| 225 | 
 | 
| 226 | 
 | 
| 227 | 
int kdT_SavePhotons (const struct PhotonMap *pmap, FILE *out) | 
| 228 | 
{ | 
| 229 | 
   unsigned long        i, j; | 
| 230 | 
   Photon               *p = (Photon*)pmap -> store.nodes; | 
| 231 | 
    | 
| 232 | 
   for (i = 0; i < pmap -> numPhotons; i++, p++) { | 
| 233 | 
      /* Write photon attributes */ | 
| 234 | 
      for (j = 0; j < 3; j++) | 
| 235 | 
         putflt(p -> pos [j], out); | 
| 236 | 
          | 
| 237 | 
      /* Bytewise dump otherwise we have portability probs */ | 
| 238 | 
      for (j = 0; j < 3; j++)  | 
| 239 | 
         putint(p -> norm [j], 1, out); | 
| 240 | 
          | 
| 241 | 
#ifdef PMAP_FLOAT_FLUX | 
| 242 | 
      for (j = 0; j < 3; j++)  | 
| 243 | 
         putflt(p -> flux [j], out); | 
| 244 | 
#else | 
| 245 | 
      for (j = 0; j < 4; j++)  | 
| 246 | 
         putint(p -> flux [j], 1, out); | 
| 247 | 
#endif | 
| 248 | 
 | 
| 249 | 
      putint(p -> primary, sizeof(p -> primary), out); | 
| 250 | 
      putint(p -> flags, 1, out); | 
| 251 | 
             | 
| 252 | 
      if (ferror(out)) | 
| 253 | 
         return -1; | 
| 254 | 
   } | 
| 255 | 
    | 
| 256 | 
   return 0; | 
| 257 | 
} | 
| 258 | 
 | 
| 259 | 
 | 
| 260 | 
 | 
| 261 | 
int kdT_LoadPhotons (struct PhotonMap *pmap, FILE *in) | 
| 262 | 
{ | 
| 263 | 
   unsigned long  i, j; | 
| 264 | 
   Photon         *p; | 
| 265 | 
    | 
| 266 | 
   /* Allocate kd-tree based on initialised pmap -> numPhotons */ | 
| 267 | 
   pmap -> store.nodes = calloc(sizeof(Photon), pmap -> numPhotons); | 
| 268 | 
   if (!pmap -> store.nodes)  | 
| 269 | 
      error(SYSTEM, "failed kd-tree allocation in kdT_LoadPhotons"); | 
| 270 | 
       | 
| 271 | 
   /* Get photon attributes */    | 
| 272 | 
   for (i = 0, p = pmap -> store.nodes; i < pmap -> numPhotons; i++, p++) { | 
| 273 | 
      for (j = 0; j < 3; j++)  | 
| 274 | 
         p -> pos [j] = getflt(in); | 
| 275 | 
          | 
| 276 | 
      /* Bytewise grab otherwise we have portability probs */ | 
| 277 | 
      for (j = 0; j < 3; j++)  | 
| 278 | 
         p -> norm [j] = getint(1, in); | 
| 279 | 
          | 
| 280 | 
#ifdef PMAP_FLOAT_FLUX | 
| 281 | 
      for (j = 0; j < 3; j++)  | 
| 282 | 
         p -> flux [j] = getflt(in); | 
| 283 | 
#else       | 
| 284 | 
      for (j = 0; j < 4; j++)  | 
| 285 | 
         p -> flux [j] = getint(1, in); | 
| 286 | 
#endif | 
| 287 | 
 | 
| 288 | 
      p -> primary = getint(sizeof(p -> primary), in); | 
| 289 | 
      p -> flags = getint(1, in); | 
| 290 | 
       | 
| 291 | 
      if (feof(in)) | 
| 292 | 
         return -1; | 
| 293 | 
   } | 
| 294 | 
    | 
| 295 | 
   return 0; | 
| 296 | 
} | 
| 297 | 
 | 
| 298 | 
 | 
| 299 | 
 | 
| 300 | 
void kdT_InitFindPhotons (struct PhotonMap *pmap) | 
| 301 | 
{ | 
| 302 | 
   pmap -> squeue.len = pmap -> maxGather + 1; | 
| 303 | 
   pmap -> squeue.node = calloc(pmap -> squeue.len, | 
| 304 | 
                                sizeof(PhotonSearchQueueNode)); | 
| 305 | 
   if (!pmap -> squeue.node)  | 
| 306 | 
      error(SYSTEM, "can't allocate photon search queue"); | 
| 307 | 
} | 
| 308 | 
 | 
| 309 | 
 | 
| 310 | 
 | 
| 311 | 
static void kdT_FindNearest (PhotonMap *pmap, const float pos [3],  | 
| 312 | 
                             const float norm [3], unsigned long node) | 
| 313 | 
/* Recursive part of kdT_FindPhotons(). Locate pmap -> squeue.len nearest | 
| 314 | 
 * neighbours to pos with similar normal and return in search queue starting | 
| 315 | 
 * at pmap -> squeue.node.  Note that all heap and queue indices are | 
| 316 | 
 * 1-based, but accesses to the arrays are 0-based!  */ | 
| 317 | 
{ | 
| 318 | 
   Photon                  *p = (Photon*)pmap -> store.nodes + node - 1; | 
| 319 | 
   unsigned                i, j; | 
| 320 | 
   /* Signed distance to current photon's splitting plane */ | 
| 321 | 
   float                   d = pos [p -> discr] - p -> pos [p -> discr],  | 
| 322 | 
                           d2 = d * d, dv [3]; | 
| 323 | 
   PhotonSearchQueueNode*  sq = pmap -> squeue.node; | 
| 324 | 
   const unsigned          sqSize = pmap -> squeue.len; | 
| 325 | 
    | 
| 326 | 
   /* Search subtree closer to pos first; exclude other subtree if the  | 
| 327 | 
      distance to the splitting plane is greater than maxDist */ | 
| 328 | 
   if (d < 0) { | 
| 329 | 
      if (node << 1 <= pmap -> numPhotons)  | 
| 330 | 
         kdT_FindNearest(pmap, pos, norm, node << 1); | 
| 331 | 
          | 
| 332 | 
      if (d2 < pmap -> maxDist2 && node << 1 < pmap -> numPhotons)  | 
| 333 | 
         kdT_FindNearest(pmap, pos, norm, (node << 1) + 1); | 
| 334 | 
   } | 
| 335 | 
   else { | 
| 336 | 
      if (node << 1 < pmap -> numPhotons)  | 
| 337 | 
         kdT_FindNearest(pmap, pos, norm, (node << 1) + 1); | 
| 338 | 
          | 
| 339 | 
      if (d2 < pmap -> maxDist2 && node << 1 <= pmap -> numPhotons)  | 
| 340 | 
         kdT_FindNearest(pmap, pos, norm, node << 1); | 
| 341 | 
   } | 
| 342 | 
 | 
| 343 | 
   /* Reject photon if normal faces away (ignored for volume photons) with | 
| 344 | 
    * tolerance to account for perturbation; note photon normal is coded | 
| 345 | 
    * in range [-127,127], hence we factor this in */ | 
| 346 | 
   if (norm && DOT(norm, p -> norm) <= PMAP_NORM_TOL * 127 * frandom()) | 
| 347 | 
      return; | 
| 348 | 
       | 
| 349 | 
   if (isContribPmap(pmap) && pmap -> srcContrib) { | 
| 350 | 
      /* Lookup in contribution photon map */ | 
| 351 | 
      OBJREC *srcMod;  | 
| 352 | 
      const int srcIdx = photonSrcIdx(pmap, p); | 
| 353 | 
       | 
| 354 | 
      if (srcIdx < 0 || srcIdx >= nsources) | 
| 355 | 
         error(INTERNAL, "invalid light source index in photon map"); | 
| 356 | 
       | 
| 357 | 
      srcMod = findmaterial(source [srcIdx].so); | 
| 358 | 
 | 
| 359 | 
      /* Reject photon if contributions from light source which emitted it | 
| 360 | 
       * are not sought */ | 
| 361 | 
      if (!lu_find(pmap -> srcContrib, srcMod -> oname) -> data) | 
| 362 | 
         return; | 
| 363 | 
 | 
| 364 | 
      /* Reject non-caustic photon if lookup for caustic contribs */ | 
| 365 | 
      if (pmap -> lookupCaustic & !p -> caustic) | 
| 366 | 
         return; | 
| 367 | 
   } | 
| 368 | 
    | 
| 369 | 
   /* Squared distance to current photon (note dist2() requires doubles) */ | 
| 370 | 
   VSUB(dv, pos, p -> pos); | 
| 371 | 
   d2 = DOT(dv, dv); | 
| 372 | 
 | 
| 373 | 
   /* Accept photon if closer than current max dist & add to priority queue */ | 
| 374 | 
   if (d2 < pmap -> maxDist2) { | 
| 375 | 
      if (pmap -> squeue.tail < sqSize) { | 
| 376 | 
         /* Priority queue not full; append photon and restore heap */ | 
| 377 | 
         i = ++pmap -> squeue.tail; | 
| 378 | 
          | 
| 379 | 
         while (i > 1 && sq [(i >> 1) - 1].dist2 <= d2) { | 
| 380 | 
            sq [i - 1].idx    = sq [(i >> 1) - 1].idx; | 
| 381 | 
            sq [i - 1].dist2  = sq [(i >> 1) - 1].dist2; | 
| 382 | 
            i >>= 1; | 
| 383 | 
         } | 
| 384 | 
          | 
| 385 | 
         sq [--i].idx = (PhotonIdx)p; | 
| 386 | 
         sq [i].dist2 = d2; | 
| 387 | 
         /* Update maxDist if we've just filled the queue */ | 
| 388 | 
         if (pmap -> squeue.tail >= pmap -> squeue.len) | 
| 389 | 
            pmap -> maxDist2 = sq [0].dist2; | 
| 390 | 
      } | 
| 391 | 
      else { | 
| 392 | 
         /* Priority queue full; replace maximum, restore heap, and  | 
| 393 | 
            update maxDist */ | 
| 394 | 
         i = 1; | 
| 395 | 
          | 
| 396 | 
         while (i <= sqSize >> 1) { | 
| 397 | 
            j = i << 1; | 
| 398 | 
            if (j < sqSize && sq [j - 1].dist2 < sq [j].dist2)  | 
| 399 | 
               j++; | 
| 400 | 
            if (d2 >= sq [j - 1].dist2)  | 
| 401 | 
               break; | 
| 402 | 
            sq [i - 1].idx    = sq [j - 1].idx; | 
| 403 | 
            sq [i - 1].dist2  = sq [j - 1].dist2; | 
| 404 | 
            i = j; | 
| 405 | 
         } | 
| 406 | 
          | 
| 407 | 
         sq [--i].idx = (PhotonIdx)p; | 
| 408 | 
         sq [i].dist2 = d2; | 
| 409 | 
         pmap -> maxDist2 = sq [0].dist2; | 
| 410 | 
      } | 
| 411 | 
   } | 
| 412 | 
} | 
| 413 | 
 | 
| 414 | 
 | 
| 415 | 
 | 
| 416 | 
void kdT_FindPhotons (struct PhotonMap *pmap, const FVECT pos,  | 
| 417 | 
                      const FVECT norm) | 
| 418 | 
{ | 
| 419 | 
   float p [3], n [3]; | 
| 420 | 
    | 
| 421 | 
   /* Photon pos & normal stored at lower precision */ | 
| 422 | 
   VCOPY(p, pos); | 
| 423 | 
   VCOPY(n, norm); | 
| 424 | 
   kdT_FindNearest(pmap, p, n, 1); | 
| 425 | 
} | 
| 426 | 
 | 
| 427 | 
 | 
| 428 | 
 | 
| 429 | 
static void kdT_Find1Nearest (PhotonMap *pmap, const float pos [3],  | 
| 430 | 
                              const float norm [3], Photon **photon,  | 
| 431 | 
                              unsigned long node) | 
| 432 | 
/* Recursive part of kdT_Find1Photon().  Locate single nearest neighbour to | 
| 433 | 
 * pos with similar normal.  Note that all heap and queue indices are | 
| 434 | 
 * 1-based, but accesses to the arrays are 0-based!  */ | 
| 435 | 
{ | 
| 436 | 
   Photon   *p = (Photon*)pmap -> store.nodes + node - 1; | 
| 437 | 
   /* Signed distance to current photon's splitting plane */ | 
| 438 | 
   float    d  = pos [p -> discr] - p -> pos [p -> discr], d2 = d * d,  | 
| 439 | 
            dv [3]; | 
| 440 | 
    | 
| 441 | 
   /* Search subtree closer to pos first; exclude other subtree if the  | 
| 442 | 
      distance to the splitting plane is greater than maxDist */ | 
| 443 | 
   if (d < 0) { | 
| 444 | 
      if (node << 1 <= pmap -> numPhotons)  | 
| 445 | 
         kdT_Find1Nearest(pmap, pos, norm, photon, node << 1); | 
| 446 | 
          | 
| 447 | 
      if (d2 < pmap -> maxDist2 && node << 1 < pmap -> numPhotons)  | 
| 448 | 
         kdT_Find1Nearest(pmap, pos, norm, photon, (node << 1) + 1); | 
| 449 | 
   } | 
| 450 | 
   else { | 
| 451 | 
      if (node << 1 < pmap -> numPhotons)  | 
| 452 | 
         kdT_Find1Nearest(pmap, pos, norm, photon, (node << 1) + 1); | 
| 453 | 
          | 
| 454 | 
      if (d2 < pmap -> maxDist2 && node << 1 <= pmap -> numPhotons)  | 
| 455 | 
         kdT_Find1Nearest(pmap, pos, norm, photon, node << 1); | 
| 456 | 
   } | 
| 457 | 
    | 
| 458 | 
   /* Squared distance to current photon */ | 
| 459 | 
   VSUB(dv, pos, p -> pos); | 
| 460 | 
   d2 = DOT(dv, dv); | 
| 461 | 
    | 
| 462 | 
   if (d2 < pmap -> maxDist2 &&  | 
| 463 | 
       DOT(norm, p -> norm) > PMAP_NORM_TOL * 127 * frandom()) { | 
| 464 | 
      /* Closest photon so far with similar normal. We allow for tolerance | 
| 465 | 
       * to account for perturbation in the latter; note the photon normal | 
| 466 | 
       * is coded in the range [-127,127], hence we factor this in  */ | 
| 467 | 
      pmap -> maxDist2 = d2; | 
| 468 | 
      *photon = p; | 
| 469 | 
   } | 
| 470 | 
} | 
| 471 | 
 | 
| 472 | 
 | 
| 473 | 
 | 
| 474 | 
void kdT_Find1Photon (struct PhotonMap *pmap, const FVECT pos,  | 
| 475 | 
                      const FVECT norm, Photon *photon) | 
| 476 | 
{ | 
| 477 | 
   float    p [3], n [3]; | 
| 478 | 
   Photon   *pnn; | 
| 479 | 
    | 
| 480 | 
   /* Photon pos & normal stored at lower precision */ | 
| 481 | 
   VCOPY(p, pos); | 
| 482 | 
   VCOPY(n, norm); | 
| 483 | 
   kdT_Find1Nearest(pmap, p, n, &pnn, 1); | 
| 484 | 
   memcpy(photon, pnn, sizeof(Photon)); | 
| 485 | 
} | 
| 486 | 
 | 
| 487 | 
 | 
| 488 | 
 | 
| 489 | 
int kdT_GetPhoton (const struct PhotonMap *pmap, PhotonIdx idx, | 
| 490 | 
                   Photon *photon) | 
| 491 | 
{ | 
| 492 | 
   memcpy(photon, idx, sizeof(Photon));    | 
| 493 | 
   return 0; | 
| 494 | 
} | 
| 495 | 
 | 
| 496 | 
 | 
| 497 | 
 | 
| 498 | 
Photon *kdT_GetNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx) | 
| 499 | 
{ | 
| 500 | 
   return idx; | 
| 501 | 
} | 
| 502 | 
 | 
| 503 | 
 | 
| 504 | 
 | 
| 505 | 
PhotonIdx kdT_FirstPhoton (const struct PhotonMap* pmap) | 
| 506 | 
{ | 
| 507 | 
   return pmap -> store.nodes; | 
| 508 | 
} | 
| 509 | 
 | 
| 510 | 
 | 
| 511 | 
 | 
| 512 | 
void kdT_Delete (PhotonKdTree *kdt) | 
| 513 | 
{ | 
| 514 | 
   free(kdt -> nodes); | 
| 515 | 
   kdt -> nodes = NULL; | 
| 516 | 
} |