sc.m

SC Display/output truecolor images with a range of colormaps

            I = I - meanCol(ones(x*y, 1, 'uint8'),:);
        end
        % Calculate top 3 principle components
        I = calc_prin_comps(I, 3);
        % Normalize each channel independently
        if isBsx
            I = bsxfun(@minus, I, min(I, [], 1));
            I = bsxfun(@times, I, 1./max(I, [], 1));
        else
            for a = 1:3
                I(:,a) = I(:,a) - min(I(:,a));
                I(:,a) = I(:,a) / max(I(:,a));
            end
        end
        % Put components in order of human eyes' response to channels
        I = I(:,[2 1 3]);
        limits = []; % This colourmap doesn't have a valid colourbar
%% Stereo (anaglyph)
    case 'stereo'
        % Convert 2 colour images to intensity images
        % Show first channel as red and second channel as cyan
        A = rgb2grey(I(:,1:floor(end/2)), limits, false);
        I = rgb2grey(I(:,floor(end/2)+1:end), limits, false);
        if reverseMap
            I(:,2:3) = A(:,1:2); % Make first image cyan
        else
            I(:,1) = A(:,1); % Make first image red
        end
        limits = []; % This colourmap doesn't have a valid colourbar
%% Coloured anaglyph
    case 'stereo_col'
        if c ~= 6
            error('''stereo_col'' requires a 6 channel image');
        end
        I = normalize(I, limits);
        % Red channel from one image, green and blue from the other
        if reverseMap
            I(:,1) = I(:,4); % Make second image red
        else
            I(:,2:3) = I(:,5:6); % Make first image red
        end
        I = I(:,1:3);
        limits = []; % This colourmap doesn't have a valid colourbar
%% None
    case 'none'
        % No colour map - just output the image
        if c ~= 3
            [I limits] = grey(I, limits, reverseMap);
        else
            I = intensity(I(:), limits, reverseMap);
            limits = [];
        end
%% Grey
    case {'gray', 'grey'}
        % Greyscale
        [I limits] = grey(I, limits, reverseMap);
%% Jet
    case 'jet'
        % Dark blue to dark red, through green
        [I limits] = jet_helper(I, limits, reverseMap);
    case 'jet2'
        % Like jet, but starts in black and goes to saturated red
        [I limits] = interp_map(I, limits, reverseMap, [0 0 0; 0.5 0 0.5; 0 0 0.9; 0 1 1; 0 1 0; 1 1 0; 1 0 0]);
%% Hot
    case 'hot'
        % Black to white through red and yellow
        [I limits] = interp_map(I, limits, reverseMap, [0 0 0 3; 1 0 0 3; 1 1 0 2; 1 1 1 0]);
    case 'hot2'
        % Like hot, but equally spaced
        [I limits] = intensity(I, limits, reverseMap); % Intensity map
        I = I * 3;
        I = [I, I-1, I-2];
        I = min(max(I, 0), 1); % Truncate
    case {'hotter', 'hot*'}
        % Converts to linear greyscale
        [I limits] = interp_map(I, limits, reverseMap, [0 0 0 299; 1 0 0 587; 1 1 0 114; 1 1 1 0]);
%% Thermal
    case 'thermal'
        % Black, purple, red, orange, yellow - typical for thermal imaging
        [I limits] = interp_map(I, limits, reverseMap, [0 0 0; 0.3 0 0.7; 1 0.2 0; 1 1 0; 1 1 1]);
    case 'thermal*'
        % Converts to linear greyscale
        [I limits] = interp_map(I, limits, reverseMap, [0 0 0 1695; 0.3 0 0.7 2469; 1 0.2 0 4696; 1 1 0 1140; 1 1 1 0]);
%% Contrast
    case 'contrast'
        % A high contrast, full-colour map that goes from black to white
        % linearly when converted to greyscale, and passes through all the
        % corners of the RGB colour cube
        [I limits] = interp_map(I, limits, reverseMap, [0 0 0 114; 0 0 1 185; 1 0 0 114; 1 0 1 174;...
                                                        0 1 0 114; 0 1 1 185; 1 1 0 114; 1 1 1 0]);
%% HSV
    case 'hsv'
        % Cycle through hues
        [I limits] = intensity(I, limits, reverseMap); % Intensity map
        I = hsv_helper(I);
%% Bone
    case 'bone'
        % Greyscale with a blue tint
        [I limits] = interp_map(I, limits, reverseMap, [0 0 0 3; 21 21 29 3; 42 50 50 2; 64 64 64 1]/64);
    case 'bone2'
        % Like bone, but equally spaced
        [I limits] = intensity(I, limits, reverseMap); % Intensity map
        J = [I-2/3, I-1/3, I];
        J = max(min(J, 1/3), 0) * (2 / 5);
        I = I * (13 / 15);
        I = J + I(:,[1 1 1]);
%% Colourcube
    case {'colorcube', 'colourcube'}
        % Psychedelic colourmap inspired by MATLAB's version
        [I limits] = intensity(I, limits, reverseMap); % Intensity map
        step = 4;
        I = I * (step * (1 - eps));
        J = I * step;
        K = floor(J);
        I = cat(3, mod(K, step)/(step-1), J - floor(K), mod(floor(I), step)/(step-1));
%% Cool
    case 'cool'
        % Cyan through to magenta
        [I limits] = intensity(I, limits, reverseMap); % Intensity map
        I = [I, 1-I, ones(size(I))];
%% Spring
    case 'spring'
        % Magenta through to yellow
        [I limits] = intensity(I, limits, reverseMap); % Intensity map
        I = [ones(size(I)), I, 1-I];
%% Summer
    case 'summer'
        % Darkish green through to pale yellow
        [I limits] = intensity(I, limits, reverseMap); % Intensity map
        I = [I, 0.5+I*0.5, 0.4*ones(size(I))];
%% Autumn
    case 'autumn'
        % Red through to yellow
        [I limits] = intensity(I, limits, reverseMap); % Intensity map
        I = [ones(size(I)), I, zeros(size(I))];
%% Winter
    case 'winter'
        % Blue through to turquoise
        [I limits] = intensity(I, limits, reverseMap); % Intensity map
        I = [zeros(size(I)), I, 1-I*0.5];
%% Copper
    case 'copper'
        % Black through to copper
        [I limits] = intensity(I, limits, reverseMap); % Intensity map
        I = [I*(1/0.8), I*0.78, I*0.5];
        I = min(max(I, 0), 1); % Truncate
    case {'copper2', 'copper*'}
        % Converts to greyscale
        [I limits] = interp_map(I, limits, reverseMap, [0 0 0; 0.2651 0.2426 0.2485; 0.666 0.4399 0.3738; 0.8118 0.7590 0.5417; 1 1 1]);
%% Pink
    case 'pink'
        % Greyscale with a pink tint
        [I limits] = intensity(I, limits, reverseMap); % Intensity map
        J = I * (2 / 3);
        I = [I, I-1/3, I-2/3];
        I = max(min(I, 1/3), 0);
        I = I + J(:,[1 1 1]);
        I = sqrt(I);
%% Sepia
    case 'sepia'
        % Greyscale with a brown (sepia) tint
        [I limits] = interp_map(I, limits, reverseMap, [0 0 0 5; 0.1 0.05 0 85; 1 0.9 0.8 10; 1 1 1 0]);
%% Bled
    case 'bled'
        % Black to red, through blue
        [I limits] = bled(I, limits, reverseMap);
%% Earth
    case 'earth'
        % High contrast, converts to linear scale in grey, strong
        % shades of green
        table = [0 0 0; 0 0.1104 0.0583; 0.1661 0.1540 0.0248; 0.1085 0.2848 0.1286;...
            0.2643 0.3339 0.0939; 0.2653 0.4381 0.1808; 0.3178 0.5053 0.3239;...
            0.4858 0.5380 0.3413; 0.6005 0.5748 0.4776; 0.5698 0.6803 0.6415;...
            0.5639 0.7929 0.7040; 0.6700 0.8626 0.6931; 0.8552 0.8967 0.6585;...
            1 0.9210 0.7803; 1 1 1];
        [I limits] = interp_map(I, limits, reverseMap, table);
%% Pinker
    case 'pinker'
        % High contrast, converts to linear scale in grey, strong
        % shades of pink
        table = [0 0 0; 0.0455 0.0635 0.1801; 0.2425 0.0873 0.1677;...
            0.2089 0.2092 0.2546; 0.3111 0.2841 0.2274; 0.4785 0.3137 0.2624;...
            0.5781 0.3580 0.3997; 0.5778 0.4510 0.5483; 0.5650 0.5682 0.6047;...
            0.6803 0.6375 0.5722; 0.8454 0.6725 0.5855; 0.9801 0.7032 0.7007;...
            1 0.7777 0.8915; 0.9645 0.8964 1; 1 1 1];
        [I limits] = interp_map(I, limits, reverseMap, table);
%% Pastel
    case 'pastel'
        % High contrast, converts to linear scale in grey, strong
        % pastel shades
        table = [0 0 0; 0.4709 0 0.018; 0 0.3557 0.6747; 0.8422 0.1356 0.8525;
            0.4688 0.6753 0.3057; 1 0.6893 0.0934; 0.9035 1 0; 1 1 1];
        [I limits] = interp_map(I, limits, reverseMap, table);
%% Bright
    case 'bright'
        % High contrast, converts to linear scale in grey, strong
        % saturated shades
        table = [0 0 0; 0.3071 0.0107 0.3925; 0.007 0.289 1; 1 0.0832 0.7084;
            1 0.4447 0.1001; 0.5776 0.8360 0.4458; 0.9035 1 0; 1 1 1];
        [I limits] = interp_map(I, limits, reverseMap, table);
%% Unknown colourmap
    otherwise
        error('Colormap ''%s'' not recognised.', map);
end
return


%% Display image
function display_image(I, map, limits, reverseMap)
% Clear the axes
cla(gca, 'reset');
% Display the image - using image() is fast
hIm = image(I);
% Get handles to the figure and axes (now, as the axes may have
% changed)
hFig = gcf; hAx = gca;
% Axes invisible and equal
set(hFig, 'Units', 'pixels');
set(hAx, 'Visible', 'off', 'DataAspectRatio', [1 1 1], 'DrawMode', 'fast');
% Make title and axis labels visible
set(get(hAx, 'XLabel'), 'Visible', 'on');
set(get(hAx, 'YLabel'), 'Visible', 'on');
set(get(hAx, 'Title'), 'Visible', 'on');
% Set data for a colorbar
if ~isempty(limits) && limits(1) ~= limits(2)
    colBar = (0:255) * ((limits(2) - limits(1)) / 255) + limits(1);
    colBar = squeeze(sc(colBar, map, limits));
    if reverseMap
        colBar = colBar(end:-1:1,:);
    end
    set(hFig, 'Colormap', colBar);
    set(hAx, 'CLim', limits);
    set(hIm, 'CDataMapping', 'scaled');
end
% Only resize image if it is alone in the figure
if numel(findobj(get(hFig, 'Children'), 'Type', 'axes')) > 1
    return
end
% Could still be the first subplot - do another check
axesPos = get(hAx, 'Position');
if isequal(axesPos, get(hFig, 'DefaultAxesPosition'))
    % Default position => not a subplot
    % Fill the window
    set(hAx, 'Units', 'normalized', 'Position', [0 0 1 1]);
    axesPos = [0 0 1 1];
end
if ~isequal(axesPos, [0 0 1 1]) || strcmp(get(hFig, 'WindowStyle'), 'docked')
    % Figure not alone, or docked. Either way, don't resize.
    return
end
% Get the size of the monitor we're on
figPosCur = get(hFig, 'Position');
MonSz = get(0, 'MonitorPositions');
MonOn = size(MonSz, 1);
if MonOn > 1
    figCenter = figPosCur(1:2) + figPosCur(3:4) / 2;
    figCenter = MonSz - repmat(figCenter, [MonOn 2]);
    MonOn = all(sign(figCenter) == repmat([-1 -1 1 1], [MonOn 1]), 2);
    MonOn(1) = MonOn(1) | ~any(MonOn);
    MonSz = MonSz(MonOn,:);
end
MonSz(3:4) = MonSz(3:4) - MonSz(1:2) + 1;
% Check if the window is maximized
% This is a hack which may only work on Windows! No matter, though.
if isequal(MonSz([1 3]), figPosCur([1 3]))
    % Leave maximized
    return
end
% Compute the size to set the window
MaxSz = MonSz(3:4) - [20 120];
ImSz = [size(I, 2) size(I, 1)];
RescaleFactor = min(MaxSz ./ ImSz);
if RescaleFactor > 1
    % Integer scale for enlarging, but don't make too big
    MaxSz = min(MaxSz, [1000 680]);
    RescaleFactor = max(floor(min(MaxSz ./ ImSz)), 1);
end
figPosNew = ceil(ImSz * RescaleFactor);
% Don't move the figure if the size isn't changing
if isequal(figPosCur(3:4), figPosNew)
    return
end
% Keep the centre of the figure stationary
figPosNew = [max(1, floor(figPosCur(1:2)+(figPosCur(3:4)-figPosNew)/2)) figPosNew];
% Ensure the figure bar is in bounds
figPosNew(1:2) = min(figPosNew(1:2), MonSz(1:2)+MonSz(3:4)-[6 101]-figPosNew(3:4));
set(hFig, 'Position', figPosNew);
return

%% Parse input variables
function [map limits mask] = parse_inputs(I, inputs, y, x)

% Check the first two arguments for the colormap and limits
ninputs = numel(inputs);
map = 'none';
limits = [];
mask = 1;
for a = 1:min(2, ninputs)
    if ischar(inputs{a}) && numel(inputs{a}) > 1
        % Name of colormap
        map = inputs{a};
    elseif isnumeric(inputs{a})
        [p q r] = size(inputs{a});
        if (p * q * r) == 2
            % Limits
            limits = double(inputs{a});
        elseif p > 1 && (q == 3 || q == 4) && r == 1
            % Table-based colormap
            map = inputs{a};
        else
            break;
        end
    else
        break;
    end
    mask = mask + 1;
end
% Check for following inputs
if mask > ninputs
    mask = cell(2, 0);
    return
end
% Following inputs must either be colour/mask pairs, or a colour for NaNs
if ninputs - mask == 0
    mask = cell(2, 1);
    mask{1} = inputs{end};
    mask{2} = ~all(isfinite(I), 3);
elseif mod(ninputs-mask, 2) == 1
    mask = reshape(inputs(mask:end), 2, []);
else
    error('Error parsing inputs');
end
% Go through pairs and generate
for a = 1:size(mask, 2)
    % Generate any masks from functions
    if isa(mask{2,a}, 'function_handle')
        mask{2,a} = mask{2,a}(I);
    end
    if ~islogical(mask{2,a})
        error('Mask is not a logical array');
    end
    if ~isequal(size(mask{2,a}), [y x])
        error('Mask does not match image size');
    end
    if ischar(mask{1,a})
        if numel(mask{1,a}) == 1
            % Generate colours from MATLAB colour strings
            mask{1,a} = rem(floor((strfind('kbgcrmyw', mask{1,a}) - 1) * [0.25 0.5 1]), 2);
        else
            % Assume it's a colormap name
            mask{1,a} = sc(I, mask{1,a});
        end
    end
    mask{1,a} = reshape(mask{1,a}, [], 3);
    if size(mask{1,a}, 1) ~= y*x && size(mask{1,a}, 1) ~= 1
        error('Replacement color/image of unexpected dimensions');
    end
    if size(mask{1,a}, 1) ~= 1
        mask{1,a} = mask{1,a}(mask{2,a},:);
    end