Color Maps (#752)

* Color Maps

* Make color maps constexpr

[ci skip]

* Add newton fractal example

[ci skip]

* Remove some unused code.

* Make the color map base class generic in size

Fix naming convention
[ci skip]

* Begin documentation.

* Move helper functions from example into tools header

[ci skip]

* Update docs and remove non-ASCII characters from example

* Add image to docs

* Reduce size of virdis_newton_fractal from 1.31MB to 131KB

[ci skip]

* Add performance file

* Don't force linear complexity and fix CI failure for old clang versions

* Convert color_maps to free functions.

* Add missing header and remove constexpr test

* Convert tabs to spaces

[ci skip]

* Fix compile tests and make static constexpr uniform across data

* Add swatches to docs.

* Fix image links in docs

[ci skip]

Co-authored-by: Nick Thompson <[email protected]>

doc/internals/color_maps.qbk

@@ -0,0 +1,98 @@
+[/
+  Copyright Nick Thompson, Matt Borland, 2022
+  Distributed under the Boost Software License, Version 1.0.
+  (See accompanying file LICENSE_1_0.txt or copy at
+  http://www.boost.org/LICENSE_1_0.txt).
+]
+
+[section:color_maps Color Maps]
+
+[h4 Synopsis]
+
+``
+#include <boost/math/tools/color_maps.hpp>
+
+namespace boost::math::tools {
+
+template<typename Real>
+std::array<Real, 3> viridis(Real x);
+
+template<typename Real>
+std::array<Real, 3> plasma(Real x);
+
+template<typename Real>
+std::array<Real, 3> black_body(Real x);
+
+template<typename Real>
+std::array<Real, 3> inferno(Real x);
+
+template<typename Real>
+std::array<Real, 3> smooth_cool_warm(Real x);
+
+template<typename Real>
+std::array<Real, 3> kindlmann(Real x);
+
+template<typename Real>
+std::array<Real, 3> extended_kindlmann(Real x);
+
+template<typename Real>
+std::array<uint8_t, 4> to_8bit_rgba(std::array<Real, 3> const & color);
+
+} // namespaces
+``
+
+[heading Description]
+
+Abstractly, a color map is any function which maps [0, 1] -> [0, 1]^3.
+As stated, this definition is too broad to be useful, so in Boost, we restrict our attention to the subset of color maps which are useful for the understanding of scientific data.
+[@https://www.kennethmoreland.com/color-advice/BadColorMaps.pdf Much research] has demonstrated that color maps differ wildly in their usefulness for interpreting quantitative data; see [@https://www.kennethmoreland.com/color-advice/ here] for details.
+In addition, different color maps are useful in different contexts.
+For example, the `smooth_cool_warm` color map is useful for examining surfaces embedded in 3-space which have scalar fields defined on them, whereas the `inferno` color map is better for understanding 2D data.
+
+Despite the fact that a color map, per our definition, has a domain of [0, 1], we nonetheless do not throw an exception if the value provided falls outside this range.
+This is for two reasons: First, visualizations are themselves amazing debuggers, and if we threw an exception the calculation would not complete and visual debugging would be inaccessible.
+Second, often small changes in floating point rounding cause the value provided to be only slightly below zero, or just slightly above 1.
+Hence, we make a call to `std::clamp` before interpolating into the color table. 
+
+For an example of how to use these facilites please refer to `example/color_maps_example.cpp`
+Note: To compile the examples directly you will need to have [@https://github.com/lvandeve/lodepng lodepng].
+An example of the viridis color map using [@https://en.wikipedia.org/wiki/Newton_fractal the newton fractal] is shown below:
+
+[$../images/viridis_newton_fractal.png]
+
+Swatches of each are listed below:
+
+[$../images/smooth_cool_warm_swatch.png]
+
+*Smooth cool warm*
+
+[$../images/viridis_swatch.png]
+
+*Viridis*
+
+[$../images/plasma_swatch.png]
+
+*Plasma*
+
+[$../images/black_body_swatch.png]
+
+*Black body*
+
+[$../images/inferno_swatch.png]
+
+*Inferno*
+
+[$../images/kindlmann_swatch.png]
+
+*Kindlmann*
+
+[$../images/extended_kindlmann_swatch.png]
+
+*Extended Kindlmann*
+
+[heading References]
+
+* Ken Moreland. ['Why We Use Bad Color Maps and What You Can Do About it] .
+
+[endsect]
+[/section:color_maps Color Maps]

doc/math.qbk

@@ -787,6 +787,7 @@ and as a CD ISBN 0-9504833-2-X  978-0-9504833-2-0, Classification 519.2-dc22.
 [include internals/minimax.qbk]
 [include internals/relative_error.qbk] [/ uncertainty of floating-point calculations.]
 [include internals/test_data.qbk]
+[include internals/color_maps.qbk]
 [endsect] [/section:internals]
 
 [endmathpart] [/mathpart internals Internal Details: Series, Rationals and Continued Fractions, Testing, and Development Tools]

example/color_maps_example.cpp

@@ -0,0 +1,221 @@
+//  (C) Copyright Nick Thompson 2021.
+//  (C) Copyright Matt Borland 2022.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <cmath>
+#include <cstdint>
+#include <array>
+#include <complex>
+#include <tuple>
+#include <iostream>
+#include <vector>
+#include <limits>
+#include <boost/math/tools/color_maps.hpp>
+
+#if !__has_include("lodepng.h")
+ #error "lodepng.h is required to run this example."
+#endif
+#include "lodepng.h"
+#include <iostream>
+#include <string>
+#include <vector>
+
+
+// In lodepng, the vector is expected to be row major, with the top row
+// specified first. Note that this is a bit confusing sometimes as it's more
+// natural to let y increase moving *up*.
+unsigned write_png(const std::string &filename,
+                   const std::vector<std::uint8_t> &img, std::size_t width,
+                   std::size_t height) {
+  unsigned error = lodepng::encode(filename, img, width, height,
+                                   LodePNGColorType::LCT_RGBA, 8);
+  if (error) {
+    std::cerr << "Error encoding png: " << lodepng_error_text(error) << "\n";
+  }
+  return error;
+}
+
+
+// Computes ab - cd.
+// See: https://pharr.org/matt/blog/2019/11/03/difference-of-floats.html
+template <typename Real>
+inline Real difference_of_products(Real a, Real b, Real c, Real d)
+{
+    Real cd = c * d;
+    Real err = std::fma(-c, d, cd);
+    Real dop = std::fma(a, b, -cd);
+    return dop + err;
+}
+
+template<typename Real>
+auto fifth_roots(std::complex<Real> z)
+{
+    std::complex<Real> v = std::pow(z,4);
+    std::complex<Real> dw = Real(5)*v;
+    std::complex<Real> w = v*z - Real(1);
+    return std::make_pair(w, dw);
+}
+
+template<typename Real>
+auto g(std::complex<Real> z)
+{
+    std::complex<Real> z2 = z*z;
+    std::complex<Real> z3 = z*z2;
+    std::complex<Real> z4 = z2*z2;
+    std::complex<Real> w = z4*(z4 + Real(15)) - Real(16);
+    std::complex<Real> dw = Real(4)*z3*(Real(2)*z4 + Real(15));
+    return std::make_pair(w, dw);
+}
+
+template<typename Real>
+std::complex<Real> complex_newton(std::function<std::pair<std::complex<Real>,std::complex<Real>>(std::complex<Real>)> f, std::complex<Real> z)
+{
+    // f(x(1+e)) = f(x) + exf'(x)
+    bool close = false;
+    do
+    {
+        auto [y, dy] = f(z);
+        z -= y/dy;
+        close = (abs(y) <= 1.4*std::numeric_limits<Real>::epsilon()*abs(z*dy));
+    } while(!close);
+    return z;
+}
+
+template<typename Real>
+class plane_pixel_map
+{
+public:
+    plane_pixel_map(int64_t image_width, int64_t image_height, Real xmin, Real ymin)
+    {
+        image_width_ = image_width;
+        image_height_ = image_height;
+        xmin_ = xmin;
+        ymin_ = ymin;
+    }
+
+    std::complex<Real> to_complex(int64_t i, int64_t j) const {
+        Real x = xmin_ + 2*abs(xmin_)*Real(i)/Real(image_width_ - 1);
+        Real y = ymin_ + 2*abs(ymin_)*Real(j)/Real(image_height_ - 1);
+        return std::complex<Real>(x,y);
+    }
+
+    std::pair<int64_t, int64_t> to_pixel(std::complex<Real> z) const {
+        Real x = z.real();
+        Real y = z.imag();
+        Real ii = (image_width_ - 1)*(x - xmin_)/(2*abs(xmin_));
+        Real jj = (image_height_ - 1)*(y - ymin_)/(2*abs(ymin_));
+
+        return std::make_pair(std::round(ii), std::round(jj));
+    }
+
+private:
+    int64_t image_width_;
+    int64_t image_height_;
+    Real xmin_;
+    Real ymin_;
+};
+
+int main(int argc, char** argv)
+{
+    using Real = double;
+    using boost::math::tools::viridis;
+    using std::sqrt;
+
+    std::function<std::array<Real, 3>(Real)> color_map = viridis<Real>;
+    std::string requested_color_map = "viridis";
+    if (argc == 2) {
+       requested_color_map = std::string(argv[1]);
+       if (requested_color_map == "smooth_cool_warm") {
+          color_map = boost::math::tools::smooth_cool_warm<Real>;
+       }
+       else if (requested_color_map == "plasma") {
+          color_map = boost::math::tools::plasma<Real>;
+       }
+       else if (requested_color_map == "black_body") {
+          color_map = boost::math::tools::black_body<Real>;
+       }
+       else if (requested_color_map == "inferno") {
+          color_map = boost::math::tools::inferno<Real>;
+       }
+       else if (requested_color_map == "kindlmann") {
+          color_map = boost::math::tools::kindlmann<Real>;
+       }
+       else if (requested_color_map == "extended_kindlmann") {
+          color_map = boost::math::tools::extended_kindlmann<Real>;
+       }
+       else {
+          std::cerr << "Could not recognize color map " << argv[1] << ".";
+          return 1;
+       }
+    }
+    constexpr int64_t image_width = 1024;
+    constexpr int64_t image_height = 1024;
+    constexpr const Real two_pi = 6.28318530718;
+
+    std::vector<std::uint8_t> img(4*image_width*image_height, 0);
+    plane_pixel_map<Real> map(image_width, image_height, Real(-2), Real(-2));
+
+    for (int64_t j = 0; j < image_height; ++j)
+    {
+        std::cout << "j = " << j << "\n";
+        for (int64_t i = 0; i < image_width; ++i)
+        {
+            std::complex<Real> z0 = map.to_complex(i,j);
+            auto rt = complex_newton<Real>(g<Real>, z0);
+            // The root is one of exp(2*pi*ij/5). Therefore, it can be classified by angle.
+            Real theta = std::atan2(rt.imag(), rt.real());
+            // Now theta in [-pi,pi]. Get it into [0,2pi]:
+            if (theta < 0) {
+                theta += two_pi;
+            }
+            theta /= two_pi;
+            if (std::isnan(theta)) {
+                std::cerr << "Theta is a nan!\n";
+            }
+            auto c = boost::math::tools::to_8bit_rgba(color_map(theta));
+            int64_t idx = 4 * image_width * (image_height - 1 - j) + 4 * i;
+            img[idx + 0] = c[0];
+            img[idx + 1] = c[1];
+            img[idx + 2] = c[2];
+            img[idx + 3] = c[3];
+        }
+    }
+
+    std::array<std::complex<Real>, 8> roots;
+    roots[0] = -Real(1);
+    roots[1] = Real(1);
+    roots[2] = {Real(0), Real(1)};
+    roots[3] = {Real(0), -Real(1)};
+    roots[4] = {sqrt(Real(2)), sqrt(Real(2))};
+    roots[5] = {sqrt(Real(2)), -sqrt(Real(2))};
+    roots[6] = {-sqrt(Real(2)), -sqrt(Real(2))};
+    roots[7] = {-sqrt(Real(2)), sqrt(Real(2))};
+
+    for (int64_t k = 0; k < 8; ++k)
+    {
+        auto [ic, jc] = map.to_pixel(roots[k]);
+
+        int64_t r = 7;
+        for (int64_t i = ic - r; i < ic + r; ++i)
+        {
+            for (int64_t j = jc - r; j < jc + r; ++j)
+            {
+                if ((i-ic)*(i-ic) + (j-jc)*(j-jc) > r*r)
+                {
+                    continue;
+                }
+                int64_t idx = 4 * image_width * (image_height - 1 - j) + 4 * i;
+                img[idx + 0] = 0;
+                img[idx + 1] = 0;
+                img[idx + 2] = 0;
+                img[idx + 3] = 0xff;
+            }
+        }
+    }
+
+    // Requires lodepng.h
+    // See: https://github.com/lvandeve/lodepng for download and compilation instructions
+    write_png(requested_color_map + "_newton_fractal.png", img, image_width, image_height);
+}