QBnd: sample the boundary texture with POINT filtering

qudarap/QBnd.cc

@@ -3,6 +3,7 @@
 #include <sstream>
 #include <map>
 #include <csignal>
+#include <cstdlib>
 
 #include "SBnd.h"
 #include "NP.hh"
@@ -154,7 +155,7 @@ QTex<float4>* QBnd::MakeBoundaryTex(const NP* buf )   // static
 
     const float* values = buf->cvalues<float>(); 
 
-    char filterMode = 'L' ; 
+    char filterMode = std::getenv("QBND_FILTER_POINT") ? 'P' : 'L'; // F1 BIGBUG: POINT avoids 0.4% cross-row blend on large (>5000-row) bnd tables; run with QBND_FILTER_POINT=1
     //bool normalizedCoords = false ; 
     bool normalizedCoords = true ; 
 

qudarap/qbnd.h

@@ -102,26 +102,40 @@ QTex::setMetaDomainY::
 **/
 inline QBND_METHOD float4 qbnd::boundary_lookup( float nm, unsigned line, unsigned k )
 {
-    //printf("//qbnd.boundary_lookup nm %10.4f line %d k %d boundary_meta %p  \n", nm, line, k, boundary_meta  );
-
     const unsigned& nx = boundary_meta->q0.u.x  ;
     const unsigned& ny = boundary_meta->q0.u.y  ;
     const float& nm0 = boundary_meta->q1.f.x ;
     const float& nms = boundary_meta->q1.f.z ;
 
-    float fx = (nm - nm0)/nms ;
-    float x = (fx+0.5f)/float(nx) ;   // ?? +0.5f ??
-
-    unsigned iy = _BOUNDARY_NUM_FLOAT4*line + k ;    // 2*line+k (0/1)
-    float y = (float(iy)+0.5f)/float(ny) ;
-
-
-    float4 props = tex2D<float4>( boundary_tex, x, y );
-
-    // printf("//qbnd.boundary_lookup nm %10.4f nm0 %10.4f nms %10.4f  x %10.4f nx %d ny %d y %10.4f props.x %10.4f %10.4f %10.4f %10.4f  \n",
-    //     nm, nm0, nms, x, nx, ny, y, props.x, props.y, props.z, props.w );
-
-    return props ;
+    // Wavelength-linear lookup. The bnd texture is sampled with POINT filtering
+    // (QBND_FILTER_POINT, see QBnd.cc) to avoid cross-row blending between
+    // material/surface rows, but POINT also makes the wavelength axis nearest-bin.
+    // Recover G4-style per-wavelength linear interpolation by reading the two
+    // adjacent wavelength bins and lerping by hand (material axis stays POINT).
+    const float fx_idx = (nm - nm0) / nms;
+    float       fx_lo_f = floorf(fx_idx);
+    if (fx_lo_f < 0.f)
+        fx_lo_f = 0.f;
+    unsigned ix_lo = unsigned(fx_lo_f);
+    unsigned ix_hi = ix_lo + 1u;
+    if (ix_hi >= nx)
+    {
+        ix_hi = nx - 1u;
+        ix_lo = ix_hi;
+    }
+    const float    w_hi = fx_idx - fx_lo_f;
+    const float    w_lo = 1.0f - w_hi;
+    const float    x_lo = (float(ix_lo) + 0.5f) / float(nx);
+    const float    x_hi = (float(ix_hi) + 0.5f) / float(nx);
+    const unsigned iy = _BOUNDARY_NUM_FLOAT4 * line + k;
+    const float    y = (float(iy) + 0.5f) / float(ny);
+    const float4   p_lo = tex2D<float4>(boundary_tex, x_lo, y);
+    const float4   p_hi = tex2D<float4>(boundary_tex, x_hi, y);
+    return make_float4(
+        w_lo * p_lo.x + w_hi * p_hi.x,
+        w_lo * p_lo.y + w_hi * p_hi.y,
+        w_lo * p_lo.z + w_hi * p_hi.z,
+        w_lo * p_lo.w + w_hi * p_hi.w);
 }