function varargout = cptcmap(varargin)
%CPTCMAP Apply a .cpt file as colormap to an axis
%
% cptcmap(name);
% cptcmap(name, ax);
% cptcmap(... param, val, ...);
% [cmap, lims, ticks, bfncol, ctable] = cptcmap(...)
%
% This function creates and applies a colormap defined in a color palette
% table (.cpt file). For a full description of the cpt file format, see
% the Generic Mapping Tools documentation (http://gmt.soest.hawaii.edu/).
% Color palette files provide more flexible colormapping than Matlab's
% default schemes, including both discrete and continuous gradients, as
% well as easier direct color mapping.
%
% Limitations: X11 color names not supported, patterns not supported, CMYK
% not supported yet
%
% Input variables:
%
% name: .cpt file name. You may either specify either the full
% path, or just the file name. In the latter case, the
% function will look for the file in the folder specified by
% the cptpath variable in the first line of code; by default
% this folder is located in the same location as cptcmap.m
% and is called cptfiles.
%
% ax: handle of axis or axes where colormap should be applied
% (colormaps will effect the entire figure(s), but axis clim
% adjustments for direct scaling will only affect the
% specified axes). If no axis is specified and no output
% variables are supplied, colormap will be applied to the
% current axis. If no axis is specified and output variables
% are supplied, the colormap will not be applied to any axes.
%
% 'showall': When this option is used, a figure is created displaying
% colorbars for all colormaps contained in the .cpt folder.
% Color limits of each colormap are listed along with the
% names of each. A small tick mark indicates the location of
% 0, where applicable. NOTE: the number of columns to use
% for display is hard-coded. As you start collecting more
% color palettes, the figure may get too cluttered and you
% may have to adjust this (variable ncol is the plotcmaps
% subfunction).
%
% Optional input variables (passed as parameter/value pairs):
%
% 'mapping': 'scaled' or 'direct'. Scaled mapping spreads the colormap
% to cover the color limits of the figure. Direct mapping
% resets the color limits of the axes so that colors are
% mapped to the levels specified by the .cpt file. ['scaled']
%
% 'ncol': number of colors in final colormap. If not included or NaN,
% this function will try to choose the fewest number of
% blocks needed to display the colormap as accurately as
% possible. I have arbitrarily chosen that it will not try to
% create more than 256 colors in the final colormap when
% using this automatic scheme. However, you can manually set
% ncol higher if necessary to resolve all sharp breaks and
% gradients in the colormap.
%
% 'flip': if true, reverse the colormap order [false]
%
% Output variables:
%
% cmap: ncol x 3 colormap array
%
% lims: 1 x 2 array holding minimum and maximum values for which
% the colormap is defined.
%
% ticks: vector of tick values specifying where colors were defined
% in the original file
%
% bfncol: 3 x 3 colormap array specifying the colors defined for
% background (values lower than lowest color limit),
% foreground (values higher than highest color limit), and
% NaN values. These do not affect the resulting colormap,
% but can be applied by the user to replicate the behavior
% seen in GMT.
%
% ctable: n x 8 color palette table, translated to Matlab color
% space. Column 1 holds the lower limit of each color cell,
% columns 2-4 the RGB values corresponding to the lower
% limit, column 5 the upper limit of the color cell, and
% columns 6-8 the RGB values of the upper limit. When the
% lower and upper colors are the same, this defines a
% solid-colored cell; when they are different, colors are
% linearly interpolated between the endpoints.
%
% Example:
%
% [lat, lon, z] = satbath(10);
% pcolor(lon, lat, z);
% shading flat;
% cptcmap('GMT_globe', 'mapping', 'direct');
% colorbar;
% Copyright 2011 Kelly Kearney
% File Exchange update: 4/12/2011
%------------------------------
% Parse input
%------------------------------
% The .cpt file folder. By default, the cptfiles folder is located in
% the same place as the cptcmap.m file. If you change this on your
% computer, change the cptpath definition to reflect the new location.
cptpath = fullfile(fileparts(which('cptcmap')), 'cptfiles');
if ~exist(cptpath, 'dir')
error('You have moved the cptfiles directory. Please modify the cptpath variable in this code to point to the directory where your.cpt files are stored');
end
% Check for 'showall' option first
if nargin == 1 & strcmp(varargin{1}, 'showall')
plotcmaps(cptpath);
return;
end
% Name of file
[blah, blah, ext] = fileparts(varargin{1});
if isempty(ext)
varargin{1} = [varargin{1} '.cpt'];
end
if exist(varargin{1}, 'file') % full filename and path given
filename = varargin{1};
else % only file name given
[blah,blah,ext] = fileparts(varargin{1});
if ~isempty(ext) % with extension
filename = fullfile(cptpath, varargin{1});
else % without extension
filename = fullfile(cptpath, [varargin{1} '.cpt']);
end
if ~exist(filename, 'file')
error('Specified .cpt file not found');
end
end
% Axes to which colormap will be applied
if nargin > 1 && isnumeric(varargin{2}) && all(ishandle(varargin{2}(:)))
ax = varargin{2};
pv = varargin(3:end);
applycmap = true;
elseif nargout == 0
ax = gca;
pv = varargin(2:end);
applycmap = true;
else
pv = varargin(2:end);
applycmap = false;
end
% Optional paramter/value pairs
p = inputParser;
p.addParamValue('mapping', 'scaled', @(x) any(strcmpi(x, {'scaled', 'direct'})));
p.addParamValue('ncol', NaN, @(x) isscalar(x) && isnumeric(x));
p.addParamValue('flip', false, @(x) isscalar(x) && islogical(x));
p.parse(pv{:});
Opt = p.Results;
%------------------------------
% Calculate colormap and apply
%------------------------------
[cmap, lims,ticks,bfncol,ctable] = cpt2cmap(filename, Opt.ncol);
if Opt.flip
if strcmp(Opt.mapping, 'direct')
warning('Flipping colormap with direct mapping may lead to odd color breaks');
end
cmap = flipud(cmap);
end
if applycmap
for iax = 1:numel(ax)
axes(ax(iax));
if strcmp(Opt.mapping, 'direct')
set(ax(iax), 'clim', lims);
end
colormap(cmap);
end
end
%------------------------------
% Output
%------------------------------
allout = {cmap, lims, ticks, bfncol, ctable};
varargout(1:nargout) = allout(1:nargout);
%------------------------------
% Subfunction: Read colormap
% from file
%------------------------------
function [cmap, lims, ticks, bfncol, ctable] = cpt2cmap(file, ncol)
% Read file
fid = fopen(file);
txt = textscan(fid, '%s', 'delimiter', '\n');
txt = txt{1};
fclose(fid);
isheader = strncmp(txt, '#', 1);
isfooter = strncmp(txt, 'B', 1) | strncmp(txt, 'F', 1) | strncmp(txt, 'N', 1);
% Extract color data, ignore labels (errors if other text found)
ctabletxt = txt(~isheader & ~isfooter);
ctable = str2num(strvcat(txt(~isheader & ~isfooter)));
if isempty(ctable)
nr = size(ctabletxt,1);
ctable = cell(nr,1);
for ir = 1:nr
ctable{ir} = str2num(strvcat(regexp(ctabletxt{ir}, '[\d\.-]*', 'match')))';
end
try
ctable = cell2mat(ctable);
catch
error('Cannot parse this format .cpt file yet');
end
end
% Determine which color model is used (RGB, HSV, CMYK, names, patterns,
% mixed)
[nr, nc] = size(ctable);
iscolmodline = cellfun(@(x) ~isempty(x), regexp(txt, 'COLOR_MODEL'));
if any(iscolmodline)
colmodel = regexprep(txt{iscolmodline}, 'COLOR_MODEL', '');
colmodel = strtrim(lower(regexprep(colmodel, '[#=]', '')));
else
if nc == 8
colmodel = 'rgb';
elseif nc == 10
colmodel = 'cmyk';
else
error('Cannot parse this format .cpt file yet');
end
end
% try
% temp = str2num(strvcat(txt(~isheader & ~isfooter)));
% if size(temp,2) == 8
% colmodel = 'rgb';
% elseif size(temp,2) == 10
% colmodel = 'cmyk';
% else % grayscale, maybe others
% error('Cannot parse this format .cpt file yet');
% end
% catch % color names, mixed formats, dash placeholders
% error('Cannot parse this format .cpt file yet');
% end
% end
%
%
% iscmod = strncmp(txt, '# COLOR_MODEL', 13);
%
%
% if ~any(iscmod)
% isrgb = true;
% else
% cmodel = strtrim(regexprep(txt{iscmod}, '# COLOR_MODEL =', ''));
% if strcmp(cmodel, 'RGB')
% isrgb = true;
% elseif strcmp(cmodel, 'HSV')
% isrgb = false;
% else
% error('Unknown color model: %s', cmodel);
% end
% end
% Reformat color table into one column of colors
cpt = zeros(nr*2, 4);
cpt(1:2:end,:) = ctable(:,1:4);
cpt(2:2:end,:) = ctable(:,5:8);
% Ticks
ticks = unique(cpt(:,1));
% Choose number of colors for output
if isnan(ncol)
endpoints = unique(cpt(:,1));
% For gradient-ed blocks, ensure at least 4 steps between endpoints
issolid = all(ctable(:,2:4) == ctable(:,6:8), 2);
for ie = 1:length(issolid)
if ~issolid(ie)
temp = linspace(endpoints(ie), endpoints(ie+1), 11)';
endpoints = [endpoints; temp(2:end-1)];
end
end
endpoints = sort(endpoints);
% Determine largest step size that resolves all endpoints
space = diff(endpoints);
space = unique(space);
space = roundn(space, -3); % To avoid floating point issues when converting to integers
space = round(space*1e3)/1e3;
nspace = length(space);
if ~isscalar(space)
fac = 1;
tol = .001;
while 1
if all(space >= 1 & (space - round(space)) < tol)
space = round(space);
break;
else
space = space * 10;
fac = fac * 10;
end
end
pairs = nchoosek(space, 2);
np = size(pairs,1);
commonsp = zeros(np,1);
for ip = 1:np
commonsp(ip) = gcd(pairs(ip,1), pairs(ip,2));
end
space = min(commonsp);
space = space/fac;
end
ncol = (max(endpoints) - min(endpoints))./space;
ncol = min(ncol, 256);
end
% Remove replicates and mimic sharp breaks
isrep = [false; ~any(diff(cpt),2)];
cpt = cpt(~isrep,:);
difc = diff(cpt(:,1));
minspace = min(difc(difc > 0));
isbreak = [false; difc == 0];
cpt(isbreak,1) = cpt(isbreak,1) + .01*minspace;
% Parse background, foreground, and nan colors
footer = txt(isfooter);
bfncol = nan(3,3);
for iline = 1:length(footer)
if strcmp(footer{iline}(1), 'B')
bfncol(1,:) = str2num(regexprep(footer{iline}, 'B', ''));
elseif strcmp(footer{iline}(1), 'F')
bfncol(2,:) = str2num(regexprep(footer{iline}, 'F', ''));
elseif strcmp(footer{iline}(1), 'N')
bfncol(3,:) = str2num(regexprep(footer{iline}, 'N', ''));
end
end
% Convert to Matlab-format colormap and color limits
lims = [min(cpt(:,1)) max(cpt(:,1))];
endpoints = linspace(lims(1), lims(2), ncol+1);
midpoints = (endpoints(1:end-1) + endpoints(2:end))/2;
cmap = interp1(cpt(:,1), cpt(:,2:4), midpoints);
switch colmodel
case 'rgb'
cmap = cmap ./ 255;
bfncol = bfncol ./ 255;
ctable(:,[2:4 6:8]) = ctable(:,[2:4 6:8]) ./ 255;
case 'hsv'
cmap(:,1) = cmap(:,1)./300;
cmap = hsv2rgb(cmap);
bfncol(:,1) = bfncol(:,1)./300;
bfncol = hsv2rgb(bfncol);
ctable(:,2) = ctable(:,2)./300;
ctable(:,6) = ctable(:,6)./300;
ctable(:,2:4) = hsv2rgb(ctable(:,2:4));
ctable(:,6:8) = hsv2rgb(ctable(:,6:8));
case 'cmyk'
error('CMYK color conversion not yet supported');
end
% Rouding issues: occasionally, the above calculations lead to values just
% above 1, which colormap doesn't like at all. This is a bit kludgy, but
% should solve those issues
isnear1 = cmap > 1 & (abs(cmap-1) < 2*eps);
cmap(isnear1) = 1;
%------------------------------
% Subfunction: Display all
% colormaps
%------------------------------
function plotcmaps(folder)
Files = dir(fullfile(folder, '*.cpt'));
nfile = length(Files);
ncol = 3;
nr = ceil(nfile/ncol);
width = (1 - .05*2)/ncol;
height = (1-.05*2)/nr;
left = .05 + (0:ncol-1)*width;
bot = .05 + (0:nr-1)*height;
[l, b] = meshgrid(left, bot);
w = width * .8;
h = height * .4;
figure('color','w');
ax = axes('position', [0 0 1 1]);
hold on;
for ifile = 1:nfile
[cmap,blah,blah,blah,ctable] = cptcmap(Files(ifile).name);
[x,y,c] = ctable2patch(ctable);
xtick = unique(x);
dx = max(x(:)) - min(x(:));
xsc = ((x-min(xtick))./dx).*w + l(ifile);
ysc = y.*h + b(ifile);
xrect = [0 1 1 0 0] .*w + l(ifile);
yrect = [1 1 0 0 1] .*h + b(ifile);
xticksc = ((xtick-min(xtick))./dx).*w + l(ifile);
x0 = interp1(xtick, xticksc, 0);
y0 = b(ifile) + [0 .2*h NaN .8*h h];
x0 = ones(size(y0))*x0;
lbl = sprintf('%s [%g, %g]', regexprep(Files(ifile).name,'\.cpt',''), min(x(:)), max(x(:)));
patch(xsc, ysc, c, 'edgecolor', 'none');
line(xrect, yrect, 'color', 'k');
line(x0, y0, 'color', 'k');
text(l(ifile), b(ifile)+h, lbl, 'interpreter', 'none', 'fontsize', 8, 'verticalalignment', 'bottom', 'horizontalalignment', 'left');
end
set(ax, 'ylim', [0 1], 'xlim', [0 1], 'visible', 'off');
% Determine patch coordinates
function [x,y,c] = ctable2patch(ctable)
np = size(ctable,1);
x = zeros(4, np);
y = zeros(4, np);
c = zeros(4, np, 3);
y(1:2,:) = 1;
for ip = 1:np
x(:,ip) = [ctable(ip,1) ctable(ip,5) ctable(ip,5) ctable(ip,1)];
c(:,ip,:) = [ctable(ip,2:4); ctable(ip,6:8); ctable(ip,6:8); ctable(ip,2:4)];
end
