n/b/patch: add blur width parameter

noise/basis/patch.hpp

@@ -17,12 +17,13 @@
 #ifndef __UTIL_NOISE_BASIS_PATCH_HPP
 #define __UTIL_NOISE_BASIS_PATCH_HPP
 
+#include "../basis.hpp"
 #include "../../point.hpp"
 
 namespace util { namespace noise { namespace basis {
     template <typename T>
     struct patch {
-        patch (seed_t);
+        patch (seed_t, T width = 0);
         
         range<T> bounds (void) const;
         T operator() (point2<T>) const;
@@ -30,10 +31,17 @@ namespace util { namespace noise { namespace basis {
         seed_t seed (void) const;
         seed_t seed (seed_t);
 
+        T width (void) const;
+        T width (T);
+
     private:
         point2<T> centroid (util::point2i) const;
         T generate (util::point2i) const;
 
+        static constexpr T THRESHOLD = 1 - T(0.999);
+
+        T m_width;
+        T m_power;
         seed_t m_seed;
     };
 } } }

noise/basis/patch.ipp

@@ -19,28 +19,37 @@
 #endif
 #define __UTIL_NOISE_BASIS_PATCH_IPP
 
+#include "../../types.hpp"
+#include "../../ray.hpp"
+#include "../../vector.hpp"
 
-///////////////////////////////////////////////////////////////////////////////
+
-template <typename T>
+namespace util { namespace noise { namespace basis {
-util::noise::basis::patch<T>::patch (seed_t _seed):
+    ///////////////////////////////////////////////////////////////////////////
+    template <typename T>
+    patch<T>::patch (seed_t _seed, T _width):
+        m_width (_width),
+        m_power (exactly_zero (_width)
+                 ? std::numeric_limits<T>::infinity ()
+                 : std::log (THRESHOLD) / std::log (1 - _width)),
         m_seed (_seed)
-{ ; }
+    { ; }
 
 
-///////////////////////////////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////////
-template <typename T>
+    template <typename T>
-util::range<T>
+    range<T>
-util::noise::basis::patch<T>::bounds (void) const
+    patch<T>::bounds (void) const
-{
+    {
         return { T{0}, T{1} };
-}
+    }
 
 
-///////////////////////////////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////////
-template <typename T>
+    template <typename T>
-T
+    T
-util::noise::basis::patch<T>::operator () (point2<T> p) const
+    patch<T>::operator () (point2<T> p) const
-{
+    {
         // extract integer and fractional parts. be careful to always round down
         // (particularly with negatives) and avoid rounding errors.
         auto p_int = p.template cast<intmax_t> ();
@@ -48,48 +57,108 @@ util::noise::basis::patch<T>::operator () (point2<T> p) const
         if (p.y < 0) p_int.y -= 1;
         auto p_rem = (p - p_int).template as<point> ();
 
-    T closest = std::numeric_limits<T>::infinity ();
+        static const util::vector2i OFFSETS[] = {
-    T value;
+                                    {  0, -2 },
+                        { -1, -1 }, {  0, -1 }, {  1, -1 },
+            { -2,  0 }, { -1,  0 }, {  0,  0 }, {  1,  0 }, {  2,  0 },
+                        { -1,  1 }, {  0,  1 }, {  1,  1 },
+                                    {  0,  2 },
+        };
 
-    for (signed y = -1; y <= 1; ++y)
+        static const size_t COUNT = elems (OFFSETS);
-        for (signed x = -1; x <= 1; ++x) {
-            util::vector2i offset {x, y};
-            auto c = centroid (p_int + offset);
-            auto d = util::distance (p_rem, c + offset);
 
-            if (d < closest) {
+        // find the distances to each neighbour's centroid
-                closest = d;
+        util::point2<T> centres[COUNT];
-                value = generate (p_int + offset);
+        for (size_t i = 0; i < COUNT; ++i)
-            }
+            centres[i] = centroid (p_int + OFFSETS[i]) + OFFSETS[i];
+
+        T distances[COUNT];
+        for (size_t i = 0; i < COUNT; ++i)
+            distances[i] = util::distance (p_rem, centres[i]);
+
+        // sort the distances, using indices so we can use 'offsets' to generate values
+        unsigned indices[COUNT];
+        std::iota (std::begin (indices), std::end (indices), 0);
+
+        std::sort (std::begin (indices),
+                   std::end (indices),
+                   [&] (auto a, auto b) {
+            return distances[a] < distances[b];
+        });
+
+        // calculate normalisation constants for the 9 nearest points. the
+        // neighbourhood size is implicitly specified by the 1.5 unit maximum
+        // distance.
+        constexpr auto MAX_DISTANCE = std::hypot (1.5f, 1.5f);
+        const auto lo = distances[indices[0]];
+        const auto hi = std::min (distances[indices[COUNT-1]], MAX_DISTANCE);
+
+        T out  = 0.f;
+        T sumw = 0.f;
+
+        // sum the weight values of each neighbour. weight by a function of
+        // the distance. we use an power function which allows a known width
+        // to blend.
+        for (size_t i = 0; i < COUNT && distances[indices[i]] <= MAX_DISTANCE; ++i)
+        {
+            auto v = generate (p_int + OFFSETS[indices[i]]);
+            auto d = (distances[indices[i]] - lo) / (hi - lo);
+            auto w = std::pow (1 - d, m_power);
+
+            sumw += w;
+            out += v * w;
         }
 
-    return value;
+        return out / sumw;
-}
+    }
 
 
-///////////////////////////////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////////
-template <typename T>
+    template <typename T>
-util::noise::seed_t
+    seed_t
-util::noise::basis::patch<T>::seed (void) const
+    patch<T>::seed (void) const
-{
+    {
         return m_seed;
-}
+    }
 
 
-//-----------------------------------------------------------------------------
+    //-------------------------------------------------------------------------
-template <typename T>
+    template <typename T>
-util::noise::seed_t
+    seed_t
-util::noise::basis::patch<T>::seed (util::noise::seed_t _seed)
+    patch<T>::seed (util::noise::seed_t _seed)
-{
+    {
         return m_seed = _seed;
-}
+    }
 
 
-///////////////////////////////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////////
-template <typename T>
+    template <typename T>
-util::point2<T>
+    T
-util::noise::basis::patch<T>::centroid (util::point2i p) const
+    patch<T>::width (void) const
-{
+    {
+        return m_width;
+    }
+
+
+    //-------------------------------------------------------------------------
+    template <typename T>
+    T
+    patch<T>::width (T _width)
+    {
+        m_width = _width;
+        m_power = exactly_zero (_width)
+                  ? std::numeric_limits<T>::infinity ()
+                  : std::log (THRESHOLD) / std::log (1 - _width);
+
+        return m_width;
+    }
+
+
+    ///////////////////////////////////////////////////////////////////////////
+    template <typename T>
+    util::point2<T>
+    patch<T>::centroid (util::point2i p) const
+    {
         using util::hash::murmur2::mix;
 
         auto u = mix (m_seed, mix (uint64_t (p.x), uint64_t (p.y)));
@@ -102,16 +171,17 @@ util::noise::basis::patch<T>::centroid (util::point2i p) const
 
         CHECK_LIMIT (r, T{0}, T{1});
         return r;
-}
+    }
 
 
-//-----------------------------------------------------------------------------
+    //-------------------------------------------------------------------------
-template <typename T>
+    template <typename T>
-T
+    T
-util::noise::basis::patch<T>::generate (util::point2i p) const
+    patch<T>::generate (util::point2i p) const
-{
+    {
         using util::hash::murmur2::mix;
 
         auto u = mix (m_seed, mix (uint64_t (p.x), uint64_t (p.y)));
         return (u & 0xffff) / T{0xffff};
-}
+    }
+} } }

tools/noise.cpp

@@ -209,6 +209,7 @@ main (int argc, char **argv)
     float scale = 1.f;
     float turbulence = 0.f;
     unsigned single = 0;
+    float width = 0;
 
     // fill variables from arguments
     util::cmdopt::parser args;
@@ -222,7 +223,8 @@ main (int argc, char **argv)
     args.add<util::cmdopt::option::count<unsigned>>  ('1', "single",        "single octave",            single);
     args.add<util::cmdopt::option::value<float>>     ('H', "hurst",         "Hurst exponent",           H);
     args.add<util::cmdopt::option::value<float>>     ('x', "scale",         "frequency multiplier",     scale);
-    args.add<util::cmdopt::option::value<float>>     ('t', "turbulence","turbulence scale",         turbulence);
+    args.add<util::cmdopt::option::value<float>>     ('t', "turbulence",    "turbulence scale",         turbulence);
+    args.add<util::cmdopt::option::value<float>>     ('W', "patch-width",   "patch blur width",         width);
 
     args.scan (argc, argv);
 
@@ -293,7 +295,7 @@ main (int argc, char **argv)
         }
 
         case PATCH: {
-            b.reset<util::noise::basis::patch<float>> (seed);
+            b.reset<util::noise::basis::patch<float>> (seed, width);
             break;
         }
 
@@ -343,6 +345,7 @@ main (int argc, char **argv)
     auto offset = *range.first;
     auto div    = *range.second - *range.first;
 
+    std::cerr << '[' << *range.first << ',' << *range.second << "]\n";
     std::transform (img.begin (), img.end (), img.begin (), [offset,div] (auto i) { return (i - offset) / div; });
 
     // write the images to disk