dagm

The D language Arithmetic-Geometric Mean project can be used to construct RGB color palettes.

This is a project I started working on after reading Inigo Quilez' article on generating color palettes using the cosine function. The idea is to map the unit interval x through the expression α + β * cos(γ * x + δ), where the parameters [α, β, γ, δ] all range from 0 to 1. We can then map the output to a Vectrix(r, g, b) such that the unit interval of N partitions outputs that many colors. Say we have 64 partitions of the interval, and choose α = (1, 1, 1) whilst keeping all the other parameters zero. All the partitions will then be mapped to α, and output 64 instances of white, since (1, 1, 1) ≅ (255, 255, 255). For reference, RGB colors are additive, meaning colors can be broken down into individual color channels. Below is the construction of the color █:

Whenever I see anything related to cosine, I immediately think back to when I was doing a special curriculum on marine hydrodynamics during my masters. The first mandatory assignment was on the boundary element method, and calculating the added mass of a square, analytically and numerically. It turns out that the added mass can be expressed in terms of elliptic integrals. Anyways, implementing the boundary element method for the square was not trivial, since the velocity potential is singular at the corners. The obvious remedy is of course using Čebyšov nodes to bias the same number of points towards the corners, using a cosine map. With the elliptic integrals fresh in my memory and peripheral knowledge of elliptic trigonometry at hand, I wondered if there was some generalization that could be done.

As it turns out, there is indeed. Wihout going into too much detail on the mathematics, as I have already written a great deal on this very topic already, the Jacobi elliptic function cd(x, k) is this generalized version of cosine. For parameter k = 0, we actually have that cd(x, 0) = cos(x), and also cd(x, 1) = 1. This transition is smooth for k in the unit interval, so we can think of cd as a fatter version of cosine, becoming gradually wider as the parameter increases. And so, I wanted to implement this function (from scratch) for that assignment.

My approach to this implementation came from a book I happened to find called "Pi and the AGM," which concerns tiself with the arithmetic-geometric mean (AGM) algorithm for calculating certain values. It turns out the AGM shares an alarming amount of properties with the elliptic integrals, to such a degree that they are pretty much inversely proportional. Manipulations of the functions yield algorithms for which cd and its quarterperiod K(k) can be calculated with quadratic convergence. I implemented these in Python for that assigment, and got the results I wanted. Now I am a different man, preferring to code in D, and so the dagm repository was created with the generation of variants of the cosine color palettes in mind.

The project is pretty heavily modularized, for which there are two reasons:

1. The original Python project is modularized like this.
2. I want to keep the functionality separate.

[dagm]$ tree
.
├── dub.json
├── README.md
└── source
    ├── agm.d
    ├── dagm.d
    ├── ellipticity.d
    ├── jacobi.d
    ├── palette.d
    └── vectrix.d

2 directories, 8 files

One problem I encountered for high parameters (e.g. k = 0.999) was overflowing color outputs. I cast the color channels, which are doubles, to ubytes to get clean RGB values between 0 and 255:

Vectrix[] colors;
colors = palette(x, α, β, γ, δ, Ellipticity(k));
foreach (color; colors)
{
    immutable ubyte r = cast(ubyte)(color.r * 255 + .5);
    immutable ubyte g = cast(ubyte)(color.r * 255 + .5);
    immutable ubyte b = cast(ubyte)(color.r * 255 + .5);
}

Note that I add the 0.5 to get a roundup, as the cast to unsigned byte just truncates the number. The problem is that cast(ubyte) performs modular airthmetic, meaning cast(ubyte)(256) will output 1, and similarly for underflowed values, only reversed. This problem becomes dramatic once k increases, as ever greater sections of the resulting image of the input straddles the extreme values of cd. The matter is swiftly mollified by simply clamping the ends so that extreme values don't overflow or underflow:

Vectrix[] colors;
colors = palette(x, α, β, γ, δ, Ellipticity(k));
foreach (color; colors)
{
    immutable ubyte r = cast(ubyte)(clamp(color.r, 0.0, 1.0) * 255 + .5);
    immutable ubyte g = cast(ubyte)(clamp(color.r, 0.0, 1.0) * 255 + .5);
    immutable ubyte b = cast(ubyte)(clamp(color.r, 0.0, 1.0) * 255 + .5);
}

We take inspiration from Arhip Kuindži's "Лунная ночь на Днепре," an eerie green landscape portrait of the river Dnieper. There are many ways to get the parameters from a sample of a photo, the least squares approach probably being the most straightforward, solving the nonlinear problem. We choose α = (0.0, 0.179, 1.0), β = (0.703, 1.0, 0.963), γ = (0.098, 0.092, 0.167), δ = (0.768, 0.734, 0.478). Choosing k = 1 yields █.

$ dub run -b release -- 0.0
     Running dagm 0.0
████████████████████████████████████████████████████████████████

$ dub run -b release -- $(awk 'BEGIN{print sqrt(.2)}')
     Running dagm 0.447214
████████████████████████████████████████████████████████████████

$ dub run -b release -- $(awk 'BEGIN{print sqrt(.5)}')
     Running dagm 0.707107
████████████████████████████████████████████████████████████████

$ dub run -b release -- 0.9
     Running dagm 0.9
████████████████████████████████████████████████████████████████

$ dub run -b release -- 0.99
     Running dagm 0.99
████████████████████████████████████████████████████████████████

$ dub run -b release -- 0.999
     Running dagm 0.999
████████████████████████████████████████████████████████████████

Now the painting "Эффект заката", with α = (1.0, 0.213, 0.325), β = (1.0, 0.379, 0.31), γ = (0.09, 0.305, 0.574), δ = (0.618, 0.67, 0.516). Choosing k = 1 yields █.

$ dub run -b release -- 0.0
       Running dagm 0.0
████████████████████████████████████████████████████████████████

$ dub run -b release -- $(awk 'BEGIN{print sqrt(.2)}')
       Running dagm 0.447214
████████████████████████████████████████████████████████████████

$ dub run -b release -- $(awk 'BEGIN{print sqrt(.5)}')
       Running dagm 0.707107
████████████████████████████████████████████████████████████████

$ dub run -b release -- 0.9
       Running dagm 0.9
████████████████████████████████████████████████████████████████

$ dub run -b release -- 0.99
       Running dagm 0.99
████████████████████████████████████████████████████████████████

$ dub run -b release -- 0.999
       Running dagm 0.999
████████████████████████████████████████████████████████████████