sc.m

SC Display/output truecolor images with a range of colormaps

function I = sc(I, varargin)
%SC  Display/output truecolor images with a range of colormaps
%
% Examples:
%   sc(image)
%   sc(image, limits)
%   sc(image, map)
%   sc(image, limits, map)
%   sc(image, map, limits)
%   sc(..., col1, mask1, col2, mask2,...)
%   out = sc(...)
%   sc
%
% Generates a truecolor RGB image based on the input values in 'image' and
% any maximum and minimum limits specified, using the colormap specified.
% The image is displayed on screen if there is no output argument.
% 
% SC has these advantages over MATLAB image rendering functions:
%   - images can be displayed or output; makes combining/overlaying images
%     simple.
%   - images are rendered/output in truecolor (RGB [0,1]); no nasty
%     discretization of the input data.
%   - many special, built-in colormaps for viewing various types of data.
%   - linearly interpolates user defined linear and non-linear colormaps.
%   - no border and automatic, integer magnification (unless figure is
%     docked or maximized) for better display.
%   - multiple images can be generated for export simultaneously.
%
% For a demonstration, simply call SC without any input arguments.
%
% IN:
%   image - MxNxCxP or 3xMxNxP image array. MxN are the dimensions of the
%           image(s), C is the number of channels, and P the number of
%           images. If P > 1, images can only be exported, not displayed.
%   limits - [min max] where values in image less than min will be set to
%            min and values greater than max will be set to max.
%   map - Kx3 or Kx4 user defined colormap matrix, where the optional 4th
%         column is the relative distance between colours along the scale,
%         or a string containing the name of the colormap to use to create
%         the output image. Default: 'none', which is RGB for 3-channel
%         images, grayscale otherwise.  Conversion of multi-channel images
%         to intensity for intensity-based colormaps is done using the L2
%         norm. Most MATLAB colormaps are supported. All named colormaps
%         can be reversed by prefixing '-' to the string. This maintains
%         integrity of the colorbar. Special, non-MATLAB colormaps are:
%      'contrast' - a high contrast colormap for intensity images that
%                   maintains intensity scale when converted to grayscale,
%                   for example when printing in black & white.
%      'prob' - first channel is plotted as hue, and the other channels
%               modulate intensity. Useful for laying probabilites over
%               images.
%      'prob_jet' - first channel is plotted as jet colormap, and the other
%                   channels modulate intensity.
%      'diff' - intensity values are marked blue for > 0 and red for < 0.
%               Darker colour means larger absolute value. For multi-
%               channel images, the L2 norm of the other channels sets
%               green level. 3 channel images are converted to YUV and
%               images with more that 3 channels are projected onto the
%               principle components first.
%      'compress' - compress many channels to RGB while maximizing
%                   variance.
%      'flow' - display two channels representing a 2d Cartesian vector as
%               hue for angle and intensity for magnitude (darker colour
%               indicates a larger magnitude).
%      'phase' - first channel is intensity, second channel is phase in
%                radians. Darker colour means greater intensity, hue
%                represents phase from 0 to 2 pi.
%      'stereo' - pair of concatenated images used to generate a red/cyan
%                 anaglyph.
%      'stereo_col' - pair of concatenated RGB images used to generate a
%                     colour anaglyph.
%      'rand' - gives an index image a random colormap. Useful for viewing
%               segmentations.
%      'rgb2gray' - converts an RGB image to grayscale in the same fashion
%                   as MATLAB's rgb2gray (in the image processing toolbox).
%   col/mask pairs - Pairs of parameters for coloring specific parts of the
%                    image differently. The first (col) parameter can be
%                    a MATLAB color specifier, e.g. 'b' or [0.5 0 1], or
%                    one of the colormaps named above, or an MxNx3 RGB
%                    image. The second (mask) paramater should be an MxN
%                    logical array indicating those pixels (true) whose
%                    color should come from the specified color parameter.
%                    If there is only one col parameter, without a mask
%                    pair, then mask = any(isnan(I, 3)), i.e. the mask is
%                    assumed to indicate the location of NaNs. Note that
%                    col/mask pairs are applied in order, painting over
%                    previous pixel values.
%
% OUT:
%   out - MxNx3xP truecolour (double) RGB image array in range [0, 1]
%
% See also IMAGE, IMAGESC, IMSHOW, COLORMAP, COLORBAR.

% $Id: sc.m,v 1.84 2009/03/17 10:38:54 ojw Exp $
% Copyright: Oliver Woodford, 2007

%% Check for arguments
if nargin == 0
    % If there are no input arguments then run the demo
    if nargout > 0
        error('Output expected from no inputs!');
    end
    demo; % Run the demo
    return
end

%% Size our image(s)
[y x c n] = size(I);
I = reshape(I, y, x, c, n);

%% Check if image is given with RGB colour along the first dimension
if y == 3 && c > 3
    % Flip colour to 3rd dimension
    I = permute(I, [2 3 1 4]);
    [y x c n] = size(I);
end

%% Don't do much if I is empty
if isempty(I)
    if nargout == 0
        % Clear the current axes if we were supposed to display the image
        cla; axis off;
    else
        % Create an empty array with the correct dimensions
        I = zeros(y, x, (c~=0)*3, n);
    end
    return
end

%% Check for multiple images
% If we have a non-singleton 4th dimension we want to display the images in
% a 3x4 grid and use buttons to cycle through them
if n > 1
    if nargout > 0
        % Return transformed images in an YxXx3xN array
        A = zeros(y, x, 3, n);
        for a = 1:n
            A(:,:,:,a) = sc(I(:,:,:,a), varargin{:});
        end
        I = A;
    else
        % Removed functionality
        fprintf([' SC no longer supports the display of multiple images. The\n'...
                 ' functionality has been incorporated into an improved version\n'...
                 ' of MONTAGE, available on the MATLAB File Exchange at:\n'...
                 '    http://www.mathworks.com/matlabcentral/fileexchange/22387\n']);
        clear I;
    end
    return
end

%% Parse the input arguments coming after I (1st input)
[map limits mask] = parse_inputs(I, varargin, y, x);

%% Call the rendering function
I = reshape(double(real(I)), y*x, c); % Only work with real doubles
if ~ischar(map)
    % Table-based colormap
    reverseMap = false;
    [I limits] = interp_map(I, limits, reverseMap, map);
else
    % If map starts with a '-' sign, invert the colourmap
    reverseMap = map(1) == '-';
    map = lower(map(reverseMap+1:end));
    
    % Predefined colormap
    [I limits] = colormap_switch(I, map, limits, reverseMap, c);
end

%% Update any masked pixels
I = reshape(I, y*x, 3);
for a = 1:size(mask, 2)
    I(mask{2,a},1) = mask{1,a}(:,1);
    I(mask{2,a},2) = mask{1,a}(:,2);
    I(mask{2,a},3) = mask{1,a}(:,3);
end
I = reshape(I, [y x 3]); % Reshape to correct size


%% Only display if the output isn't used
if nargout == 0
    display_image(I, map, limits, reverseMap);
    % Don't print out the matrix if we've forgotten the ";"
    clear I
end
return

%% Colormap switch
function [I limits] = colormap_switch(I, map, limits, reverseMap, c)
% Large switch statement for all the colourmaps
switch map
%% Prism
    case 'prism'
        % Similar to the MATLAB internal prism colormap, but only works on
        % index images, assigning each index (or rounded float) to a
        % different colour
        [I limits] = index_im(I);
        % Generate prism colourmap
        map = prism(6);
        if reverseMap
            map = map(end:-1:1,:); % Reverse the map
        end
        % Lookup the colours
        I = mod(I, 6) + 1;
        I = map(I,:);
%% Rand
    case 'rand'
        % Assigns a random colour to each index
        [I limits num_vals] = index_im(I);
        % Generate random colourmap
        map = rand(num_vals, 3);
        % Lookup the colours
        I = map(I,:);
%% Diff
    case 'diff'
        % Show positive as blue and negative as red, white is 0
        switch c
            case 1
                I(:,2:3) = 0;
            case 2
                % Second channel can only have absolute value
                I(:,3) = abs(I(:,2));
            case 3
                % Diff of RGB images - convert to YUV first
                I = rgb2yuv(I);
                I(:,3) = sqrt(sum(I(:,2:end) .^ 2, 2)) ./ sqrt(2);
            otherwise
                % Use difference along principle component, and other
                % channels to modulate second channel
                I = calc_prin_comps(I);
                I(:,3) = sqrt(sum(I(:,2:end) .^ 2, 2)) ./ sqrt(c - 1);
                I(:,4:end) = [];
        end
        % Generate limits
        if isempty(limits)
            limits = [min(I(:,1)) max(I(:,1))];
        end
        limits = max(abs(limits));
        if limits
            % Scale
            if c > 1
                I(:,[1 3]) = I(:,[1 3]) / limits;
            else
                I = I / (limits * 0.5);
            end
        end
        % Colour
        M = I(:,1) > 0;
        I(:,2) = -I(:,1) .* ~M;
        I(:,1) = I(:,1) .* M;
        if reverseMap
            % Swap first two channels
            I = I(:,[2 1 3]);
        end
        %I = 1 - I * [1 0.4 1; 0.4 1 1; 1 1 0.4]; % (Green/Red)
        I = 1 - I * [1 1 0.4; 0.4 1 1; 1 0.4 1]; % (Blue/Red)
        I = min(max(I, 0), 1);
        limits = [-limits limits]; % For colourbar
%% Flow
    case 'flow'
        % Calculate amplitude and phase, and use 'phase'
        if c ~= 2
            error('''flow'' requires two channels');
        end
        A = sqrt(sum(I .^ 2, 2));
        if isempty(limits)
            limits = [min(A) max(A)*2];
        else
            limits = [0 max(abs(limits)*sqrt(2))*2];
        end
        I(:,1) = atan2(I(:,2), I(:,1));
        I(:,2) = A;
        if reverseMap
            % Invert the amplitude
            I(:,2) = -I(:,2);
            limits = -limits([2 1]);
        end
        I = phase_helper(I, limits, 2); % Last parameter tunes how saturated colors can get
        % Set NaNs (unknown flow) to 0
        I(isnan(I)) = reverseMap;
        limits = []; % This colourmap doesn't have a valid colourbar
%% Phase
    case 'phase'
        % Plot amplitude as intensity and angle as hue
        if c < 2
            error('''phase'' requires two channels');
        end
        if isempty(limits)
            limits = [min(I(:,1)) max(I(:,1))];
        end
        if reverseMap
            % Invert the phase
            I(:,2) = -I(:,2);
        end
        I = I(:,[2 1]);
        if diff(limits)
            I = phase_helper(I, limits, 1.3); % Last parameter tunes how saturated colors can get
        else
            % No intensity - just cycle hsv
            I = hsv_helper(mod(I(:,1) / (2 * pi), 1));
        end
        limits = []; % This colourmap doesn't have a valid colourbar
%% RGB2Grey
    case {'rgb2grey', 'rgb2gray'}
        % Compress RGB to greyscale
        [I limits] = rgb2grey(I, limits, reverseMap);
%% RGB2YUV
    case 'rgb2yuv'
        % Convert RGB to YUV - not for displaying or saving to disk!
        [I limits] = rgb2yuv(I);
%% YUV2RGB
    case 'yuv2rgb'
        % Convert YUV to RGB - undo conversion of rgb2yuv
        if c ~= 3
            error('''yuv2rgb'' requires a 3 channel image');
        end
        I = reshape(I, y*x, 3);
        I = I * [1 1 1; 0, -0.39465, 2.03211; 1.13983, -0.58060  0];
        I = reshape(I, y, x, 3);
        I = sc(I, limits);
        limits = []; % This colourmap doesn't have a valid colourbar
%% Prob
    case 'prob'
        % Plot first channel as grey variation of 'bled' and modulate
        % according to other channels
        if c > 1
            A = rgb2grey(I(:,2:end), [], false);
            I = I(:,1);
        else
            A = 0.5;
        end
        [I limits] = bled(I, limits, reverseMap);
        I = normalize(A + I, [-0.1 1.3]);
%% Prob_jet
    case 'prob_jet'
        % Plot first channel as 'jet' and modulate according to other
        % channels
        if c > 1
            A = rgb2grey(I(:,2:end), [], false);
            I = I(:,1);
        else
            A = 0.5;
        end
        [I limits] = jet_helper(I, limits, reverseMap);
        I = normalize(A + I, [0.2 1.8]);
%% Compress
    case 'compress'
        % Compress to RGB, maximizing variance
        % Determine and scale to limits
        I = normalize(I, limits);
        if reverseMap
            % Invert after everything
            I = 1 - I;
        end
        % Zero mean
        meanCol = mean(I, 1);
        isBsx = exist('bsxfun', 'builtin');
        if isBsx
            I = bsxfun(@minus, I, meanCol);
        else