Stochastic mipmap interpolation (AcademySoftwareFoundation#3127)

Suggested by Luke Emrose, the idea here is that we introduce two new
MIP interpolation modes: StochasticTrilinear and StochasticAniso.
These work similarly to Trilinear and Aniso, respectively, except that
instead of sampling both of the bracketing MIP levels and blending
between them, in stochastic mode we use a random variate to select ONE
of the two MIP levels with probability proportional to the computed
blending weight.

If these mip modes are used, the caller needs to pass a
well-stratified value in the new "rnd" field of TextureOpt. Also note
that the method only really makes sense in the context of a renderer
that is using many samples per pixel. As a single texture lookup, it
looks horrible and noisy, as you can see by the testshade-based tests
we do of this, which are 1spp and looks pretty rough.

My benchmarks (in a hard stress test of the texture system using
`testtex --threadtimes 8`) show that this approach has speedups of
35-40% when using between 1-4 threads, dropping to more like 25-30%
speedup as threads are in the 12-32 range (that makes sense, it's
doing the same amount of total I/O).

In this first implementation, in order to preserve ABI compatibility
if we decide to backport this feature to the release branch,
TextureOpt::rnd is added as a union co-member with the unused "bias"
field. For TextureOptBatch, we had to add it outright, but it's
guarded by an `#if` and will not be backported to 2.3. (We may change
our minds and not backport it at all, I just wanted to leave the
option open with this initial patch. Also, there will probably be
additional further ABI-breaking changes coming later that will
definitely not be backported.)

In this first phase, we only use the stochastic component to turn the
2-MIP-level lookup into one level lookup.  It has been further
suggested that we could apply stochastic techniques to anisotropic
filtering (instead of many samples spaced across the major axis of the
filter footprint, take one sample distributed across the axis), and
even for the ordinary bilinear or bicubic within each level (instead
of 4 or 16 weighted texels, just retrieve one texel chosen with the
weights as probabilities). But these deserve to be benchmarked and
have their quality evaluated separtely and deliberately, and so will
be relegated to a later investigation and implementation.

Author: lgritz

src/cmake/testing.cmake

@@ -161,6 +161,7 @@ macro (oiio_add_all_tests)
                     misnamed-file
                     texture-crop texture-cropover
                     texture-filtersize
+                    texture-filtersize-stochastic
                     texture-overscan
                     texture-wrapfill
                     texture-res texture-maxres
@@ -191,6 +192,8 @@ macro (oiio_add_all_tests)
                     texture-interp-closest
                     texture-mip-nomip texture-mip-onelevel
                     texture-mip-trilinear
+                    texture-mip-stochastictrilinear
+                    texture-mip-stochasticaniso
                     texture-missing
                     texture-pointsample
                     texture-udim texture-udim2

src/include/OpenImageIO/texture.h

@@ -20,6 +20,7 @@
 // features are supported.
 #define OIIO_TEXTURESYSTEM_SUPPORTS_CLOSE 1
 
+#define OIIO_TEXTURESYSTEM_SUPPORTS_STOCHASTIC 1
 
 
 OIIO_NAMESPACE_BEGIN
@@ -93,7 +94,9 @@ enum class MipMode {
     NoMIP,      ///< Just use highest-res image, no MIP mapping
     OneLevel,   ///< Use just one mipmap level
     Trilinear,  ///< Use two MIPmap levels (trilinear)
-    Aniso       ///< Use two MIPmap levels w/ anisotropic
+    Aniso,      ///< Use two MIPmap levels w/ anisotropic
+    StochasticTrilinear, ///< Stochastic trilinear
+    StochasticAniso, ///< Stochastic anisotropic
 };
 
 /// Interp mode determines how we sample within a mipmap level
@@ -102,7 +105,7 @@ enum class InterpMode {
     Closest,      ///< Force closest texel
     Bilinear,     ///< Force bilinear lookup within a mip level
     Bicubic,      ///< Force cubic lookup within a mip level
-    SmartBicubic  ///< Bicubic when maxifying, else bilinear
+    SmartBicubic  ///< Bicubic when magnifying, else bilinear
 };
 
 
@@ -190,7 +193,9 @@ class OIIO_API TextureOpt {
         MipModeNoMIP,      ///< Just use highest-res image, no MIP mapping
         MipModeOneLevel,   ///< Use just one mipmap level
         MipModeTrilinear,  ///< Use two MIPmap levels (trilinear)
-        MipModeAniso       ///< Use two MIPmap levels w/ anisotropic
+        MipModeAniso,      ///< Use two MIPmap levels w/ anisotropic
+        MipModeStochasticTrilinear, ///< Stochastic trilinear
+        MipModeStochasticAniso, ///< Stochastic anisotropic
     };
 
     /// Interp mode determines how we sample within a mipmap level
@@ -199,7 +204,7 @@ class OIIO_API TextureOpt {
         InterpClosest,      ///< Force closest texel
         InterpBilinear,     ///< Force bilinear lookup within a mip level
         InterpBicubic,      ///< Force cubic lookup within a mip level
-        InterpSmartBicubic  ///< Bicubic when maxifying, else bilinear
+        InterpSmartBicubic  ///< Bicubic when magnifying, else bilinear
     };
 
 
@@ -212,10 +217,8 @@ class OIIO_API TextureOpt {
         anisotropic(32), conservative_filter(true),
         sblur(0.0f), tblur(0.0f), swidth(1.0f), twidth(1.0f),
         fill(0.0f), missingcolor(nullptr),
-        // dresultds(nullptr), dresultdt(nullptr),
-        time(0.0f), bias(0.0f), samples(1),
-        rwrap(WrapDefault), rblur(0.0f), rwidth(1.0f), // dresultdr(nullptr),
-        // actualchannels(0),
+        time(0.0f), rnd(0.0f), samples(1),
+        rwrap(WrapDefault), rblur(0.0f), rwidth(1.0f),
         envlayout(0)
     { }
 
@@ -237,7 +240,10 @@ class OIIO_API TextureOpt {
     float fill;                 ///< Fill value for missing channels
     const float* missingcolor;  ///< Color for missing texture
     float time;                 ///< Time (for time-dependent texture lookups)
-    float bias;                 ///< Bias for shadows
+    union {
+        float bias;                 ///< Bias for shadows (DEPRECATED?)
+        float rnd;                  ///< Stratified sample value
+    };
     int samples;                ///< Number of samples for shadows
 
     // For 3D volume texture lookups only:
@@ -291,6 +297,9 @@ class OIIO_API TextureOptBatch {
     alignas(Tex::BatchAlign) float twidth[Tex::BatchWidth];
     alignas(Tex::BatchAlign) float rwidth[Tex::BatchWidth];
     // Note: rblur,rwidth only used for volumetric lookups
+#if OIIO_VERSION_GREATER_EQUAL(2,4,0)
+    alignas(Tex::BatchAlign) float rnd[Tex::BatchWidth];
+#endif
 
     // Options that must be the same for all points we're texturing at once
     int firstchannel = 0;                 ///< First channel of the lookup
@@ -354,15 +363,17 @@ class OIIO_API TextureOptions {
         InterpClosest,      ///< Force closest texel
         InterpBilinear,     ///< Force bilinear lookup within a mip level
         InterpBicubic,      ///< Force cubic lookup within a mip level
-        InterpSmartBicubic  ///< Bicubic when maxifying, else bilinear
+        InterpSmartBicubic  ///< Bicubic when magnifying, else bilinear
     };
 
     /// Create a TextureOptions with all fields initialized to reasonable
     /// defaults.
+    OIIO_DEPRECATED("no longer used since OIIO 1.8")
     TextureOptions();
 
     /// Convert a TextureOpt for one point into a TextureOptions with
     /// uniform values.
+    OIIO_DEPRECATED("no longer used since OIIO 1.8")
     TextureOptions(const TextureOpt& opt);
 
     // Options that must be the same for all points we're texturing at once

src/libtexture/environment.cpp

@@ -447,7 +447,8 @@ TextureSystemImpl::environment(TextureHandle* texture_handle_,
 
     TextureOpt::MipMode mipmode = options.mipmode;
     bool aniso                  = (mipmode == TextureOpt::MipModeDefault
-                  || mipmode == TextureOpt::MipModeAniso);
+                  || mipmode == TextureOpt::MipModeAniso
+                  || mipmode == TextureOpt::MipModeStochasticAniso);
 
     float aspect, trueaspect, filtwidth;
     int nsamples;
@@ -612,6 +613,7 @@ TextureSystemImpl::environment(TextureHandle* texture_handle,
             opt.tblur  = options.tblur[i];
             opt.swidth = options.swidth[i];
             opt.twidth = options.twidth[i];
+            opt.rnd    = options.rnd[i];
             Imath::V3f R_(R[i], R[i + Tex::BatchWidth],
                           R[i + 2 * Tex::BatchWidth]);
             Imath::V3f dRdx_(dRdx[i], dRdx[i + Tex::BatchWidth],

src/libtexture/texture3d.cpp

@@ -112,6 +112,8 @@ TextureSystemImpl::texture3d(TextureHandle* texture_handle_,
         &TextureSystemImpl::texture3d_lookup_nomip,
         &TextureSystemImpl::texture3d_lookup_trilinear_mipmap,
         &TextureSystemImpl::texture3d_lookup_trilinear_mipmap,
+        &TextureSystemImpl::texture3d_lookup,
+        &TextureSystemImpl::texture3d_lookup_trilinear_mipmap,
         &TextureSystemImpl::texture3d_lookup
     };
     texture3d_lookup_prototype lookup = lookup_functions[(int)options.mipmode];

src/libtexture/texturesys.cpp

@@ -1085,6 +1085,8 @@ TextureSystemImpl::texture(TextureHandle* texture_handle_,
         &TextureSystemImpl::texture_lookup_nomip,
         &TextureSystemImpl::texture_lookup_trilinear_mipmap,
         &TextureSystemImpl::texture_lookup_trilinear_mipmap,
+        &TextureSystemImpl::texture_lookup,
+        &TextureSystemImpl::texture_lookup_trilinear_mipmap,
         &TextureSystemImpl::texture_lookup
     };
     texture_lookup_prototype lookup = lookup_functions[(int)options.mipmode];
@@ -1273,6 +1275,7 @@ TextureSystemImpl::texture(TextureHandle* texture_handle,
             opt.tblur  = options.tblur[i];
             opt.swidth = options.swidth[i];
             opt.twidth = options.twidth[i];
+            opt.rnd    = options.rnd[i];
             // rblur, rwidth not needed for 2D texture
             if (dresultds) {
                 ok &= texture(texture_handle, thread_info, opt, s[i], t[i],
@@ -1482,6 +1485,10 @@ compute_miplevels(TextureSystemImpl::TextureFile& texturefile,
         float filtwidth_ras = minorlength
                               * std::min(subinfo.spec(m).width,
                                          subinfo.spec(m).height);
+        // FIXME: We should store the min(width,height) of each level directly
+        // in an array in subinfo, so we're not rifling through the specs
+        // and taking mins in this loop every single time.
+
         // Once the filter width is smaller than one texel at this level,
         // we've gone too far, so we know that we want to interpolate the
         // previous level and the current level.  Note that filtwidth_ras
@@ -1521,6 +1528,18 @@ compute_miplevels(TextureSystemImpl::TextureFile& texturefile,
         miplevel[0] = miplevel[1];
         levelblend  = 0;
     }
+    if (options.mipmode == TextureOpt::MipModeStochasticTrilinear
+        || options.mipmode == TextureOpt::MipModeStochasticAniso) {
+        // If using stochastic sampling, the random deviate is a threshold
+        // versus the levelblend to determine which ONE of the two MIP
+        // levels to use.
+        if (options.rnd > levelblend) {
+            miplevel[1] = miplevel[0];
+        } else {
+            miplevel[0] = miplevel[1];
+        }
+        levelblend = 0;
+    }
     levelweight[0] = 1.0f - levelblend;
     levelweight[1] = levelblend;
 }

src/testtex/testtex.cpp

@@ -42,6 +42,7 @@ static std::string dataformatname = "half";
 static float sscale = 1, tscale = 1;
 static float sblur = 0, tblur = -1;
 static float width       = 1;
+static float widthramp   = 0;
 static float anisoaspect = 1.0;  // anisotropic aspect ratio
 static std::string wrapmodes("periodic");
 static int anisomax           = TextureOpt().anisotropic;
@@ -138,6 +139,8 @@ getargs(int argc, const char* argv[])
       .help("Add blur (s, t) to texture lookup");
     ap.arg("--width %f:WIDTH", &width)
       .help("Multiply filter width of texture lookup");
+    ap.arg("--widthramp %f:WIDTH", &widthramp)
+      .help("If set, ramp to this width on the right side");
     ap.arg("--fill %f:FILLVAL", &fill)
       .help("Set fill value for missing channels");
     ap.arg("--wrap %s:MODE", &wrapmodes)
@@ -578,6 +581,20 @@ plain_tex_region(ImageBuf& image, ustring filename, Mapping2D mapping,
         float s, t, dsdx, dtdx, dsdy, dtdy;
         mapping(p.x(), p.y(), s, t, dsdx, dtdx, dsdy, dtdy);
 
+        if (widthramp != 0.0f) {
+            // If widthramp is set, we want to blend between the set width
+            // and the ramp width from left to right.
+            opt.swidth = OIIO::lerp(width, widthramp, s);
+            opt.twidth = opt.swidth;
+        }
+        if (mipmode == TextureOpt::MipModeStochasticTrilinear
+            || mipmode == TextureOpt::MipModeStochasticAniso) {
+            // Hash the pixel coords to get a pseudo-random variant
+            constexpr float inv = 1.0f
+                                  / float(std::numeric_limits<uint32_t>::max());
+            opt.rnd = bjhash::bjfinal(p.x(), p.y()) * inv;
+        }
+
         // Call the texture system to do the filtering.
         bool ok;
         if (use_handle)
@@ -694,6 +711,26 @@ plain_tex_region_batch(ImageBuf& image, ustring filename, Mapping2DWide mapping,
         for (int x = roi.xbegin; x < roi.xend; x += BatchWidth) {
             FloatWide s, t, dsdx, dtdx, dsdy, dtdy;
             mapping(IntWide::Iota(x), y, s, t, dsdx, dtdx, dsdy, dtdy);
+
+            if (widthramp != 0.0f) {
+                // If widthramp is set, we want to blend between the set width
+                // and the ramp width from left to right.
+                for (int i = 0; i < BatchWidth; ++i) {
+                    opt.swidth[i] = OIIO::lerp(width, widthramp, s[i]);
+                    opt.twidth[i] = opt.swidth[i];
+                }
+            }
+            if (mipmode == TextureOpt::MipModeStochasticTrilinear
+                || mipmode == TextureOpt::MipModeStochasticAniso) {
+                // Hash the pixel coords to get a pseudo-random variant
+#if OIIO_VERSION_GREATER_EQUAL(2, 4, 0)
+                constexpr float inv
+                    = 1.0f / float(std::numeric_limits<uint32_t>::max());
+                for (int i = 0; i < BatchWidth; ++i)
+                    opt.rnd[i] = bjhash::bjfinal(x + i, y) * inv;
+#endif
+            }
+
             int npoints  = std::min(BatchWidth, roi.xend - x);
             RunMask mask = RunMaskOn >> (BatchWidth - npoints);
             // Call the texture system to do the filtering.

testsuite/texture-filtersize-stochastic/ref/out.tif

@@ -0,0 +1,12 @@
+#!/usr/bin/env python
+
+# This test just views the "miplevels" texture straight on, but it uses
+# -widthramp to smoothly blend between mipmap levels from left to right
+# (wanting the 256^2 level at the left and the 64^2 level at the right).
+# Using MIPmode StochasticAniso and pseudo-random variate, this tests that
+# we are blending between the levels correctly.
+
+
+command = testtex_command (OIIO_TESTSUITE_IMAGEDIR + "/miplevels.tx",
+                           " -nowarp -res 256 256 -mipmode 6 -widthramp 4 -d uint8 -o out.tif")
+outputs = [ "out.tif" ]

testsuite/texture-filtersize-stochastic/run.py

@@ -0,0 +1,5 @@
+#!/usr/bin/env python
+
+command = testtex_command ("../common/textures/grid.tx",
+                           extraargs = "-mipmode 6 -d uint8 -o out.tif")
+outputs = [ "out.tif" ]

testsuite/texture-mip-stochasticaniso/ref/out.tif

@@ -0,0 +1,5 @@
+#!/usr/bin/env python
+
+command = testtex_command ("../common/textures/grid.tx",
+                           extraargs = "-mipmode 5 -d uint8 -o out.tif")
+outputs = [ "out.tif" ]