Logarithmic spiral, nautilus, and rainbow

The other day I stumbled into an interesting article on The Guardian online: The medieval bishop who helped to unweave the rainbow. In the article I learned for the first time of Robert Grosseteste, a 13th century British scholar (with an interesting Italian last name: Grosse teste = big heads) who was also the Bishop of Lincoln.

The Bishops’ interests and investigations covered diverse topics, making him a pre-renaissance polymath; however, it is his 1225 treatise on colour, the De Colore, that is receiving much attention.

In a recent commentary on Nature Physics (All the colours of the rainbow), and reference therein (A three-dimensional color space from the 13th century),  Smithson et al. (who also recently published a new critical edition/translation of the treatise with analysis and critical commentaries) analyze the 3D colorspace devised by Grosseteste, who claimed it allows the generation of all possible colours and to describe the variations of colours among different rainbows.

As we learn from Smithson et al., Grosseteste’s colorpsace had three dimensions, quantified by physical properties of the incident light and the medium: these are the scattering angle (which produces variation of hue within a rainbow), the purity of the scattering medium (which produces variation between different rainbows and is linked to the size of the water droplets in the rainbow), and the altitude of the sun (which produces variation in the light incident on a rainbow). The authors were able to model this colorspace and also to show that the locus of rainbow colours generated in that colorspace forms a spiral surface (a family of spiral curves, each form a specific rainbow) in  the perceptual CIELab colorspace.

I found this not only fascinating – a three-dimensional, perceptual colorspace from the 13th century!! – but also a source of renewed interest in creating the perfect perceptual colormaps by spiralling through CIELab.

My first attempt of colormap spiralling in CIELab, CubicYF, came to life by selecting hand-picked colours on CIELab colour charts at fixed lightness values (found in this document by Gernot Hoffmann). The process was described in this post, and you can see an animation of the spiral curve in CIELab space (created with the 3D color inspector plugin in ImageJ) in the video below:

Some time later, after reading this post by Rob Simmon (in particular the section on the NASA Ames Color Tool), and after an email exchange with Rob, I started tinkering with the idea of creating perceptual rainbow colormaps in CIELab programmatically, by using a helix curve or an Archimedean spiral, but reading Smithson et al. got me to try the logarithmic spiral.

So I started my experiments with a warm-up and tried to replicate a Nautilus using a logarithmic spiral with a growth ratio equal to 0.1759. You may have read that the rate at which a Nautilus shell grows can be described by the golden ratio phi, but in fact the golden spiral constructed from a golden rectangle is not a Nautilus Spiral (as an aside, as I was playing with the code I recalled reading some time ago Golden spiral, a nice blog post (with lots of code) by Cleve Moler, creator of the first version of Matlab,  who simulated a golden spiral using a continuously expanding sequence of golden rectangles and inscribed quarter circles).

My nautilus-like spiral, plotted in Figure 1, has a growth ratio of 0.1759 instead of the golden ratio of 1.618.

nautilus logarithmic spiral with growth ratio = 0.1759

Figure 1: nautilus-like spiral with growth ratio = 0.1759

And here’s the colormap (I called it logspiral) I came up with after a couple of hours of hacking: as hue cycles from 360 to 90 degrees, chroma spirals outwardly (I used a logarithmic spiral with polar equation c1*exp(c2*h) with a growth ratio c2 of 0.3 and a constant c1 of 20), and lightness increases linearly from 30 to 90.

Figure 2 shows the trajectory in the 2D CIELab a-b plane; the colours shown are the final RGB colours. In Figure 3 the trajectory is shown in 3D CIELab space. The coloured lightness profiles were made using the Colormapline submission from the Matlab File Exchange.

2D logspiral colormap in CIELab a-b plane

Figure 2: logspiral colormap trajectory in CIELab a-b plane

logspiral colormap trajectory in CIELab 3D space

Figure 3: logspiral colormap in CIELab 3D space

 

N.B. In creating logspiral, I was inspired by Figure 2 in the Nature Physics paper, but there are important differences in terms of colorspace, lightness profile and perception: I am not certain their polar coordinates are equivalent to Lightness, Chroma, and Hue, although they could; and, more importantly, the three-dimensional spirals based on Grosseteste’s colorpsace go from low lightness at low scattering angles to much higher values at mid scattering angles, and then drop again at high scattering angle (remember that these spirals describe real world rainbows), whereas lightness in logspiral lightness is strictly monotonically increasing.

In my next post I will share the Matlab code to generate a full set of logspiral colormaps sweeping the hue circle from different staring colours (and end colors) and also the slower-growing logarithmic spirals to make a set of monochromatic colormaps (similar to those in Figure 2 in the Nature Physics paper).

 

Convert Matlab colormap to Surfer colormap

In the comment section of my last post, Steve asked if I had code to generate a Surfer.clr file from my Matlab colormaps.

Some time ago I did write a simple Matlab .m file to write a colormap to a variable with the correct Surfer format, but at the time I was content to have the variable output to a .txt file, which I would then open in a text editor to add a 1-line header and change the file extension to .clr.

I went back, cleaned up the script, and automated all the formatting. This is my revised code (you may need to change the target directory c:\My Documents\MATLAB):

%% Make a Matlab colormap
% one of the colormaps from my function, Perceptually improved colormaps
sawtooth=pmkmp(256,'swtth');
%% Initialize variable for Surfer colormap
%reduce to 101 samples
sawtooth_surfer=zeros(101,5);
%% Make Surfer colormap
% add R,G,B columnns
for i=1:3 
sawtooth_surfer(:,i+1)=round(interp1([1:1:256]',sawtooth(:,i),(linspace(1,256,101))')*255);
end
% add counter and alpha (opacity) columns
sawtooth_surfer(:,1)=linspace(0,100,101);
sawtooth_surfer(:,5)=ones(101,1)*255;
%% Create output file 
filename='c:\My Documents\MATLAB\stth4surf.clr'; 
fileID=fopen(filename,'wt');
%%  Write Surfer colormap header
fprintf(fileID,'ColorMap 2 1\n');
fclose(fileID);
 
%% Add the colormap:
dlmwrite('c:\My Documents\MATLAB\stth4surf.clr',...
    sawtooth_surfer,'precision',5,'delimiter','\t', '-append');

Spectral lightness rainbow colormap

Spectral lightness rainbow

Quick post to share my replica of Ethan Montag ‘s Spectral lightness colormap from this paper. My version has a linear Lightness profile (Figure 1) that increases monotonically from 0 (black) to 100 (white) while Hue cycles from 180 to 360 degrees and then wraps around and continues from 0 to 270.

You can download the colormap files (comma separated in MS Word .doc format) with RGB values in 0-1 range and 0-255 range.

L-plot

Figure 1

Spectral-L

Figure 2

Matlab code

To run this code you will need Colorspace, a free function from Matlab File Exchange, for the color space transformations.

%% generate Chroma, Lightness, Hue components
h2=linspace(0,270,60*2.56);
h1=linspace(180,360,103);
h=horzcat(h1,h2);
c=ones(1,256)*25;
l=linspace(0,100,256);
%% merge together
lch = vertcat(l,c,h)';
%% convert to RGB
rgb = colorspace('LCH->rgb',lch);
%% Plot color-coded lightness profile
figure;
h=colormapline(1:1:256,lch(:,1),[],rgb);
set(h,'linewidth',2.8);
title('Montag Spectral L* lightness plot ','Color','k','FontSize',12,'FontWeight','demi');
%%  Pyramid test data
PY=zeros(241,241);
for i = 1:241
    temp=(0:1:i-1);
      PY(i,1:i)=temp(:);
end
test=PY.';
test=test-1;
test(test==-1)=0;
test1=test([2:1:end,1],:);
PY1=PY+test1;
PY2=fliplr(PY1);
PY3= flipud(PY1);
PY4=fliplr(PY3);
GIZA=[PY1,PY2;PY3,PY4].*2;
x=linspace(1,756,size(GIZA,1));
y=x;
[X,Y]=meshgrid(x,y);
clear i test test1 PY PY1 PY2 PY3 PY4 temp;

%% display Pyramid surface
fig1=figure;
surf(X,Y,GIZA,GIZA,'FaceColor','interp','EdgeColor','none');
view(-35,70);
colormap(rgb);
axis off;
grid off;
colorbar;
% set(fig1,'Position',[720 400 950 708]);
% set(fig1,'OuterPosition',[716 396 958 790]);
title('Montag Spectral L*','Color','k','FontSize',12,'FontWeight','demi');

Aknowledgements

The coloured lightness profiles were made using the Colormapline submission from the Matlab File Exchange.

New rainbow colormap: sawthoot-shaped lightness profile

Why another rainbow

In the comment section of my last post Steve Eddins from Mathworks reported that some Matlab users prefer Jet to Parula, the new default perceptual colormap in Matlab, because within certain ranges Jet affords a greater contrast, intended as the rate of change in lightness.

My counter-argument to that is that yes, some data may benefit from being displayed using Jet (in terms of contrast, and hence the power to resolve smaller anomalies) because of those areas of very steep rate of change of lightness, like the blue to cyan and yellow to red portions (see Figure 1). But the price one has to pay is that there is an area of very low gradient (a greenish band between cyan and yellow) where there’s nearly no contrast, which would obfuscate subtle anomalies in the data. On top of that there’s no control of where each of those areas are located, so a lot of effort has to go into trying to fit those regions of artificially high contrast to the portion of data of interest.

L_profile_jet_cl

Figure 1

Because of their high lightness, the yellow and cyan artificial edges also cause problems. In his latest blog post Steve uses a test pattern do demonstrate how they make the interpretation of trivial structures more difficult. He also explains why they occurr in some locations and not others in the first place. I wonder if the resulting regions of high lightness juxtaposed to regions of low lightness could be chromatic Mach bands.

Additionally, as Steve points out, the low-contrast juxtaposition of dark red and dark blue bands creates the visual illusion of depth (Chromostereopsis) in other positions of the test pattern, creating further confusion.

But I have some good news for the hardcore fans of Jet, and rainbow colormaps in general. I created a rainbow with a sawtooth-shaped lightness profile made up of 5 ramps, each with the same  rate of change in lightness and total lightness change of 60, and alternatively negative and positive signs. This is shown in Figure 2, and replaces the lightness profile of a basic 6-color rainbow (magenta-blue-cyan-green-yellow-red) shown in Figure 3.

Figure 2

Figure 2

Figure 2

Figure 3

With this rainbow users have the ability to apply greater contrast to their data to boost small anomalies, but in a more controlled way. The colormap is available with my File Exchange function, Perceptually improved colormaps. Below is the Matlab code I used to generate the new rainbow.

Matlab code

To run this code you will need Colorspace, a free function from Matlab File Exchange, for the color space transformations.

%% basic 6-colour rainbow
% create RGB components
m = [1, 0, 1]; % magenta
b = [0, 0, 1]; % blue
c = [0, 1, 1]; % cyan
g = [0, 1, 0]; % green
y = [1, 1, 0]; % yellow
r = [1, 0, 0]; % red
% concatenate components
rgb = vertcat(m,b,c,g,y,r);
% interpolate to 256 colours
rainbow=interp1(linspace(1, 256, 6),rgb,[1:1:256]);
%% calculate Lab components
% convert from RGB to Lab colour space
% requires this function: Colorspace transforamtions
% www.mathworks.com/matlabcentral/fileexchange/28790-colorspace-transformations
lab = colorspace('RGB->Lab',rainbow);
%% replace random lightness profile with sawtooth-shaped profile
% contrast (magnitude of lightness change) between
% each pair of adjeacent colors set to 60
L1 = [90, 30, 90, 30, 90, 30];
% interpolate to 256 lightness values
L1int = interp1(linspace(1, 256, 6),L1,[1:1:256])';
% replace
lab1 = horzcat(L1int,lab(:,2),lab(:,3));
%% new rainbow
% convert back from Lab to RGB colour space
swtth = colorspace('RGB<-Lab',lab1);

Test results

Figures 4, 5, and 6 show the three colormaps used with my Pyramid test surface (notice in Figure 5 that the green band artifact with this rainbow is even more pronounced than with jet). I welcome feedback.

Figure 4

Figure 4

Pyramid_basic_rainbow

Figure 5

Figure 4

Figure 6

Aknowledgements

The coloured lightness profiles were made using the Colormapline submission from the Matlab File Exchange.

 

What your brain does with colours when you are not “looking” – part 2

In What your brain does with colours when you are not “looking”, part 1, I displayed some audio spectrogram data (courtesy of Giuliano Bernardi at the University of Leuven) using 5 different colormaps to render the amplitude values: Jet (until recently Matlab’s standard colormap), grayscale, linear lightness rainbow, modified heated body, and cube lightness rainbow. I then asked readers to cast a vote for what they thought was the best colormap to visualize this dataset.

I was curious to see how all these colormaps fared, but my expectation was that Jet would sink to the bottom.  I was really surprised to see it came on top, one vote ahead of the linear lightness rainbow (21 and 20 votes out of 62, respectively). The modified heated body followed with 11 votes.

My surprise comes from the fact that Jet carries perceptual artifacts within the progression of colours (see for example this post). One way to demonstrate these artifacts is to convert the 2D map into a 3D surface where again we use Jet to colour amplitude values, but we use the intensities from the 2D map for the elevation. This can be done for example using the Interactive 3D Surface Plot plugin for ImageJ (as in my previous post ). The resulting surface is shown in Figure 1. This is almost exactly what your brain would do when you look at the 2D map colored with Jet in the previous post.

Surface_Plot_of_spectrogram_jet

Figure 1

In Figure 2 the same data is now displayed as a surface where amplitude values were used for the elevation, with a very light sun shading to help a bit with the perception of relief, but no colormap at all. to When comparing Figure 1 with Figure 2 one of the artifacts is immediately recognized: the highest values in Figure 2, which honours the data, become a relative low in Figure 1. This is because red has lower intensity than yellow and therefore data colored in red in 2D are plotted at a lower elevation than data colored in yellow, even though the amplitudes of the latter were lowest.

spectrogram_surf

Figure 2

For these reasons, I did not expect Jet to be the top pick. On the other hand, I think Jet is perhaps favoured because with consistent use, our brain, learns in part to accommodate for these non-perceptual artifacts in 2D maps, and because it has at least two regions of higher contrast (higher magnitude gradient) than other colormaps. Unfortunately, as I wrote in a recently published tutorial, these regions are randomly placed in the colormap, and the gradients are variable, so we gain on contrast but lose on faithfulness in representing the data structure.

Matt Hall wrote a great comment following the previous post, really making an argument for switching between multiple colormaps in the interpretation stage to explore  and highlight features in both the signal and the noise in the data, and that perhaps no single colormap is best overall. I agree 100% on almost everything Matt said, except perhaps on the best overall: looking at the 2D maps, at least with this dataset, I feel the heated body could be the best overall colormap, even if marginally. In Figure 3, Figure 4, Figure 5, and Figure 6 I show the 3D displays obtained by converting the 2D grayscale, linear lightness rainbow, modified heated body, and cube llightness rainbow, respectively. Looking at the 3D displays altogether gives me a confirmation of that feeling.

What do you think?

Surface_Plot_of_spectrogram_gray

Figure 3

Surface_Plot_of_spectrogram_lin_L_rainbow

Figure 4

Surface_Plot_of_spectrogram_mod_heated_body

Figure 5

Surface_Plot_of_spectrogram_CubicYF

Figure 6

New Matlab isoluminant colormap for azimuth data

I recently added to my Matlab File Exchange function, Perceptually improved colormaps, a colormap for periodic data like azimuth or phase. I am going to briefly showcase it using data from my degree thesis in geology, which I used before, for example in Visualization tips for geoscientists – Matlab. Figure 1, from that post, shows residual gravity anomalies in milligals.

data cube1_final_shading_slope

Figure 1

Often we’re interested in characterizing these anomalies by calculating the direction of maximum dip at each point on the surface, and for that direction display the azimuth, or dip azimuth.  I’ve done this for the surface of residual anomalies from Figure 1 and displayed the azimuth in Figure 2. Azimuth from 0 to 360 degrees are color-coded using Jet, Matlab’s standard colormap (until recently). Typically I do not trust azimuth values when the dip is close to zero because it is often contaminated by noise so I would use shading to de-saturate the colors where dip has the lowest values, but for ease of discussion I haven’t done so in this case.

Figure 2. Azimuth values color-coded with Jet.

Figure 2. Azimuth values color-coded with Jet.

There are two problems with Figure 2. First, the well-known problems with the jet colormap. For example, blue is too dark and blue areas appear as bands of constant colour. Yellow is much lighter than any other colour so we see artificial yellow edges that are not really present in the data. But there is an additional issue in Figure 2 because azimuths close in value to 0 and 360 degrees are colored with blue and red, respectively, instead of a single color as they should, causing an additional artificial edge.

In Figure 3 I recolored the map using a colormap that replicates those used in many geophysical software tools to display azimuth or phase data. This is better because it wraps around at 360 degrees but the perceptual issues are unresolved: in this case red, yellow and blue all appear as sharp perceptual edges.

standardAZ

Figure 3. Azimuth values color-coded with generic azimuth colormap.

 

isoAZ

Figure 4. Azimuth values color-coded with isoluminant azimuth colormap.

 

In Figure 4 I used my new colormap, called isoAZ (for isoluminant azimuth). This colormap is much better because not only does it wraps around at 360 degrees, but also lightness is held constant for all colors, which eliminates the perceptual anomalies. All the artificial yellow, red, and blue edges are gone, only real edges are left. This can be more easily appreciated in the figure below: if you hover with your mouse over it you are able to switch back and forth between Figure 3 and Figure 4.

isoAZstandardAZ

From an interpretation point of view, azimuths 180 degrees apart are of opposing colours, which is ideal for dip azimuth data because it allows us to easily recognize folds where dips of opposite direction are juxtaposed at an edge. One example is the sharp edge in the northwest quadrant of Figure 4, where magenta is juxtaposed to green. If you look at Figure 1 you see that there’s a relative high in this area (the edge in Figure 4) with dips of opposite direction on either side (East and West, or 0 and 360 degrees).

The colormap was created in the Lightness-Chroma-Hue color space, a polar transform of the Lab color space, where lightness is the vertical axis and at each value of lightness, chroma is the radial coordinate and hue the polar angle. One limitation of this approach is that due to theirregular  shape of the color gamut section at each lightness value, we can never exceed  chroma values  of about 38-40 (at lightness = 65 in Matlab; in Python, with extensive trial and error, I have not been able to go past 36 using the Scikit-image Color module), which make the resulting colors pale, pastely.

it creates For those that want to experiment with it further, I used just a few lines of code similar to the ones below:

radius = 38; % chroma
theta = linspace(0, 2*pi, 256)'; % hue
a = radius * cos(theta);
b = radius * sin(theta);
L = (ones(1, 256)*65)'; % lightness
Lab = [L, a, b];
RGB=colorspace('RGB<-Lab',Lab(end:-1:1,:));

This code is a modification from an example by Steve Eddins on a post on his Matlab Central blog. In Steve’s example the colormap cycles through the hues as lightness increases monotonically (which by the way is an excellent way to generate a perceptual rainbow). In this case lightness is kept constant and hue cycles through the entire 360 degrees and wraps around. Also, instead of using the Image Processing Toolbox, I used  Colorspace, a free function from Matlab File Exchange, for the color space transformations.

For data like fracture orientation where azimuths 180 degrees apart are equivalent it is better to stack two of these isoluminant colormaps in a row. In this way we place opposing colors 90 degrees apart, whereas color 180 degrees apart are the same. You can do it using Matlab commands repmat or vertcat, as below:

radius = 38; % chroma
theta = linspace(0, 2*pi, 128)'; % hue
a = radius * cos(theta);
b = radius * sin(theta);
L = (ones(1, 128)*65)'; % lightness
Lab = [L, a, b];
rgb=colorspace('RGB<-Lab',Lab(end:-1:1,:));
RGB=vertcat(rgb,rgb);

Colormap compromise

At the end the series The rainbow is dead…long live the rainbow! I introduced CubicYF, a rainbow-like color palette with a slightly compressive, monotonic increasing lightness profile. The red color is missing from the palette because green and red at people with deuteranopia, the most common type of color vision deficiency, confuse red and green, so I took one out. An alternative solution to taking the red out is to have a yellow-like color with lower lightness followed by a red-like color of higher lightness, like in the LinearL palette described The rainbow is dead…long live the rainbow! – Part 5. A third solution, which is a compromise in perceptual terms, is CubicL (sometimes I called it cube1), the color palette I used with the gravity data from my degree thesis in geology in the series Visualization tips for geoscientists. An example map from the series in which I used CubicL is in Figure 1 below:

data cube1_final_shading_slope

Figure 1

To generate CubicL  I used a compressive function for the increase in lightness very much like the one for CubicYF, except in this case I allowed lightness to decrease from about 90 in the yellow portion of the colormap to about 80 to get to the red. As a result, there is an inversion in the lightness trend, indicated by the arrow in the bottom row in Figure 2. However, I believe this is an acceptable compromise, and one that is pleasing for people accustomed to rainbow (it has red), because the inversion is smooth, with a gentle downward-facing concavity. I also believe this colormap is still a perceptual improvement over the rainbow or spectrum so commonly used in seismic visualization, which typically have at least 3 sharp lightness inversions (indicated by arrows in the top row of Figure 2 below).

spectrum_vs_cubeY

Figure 2

CubicL is one of the colormaps available in my Matlab File Exchange function Perceptually Improved Colormaps.

Parula: a new Matlab colormap

Steve Eddins of the Matwork just published a post announcing a new Matlab colormap to replace Jet. It is called Parula (more to come on his blog about this intriguing name).

parula

First impression: Parula looks good.

And while I haven’t had time to take it into Python to run a full perceptual test and into ImageJ for a colour blindness test, as a preliminary test I did convert it to grayscale with an online picture converting tool that uses the lightness information to perform the conversion (instad of just desaturating the colors) and the result shows monotonic changes in gray.

parula_compare

Looks promising… Full test to come.