Color Maps

[mborland]

@NAThompson Here's my cut at implementing color tables without as much code duplication. I also think this is a good opportunity to start using C++20 concepts. Compiler support seems pretty good.

[NAThompson]

@mborland : I like the concepts.

Two comments:

What about:

boost::math::tools::kindlemann

rather than

boost::math::tools::kindlemann_color_map

?

I wonder if the latter is too verbose.

Next, will we get any performance wins if we use:

std::array<std::array<Real, 3>, N>

rather than:

std::array<Real, N> r;
std::array<Real, N> g;
std::array<Real, N> b;

?

I also wonder if we should be using the move constructor, since these color tables are generally known at compile time. I think the tables can be static constexpr const, which makes moving them a bit awkward . . .

Is there a concept for a RandomAccessContainer which has elements which live in R^3?

[mborland]

1. I think using the full name describes the intent better, but agree it is verbose. Since it's an object type the user will only have to type that full thing once.
2. I don't think in this case we can gain anything because we can't convert this into a compressed matrix (that I know of at least)
3. Agreed this should be constexpr so now the entire class is. The test file validates use in constexpr contexts.
4. Not out of the box. You can write your own concepts which as we move forward will need to be done (e.g. to support boost.multiprecision types).

[NAThompson]

I think using the full name describes the intent better, but agree it is verbose. Since it's an object type the user will only have to type that full thing once.

So you're thinking the users will do:

auto kindlemann = boost::math::tools::kindlemann_color_map()

and the less verbose syntax will occur in user code, rather than in ours?

I don't think in this case we can gain anything because we can't convert this into a compressed matrix (that I know of at least)

Well, you might get gains from having the referenced elements being contiguous in memory, but of course it won't change any fundamental complexity . . .

... If these are all constexpr do they need to be classes? I feel like we may be able to get them to be straight function calls.

[mborland]

So you're thinking the users will do:

auto kindlemann = boost::math::tools::kindlemann_color_map()

and the less verbose syntax will occur in user code, rather than in ours?

Yes, then all the user has to do is call kindlemann(scalar)

I don't think in this case we can gain anything because we can't convert this into a compressed matrix (that I know of at least)

Well, you might get gains from having the referenced elements being contiguous in memory, but of course it won't change any fundamental complexity . . .

... If these are all constexpr do they need to be classes? I feel like we may be able to get them to be straight function calls.

To me derived classes was the easiest way to eliminate code repetition. It also allowed for an easy implementation of custom_color_map where the user can define whatever values they want.

You mentioned using boost.GIL for an example. Do you have one written already I could edit to make this work? I have never used GIL before.

[NAThompson]

You mentioned using boost.GIL for an example. Do you have one written already I could edit to make this work? I have never used GIL before.

Actually I've never used it either. We can even write a .pnm to a textfile and then screenshot it for the docs; another option is to use lodepng; that's what I did here.

Let me think of a 2D function which would nicely get plotted with these color maps. This one might be good to show off viridis; the code to reproduce it is here.

@jzmaddock , @ckormanyos : Can you think of anything that could be displayed really nicely using a color map?

[mborland]

@NAThompson I ended up modifying your example from DRT, and get the following output:

Also here is where the color table object is constructed, and here is where it is used. I don't think the calls are overly verbose.

[NAThompson]

@mborland : That looks awesome! Looks like it wasn't too much work to get there either . . .

I agree that the calls aren't too verbose, however I do wonder who will actually care to provide their own color tables. These color tables have been through extensive vetting by psychologists who measure their interpretive power. I do worry that if people do provide their own color tables, they'll use something like rainbow which is known to be terrible for scientific interpretation.

[NAThompson]

@mborland : It occurs to me that the user might want to try multiple color tables without recompilation, say, via a command line flag. Here's my first attempt:

diff --git a/example/color_maps_example.cpp b/example/color_maps_example.cpp
index 23b3a3782..b46d84a0d 100644
--- a/example/color_maps_example.cpp
+++ b/example/color_maps_example.cpp
@@ -130,9 +130,21 @@ private:
     Real ymin_;
 };
 
-int main()
+int main(int argc, char** argv)
 {
-    using Real = long double;
+    using Real = double;
+    if (argc <= 1) {
+        std::cerr << "Usage: ./color_map {color_table}, where {color_table} is one of viridis, kindlemann, plasma, or ..\n";
+	return 1;
+    }
+    std::function<std::array<Real, 3>(Real)> color_map;
+    if (argv[1] == "viridis") {
+        color_map = boost::math::tools::viridis_color_map<Real>;
+    }
+    else {
+        std::cerr << "Could not recognize color map " << argv[1] << " options are viridis, kindlemann, plasma . . .\n";
+	return 1;
+    }
     constexpr int64_t image_width = 4096;
     constexpr int64_t image_height = 4096;
     std::vector<uint8_t> img(4*image_width*image_height, 0);

but this doesn't really work . . .

[NAThompson]

Now that I think of it, I doubt we're even coming close to achieving double precision resolution. Would it make sense to replace these character by a 32 bit hexfloat representation?

[mborland]

We could. I don't think there are many displays in circulation that have a 32 bit color depth.

[NAThompson]

@mborland : I personally like huge images but the maintainers of the global boost package will complain if we add a 1.3MB file to the distribution. Could we make the image smaller?

Also could you drop a benchmark into reporting/performance?

[ckormanyos]

Can the color maps be used on histogram data from fractal iterations such as Mandelbrot set?

[ckormanyos]

You mentioned using boost.GIL for an example. Do you have one written already I could edit to make this work? I have never used GIL before.

Actually I've never used it either.

I have a non-published project that uses Boost.Gil in combination with a jpeg library that I made mostly to be portable as well. I can't find the origin of the original source code of my ported jpeg library any more. It seems to come from last century. I don't know if it can be redistributed.

Aside from the jpeg LIB potential availability/portability problem, I ended up going from boost::gil::rgb8_view_t directly to jpeg via boost::gil::jpeg_write_view. But this does mandate a jpeg eingine (LIB) in the project.

[NAThompson]

Can the color maps be used on histogram data from fractal iterations such as Mandelbrot set?

Of course!

[mborland]

@NAThompson Here is the performance data

Unable to determine clock rate from sysctl: hw.cpufrequency: No such file or directory
2022-02-01T22:19:08+02:00
Running ./color_tables_performance
Run on (10 X 24.0833 MHz CPU s)
CPU Caches:
  L1 Data 64 KiB (x10)
  L1 Instruction 128 KiB (x10)
  L2 Unified 4096 KiB (x5)
Load Average: 1.51, 1.51, 1.58
-----------------------------------------------------------------------------------------------
Benchmark                                                     Time             CPU   Iterations
-----------------------------------------------------------------------------------------------
color_table_benchmark<float>/64/real_time                  4144 ns         4144 ns       125506
color_table_benchmark<float>/128/real_time                 8230 ns         8230 ns        84417
color_table_benchmark<float>/256/real_time                16417 ns        16417 ns        42445
color_table_benchmark<float>/512/real_time                32776 ns        32776 ns        21361
color_table_benchmark<float>/1024/real_time               65499 ns        65499 ns        10596
color_table_benchmark<float>/2048/real_time              130932 ns       130932 ns         5319
color_table_benchmark<float>/4096/real_time              261722 ns       261716 ns         2671
color_table_benchmark<float>/8192/real_time              523603 ns       523604 ns         1322
color_table_benchmark<float>/16384/real_time            1046917 ns      1046917 ns          663
color_table_benchmark<float>/32768/real_time            2093535 ns      2093536 ns          334
color_table_benchmark<float>/65536/real_time            4185543 ns      4185545 ns          167
color_table_benchmark<float>/131072/real_time           8373275 ns      8373277 ns           83
color_table_benchmark<float>/262144/real_time          16747412 ns     16747405 ns           42
color_table_benchmark<float>/524288/real_time          33508554 ns     33508524 ns           21
color_table_benchmark<float>/1048576/real_time         66997783 ns     66997800 ns           10
color_table_benchmark<float>/real_time_BigO               63.90 N         63.90 N
color_table_benchmark<float>/real_time_RMS                    0 %             0 %
color_table_benchmark<double>/64/real_time                 4112 ns         4112 ns       170016
color_table_benchmark<double>/128/real_time                8185 ns         8185 ns        84702
color_table_benchmark<double>/256/real_time               16317 ns        16317 ns        42690
color_table_benchmark<double>/512/real_time               32595 ns        32595 ns        21467
color_table_benchmark<double>/1024/real_time              65187 ns        65187 ns        10642
color_table_benchmark<double>/2048/real_time             131284 ns       131265 ns         5228
color_table_benchmark<double>/4096/real_time             267619 ns       267287 ns         2637
color_table_benchmark<double>/8192/real_time             536032 ns       535327 ns         1286
color_table_benchmark<double>/16384/real_time           1059829 ns      1059212 ns          646
color_table_benchmark<double>/32768/real_time           2108699 ns      2107826 ns          327
color_table_benchmark<double>/65536/real_time           4198880 ns      4197747 ns          166
color_table_benchmark<double>/131072/real_time          8501609 ns      8492410 ns           83
color_table_benchmark<double>/262144/real_time         16936271 ns     16925381 ns           42
color_table_benchmark<double>/524288/real_time         33957917 ns     33924333 ns           21
color_table_benchmark<double>/1048576/real_time        67226867 ns     67208000 ns           10
color_table_benchmark<double>/real_time_BigO              64.27 N         64.24 N
color_table_benchmark<double>/real_time_RMS                   1 %             1 %
color_table_benchmark<long double>/64/real_time            4113 ns         4113 ns       170037
color_table_benchmark<long double>/128/real_time           8183 ns         8183 ns        84687
color_table_benchmark<long double>/256/real_time          16403 ns        16402 ns        42749
color_table_benchmark<long double>/512/real_time          32592 ns        32592 ns        21458
color_table_benchmark<long double>/1024/real_time         65143 ns        65143 ns        10656
color_table_benchmark<long double>/2048/real_time        130240 ns       130240 ns         5349
color_table_benchmark<long double>/4096/real_time        263375 ns       263214 ns         2681
color_table_benchmark<long double>/8192/real_time        522375 ns       522349 ns         1329
color_table_benchmark<long double>/16384/real_time      1041733 ns      1041734 ns          669
color_table_benchmark<long double>/32768/real_time      2082051 ns      2082048 ns          335
color_table_benchmark<long double>/65536/real_time      4164367 ns      4164363 ns          168
color_table_benchmark<long double>/131072/real_time     8332090 ns      8332096 ns           83
color_table_benchmark<long double>/262144/real_time    16670387 ns     16670262 ns           42
color_table_benchmark<long double>/524288/real_time    33341579 ns     33341571 ns           21
color_table_benchmark<long double>/1048576/real_time   66667208 ns     66667200 ns           10
color_table_benchmark<long double>/real_time_BigO         63.58 N         63.58 N
color_table_benchmark<long double>/real_time_RMS              0 %             0 %

[NAThompson]

@mborland : Shouldn't the call be constant time?

[mborland]

@NAThompson This is now clean with CI. Any other changes, or good to merge?

[NAThompson]

I still want to check if we can make the color data an array of structs, rather than three arrays, and see how it affects the performance. Gimme a few days . . .

[NAThompson]

Ok I managed to hack up a script to get the color tables into an array of structs:

using CSV
using DataFrames
csv_reader = CSV.File("plasma-table-float-0256.csv")
print("std::array<std::array<double, 3>, 256> plasma_data_ = {{")
for row in csv_reader
    println("{$(row.RGB_r), $(row.RGB_g), $(row.RGB_b)}, ")
end
print("}};")

@mborland : Accidently amend committed on a shared branch; my bad. Could you rerun the performance benchmark on your machine? I'm getting:

./a.out                                                                                           color_tables ✱ ◼
2022-02-06T10:52:23-08:00
Running ./a.out
Run on (4 X 2700 MHz CPU s)
CPU Caches:
  L1 Data 32 KiB (x2)
  L1 Instruction 32 KiB (x2)
  L2 Unified 256 KiB (x2)
  L3 Unified 3072 KiB (x1)
Load Average: 2.25, 1.68, 1.69
------------------------------------------------------------------------
Benchmark                              Time             CPU   Iterations
------------------------------------------------------------------------
color_table_benchmark<float>        2.30 ns         2.23 ns    303114574
color_table_benchmark<double>       6.56 ns         6.39 ns    108312187

Using the array of structs, I can get:

color_table_benchmark<float>        2.02 ns         2.01 ns   1000000000
color_table_benchmark<double>       6.07 ns         6.06 ns    693552276
array_of_structs_bm<float>          2.02 ns         2.01 ns   1000000000
array_of_structs_bm<double>         5.60 ns         5.59 ns    747409426

This is consistent with the amount of memory pressure being larger for doubles. It also makes it so we can dynamically bind a std::function to the color table allowing runtime selection of the map.

I had to remove the lodepng.h, because if a user just happens to have it on their syspath but hasn't linked it (the lodepng link is super weird) it will break their compile.

Nonetheless, here's an example that pumps out swatches of the color tables:

//  (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 <array>
#include <boost/math/tools/color_maps.hpp>
#include <cmath>
#include <complex>
#include <cstdint>
#include <iostream>
#include <limits>
#include <tuple>
#include <vector>

#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;
}

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 = 10;
  std::vector<std::uint8_t> img(4 * image_width * image_height, 0);

  for (int64_t j = 0; j < image_height; ++j) {
    std::cout << "j = " << j << "\n";
    for (int64_t i = 0; i < image_width; ++i) {
      double t = i / (image_width - 1.0);
      auto c = boost::math::tools::to_8bit_rgba(color_map(t));
      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];
    }
  }
  // Requires lodepng.h
  // See: https://github.com/lvandeve/lodepng for download and compilation
  // instructions
  write_png(requested_color_map + "_swatch.png", img, image_width,
            image_height);
}

[NAThompson]

@mborland : Nice diffs. I think I'll put swatches for all the color maps in the docs and then we'll be done.

[NAThompson]

@mborland : Could you put eyes on the docs for me then squash merge once you're happy?

[mborland]

Docs look good to me now so merging.