--- ray/src/rt/pmutil.c	2018/01/24 19:39:05	2.2
+++ ray/src/rt/pmutil.c	2020/04/08 15:14:21	2.5
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: pmutil.c,v 2.2 2018/01/24 19:39:05 rschregle Exp $";
+static const char RCSid[] = "$Id: pmutil.c,v 2.5 2020/04/08 15:14:21 rschregle Exp $";
 #endif
 
 /* 
@@ -12,7 +12,7 @@ static const char RCSid[] = "$Id: pmutil.c,v 2.2 2018/
        supported by the Swiss National Science Foundation (SNSF, #147053)
    ======================================================================
    
-   $Id: pmutil.c,v 2.2 2018/01/24 19:39:05 rschregle Exp $
+   $Id: pmutil.c,v 2.5 2020/04/08 15:14:21 rschregle Exp $
 */
 
 #include "pmap.h"
@@ -75,16 +75,38 @@ 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;
-         }
+         }            
       }
 }
 
@@ -109,7 +131,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 +163,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 +173,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 +203,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 +217,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 +242,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 +266,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