% EQUALISECOLOURMAP - Equalise colour contrast over a colourmap
%
% Usage: newrgbmap = equalisecolourmap(rgblab, map, formula, W, sigma, diagnostics)
%
% Arguments:     rgblab - String 'RGB' or 'LAB' indicating the type of data
%                         in map.
%                   map - A Nx3 RGB or CIELAB colour map
%               formula - String 'CIE76' or 'CIEDE2000'
%                     W - A 3-vector of weights to be applied to the
%                         lightness, chroma and hue components of the
%                         difference equation. Use [1 0 0] to only take into
%                         account lightness, [1 1 1] for the full formula.
%                         See note below.
%                 sigma - Optional Gaussian smoothing parameter, see
%                         explanation below.
%                cyclic - Flag 1/0 indicating whether the colour map is
%                         cyclic. This affects how smoothing is applied at
%                         the end points.
%           diagnostics - Optional flag 0/1 indicating whether diagnostic
%                         plots should be displayed.  Defaults to 0.
%
% Returns:    newrgbmap - RGB colour map adjusted so that the perceptual
%                         contrast of colours along the colour map is constant.
%
% Suggested parameters:
%
% The CIE76 and CIEDE2000 colour difference formulas were developed for much
% lower spatial frequencies than we are typically interested in.  Neither is
% ideal for our application.  While the CIEDE2000 lightness correction is
% probably useful the difference is barely noticeable. The CIEDE2000 chroma
% correction is definitely *not* suitable for our needs. 
%
% For colour maps with a significant range of lightness use:
%                        formula = 'CIE76' or 'CIEDE2000'
%                              W = [1 0 0]  (Only correct for lightness)
%                          sigma = 5 - 7
%
% For isoluminant or low lightness gradient colour maps use:  
%                        formula = 'CIE76'
%                              W = [1 1 1]  (Correct for colour and lightness)
%                          sigma = 5 - 7
%
% Ideally, for a colour map to be effective the perceptual contrast along the
% colour map should be constant.  Many colour maps are very poor in this regard.
% Try testing your favourite colour map on the SINERAMP test image.  The
% perceptual contrast is very much dominated by the contrast in colour lightness
% values along the map.  This function attempts to equalise the chosen
% perceptual contrast measure along a colour map by stretching and/or
% compressing sections of the colour map.
%
% This function's primary use is for the correction of colour maps generated by
% CMAP however it can be applied to any colour map.  There are limitations to
% what this function can correct.  When applied to some of MATLAB's colour maps
% such as 'jet', 'hsv' and 'cool' you get colour discontinuity artifacts because
% these colour maps have segments that are nearly constant in lightness.
% However, it does a nice job of fixing up MATLAB's 'hot', 'winter', 'spring'
% and 'autumn' colour maps.  If you do see colour discontinuities in the
% resulting colour map try changing W from [1 0 0] to [1 1 1], or some
% intermediate weighting of [1 0.5 0.5], say.
%
% Difference formula: Neither CIE76 or CIEDE2000 difference measures are ideal
% for the high spatial frequencies that we are interested in.  Empirically I
% find that CIEDE2000 seems to give slightly better results on colour maps where
% there is a significant lightness gradient (this applies to most colour maps).
% In this case you would be using a weighting vector W = [1 0 0].  For
% isoluminant, or low lightness gradient colour maps where one is using a
% weighting vector W = [1 1 1] CIE76 should be used as the CIEDE2000 chroma
% correction is inapropriate for the spatial frequencies we are interested in.
%
% Weighting vetor W: The CIEDE2000 colour difference formula incorporates the
% scaling parameters kL, kC, kH in the demonimator of the lightness, chroma, and
% hue difference components respectively.  The 3 components of W correspond to
% the reciprocal of these 3 parameters.  (I do not know why they chose to put
% kL, kC, kH in the denominator. If you wanted to ignore, say, the chroma
% component you would have to set kC to Inf, rather than setting W(2) to 0 which
% seems more sensible to me).  If you are using CIE76 then W(2) amd W(3) are
% applied to the differences in a and b.  In this case you should ensure W(2) =
% W(3).  In general, for the spatial frequencies of interest to us, lightness
% differences are overwhelmingly more important than chroma or hue and W shoud
% be set to [1 0 0]
%
% Smoothing parameter sigma: 
% The output colour map will have lightness values of constant slope magnitude.
% However, it is possible that the sign of the slope may change, for example at
% the mid point of a bilateral colour map.  This slope discontinuity of lightness
% can induce a false apparent feature in the colour map.  A smaller effect is
% also occurs for slope discontinuities in a and b.  For such colour maps it can
% be useful to introduce a small amount of smoothing of the Lab values to soften
% the transition of sign in the slope to remove this apparent feature.  However
% in doing this one creates a small region of suppressed luminance contrast in
% the colour map which induces a 'blind spot' that compromises the visibility of
% features should they fall in that data range.  Accordingly the smoothing
% should be kept to a minimum.  A value of sigma in the range 5 to 7 in a 256
% element colour map seems about right.  As a guideline sigma should not be more
% than about 1/25 of the number of entries in the colour map, preferably less.
%
% If you use the CIEDE2000 formula option this function requires Gaurav Sharma's
% MATLAB implementation of the CIEDE2000 color difference formula deltaE2000.m
% available at: http://www.ece.rochester.edu/~gsharma/ciede2000/
%
% See also: CMAP, SINERAMP, COLOURMAPPATH, ISOCOLOURMAPPATH
%
% Reference:
% Peter Kovesi. Good Colour Maps: How to Design Them.
% arXiv:1509.03700 [cs.GR] 2015
% https://arxiv.org/abs/1509.03700

% Copyright (C) 2013-2014 Peter Kovesi
% Centre for Exploration Targeting
% The University of Western Australia
% peter.kovesi at uwa edu au
% 
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
% 
% The above copyright notice and this permission notice shall be included in 
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.
%
% December 2013  - Original version
% June     2014  - Generalised and cleaned up

function newrgbmap = equalisecolourmap(rgblab, map, formula, W, sigma, cyclic, diagnostics)
    
    if ~exist('formula', 'var'), formula = 'CIE76'; end
    if ~exist('W', 'var'), W = [1 0 0]; end
    if ~exist('sigma', 'var'), sigma = 0; end
    if ~exist('cyclic', 'var'), cyclic = 0; end
    if ~exist('diagnostics', 'var'), diagnostics = 0; end
    
    N = size(map,1);  % No of colour map entries

    if N/sigma < 25
        warning(['It is not recommended that sigma be larger than 1/25 of' ...
                 ' colour map length'])
    end

    if strcmpi(rgblab, 'RGB') && (max(map(:)) > 1.01  || min(map(:) < -0.01))
        error('If map is RGB values should be in the range 0-1')
    elseif strcmpi(rgblab, 'LAB') && max(abs(map(:))) < 10
        error('If map is LAB magnitude of values are expected to be > 10')
    end
    
    % Use D65 whitepoint to match typical monitors.
    wp = whitepoint('D65');
%    wp = whitepoint('D50');
    
    % If input is RGB convert colour map to Lab. 
    if strcmpi(rgblab, 'RGB')
        rgbmap = map;
        labmap = applycform(map, makecform('srgb2lab', 'AdaptedWhitePoint', wp));
        L = labmap(:,1);
        a = labmap(:,2);
        b = labmap(:,3);
    elseif strcmpi(rgblab, 'LAB')
        labmap = map;
        rgbmap = applycform(map, makecform('lab2srgb', 'AdaptedWhitePoint', wp));
        L = map(:,1);
        a = map(:,2);
        b = map(:,3);
    else
        error('Input must be RGB or LAB')
    end

    % The following section of code computes the locations to interpolate into
    % the colour map in order to achieve equal steps of perceptual contrast.
    % The process is repeated recursively on its own output. This helps overcome
    % the approximations induced by using linear interpolation to estimate the
    % locations of equal perceptual contrast. This is mainly an issue for
    % colour maps with only a few entries.

    for iter = 1:3
        % Compute perceptual colour difference values along the colour map using
        % the chosen formula and weighting vector.
        if strcmpi(formula, 'CIE76')
            deltaE = cie76(L, a, b, W);
        elseif strcmpi(formula, 'CIEDE2000')
            if ~exist('deltaE2000','file')
                fprintf('You need Gaurav Sharma''s function deltaE2000.m\n');
                fprintf('available at:\n');
                fprintf('http://www.ece.rochester.edu/~gsharma/ciede2000/\n');
                newrgbmap = [];
                return
            end
                
            deltaE = ciede2000(L, a, b, W);
        else
            error('Unknown colour difference formula')
        end
        
        % Form cumulative sum of of delta E values.  However, first ensure all
        % values are larger than 0.001 to ensure the cumulative sum always
        % increases.
        deltaE(deltaE < 0.001) = 0.001;
        cumdE = cumsum(deltaE);
        
        % Form an array of equal steps in cumulative contrast change.
        equicumdE =  (0:(N-1))'/(N-1) * (cumdE(end)-cumdE(1)) + cumdE(1);
        
        % Solve for the locations that would give equal Delta E values.
        method = 'linear';
        newN = interp1(cumdE, (1:N)', equicumdE, method, 'extrap');
        
        % newN now represents the locations where we want to interpolate into the
        % colour map to obtain constant perceptual contrast
        L = interp1((1:N)', L, newN, method, 'extrap');  
        a = interp1((1:N)', a, newN, method, 'extrap');
        b = interp1((1:N)', b, newN, method, 'extrap');
    
        % Record initial colour differences for evaluation at the end
        if iter == 1  
            initialdeltaE = deltaE;
            initialcumdE = cumdE;
            initialequicumdE = equicumdE;
            initialnewN = newN;
        end
    end
    
    % Apply smoothing of the path in CIELAB space if requested.  The aim is to
    % smooth out sharp lightness/colour changes that might induce the perception
    % of false features.  In doing this there will be some cost to the
    % perceptual contrast at these points.
    if sigma
        L = smooth(L, sigma, cyclic);
        a = smooth(a, sigma, cyclic);
        b = smooth(b, sigma, cyclic);
    end

    % Convert map back to RGB
    newlabmap = [L a b];
    newrgbmap = applycform(newlabmap, makecform('lab2srgb', 'AdaptedWhitePoint', wp));
    
    if diagnostics
        % Compute actual perceptual contrast values achieved
        if strcmpi(formula, 'CIE76')
            newdeltaE = cie76(L, a, b, W);
        elseif strcmpi(formula, 'CIEDE2000')
            newdeltaE = ciede2000(L, a, b, W);
        else
            error('Unknown colour difference formula')
        end
    
        show(sineramp,1,'Unequalised input colour map'), colormap(rgbmap);
        show(sineramp,2,'Equalised colour map'), colormap(newrgbmap);

        figure(3), clf
        subplot(2,1,1)
        plot(initialcumdE)
        title(sprintf('Cumulative change in %s of input colour map', formula));
        axis([1 N 0 1.05*max(cumdE(:))]);
        xlabel('Colour map index');

        subplot(2,1,2)
        plot(1:N, initialdeltaE, 1:N, newdeltaE)
        axis([1 N 0 1.05*max(initialdeltaE(:))]);
        legend('Original colour map', ...
               'Adjusted colour map', 'location', 'NorthWest');
        title(sprintf('Magnitude of %s differences along colour map', formula));
        xlabel('Colour map index');
        ylabel('dE')
    
        % Convert newmap back to Lab to check for gamut clipping
        labmap = applycform(newrgbmap,...
                            makecform('srgb2lab', 'AdaptedWhitePoint', wp));
        figure(4), clf
        subplot(3,1,3)
        plot(1:N, L-labmap(:,1),...
             1:N, a-labmap(:,2),...
             1:N, b-labmap(:,3))
        legend('Lightness', 'a', 'b', 'Location', 'NorthWest')
        maxe = max(abs([L-labmap(:,1); a-labmap(:,2); b-labmap(:,3)]));
        axis([1 N -maxe maxe]);
        title('Difference between desired and achieved L a b values (gamut clipping)')
        
        % Plot RGB values
        subplot(3,1,1)
        plot(1:N, newrgbmap(:,1),'r-',...
             1:N, newrgbmap(:,2),'g-',...
             1:N, newrgbmap(:,3),'b-')
    
        legend('Red', 'Green', 'Blue', 'Location', 'NorthWest')
        axis([1 N 0 1.1]);        
        title('RGB values along colour map')
        
        % Plot Lab values
        subplot(3,1,2)
        plot(1:N, L, 1:N, a, 1:N, b)
        legend('Lightness', 'a', 'b', 'Location', 'NorthWest')
        axis([1 N -100 100]);                
        title('L a b values along colour map')
    
        
        %% Extra plots to generate figures.
        % The following plots were developed for illustrating the
        % equalization of MATLAB's hot(64) colour map.  Run using:
        % >> hoteq = equalisecolourmap('RGB', hot, 'CIE76', [1 0 0], 0, 0, 1);

        fs = 12; % fontsize
        % Colour difference values
        figure(10), clf
        plot(1:N, initialdeltaE, 'linewidth', 2);
        axis([1 N 0 1.05*max(initialdeltaE(:))]);
        title(sprintf('Magnitude of %s lightness differences along colour map', ...
                      formula), 'fontweight', 'bold', 'fontsize', fs);
        
        xlabel('Colour map index', 'fontweight', 'bold', 'fontsize', fs);
        ylabel('dE', 'fontweight', 'bold', 'fontsize', fs);

        % Cumulative difference plot showing mapping of equicontrast colours
        figure(11), clf
        plot(initialcumdE, 'linewidth', 2), hold on
        s = 7;
        s = 10;
        for n = s:s:N
            line([0 initialnewN(n) initialnewN(n)],...
                 [initialequicumdE(n) initialequicumdE(n) 0], ...
                 'linestyle', '--', 'color', [0 0 0], 'linewidth', 1.5)
 
            ah = s/2; % Arrow height and width
            aw = ah/5;
            plot([initialnewN(n)-aw initialnewN(n) initialnewN(n)+aw],...
                 [ah 0 ah],'k-', 'linewidth', 1.5)            
        end
        
        title(sprintf('Cumulative change in %s lightness of colour map', formula), ...
              'fontweight', 'bold', 'fontsize', fs);
        axis([1 N+6 0 1.05*max(cumdE(:))]);
        xlabel('Colour map index', 'fontweight', 'bold', 'fontsize', fs);
        ylabel('Equispaced cumulative contrast', 'fontweight', 'bold', 'fontsize', fs);
        hold off
        
        
    end

%----------------------------------------------------------------------------
%
% Function to smooth an array of values but also ensure end values are not
% altered or, if the map is cyclic, ensures smoothing is applied across the end
% points in a cyclic manner.  Assumed that the input data is a column vector

function Ls = smooth(L, sigma, cyclic)
    
    sze = ceil(6*sigma);  if ~mod(sze,2), sze = sze+1; end
    G = fspecial('gaussian', [sze 1], sigma);
    
    if cyclic 
        Le = [L(:); L(:); L(:)]; % Form a concatenation of 3 repetitions of the array. 
        Ls = filter2(G, Le);     % Apply smoothing filter
                                 % and then return the center section    
        Ls = Ls(length(L)+1: length(L)+length(L));
    
    else  % Non-cyclic colour map: Pad out input array L at both ends by 3*sigma
          % with additional values at the same slope.  The aim is to eliminate
          % edge effects in the filtering
        extension = (1:ceil(3*sigma))';
        
        dL1 = L(2)-L(1);
        dL2 = L(end)-L(end-1);
        Le = [-flipud(dL1*extension)+L(1); L;  dL2*extension+L(end)];
        
        Ls = filter2(G, Le);  % Apply smoothing filter
        
        % Trim off extensions
        Ls = Ls(length(extension)+1 : length(extension)+length(L));
    end
    
%----------------------------------------------------------------------------
% Delta E using the CIE76 formula + weighting

function deltaE = cie76(L, a, b, W)
    
    N = length(L);
    
    % Compute central differences 
    dL = zeros(size(L));
    da = zeros(size(a));
    db = zeros(size(b));

    dL(2:end-1) = (L(3:end) - L(1:end-2))/2;
    da(2:end-1) = (a(3:end) - a(1:end-2))/2;
    db(2:end-1) = (b(3:end) - b(1:end-2))/2;
    
    % Differences at end points
    dL(1) = L(2) - L(1);  dL(end) = L(end) - L(end-1);
    da(1) = a(2) - a(1);  da(end) = a(end) - a(end-1);
    db(1) = b(2) - b(1);  db(end) = b(end) - b(end-1);
    
    deltaE = sqrt(W(1)*dL.^2 + W(2)*da.^2 + W(3)*db.^2);
    
%----------------------------------------------------------------------------
% Delta E using the CIEDE2000 formula + weighting
%
% This function requires Gaurav Sharma's MATLAB implementation of the CIEDE2000
% color difference formula deltaE2000.m available at:
% http://www.ece.rochester.edu/~gsharma/ciede2000/

function deltaE = ciede2000(L, a, b, W)

    N = length(L);
    Lab = [L(:) a(:) b(:)];
    KLCH = 1./W;
    
    % Compute deltaE using central differences
    deltaE = zeros(N, 1);
    deltaE(2:end-1) = deltaE2000(Lab(1:end-2,:), Lab(3:end,:), KLCH)/2; 
   
    % Differences at end points    
    deltaE(1) = deltaE2000(Lab(1,:), Lab(2,:), KLCH);
    deltaE(end) = deltaE2000(Lab(end-1,:), Lab(end,:), KLCH);
    
