sc.m

SC Display/output truecolor images with a range of colormaps

end
return

%% Grey
function [I limits] = grey(I, limits, reverseMap)
% Greyscale
[I limits] = intensity(I, limits, reverseMap);
I = I(:,[1 1 1]);
return

%% RGB2grey
function [I limits] = rgb2grey(I, limits, reverseMap)
% Compress RGB to greyscale
if size(I, 2) == 3
    I = I * [0.299; 0.587; 0.114];
end
[I limits] = grey(I, limits, reverseMap);
return
        
%% RGB2YUV
function [I limits] = rgb2yuv(I)
% Convert RGB to YUV - not for displaying or saving to disk!
if size(I, 2) ~= 3
    error('rgb2yuv requires a 3 channel image');
end
I = I * [0.299, -0.14713, 0.615; 0.587, -0.28886, -0.51498; 0.114, 0.436, -0.10001];
limits = []; % This colourmap doesn't have a valid colourbar
return

%% Phase helper
function I = phase_helper(I, limits, n)
I(:,1) = mod(I(:,1)/(2*pi), 1);
I(:,2) = I(:,2) - limits(1);
I(:,2) = I(:,2) * (n / (limits(2) - limits(1)));
I(:,3) = n - I(:,2);
I(:,[2 3]) = min(max(I(:,[2 3]), 0), 1);
I = hsv2rgb(reshape(I, [], 1, 3));
return

%% Jet helper        
function [I limits] = jet_helper(I, limits, reverseMap)
% Dark blue to dark red, through green
[I limits] = intensity(I, limits, reverseMap);
I = I * 4;
I = [I-3, I-2, I-1];
I = 1.5 - abs(I);
I = min(max(I, 0), 1);
return

%% HSV helper
function I = hsv_helper(I)
I = I * 6;
I = abs([I-3, I-2, I-4]);
I(:,1) = I(:,1) - 1;
I(:,2:3) = 2 - I(:,2:3);
I = min(max(I, 0), 1);
return

%% Bled
function [I limits] = bled(I, limits, reverseMap)
% Black to red through blue
[I limits] = intensity(I, limits, reverseMap);
J = reshape(hsv_helper(I), [], 3);
if exist('bsxfun', 'builtin') 
    I = bsxfun(@times, I, J);
else
    I = J .* I(:,[1 1 1]);
end
return

%% Normalize
function [I limits] = normalize(I, limits)
if isempty(limits)
    limits = isfinite(I);
    if ~any(reshape(limits, numel(limits), 1))
        % All NaNs, Infs or -Infs
        I = double(I > 0);
        limits = [0 1];
        return
    end
    limits = [min(I(limits)) max(I(limits))];
    I = I - limits(1);
    if limits(2) ~= limits(1)
        I = I * (1 / (limits(2) - limits(1)));
    end
else
    I = I - limits(1);
    if limits(2) ~= limits(1)
        I = I * (1 / (limits(2) - limits(1)));
    end
    I = min(max(I, 0), 1);
end
return

%% Intensity maps
function [I limits] = intensity(I, limits, reverseMap)
% Squash to 1d using L2 norm
if size(I, 2) > 1
    I = sqrt(sum(I .^ 2, 2));
end
% Determine and scale to limits
[I limits] = normalize(I, limits);
if reverseMap
    % Invert after everything
    I = 1 - I;
end
return

%% Interpolate table-based map
function [I limits] = interp_map(I, limits, reverseMap, map)
% Convert to intensity
[I limits] = intensity(I, limits, reverseMap);
% Compute indices and offsets
if size(map, 2) == 4
    bins = map(1:end-1,4);
    cbins = cumsum(bins);
    bins = bins ./ cbins(end);
    cbins = cbins(1:end-1) ./ cbins(end);
    if exist('bsxfun', 'builtin')
        ind = bsxfun(@gt, I(:)', cbins(:));
    else
        ind = repmat(I(:)', [numel(cbins) 1]) > repmat(cbins(:), [1 numel(I)]);
    end
    ind = min(sum(ind), size(map, 1) - 2) + 1;
    bins = 1 ./ bins;
    cbins = [0; cbins];
    I = (I - cbins(ind)) .* bins(ind);
else
    n = size(map, 1) - 1;
    I = I(:) * n;
    ind = min(floor(I), n-1);
    I = I - ind;
    ind = ind + 1;
end
if exist('bsxfun', 'builtin')
    I = bsxfun(@times, map(ind,1:3), 1-I) + bsxfun(@times, map(ind+1,1:3), I);
else
    I = map(ind,1:3) .* repmat(1-I, [1 3]) + map(ind+1,1:3) .* repmat(I, [1 3]);
end
I = min(max(I, 0), 1); % Rounding errors can make values slip outside bounds 
return

%% Index images
function [J limits num_vals] = index_im(I)
% Returns an index image
if size(I, 2) ~= 1
    error('Index maps only work on single channel images');
end
J = round(I);
rescaled = any(abs(I - J) > 0.01);
if rescaled
    % Appears not to be an index image. Rescale over 256 indices
    m = min(I);
    m = m * (1 - sign(m) * eps);
    I = I - m;
    I = I * (256 / max(I(:)));
    J = ceil(I);
    num_vals = 256;
elseif nargout > 2
    % Output the number of values
    J = J - (min(J) - 1);
    num_vals = max(J);
end
% These colourmaps don't have valid colourbars
limits = [];    
return

%% Calculate principle components
function I = calc_prin_comps(I, numComps)
if nargin < 2
    numComps = size(I, 2);
end
% Do SVD
[I S] = svd(I, 0);
% Calculate projection of data onto components
S = diag(S(1:numComps,1:numComps))';
if exist('bsxfun', 'builtin')
    I = bsxfun(@times, I(:,1:numComps), S);
else
    I = I(:,1:numComps) .* S(ones(size(I, 1), 1, 'uint8'),:);
end
return

%% Demo function to show capabilities of sc
function demo
%% Demo gray & lack of border
figure; fig = gcf; Z = peaks(256); sc(Z);
display_text([...
' Lets take a standard, MATLAB, real-valued function:\n\n    peaks(256)\n\n'...
' Calling:\n\n    figure\n    Z = peaks(256);\n    sc(Z)\n\n'...
' gives (see figure). SC automatically scales intensity to fill the\n'...
' truecolor range of [0 1].\n\n'...
' If your figure isn''t docked, then the image will have no border, and\n'...
' will be magnified by an integer factor (in this case, 2) so that the\n'...
' image is a reasonable size.']);

%% Demo colour image display 
figure(fig); clf;
load mandrill; mandrill = ind2rgb(X, map); sc(mandrill);
display_text([...
' That wasn''t so interesting. The default colormap is ''none'', which\n'...
' produces RGB images given a 3-channel input image, otherwise it produces\n'...
' a grayscale image. So calling:\n\n    load mandrill\n'...
'    mandrill = ind2rgb(X, map);\n    sc(mandrill)\n\n gives (see figure).']);

%% Demo discretization
figure(fig); clf;
subplot(121); sc(Z, 'jet'); label(Z, 'sc(Z, ''jet'')');
subplot(122); imagesc(Z); axis image off; colormap(jet(64)); % Fix the fact we change the default depth
label(Z, 'imagesc(Z); axis image off; colormap(''jet'');');
display_text([...
' However, if we want to display intensity images in color we can use any\n'...
' of the MATLAB colormaps implemented (most of them) to give truecolor\n'...
' images. For example, to use ''jet'' simply call:\n\n'...
'    sc(Z, ''jet'')\n\n'...
' The MATLAB alternative, shown on the right, is:\n\n'...
'    imagesc(Z)\n    axis equal off\n    colormap(jet)\n\n'...
' which generates noticeable discretization artifacts.']);

%% Demo intensity colourmaps
figure(fig); clf;
subplot(221); sc(Z, 'hsv'); label(Z, 'sc(Z, ''hsv'')');
subplot(222); sc(Z, 'colorcube'); label(Z, 'sc(Z, ''colorcube'')');
subplot(223); sc(Z, 'contrast'); label(Z, 'sc(Z, ''contrast'')');
subplot(224); sc(Z-round(Z), 'diff'); label(Z, 'sc(Z-round(Z), ''diff'')');
display_text([...
' There are several other intensity colormaps to choose from. Calling:\n\n'...
'    help sc\n\n'...
' will give you a list of them. Here are several others demonstrated.']);

%% Demo saturation limits & colourmap reversal
figure(fig); clf;
subplot(121); sc(Z, [0 max(Z(:))], '-hot'); label(Z, 'sc(Z, [0 max(Z(:))], ''-hot'')');
subplot(122); sc(mandrill, [-0.5 0.5]); label(mandrill, 'sc(mandrill, [-0.5 0.5])');
display_text([...
' SC can also rescale intensity, given an upper and lower bound provided\n'...
' by the user, and invert most colormaps simply by prefixing a ''-'' to the\n'...
' colormap name. For example:\n\n'...
'    sc(Z, [0 max(Z(:))], ''-hot'');\n'...
'    sc(mandrill, [-0.5 0.5]);\n\n'...
' Note that the order of the colormap and limit arguments are\n'...
' interchangable.']);

%% Demo prob
load gatlin;
gatlin = X;
figure(fig); clf; im = cat(3, abs(Z)', gatlin(1:256,end-255:end)); sc(im, 'prob');
label(im, 'sc(cat(3, prob, gatlin), ''prob'')');
display_text([...
' SC outputs the recolored data as a truecolor RGB image. This makes it\n'...
' easy to combine colormaps, either arithmetically, or by masking regions.\n'...
' For example, we could combine an image and a probability map\n'...
' arithmetically as follows:\n\n'...
'    load gatlin\n'...
'    gatlin = X(1:256,end-255:end);\n'...
'    prob = abs(Z)'';\n'...
'    im = sc(prob, ''hsv'') .* sc(prob, ''gray'') + sc(gatlin, ''rgb2gray'');\n'...
'    sc(im, [-0.1 1.3]);\n\n'...
' In fact, that particular colormap has already been implemented in SC.\n'...
' Simply call:\n\n'...
'    sc(cat(3, prob, gatlin), ''prob'');']);

%% Demo colorbar
colorbar;
display_text([...
' SC also makes possible the generation of a colorbar in the normal way, \n'...
' with all the colours and data values correct. Simply call:\n\n'...
'    colorbar\n\n'...
' The colorbar doesn''t work with all colormaps, but when it does,\n'...
' inverting the colormap (using ''-map'') maintains the integrity of the\n'...
' colorbar (i.e. it works correctly) - unlike if you invert the input data.']);

%% Demo combine by masking
figure(fig); clf;
sc(Z, [0 max(Z(:))], '-hot', sc(Z-round(Z), 'diff'), Z < 0);
display_text([...
' It''s just as easy to combine generated images by masking too. Here''s an\n'...
' example:\n\n'...
'    im = cat(4, sc(Z, [0 max(Z(:))], ''-hot''), sc(Z-round(Z), ''diff''));\n'...
'    mask = repmat(Z < 0, [1 1 3]);\n'...
'    mask = cat(4, mask, ~mask);\n'...
'    im = sum(im .* mask, 4);\n'...
'    sc(im)\n\n'...
' In fact, SC can also do this for you, by adding image/colormap and mask\n'...
' pairs to the end of the argument list, as follows:\n\n'...
'    sc(Z, [0 max(Z(:))], ''-hot'', sc(Z-round(Z), ''diff''), Z < 0);\n\n'...
' A benefit of the latter approach is that you can still display a\n'...
' colorbar for the first colormap.']);

%% Demo texture map
figure(fig); clf;
surf(Z, sc(Z, 'contrast'), 'edgecolor', 'none');
display_text([...
' Other benefits of SC outputting the image as an array are that the image\n'...
' can be saved straight to disk using imwrite() (if you have the image\n'...
' processing toolbox), or can be used to texture map a surface, thus:\n\n'...
'    tex = sc(Z, ''contrast'');\n'...
'    surf(Z, tex, ''edgecolor'', ''none'');']);

%% Demo compress
load mri;
mri = D;
close(fig); % Only way to get round loss of focus (bug?)
figure(fig); clf;
sc(squeeze(mri(:,:,:,1:6)), 'compress');
display_text([...
' For images with more than 3 channels, SC can compress these images to RGB\n'...
' while maintaining the maximum amount of variance in the data. For\n'...
' example, this 6 channel image:\n\n'...
'    load mri\n    mri = D;\n    sc(squeeze(mri(:,:,:,1:6), ''compress'')']);

%% Demo multiple images
figure(fig); clf; im = sc(mri, 'bone');
for a = 1:12
    subplot(3, 4, a);
    sc(im(:,:,:,a));
end
display_text([...
' SC can process multiple images for export when passed in as a 4d array.\n'...
' For example:\n\n'...
'    im = sc(mri, ''bone'')\n'...
'    for a = 1:12\n'...
'       subplot(3, 4, a);\n'...
'       sc(im(:,:,:,a));\n'...
'    end']);

%% Demo user defined colormap
figure(fig); clf; sc(abs(Z), rand(10, 3)); colorbar;
display_text([...
' Finally, SC can use user defined colormaps to display indexed images.\n'...
' These can be defined as a linear colormap. For example:\n\n'...
'    sc(abs(Z), rand(10, 3))\n    colorbar;\n\n'...
' Note that the colormap is automatically linearly interpolated.']);

%% Demo non-linear user defined colormap
figure(fig); clf; sc(abs(Z), [rand(10, 3) exp((1:10)/2)']); colorbar;
display_text([...
' Non-linear colormaps can also be defined by the user, by including the\n'...
' relative distance between the given colormap points on the colormap\n'...
' scale in the fourth column of the colormap matrix. For example:\n\n'...
'    sc(abs(Z), [rand(10, 3) exp((1:10)/2)''])\n    colorbar;\n\n'...
' Note that the colormap is still linearly interpolated between points.']);

clc; fprintf('End of demo.\n');
return

%% Some helper functions for the demo
function display_text(str)
clc;
fprintf([str '\n\n']);
fprintf('Press a key to go on.\n');
figure(gcf);
waitforbuttonpress;
return

function label(im, str)
text(size(im, 2)/2, size(im, 1)+12, str,...
    'Interpreter', 'none', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle');
return