--- ray/src/rt/pmutil.c	2016/09/26 20:19:30	2.1
+++ ray/src/rt/pmutil.c	2025/01/10 21:43:22	2.7
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: pmutil.c,v 2.1 2016/09/26 20:19:30 greg Exp $";
+static const char RCSid[] = "$Id: pmutil.c,v 2.7 2025/01/10 21:43:22 greg Exp $";
 #endif
 
 /* 
@@ -12,7 +12,7 @@ static const char RCSid[] = "$Id: pmutil.c,v 2.1 2016/
        supported by the Swiss National Science Foundation (SNSF, #147053)
    ======================================================================
    
-   $Id: pmutil.c,v 2.1 2016/09/26 20:19:30 greg Exp $
+   $Id: pmutil.c,v 2.7 2025/01/10 21:43:22 greg Exp $
 */
 
 #include "pmap.h"
@@ -28,7 +28,7 @@ extern char *octname;
 /* Photon map lookup functions per type */
 void (*pmapLookup [NUM_PMAP_TYPES])(PhotonMap*, RAY*, COLOR) = {
    photonDensity, photonPreCompDensity, photonDensity, volumePhotonDensity,
-   photonDensity, NULL
+   photonDensity, photonDensity
 };
 
 
@@ -56,9 +56,9 @@ void loadPmaps (PhotonMap **pmaps, const PhotonMapPara
    struct stat octstat, pmstat;
    PhotonMap *pm;
    PhotonMapType type;
-   
+
    for (t = 0; t < NUM_PMAP_TYPES; t++)
-      if (setPmapParam(&pm, parm + t)) {         
+      if (setPmapParam(&pm, parm + t)) {
          /* Check if photon map newer than octree */
          if (pm -> fileName && octname &&
              !stat(pm -> fileName, &pmstat) && !stat(octname, &octstat) && 
@@ -67,7 +67,7 @@ void loadPmaps (PhotonMap **pmaps, const PhotonMapPara
                     pm -> fileName);
             error(USER, errmsg);
          }
-         
+
          /* Load photon map from file and get its type */
          if ((type = loadPhotonMap(pm, pm -> fileName)) == PMAP_TYPE_NONE)
             error(USER, "failed loading photon map");
@@ -75,15 +75,44 @@ void loadPmaps (PhotonMap **pmaps, const PhotonMapPara
          /* Assign to appropriate photon map type (deleting previously
           * loaded photon map of same type if necessary) */
          if (pmaps [type]) {
+            sprintf(errmsg, "multiple %s photon maps, dropping previous",
+                    pmapName [type]);
+            error(WARNING, errmsg);
             deletePhotons(pmaps [type]);
             free(pmaps [type]);
          }
          pmaps [type] = pm;
-         
-         /* Check for invalid density estimate bandwidth */                            
+
+         /* Check for valid density estimate bandwidths */
+         if ((pm -> minGather > 1 || pm -> maxGather > 1) &&
+             (type == PMAP_TYPE_PRECOMP)) {
+            /* Force bwidth to 1 for precomputed pmap */
+            error(WARNING, "ignoring bandwidth for precomp photon map");
+            pm -> minGather = pm -> maxGather = 1;
+         }
+
+         if ((pm -> maxGather > pm -> minGather) && 
+             (type == PMAP_TYPE_VOLUME)) {            
+            /* Biascomp for volume pmaps (see volumePhotonDensity() below) 
+               is considered redundant, and there's probably no point in 
+               recovering by using the lower bandwidth, since it's probably
+               not what the user wants, so bail out. */
+            sprintf(errmsg, 
+                    "bias compensation is not available with %s photon maps",
+                    pmapName [type]);
+            error(USER, errmsg);
+         }
+
          if (pm -> maxGather > pm -> numPhotons) {
-            error(WARNING, "adjusting density estimate bandwidth");
-            pm -> minGather = pm -> maxGather = pm -> numPhotons;
+            /* Clamp lookup bandwidth to total number of photons (minus one,
+               since density estimate gets extra photon to obtain averaged 
+               radius) */
+            sprintf(
+               errmsg, "clamping density estimate bandwidth to %ld",
+               pm -> numPhotons
+            );
+            error(WARNING, errmsg);
+            pm -> minGather = pm -> maxGather = pm -> numPhotons - 1;
          }
       }
 }
@@ -98,6 +127,7 @@ void cleanUpPmaps (PhotonMap **pmaps)
       if (pmaps [t]) {
          deletePhotons(pmaps [t]);
          free(pmaps [t]);
+         pmaps [t] = NULL;
       }
    }
 }
@@ -109,7 +139,7 @@ void photonDensity (PhotonMap *pmap, RAY *ray, COLOR i
 /* Photon density estimate. Returns irradiance at ray -> rop. */
 {
    unsigned                      i;
-   float                         r;
+   float                         r2;
    COLOR                         flux;
    Photon                        *photon;
    const PhotonSearchQueueNode   *sqn;
@@ -141,9 +171,9 @@ void photonDensity (PhotonMap *pmap, RAY *ray, COLOR i
       /* No bias compensation. Just do a plain vanilla estimate */
       sqn = pmap -> squeue.node + 1;
       
-      /* Average radius between furthest two photons to improve accuracy */      
-      r = max(sqn -> dist2, (sqn + 1) -> dist2);
-      r = 0.25 * (pmap -> maxDist2 + r + 2 * sqrt(pmap -> maxDist2 * r));   
+      /* Average radius^2 between furthest two photons to improve accuracy */
+      r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
+      r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));
       
       /* Skip the extra photon */
       for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) {
@@ -151,14 +181,14 @@ void photonDensity (PhotonMap *pmap, RAY *ray, COLOR i
          getPhotonFlux(photon, flux);         
 #ifdef PMAP_EPANECHNIKOV
          /* Apply Epanechnikov kernel to photon flux based on photon dist */
-         scalecolor(flux, 2 * (1 - sqn -> dist2 / r));
+         scalecolor(flux, 2 * (1 - sqn -> dist2 / r2));
 #endif   
          addcolor(irrad, flux);
       }
       
       /* Divide by search area PI * r^2, 1 / PI required as ambient 
          normalisation factor */         
-      scalecolor(irrad, 1 / (PI * PI * r)); 
+      scalecolor(irrad, 1 / (PI * PI * r2)); 
       
       return;
    }
@@ -181,8 +211,11 @@ void photonPreCompDensity (PhotonMap *pmap, RAY *r, CO
    if (r -> ro && islight(objptr(r -> ro -> omod) -> otype)) 
       return;
       
-   find1Photon(preCompPmap, r, &p);
-   getPhotonFlux(&p, irrad);
+   if (find1Photon(preCompPmap, r, &p))
+      /* p contains a found photon, so get its irradiance, otherwise it
+       * remains zero under the assumption all photons are too distant
+       * to contribute significantly */ 
+      getPhotonFlux(&p, irrad);
 }
 
 
@@ -192,7 +225,7 @@ void volumePhotonDensity (PhotonMap *pmap, RAY *ray, C
 /* Photon volume density estimate. Returns irradiance at ray -> rop. */
 {
    unsigned                      i;
-   float                         r, gecc2, ph;
+   float                         r2, gecc2, ph;
    COLOR                         flux;
    Photon                        *photon;
    const PhotonSearchQueueNode   *sqn;
@@ -217,9 +250,9 @@ void volumePhotonDensity (PhotonMap *pmap, RAY *ray, C
       gecc2 = ray -> gecc * ray -> gecc;
       sqn = pmap -> squeue.node + 1;
       
-      /* Average radius between furthest two photons to improve accuracy */      
-      r = max(sqn -> dist2, (sqn + 1) -> dist2);
-      r = 0.25 * (pmap -> maxDist2 + r + 2 * sqrt(pmap -> maxDist2 * r));   
+      /* Average radius^2 between furthest two photons to improve accuracy */      
+      r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
+      r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));
       
       /* Skip the extra photon */
       for (i = 1; i < pmap -> squeue.tail; i++, sqn++) {
@@ -241,7 +274,7 @@ void volumePhotonDensity (PhotonMap *pmap, RAY *ray, C
       
       /* Divide by search volume 4 / 3 * PI * r^3 and phase function
          normalization factor 1 / (4 * PI) */
-      scalecolor(irrad, 3 / (16 * PI * PI * r * sqrt(r)));
+      scalecolor(irrad, 3 / (16 * PI * PI * r2 * sqrt(r2)));
       return;
    }
 #if 0