Oh, the new, perceptual MatplotLib colormaps…..
Here’s one stunning, recent Truecolor image of Pluto from the New Horizons mission:
Original image: The Rich Color Variations of Pluto. Credit: NASA/JHUAPL/SwRI. Click on the image to view the full feature on New Horizon’s site
Below, I recolored using two of the new colormaps:
Recolored images: I like Viridis, by it is Inferno that really brings to life this image, because of its wider hue and lightness range!
I started this blog in 2012; in these 3 1/2 years it has been a wonderful way to channel some of my interests in image processing, geophysics, and visualization (in particular colour), and more recently Python.
During this time, among other things, I learned how to build and maintain a blog, I packaged a popular Matlab function, wrote an essay for Agile Geoscience’s first book on Geophysics, presented at the 2012 CSEG Geoconvention, and wrote two tutorials for The Leading Edge. Last, but not least, I made many new friends and professional connections.
Starting with 2016 I would like to concentrate my efforts on building useful (hopefully) and fun (for me at least) open source (this one is for sure) tools in Python. This is going to be my modus operandi:
Here are a couple of examples of ongoing projects:
The idea for this project was inspired by Matt Hall of Agile Geoscience. The plan is to eventually build a web app that will:
This is a project started at the 2015 Calgary Geoscience Hackathon organized by Agile Geoscience with Elwyn Galloway, Evan Saltman, and Ben Bougher. The original idea, proposed by Elwyn at the Hackathon, was to make an app that would turn an image of geological sketch into a model, as in the figure below.
The implementation of the finished app involves using morphological filtering and other image processing methods to enhance the sketch image and convert it into a model with discrete bodies, then pass it on to Agile’s modelr.io to create a synthetic.
The remote sensing used to study the human body is very similar to the remote sensing used to study the subsurface. Apart from a scaling factor (due to the different frequencies of the signals used) the only major difference between the two methods of investigation is in that radiologists and doctors looking at an x-ray, ultrasound, or CAT scan image know what to look for in those images, as bones, tissues, and anomalies, have known characteristics, whereas the subsurface is always to a large extent unknown.
In this short visual post I am going to use a CAT scan of King Tut’s skull to explore the effect on the image quality of progressive decimation of the data followed by upsampling it back to the initial size. I will also look at the effect of these manipulations on the results of edge detection.
With this I want to simulate the progressive reduction in imaging quality that happens when going from high density 3D seismic acquisition to medium density 3D seismic to high quality, but sparse 2D seismic lines.
Here’s the input image in Figure 1.
Figure 1. CAT scan of King Tut’s skull – Supreme Council of Antiquities. guardians.net/hawass/press_release_tutankhamun_ct_scan_results.htm
In Figure 2 I am showing the image after import into a Jupiter Notebook and conversion to grayscale, and the result of edge detection using the Sobel filter. Notice the excellent quality of the edge detection result.
Figure 2. Original image, or ground truth for the experiment, and edge detection result.
To simulate a high-resolution 3D seismic acquisition I decimated the original image by a factor of 4 in both directions. The resulting image (no interpolation) is, shown in Figure 3, is of good quality, and so is the edge detection result.
Figure 3. Simulated high-resolution 3D survey and edge detection result.
The image in Figure 4 results from a further decimation by a factor of 2 of the image in Figure 3, then interpolation to upsample to the same size as the image in Figure 4. The image and the edge detection are still of fair quality overall, but some of the smaller features have either disappeared, merged, or faded.
Figure 4. Simulated medium resolution 3D survey and edge detection result.
Now look at Figure 5: this is the equivalent of a high quality (in one direction) 2D dataset. Although we can still guess at what this represents, I would argue this is a result of our a priori knowledge of what it is supposed to represent – a human skull; and yet I don’t think anybody would want their doctor to make a diagnosis based on this image.
Figure 5. Simulated set of very high-resolution 2D lines.
The image in Figure 6 results from 2D interpolation (my intention is to simulate the result we would get by gridding 2D data to get a continuous image. We can now definitely interpret this as a skull, but the edge detection result is very unsatisfactory.
Figure 6. Simulated interpolation of 2D lines.
In future post we will explore the effects of adding periodic noise (similar to seismic acquisition footprint) on these images and on the edge detection results. I will also show you how to remove it using 2D FFT filters, as promised (now more than a year ago) in my post Moiré Patterns.
In Visualization tips for geoscientists: MATLAB, Part III I showed there’s a qualitative correlation between occurrences of antimony mineralization in the southern Tuscany mining district and the distance from lineaments derived from the total horizontal derivative (also called maximum horizontal gradient).
Let’s take a look at the it below (distance from lineaments increases as the color goes from blue to green, yellow, and then red).
However, in a different map in the same post I showed that lineaments derived using the maxima of the hyperbolic tilt angle (Cooper and Cowan, 2006, Enhancing potential field data using filters based on the local phase) are offset systematically from those derived using the total horizontal derivative.
Let’s take a look at the it below: in this case Bouguer gravity values increase as the color goes from blue to green, yellow, and then red; white polygons are basement outcrops.
The lineaments from the total horizontal derivative are in black, those from the maxima of hyperbolic tilt angle are in gray. Which lineaments should be used?
The ideal way to map the location of density contrast edges (as a proxy for geological contacts) would be to gravity forward models, or even 3D gravity inversion, ideally constrained by all available independent data sources (magnetic or induced-polarization profiles, exploratory drilling data, reflection seismic interpretations, and so on).
The next best approach is to map edges using a number of independent gravity data enhancements, and then only use those that collocate.
Cooper and Cowan (same 2006 paper) demonstrate that no single-edge detector method is a perfect geologic-contact mapper. Citing Pilkington and Keating (2004, Contact mapping from gridded magnetic data – A comparison of techniques) they conclude that the best approach is to use “collocated solutions from different methods providing increased confidence in the reliability of a given contact location”.
I show an example of such a workflow in the image below. In the first column from the left is a map of the residual Bouguer gravity from a smaller area of interest in the southern Tuscany mining district (where measurements were made on a denser station grid). In the second column from the left are the lineaments extracted using three different (and independent) derivative-based data enhancements followed by skeletonization. The same lineaments are superimposed on the original data in the third column from the left. Finally, in the last column, the lineaments are combined into a single collocation map to increase confidence in the edge locations (I applied a mask so as to display edges only where at least two methods collocate).
If you want to learn more about this method, please read my note in the Geophysical tutorial column of The Leading Edge, which is available with open access here.
To run the open source Python code, download the iPython/Jupyter Notebook from GitHub.
With this notebook you will be able to:
1) create a full suite of derivative-based enhanced gravity maps;
2) extract and refine lineaments to map edges;
3) create a collocation map.
These technique can be easily adapted to collocate lineaments derived from seismic data, with which the same derivative-based enhancements are showing promising results (Russell and Ribordy, 2014, New edge detection methods for seismic interpretation.)
In the first post of this series I argued that we should not build colormaps for azimuth (or phase) data by interpolating linearly between fully saturated hues in RGB or HSL space.
A first step towards the ideal colormap for azimuth data would be to interpolate between isoluminant colours instead. Kindlmann et al. (2002) published isoluminant RGB values for red, yellow, green, cyan, blue, and magenta based on a user study. The code in the next block show how to interpolate between those published colours to get 256-sample R, G, and B arrays (with magenta repeated at both ends), which can then be combined into a isoluminant colormap for azimuth data.
01 r = np.array([0.718, 0.847, 0.527, 0.000, 0.000, 0.316, 0.718]) 02 g = np.array([0.000, 0.057, 0.527, 0.592, 0.559, 0.316, 0.000]) 03 b = np.array([0.718, 0.057, 0.000, 0.000, 0.559, 0.991, 0.718]) 04 x = np.linspace(0,256,7) 05 xnew = np.arange(256) 06 r256 = np.interp(xnew, x, r) 07 g256 = np.interp(xnew, x, g) 08 b256 = np.interp(xnew, x, b)
This is a good example in general of how to interpolate to a finer sampling one or more sequence of values using the interp function from the Numpy library. In line 04 we define 7 evenly spaced points between 0 and 255; this will be the sample coordinate for the r, g, and b colours created in lines 01-03. In line 05 we create the new coordinates we will be interpolating r, g, and b values at in lines 06-08 (all integers between 0 and 256). The full code will come in the Notebook accompanying the last post in this series.
This new colormap is used in the bottom map of the figure below, whereas in the top map we used a conventional HSV azimuth colormap (both maps show the dip azimuth calculated on the Penobscot horizon). The differences are subtle, but with the isoluminant colormap we are guaranteed there are no perceptual artifacts due to the random variations in lightness of the fully saturated HSV colors.
Another possible strategy to create a perceptual colormap for azimuth data would be to set lightness and chroma to constant values in LCH space and interpolate between hues. This is the Matlab colormap I previously created, and shown in Figure 4 of New Matlab isoluminant colormap for azimuth data. In the next post, I will show how to do this in Python.
The last two posts on Agile show you how to corender seismic amplitude and continuity from a time slice using a 2D colormap, and then how to corender 3 attributes from a horizon slice.
Kindlmann, E. et al. (2002). Face-based Luminance Matching for Perceptual Colour map Generation – Proceedings of the IEEE conference on Visualization.
In New Matlab isoluminant colormap for azimuth data I showcased a Matlab colormap that I believe is perceptually superior to the conventional, HSV-based colormaps for azimuth data, in that it does not superimposes on the data the color artifacts that plague all rainbows. However, it still has a limitation, which is that the main colours do not correspond exactly to the four compass directions N, E, W, and S.
My intention with this series is to go back to square one, deconstruct the conventional colormaps for azimuth, and build a new one that has all the desired properties of both perceptual linearity, and correct location of the main colors. All reproducible in Python.
If we wanted to build from scratch a colormap for azimuth (or phase) data the main tasks would be to generate a sequence of distinguishable colours at opposite quadrants, or compass directions (like 0 and 180 degrees, or N and S), and to wrap around the sequence with the same colour at the two ends.
But to do that, we should avoid interpolating linearly between fully saturated hues in RGB or HSL space.
To illustrate why, it is useful to look at the figure below. On the left is a hue circle with primary, secondary, and tertiary colours in a counter-clockwise sequence: red, rose, magenta, violet, blue, azure, cyan, aquamarine, electric green, chartreuse, yellow, and orange. The colour chips are placed at evenly spaced angular distances according to their hue (in radians).
Left, primary, secondary, and tertiary colour chips arranged using hue for angular distance; right, the same colour chips arranged using intensity for angular distance.
This looks familiar and seems like a natural ordering of colors, so we may be tempted in building a colormap, to just take that sequence, wrap it around at the red (or the magenta) and linearly interpolate to 256 colours to get a continuous colormap [1], and use it for azimuth data, which is how usually the conventional azimuth colormaps are built.
On the right side in the figure the chips have been rearranged according to their intensity on a counter-clockwise sequence from 0 to 255 with 0 at three hours; so, for example blue, which is the darkest colour with an intensity of 29, is close to the beginning of the sequence, and yellow, the brightest with an intensity of 225, is close to the end. Notice that the chips are no longer equidistant.
The most striking is that the blue and the yellow chips are more separated than the other chips, and for this reason blue and yellow features seem to stand out a lot more in a map when using this color sequence, which can be both distracting and confusing. A good example is Figure 3 in New Matlab isoluminant colormap for azimuth data.
Also, yellow and red, being two chips apart in the left circle in the figure above, are used to colour azimuths 60 degrees apart, and so do cyan and green. However, if we look at the right circle, we realize that the yellow and red chips are much further apart than the cyan and green chips [2] in the perceptual dimension of intensity; therefore, features colored in yellow and red could be perceived as much further apart (in azimuth) than cyan and green.
These differences may be subtle, but in my opinion they become important when dip azimuth is combined with other attributes, perhaps using a 3D colormap, and the resulting map is used for detailed structural interpretation. There is a really good example of this type of 3D colormap in Chopra and Marfurt (2007), where dip azimuth is rendered with hue modulation, dip magnitude with saturation modulation, and coherence with lightness modulation.
A code snippet with the main Python commands to generate the two polar scatterplots in the figure is listed, and explained below. The full code can be found in this Jupiter Notebook.
01 import matplotlib.colors as clr 02 keys=['red', '#FF007F', 'magenta', '#7F00FF', 'blue', '#0080FF','cyan', '#00FF80', '#00FF00', '#7FFF00', 'yellow', '#FF7F00'] 03 my_cmap = clr.ListedColormap(keys) 04 x = np.arange(12) 05 color = my_cmap(x) 06 n = 12 07 theta = 2*np.pi*(np.linspace(0,1,13)) 08 r = np.ones(13)*2.5 09 area = 200*r**2 # size of color chips 10 c = plt.scatter(theta, r, c=color, s=area) 11 theta_i = 2*np.pi*(sorted_intensity/255.0) 12 colors = my_sorted_cmap(np.arange(12)) 13 c = plt.scatter(theta_i, r, c=colors, s=area)
In line 01 we import the Colors module from the Matplotlib library, then line 02 creates the desired sequence of colours (red, rose, magenta, violet, blue, azure, cyan, aquamarine, electric green, chartreuse, yellow, and orange) using either the name or Hex code, and line 03 generates the colormap. Then we use lines 04 and 05 to assign colours to the chips in the first scatterplot (left), and lines 06, 07, and 09 to specify the number of chips, the angular distances between chips, and the area of the chips, respectively. Line 10 generates the plot. The modifications in lines 11-14 will result in the scatterplot on the right side in the figure (the sorted intensity is calculated in much the same way as in my Geophysical tutorial – How to evaluate and compare colormaps in Python).
[1] Or, perhaps, just create 12 discrete colour classes to group azimuth values in bins of pi/6 (30 degrees) each, and wrap around again at the magenta, to generate a discrete colormap.
[2] The green chip is almost completely covered by the orange chip.
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.
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.
Figure 2: logspiral colormap trajectory in CIELab a-b plane
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).
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');