Code covered by the BSD License  

Highlights from
simulate motion in Cartesian MRI

image thumbnail

simulate motion in Cartesian MRI

by

 

24 May 2013 (Updated )

Simulates Point Spread Function (PSF) and k-space data for motion in Cartesian sequences.

simCartesianMRI( varargin )
function [PSF, ksp] = simCartesianMRI( varargin )
% simCartesianMRI simulate Point Spread Function (PSF) for Cartesian FLASH data with modelled motion
%
% [PSF, ksp] = sim_motion_cartesian( 'param1', value1, 'param2', value2, ... )
%
% All input parameters are optional. Reasonable defaults are assumed if undefined.
%
% Paramters
% =========
%       TR          repetition time of the simulated FLASH sequence
%                   Unit is seconds!
%                   Default: 5e-3
%
%       matrix      matrix size of the k-space grid:
%                    *  matrix(1) == x-size or read points
%                    *  matrix(2) == y-size or phase enc. steps
%                    *  matrix(3) == z-size or partitions (optional)
%                   Default: [512 512 1]
%
%       mask        Vector of bins at which data are acquired
%                    *  [ 1 1 1 1 1 1 1 1 1 ] --> 10 samples at  (1...10) x TR
%                    *  true(256,1)  --> 256 samples at (1...256) x TR
%                    *  logical(mod(1:256,3)) --> skip every 3rd line
%                   Default: true( matrix(2)*matrix(3), 1 )
%
%       os          Read oversampling factor
%                   Default: 1
%
%   Additionally, the parameters of the motion model can be specified, for instance
%   simCartesianMRI( ..., 'sinXampl',  8.5, 'sinXfreq', 5 )
%   will simulate sinusoidal motion in x-direction with 8.5 Px amplitude and 5 Hz frequency.
%   All available motion parameters can be found by simply calling getMotion() without argument.
%
% Output
% ======
%       PSF         resulting PSFs for all sampling patterns
%                       size( PSF ) ==  Ny   Nread   Nz   Nsamplings
%
%       ksp         k-spaces including all phase "errors" and averages
%
%
% == NOTE ==
%
%       Further motion models can be incorporated quite easily
%       --> look inside getMotion.m
%
%       Reordering methods are (to be) implemented in getReorderLinesCartesian.m
%
%
% michael.voelker@mr-bavaria.de, 2013

    % which sampling patterns to simulate
    % (to be) implemented in getReorderLinesCartesian.m
    samplingType = { 'linear'; 'shuffle'; 'random' };

    % default values (SI units!):
    defMask = ':';              % acquire one line per matrix line
    defTR   = 5e-3;             % 5ms
    defSize = [512 512 1];      % matrix size
    defOs   = 1;                % read oversampling factor

    % parse input as key-value pairs
    p = inputParser;
    p.KeepUnmatched = true;     % unmatched input is interpreted as a type of motion and passed to getMotion()

    % addParamValue( 'varName'    , default, checkFunction   );
    p.addParamValue( 'mask'       , defMask, @checkMask );
    p.addParamValue( 'TR'         , defTR  , @(t) validateattributes( t, {'numeric'}, {'nonempty', 'scalar', 'finite', 'real', 'positive'} ) );
    p.addParamValue( 'matrix'     , defSize, @(m) validateattributes( m, {'numeric'}, {'nonempty', 'vector', 'finite', 'positive', 'integer' } ) );
    p.addParamValue( 'os'         , defOs  , @(n) validateattributes( m, {'numeric'}, {'nonempty', 'scalar', 'finite', 'positive', 'integer' } ) );
    p.parse( varargin{:} )

    mask    = p.Results.mask;
    TR      = p.Results.TR;
    os      = p.Results.os;

    switch numel( p.Results.matrix )
        case 2                          % 2D
            szZ = 1;
        case 3                          % 3D
            szZ = p.Results.matrix(3);
        otherwise
            error( 'simCartesianMRI:MatrixSize', '''matrix'' must be a vector with 2 or 3 elements.')
    end
    szRead  = p.Results.matrix(1);
    szY     = p.Results.matrix(2); 

    % Simple parser to extract the desired motion model
    toRemove = [];
    for n = 1:length( varargin )
        if strcmpi( varargin{n}, 'mask' ) || strcmpi(varargin{n}, 'TR') || strcmpi(varargin{n}, 'matrix') || strcmpi(varargin{n}, 'os')

            toRemove = [ toRemove  n   ];
            if isnumeric(varargin{n+1}) || islogical( varargin{n+1} ) || isequal( varargin{n+1}, ':' )
                toRemove = [ toRemove  n+1 ];
            else
                error( 'simCartesianMRI:ValMissing', 'Missing value after "%s"', varargin{n} )
            end
        end
    end

    motionType = varargin;
    motionType(toRemove) = [];

    if isequal( mask, ':' )
        Nexc = szY * szZ;
        NUsedLines = Nexc;
    else
        mask = logical( mask );
        Nexc = numel( mask );
        NUsedLines = nnz( mask );       % [N]o. of [n]on-[z]ero elements
    end


    % prettyprint parameters
    %
    % --> http://www.mathworks.com/matlabcentral/fileexchange/33717
    %
    fprintf('\n______________________________________________________________________________________________\n')
    displaytable( [ szRead      ,   szY     ,  szZ      ,  1e3*TR   ,      Nexc      ,  NUsedLines   ,   os         ] ,     ...
                {'read points'  ,  'y-size' , 'z-size'  ,  'TR/ms'  , '#excitations' , '#kept Lines' , 'Oversampling'} )
    fprintf('______________________________________________________________________________________________\n\n')


    Nlines  = szY * szZ;
    Nread   = szRead * os;

    kR      = ( (1:Nread).'  -  floor(Nread/2 + 1) ) ./ Nread;      % column vector with k-space read coordinates

    % all times points of the sequence (including the ones that we ignore with the mask)
    t  = (1:Nexc).' .* TR;

    % generate motion
    pos = getMotion( t, motionType{:} );       % size(pos) = [ length(t)   3 ]

    % Handle the mask
    if NUsedLines ~= Nexc
        t   = t( mask );        % time points to keep
        pos = pos( mask, : );   % relevant positions
    end

    NumSamplings = size( samplingType, 1);

    % preallocation
    ksp = complex( zeros( szRead, szY, szZ, NumSamplings, 'single' ) );

    % process the sampling types
    tic
    for Nsmpl = 1:NumSamplings      % parfor does not seem to be helpful.
                                    % The code is vectorized anyway.

        currSampling = samplingType{Nsmpl,1};

        [ idxY, idxZ ] = getReorderLinesCartesian( 'reorderString', currSampling, 'szY', szY, 'szZ', szZ, 'Nexc', Nexc );

        if NUsedLines ~= Nexc
            idxY = idxY( mask );      % indices to keep
            idxZ = idxZ( mask );
        end

        kY = (idxY.' - floor(szY/2 + 1) ) ./  szY;      %  row vectors with k-space "trajectory"
        kZ = (idxZ.' - floor(szZ/2 + 1) ) ./  szZ;      %


        % ( ===================================================================
        % Apply Fourier-Shift-Theorem for the 3 axes:
        %
        %   A displacement is assumed to be constant during readout.
        %
            data = bsxfun( @plus,                   ...
                            kR * pos(:,1).',        ...     % Read (freq. encoding)
                            (   kY .*  pos(:,2).'   ...     % Y (1st phase encoding direction)
                              + kZ .*  pos(:,3).'   ...     % Z (2nd phase encoding direction)
                            )                       ...
                          );

            data = exp( -2i * pi * data );      %  Nread  x  NUsedLines
        %
        % ) ===================================================================

        clear  kY  kZ

        % eventually add some noise...
        % sigma = 0.1 * sqrt( szRead * Nlines );
        % data = data  +  sigma * complex( randn(size(data),class(data)),  randn(size(data),class(data)) );


        % ( ===================================================================
        % reordering and averaging if necessary
        %
        N = 1;
        linIdx = sub2ind( [szY  szZ], idxY, idxZ );     % 1D index for the lines (2 dimensions "linearized")

        if isequal( sort(linIdx), unique(linIdx) )  % determine if indices are unique (no repetition)

            toKeep = any( data, 1 );
            notIdx = setdiff( 1:Nlines, linIdx );   % the lines that should *not* be assigned
            data(:,linIdx) = data(:,toKeep);
            data(:,notIdx) = 0;

            notSampled = sum(data(:) == 0) ./ Nread;
            overSampled = 0;

        else    % Averaging is necessary

            nPtsMax = 20e6;
            
            if Nread * NUsedLines < nPtsMax
                % Find index for each and every data sample point (linear index for 3D)
                linIdx = bsxfun( @plus, (1:Nread).', Nread * (linIdx(:) - 1).' );

                data = accumarray( linIdx(:),   data(:), [Nread*Nlines   1], @sum );    % Move data to their intended postions on a Cartesian grid and eventually sum their signals.
                N    = accumarray( linIdx(:), single(1), [Nread*Nlines   1], @sum );    % How many source data points for each position on the grid?

            else
                % Fallback for very large problems (to avoid out-of-memory)
                %
                numPack = round(Nread * NUsedLines / nPtsMax);      % divide the labour in that many packages
                p = 1:ceil(NUsedLines/numPack);
                readPts = (1:Nread).';
                linIdx = linIdx(:).';

                tmp = single(0);
                N   = single(0);

                for n = 1:numPack

                    if n == numPack
                        p( p > numel(linIdx) ) = [];
                        if numel(p) < 1
                            break
                        end
                    end

                    idx = bsxfun( @plus, readPts, Nread * (linIdx(p) - 1) );
                    linIdx(p) = [];

                    tmp = tmp + accumarray( idx(:), makeCol(data(:,p)), [Nread*Nlines   1], @sum );
                    N   = N   + accumarray( idx(:),     single(1)     , [Nread*Nlines   1], @sum );
                end
                data = tmp;
                clear  tmp
            end

            data = data ./ max(N, 1);                   % average

            overSampled = nnz(N) ./ Nread;
            notSampled  = Nlines - overSampled;
        end

        data = reshape( data, Nread, szY, szZ );    % make it 3D again
        %
        % ) ===================================================================

        clear  linIdx  notIdx  toKeep  idxY  idxZ

        % print some useful info for non-uniform sampling patterns
        %
        if notSampled ~= 0  ||  max(N) ~= min(N)
            fprintf(['In sampling type "%s":\n'    ...
                     '\t %g\tlines were not acquired (%0.3g%% of all lines)\n'  ...
                     '\t %g\tlines were sampled more than once (%0.3g%% of all lines)\n'    ...
                     '\t %g\t\twas the maximum No. of repeated samplings.\n\n'],    ...
                            currSampling, notSampled, 100*notSampled/Nlines, overSampled, 100*overSampled/Nlines , max(N) )
        end

        clear  N

        % remove read oversampling:
        %
        if os > 1
            osIdx = [ 1 : floor(Nread/2)*(1/os)   ,  1+ceil(Nread*(1-1/os/2)) : Nread ];
            data =  fft( data, [], 1);                      % go into spatial domain; phase issues due to missing fftshift() and fft() instead of ifft(), but they will be fixed in next line
            data = ifft( data(osIdx,:,:), [],1) ./ os;      % cut out what's interesting and go back to pure k-space
        end                                                 %  ifft() has a scaling factor 1/N, but we reduced N: N = N/os

        % Save normalized k-spaces
        %
        ksp(:,:,:,Nsmpl) = data ./ max(abs(data(:)));

        clear  data

    end % for smpl = 1:NumSamplings

    clear   mask  kR  % pos

    % follow regular matrix convention: Y x X x Z
    ksp = permute( ksp, [ 2 1 3 4 ] );

    % Calculate and normalize the Point Spread Functions
    %
    shifts = [szY  szRead  szZ  0 ] ./ 2;
    PSF = abs( ifft(ifft(ifft(  cmshiftnd(ksp,shifts), [],1),[],2),[],3) );
    PSF = cmshiftnd( PSF, shifts );
    PSF = bsxfun( @rdivide, PSF, max(max(max( PSF,[],1),[],2),[],3) );

    fprintf('Finished simulations after %g seconds.\n', toc )

    maxVal = 1.0;   % if maxVal < 1: cut off the PSF at maxVal, so we can better see the interesting, non-zero parts


    % free as much memory as possible as early as possible
    switch nargout
        case 0
            clear  ksp  PSF
        case 1
            clear  ksp
    end

    myColormap = morgenstemning(256);
    
    if szZ == 1
      % 2D:

        figure
        title(sprintf('TR = %g ms,  %d lines, %g s measurment time', 1e3*TR, NUsedLines, max(t) ))
        for n = 1:NumSamplings
            subplot( 2, NumSamplings, n )
            imagesc( abs(ksp(:,:,1,n)), [0  maxVal] )
            % imagesc(angle(ksp(:,:,1,n)) ./ pi)
            title( ['k-space, \bf' samplingType{n,1}] )
            axis image
            axis off
            colormap( myColormap )
            
            subplot( 2, NumSamplings, n + NumSamplings )
            imagesc( abs(PSF(:,:,1,n)), [0  maxVal] )
            title( ['PSF, \bf' samplingType{n,1}] )
            axis image
            axis off
            colormap( myColormap )
        end

    else
      % 3D:

        % Define what we want to show.
        % We start with "ksp", and append "PSF", afterwards.
        %
        toPlot = abs( ksp );
        %toPlot = angle( ksp )  ./ pi;

        toPlot = cat( 1, maxVal * toPlot,  PSF );   % catenate k-spaces and PSFs along vertical direction
        toPlot = permute( toPlot, [1 2 4 3] );
        toPlot = reshape( toPlot, [], NumSamplings*szRead, szZ );

        titleString = sprintf('TR = %g ms,  %d lines, %g s measurment time\n', 1e3*TR, NUsedLines, max(t) );

        titleString = [ titleString  sprintf( 'top: k-space magnitude    /    bottom: PSFs\n\\bf' ) ];
        for n = 1:NumSamplings
            titleString = [ titleString   samplingType{n,1}  '                  ' ];
        end
      
        disp( '===================================================================' )
        disp( 'In the window, navigate in partition direction using up/down arrow.' )
        disp( '===================================================================' )

        toPlot = reshape( toPlot, [], szRead * NumSamplings, 1, szZ );
        figure
        imdisp_pe( toPlot, 'Size', 1, 'DisplayRange', [0  maxVal] )
        title( titleString )
        colorbar
        colormap( myColormap )
    end

    plotMotion = true;
    
    if plotMotion
        figure
        hold on
        plot( t, pos(:,1), '+b')
        plot( t, pos(:,2), 'xr')
        plot( t, pos(:,3), 'sk')
        hold off
        xlabel('t / s')
        ylabel('pixels')
        legend('\Delta X', '\Delta Y', '\Delta Z')
    end

end % of sim_motion_cartesian()



% ====== INPUT VALIDATION HELPER(S) ======

function isValid = checkMask( mask )
    isValid = false;
    helpStr = ['Input must be a binary vector (only 1/0 or true/false values)\n'     ...
               'indicating the acquired points in time, or a single '':'', to index every line.'];

    try
        validateattributes( mask, {'numeric', 'logical'}, {'nonempty', '2d', 'binary'} )
    catch
        if ~isequal(mask, ':')   % the only non-scalar input allowed is ':', to index everything
            error( 'checkMask:BadScalar', helpStr )
        end
    end
    isValid = true;

end % of checkMask()



%###########################################################################################################################################
%#                                                                                                                                        ##
%#                                                                                                                                        ##
%#  ################                       ##        ##                                 ##               ##                               ##
%#       ##         ##                     ##        #                                  ##               ##                               ##
%#       ##         ##                     ##                                           ##               ##                               ##
%#       ##         ## ###    ####         ##        ##  ## ###    ####    #####        ## ###    ####   ##         ####   ##  ##  ##     ##
%#       ##         ###  ##  ##  ##        ##        ##  ###  ##  ##  ##  ##            ###  ##  ##  ##  ##        ##  ##  ##  ##  ##     ##
%#       ##         ##   ##  #####         ##        ##  ##   ##  #####    ####         ##   ##  #####   ##        #    #  ##  ##  ##     ##
%#       ##         ##   ##  ##   #         ##   ##  ##  ##   ##  ##   #      ##        ###  ##  ##   #   ##   ##  ##  ##  ##  ##  ##     ##
%#       ##         ##   ##   ####           ####    ##  ##   ##   ####   #####         ## ###    ####     ####     ####    ###  ###      ##
%#                                                                                                                                        ##
%#                                                                                                                                        ##
%#                                 ##                     ##                                                 ##   #######                 ##
%#                                 ##                     #                                                  ##   ##   ##                 ##
%#                                 ##                                                                        ##    ## ##                  ##
%#   ####   ### ##  ## ###         ## ###    ####         ##   ### ##  ## ###    ####   ## ###   ####    ### ##     ###                   ##
%#  ##     ##  ###  ###  ##        ###  ##  ##  ##        ##  ##  ###  ###  ##  ##  ##  ###     ##  ##  ##  ###      #                    ##
%#  #      #    ##  ##   ##        ##   ##  #####         ##  #    ##  ##   ##  #    #  ##      #####   ##   ##                           ##
%#  ##     ##  ###  ##   ##        ###  ##  ##   #        ##  ##  ###  ##   ##  ##  ##  ##      ##   #  ##  ###      #                    ##
%#   ####   ### ##  ##   ##        ## ###    ####         ##   ### ##  ##   ##   ####   ##       ####    ### ##     ###                   ##
%#                                                                 ##                                                                     ##
%#                                                            ##   ##                                                                     ##
%#                                                             #####                                                                      ##
%#                                                                                                                                        ##
%#                                                                                                                                        ##
%#   (They contain only dependencies)                                                                                                     ##
%#                                                                                                                                        ##
%###########################################################################################################################################

function x = makeCol( x )
    x = x(:);
end
function x = makeRow( x )
    x = x(:).';
end

function x = cmshiftnd( x, shifts)

%Function to circularly shift N-D arrays
%M.A. Griswold 9/11/98

% minor tweaks:
%   *   x = f(x) may be faster / more RAM-economic
%   *   return immediately if nargin < 2 or shifts is all zeros
% M. Völker 03.12.2012

    if nargin < 2 || all(shifts(:) == 0)
       return                       % no shift
    end

    sz      = size( x );
    numDims = ndims(x);             % number of dimensions
    idx = cell(1, numDims);         % creates cell array of empty matrices,
                                    % one cell for each dimension

    for k = 1:numDims               %#ok <--- no parfor-hint

        m = sz(k);
        p = ceil(shifts(k));

        if p < 0
            p = m + p;
        end

        idx{k} = [p+1:m  1:p];
    end

    % Use comma-separated list syntax for N-D indexing.
    x = x(idx{:});

end % of cmshiftnd()


function hIm = imdisp_pe(I, varargin)
%IMDISP  Display one or more images nicely
%
% Examples:
%   imdisp
%   imdisp(I)
%   imdisp(I, map)
%   imdisp(..., param1, value1, param2, value2, ...)
%   h = imdisp(...)
%
% This function displays one or more images nicely. Images can be defined
% by arrays or filenames. Multiple images can be input in a cell array or
% stacked along the fourth dimension, and are displayed as a grid of
% subplots (an improvement over MONTAGE). The size of grid is calculated or
% user defined. The figure size is set so that images are magnified by an
% integer value.
%
% If the image grid size is user defined, images not fitting in the grid
% can be scrolled through using the following key presses:
%    Up - Back a row.
%    Down - Forward a row.
%    Left - Back a page (or column if there is only one row).
%    Right - Forward a page (or column if there is only one row).
%    Shift - 2 x speed.
%    Ctrl - 4 x speed.
%    Shift + Ctrl - 8 x speed.
%
% This allows fast scrolling through a movie or image stack, e.g.
%    imdisp(imstack, 'Size', 1)
% The function can be used as a visual DIR, e.g.
%    imdisp()
% to display all images in the current directory on a grid, or
%    imdisp({}, 'Size', 1)
% to scroll through them one at a time.
%
% IN:
%   I - MxNxCxP array of images, or 1xP cell array. C is 1 for indexed
%       images or 3 for RGB images. P is the number of images. If I is a
%       cell array then each cell must contain an image. Images can equally
%       be defined by filenames. If I is an empty cell array then all the
%       images in the current directory are used. Default: {}.
%   map - Kx3 colormap to be used with indexed images. Default: gray(256).
%   Optional parameters - name, value parameter pairs for the following:
%      'Size' - [H W] size of grid to display image on. If only H is given
%               then W = H. If either H or W is NaN then the number of rows
%               or columns is chosen such that all images fit. If both H
%               and W are NaN or the array is empty then the size of grid
%               is chosen to fit all images in as large as possible.
%               Default: [].
%      'Indices' - 1xL list of indices of images to display. Default: 1:P.
%      'Border' - [B R] borders to give each image top and bottom (B) and
%                 left and right (R), to space out images. Borders are
%                 normalized to the subplot size, i.e. B = 0.01 gives a
%                 border 1% of the height of each subplot. If only B is
%                 given, R = B. Default: 0.01.
%      'DisplayRange' - [LOW HIGH] display range for indexed images.
%                       Default: [min(I(:)) max(I(:))].
%      'Map' - Kx3 colormap or (additionally from above) name of MATLAB
%              colormap, for use with indexed images. Default: gray(256).
%
%%%%% ADDED BY P. EHSES %%%%%
%      'AspectRatio' - [scaleX scaleY scaleZ] change aspect ratio.
%                      since output is 2D, [scaleX scaleY] is sufficient
%                      (scaleZ only for consistency with other matlab fcts)
%                      Default: [1 1 1]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% OUT:
%   h - HxW array of handles to images.
%
%   See also IMAGE, IMAGESC, IMSHOW, MONTAGE.

% $Id: imdisp.m,v 1.5 2009/03/31 17:25:57 ojw Exp $

    % Parse inputs
    [map layout gap indices lims aspectratio] = parse_inputs(varargin);

    if nargin == 0 || (iscell(I) && isempty(I))
        % Read in all the images in the directory
        I = get_im_names;
        if isempty(I)
            % No images found
            if nargout > 0
                hIm = [];
            end
            return
        end
    end

    % Check if input is filenames
    if ischar(I)
        [x y c] = size(I);
        if (x > 1 && y > 1) || c > 1
            I = num2cell(I, 2);
        else
            I = {I(:)'};
        end
    end

    % Get limits, etc.
    if isnumeric(I)
        [y x c n] = size(I);
        if isempty(lims)
            lims = min_max(I);
        elseif isequal(0, lims)
            lims = default_limits(I);
        end
        if isfloat(I) && c == 3 && n > 1
            I = uint8(I * 256 - 0.5);
            lims = round(lims * 256 - 0.5);
        end
        I = squeeze(num2cell(I, [1 2 3]));
    elseif iscell(I)
        A = I{1};
        if ischar(A)
            A = imread_rgb(A);
            I{1} = A;
        end
        % Assume all images are the same size and type as the first
        [y x c] = size(A);
        if isempty(lims) || isequal(0, lims)
            lims = default_limits(A);
        end
    else
        error('I not of recognized type.');
    end

    % Select indexed images
    if ~isequal(indices, -1)
        I = I(indices);
    end
    n = numel(I);

    % Get the current figure
    hFig = get(0, 'CurrentFigure');
    if isempty(hFig)
        % Create a new figure
        hFig = figure;
    elseif n > 1
        % Clear the figure
        clf(hFig, 'reset');
    end

    % Set the colormap
    set(hFig, 'Colormap', map);

    % Display the image(s)
    if n == 0
        hIm = display_image([], gca, [0 1],aspectratio);

        if nargout == 0
            clear hIm % Avoid printing this out
        end
        return
    elseif n == 1
        % IMSHOW mode
        % Display the single image
        hAx = gca;
        hIm = display_image(I{1}, hAx, lims,aspectratio);

        if nargout == 0
            clear hIm % Avoid printing this out
        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');
        newAxesPos = [gap(1) gap(end) 1-2*gap(1) 1-2*gap(end)];
        if isequal(axesPos, get(hFig, 'DefaultAxesPosition'))
            % Default position => not a subplot
            % Fill the window
            set(hAx, 'Units', 'normalized', 'Position', newAxesPos);
            axesPos = newAxesPos;
        end
        if ~isequal(axesPos, newAxesPos) || strcmp(get(hFig, 'WindowStyle'), 'docked')
            % Figure not alone, or docked. Either way, don't resize.
            return
        end
        layout = [1 1];
    else
        % MONTAGE mode
        % Compute a good layout
        layout = choose_layout(n, y, x, layout,aspectratio);

        % Create a data structure to store the data in
        num = prod(layout);
        state.n = num * ceil(n / num);
        hIm = zeros(layout);
        hAx = zeros(layout);
        I = [I(:); cell(state.n-n, 1)];

        % Set the first lot of images
        index = mod(0:num-1, state.n) + 1;
        hw = 1 ./ layout;
        gap = gap ./ layout;
        dims = hw - 2 * gap;
        dims = dims([2 1]);
        for a = 1:layout(1)
            for b = 1:layout(2)
                c = index(b + (layout(1) - a) * layout(2));
                A = I{c};
                if ischar(A)
                    A = imread_rgb(A);
                    I{c} = A;
                end
                hAx(a,b) = axes('Position', [(b-1)*hw(2)+gap(2) (a-1)*hw(1)+gap(1) dims], 'Units', 'normalized');
                hIm(a,b) = display_image(A, hAx(a,b), lims,aspectratio);
            end
        end

        % Check if we need to be able to scroll through images
        if n > num
            % Intialize rest of data structure
            state.hIm = hIm;
            state.hAx = hAx;
            state.index = 1;
            state.layout = layout;
            state.I = I;
            % Set the callback for image navigation, and save the image data in the figure
            set(hFig, 'KeyPressFcn', @keypress_callback, 'Interruptible', 'off', 'UserData', state);
        end

        if nargout == 0
            clear hIm % Avoid printing this out
        end
    end

    % Set the figure size well
    % Compute the image size
    % ImSz = layout([2 1]) .* [x y] ./ (1 - 2 * gap([end 1]));

    % pehses:
    ImSz = layout([2 1]) .* [x y] ./ aspectratio(1:2) ./ (1 - 2 * gap([end 1]));

    % 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];
    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 the figure size and position
    set(hFig, 'Position', figPosNew);
    return
end

% Keypress callback
% The function which does all the display stuff
function keypress_callback(fig, event_data)
    % Check what key was pressed and update the image index as necessary
    switch event_data.Character
        case 28 % Left
            up = -1; % Back a page
        case 29 % Right
            up = 1; % Forward a page
        case 30 % Up
            up = -0.1; % Back a row
        case 31 % Down
            up = 0.1; % Forward a row
        otherwise
            % Another key was pressed - ignore it
            return
    end
    % Use control and shift for faster scrolling
    if ~isempty(event_data.Modifier)
        up = up * (2 ^ (strcmpi(event_data.Modifier, {'shift', 'control'}) * [1; 2]));
    end
    % Get the state data, if not given
    state = get(fig, 'UserData');
    % Get the current index
    index = state.index;
    % Get number of images
    n = prod(state.layout);
    % Generate 12 valid indices
    if abs(up) < 1
        % Increment by row
        index = index + state.layout(2) * (up * 10) - 1;
    else
        if state.layout(1) == 1
            % Increment by column
            index = index + up - 1;
        else
            % Increment by page
            index = index + n * up - 1;
        end
    end
    index = mod(index:index+n, state.n) + 1;
    % Plot the images
    figure(fig);
    for a = 1:state.layout(1)
        for b = 1:state.layout(2)
            % Get the image
            c = index(b + (state.layout(1) - a) * state.layout(2));
            A = state.I{c};
            if ischar(A)
                % Filename - read the image from disk
                A = imread_rgb(A);
                state.I{c} = A;
            end
            % Set the image data
            set(state.hIm(a,b), 'CData', A);
            % Reset the axes limits
            if ~isempty(A)
                set(state.hAx(a,b), 'XLim', [0.5 size(A, 2)+0.5], 'YLim', [0.5 size(A, 1)+0.5]);
            end
        end
    end
    drawnow;
    % Save the current index
    state.index = index(1);
    set(fig, 'UserData', state);
    return
end

% Display the image
function hIm = display_image(A, hAx, lims, aspectratio)
    if isempty(A)
        hIm = image(zeros(1, 1, 3));
        set(hIm, 'CData', []);
    else
        hIm = image(A);
    end
    set(hAx, 'Visible', 'off', 'DataAspectRatio', aspectratio, 'DrawMode', 'fast', 'CLim', lims);
    set(get(hAx, 'XLabel'), 'Visible', 'on');
    set(get(hAx, 'YLabel'), 'Visible', 'on');
    set(get(hAx, 'Title'), 'Visible', 'on');
    set(hIm, 'CDataMapping', 'scaled');
    return
end

% Choose a good layout for the images
function layout = choose_layout(n, y, x, layout, aspectratio) % pehses: wip: have to include aspectratio in calculation
    layout = reshape(layout, 1, min(numel(layout), 2));
    v = numel(layout);
    N = isnan(layout);
    if v == 0 || all(N)
        sz = get(0, 'ScreenSize');
        sz = sz(3:4) ./ [x y];
        layout = ceil(sz([2 1]) ./ sqrt(prod(sz) / n));
        switch ([prod(layout - [1 0]) prod(layout - [0 1])] >= n) * [2; 1]
            case 0
            case 1
                layout = layout - [0 1];
            case 2
                layout = layout - [1 0];
            case 3
                if min(sz .* (layout - [0 1])) > min(sz .* (layout - [1 0]))
                    layout = layout - [0 1];
                else
                    layout = layout - [1 0];
                end
        end
    elseif v == 1
        layout = layout([1 1]);
    elseif any(N)
        layout(N) = ceil(n / layout(~N));
    end
    return
end

% Read image to uint8 rgb array
function A = imread_rgb(name)
    try
        [A map] = imread(name);
    catch
        % Format not recognized by imread, so create a red cross (along diagonals)
        A = eye(101) | diag(ones(100, 1), 1) | diag(ones(100, 1), -1);
        A = uint8(255 * (1 - (A | flipud(A))));
        A = cat(3, zeros(size(A), 'uint8')+uint8(255), A, A);
        return
    end
    if ~isempty(map)
        map = uint8(map * 256 - 0.5);
        A = reshape(map(A,:), [size(A) size(map, 2)]);
    elseif size(A, 3) == 4
        ll = lower(name(end));
        if ll == 'f'
            % TIFF in CMYK colourspace - convert to RGB
            error('CMYK image files not yet supported - please fix.');
        elseif ll == 's'
            % RAS in RGBA colourspace - convert to RGB
            error('RGBA image files not yet supported - please fix.');
        end
    end
    return
end

% Get the names of all images in a directory
function L = get_im_names
    D = dir;
    n = 0;
    L = cell(size(D));
    % Go through the directory list
    for a = 1:numel(D)
        % Check if file is a supported image type
        if numel(D(a).name) > 4 && ~D(a).isdir && (any(strcmpi(D(a).name(end-3:end), {'.png', '.tif', '.jpg', '.bmp', '.ppm', '.pgm', '.pbm', '.gif', '.ras'})) || any(strcmpi(D(a).name(end-4:end), {'.tiff', '.jpeg'})))
            n = n + 1;
            L{n} = D(a).name;
        end
    end
    L = L(1:n);
    return
end

% Parse inputs
function [map layout gap indices lims aspectratio] = parse_inputs(inputs)

    % Set defaults
    map = [];
    layout = [];
    gap = 0;
    indices = -1;
    lims = 0;
    aspectratio = [1 1 1];

    % Check for map
    if numel(inputs) && isnumeric(inputs{1}) && size(inputs{1}, 2) == 3
        map = inputs{1};
        inputs = inputs(2:end);
    end

    % Go through option pairs
    for a = 1:2:numel(inputs)
        switch lower(inputs{a})
            case 'map'
                map = inputs{a+1};
                if ischar(map)
                    map = feval(map, 256);
                end
            case {'size', 'grid'}
                layout = inputs{a+1};
            case {'gap', 'border'}
                gap = inputs{a+1};
            case 'indices'
                indices = inputs{a+1};
            case {'lims', 'displayrange'}
                lims = inputs{a+1};
            case {'aspectratio'} % pehses
                aspectratio = inputs{a+1};

                % normalize aspectratio:
                aspectratio = aspectratio/min(aspectratio);

                % for 2D images, only [x y] is really necessary
                % but the definition is [x y z], so fix that
                if length(aspectratio) == 2
                    aspectratio = [aspectratio 1];
                end

            otherwise
                error('Input option %s not recognized', inputs{a});
        end
    end

    if isempty(map)
       map = gray(256);
    end
    return
end

% Return default limits for the image type
function lims = default_limits(A)
    if size(A, 3) == 1
        lims = min_max(A);
    else
        lims = [0 1];
        if ~isfloat(A)
            lims = lims * double(intmax(class(A)));
        end
    end
    return
end

% Return minimum and maximum values
function lims = min_max(A)
    M = isfinite(A);
    lims = [min(A(M)) max(A(M))];
    return
end

function displaytable(data, colheadings, wid, fms, rowheadings, fid, colsep, rowending)
% Prints formatted matrix of numerical data with headings
%
% Syntax
%
% displaytable(data, colheadings, wid, fms, rowheadings, fid, colsep, rowending)
%
% Input
%
% data: a matrix or cell array, containing the data to be put in the table.
% If a matrix the numbers in the table will be printed using the default
% format specifier (f), or the specifier(s) supplied in fms. data can also
% be a cell array containing a mixture of strings and numbers. each cell in
% this case can contain only a single scalar value. In this case numbers
% are printed as in the matrix case, while strings are printed up to the
% maximum width of the column.
%
% colheadings: a cell array of strings for the headings of each column. Can
% be an empty cell array if no headings are required. If no column widths
% are supplied (see below), the columns will be made wide enough to
% accomdate the headings.
%
% wid: (optional) scalar or vector of column widths to use for the table.
% If scalar, every column will have the same width. If a vector it must be
% of the same length as the number of columns of data. If not supplied, and
% column headers are supplied, the width of the widest column header will
% be used. If not supplied and column headers are not supplied, a default
% with of 16 characters is used.
%
% fms: (optional) a string, or cell array of strings containing format
% specifiers for formatting the numerical output in each column. If a
% single string, the same specifier is used for every column. If not
% supplied, the 'g' specifier is used for every column.
%
% rowheadings: (optional) a cell array of strings for the start of each
% row. Can be an empty cell array if no row headings are required. If row
% headings are supplied, the first column will be made wide enough to
% accomodate all the headings.
%
% fid: (optional) the file id to print to. Use 1 for stdout (to print to
% the command line).
%
% colsep: (optional) A string or character to insert between every column.
% The default separation string is ' | ', i.e. a space followed by a
% vertical bar, followed by a space. A table suitible for inclusion in a
% LaTeX document can be created using the ' & ' string, for example.
%
% rowending: (optional) An optional string or character to be appended at
% the end of every row. Default is an empty string.
%
% Example 1 - Basic useage
%
% colheadings = {'number of projects','sales','profit'};
% rowheadings = {'Jimmy Slick', 'Norman Noob'}
% data = [3 rand(1) rand(1); 1 rand(1) rand(1)];
%
% To format the first number in each row as a decimal (%d), the second
% number %16.4f, and the third as %16.5E do the following:
%
% wid = 16;
% fms = {'d','.4f','.5E'};
%
% In this case 16 will be the field width, and '.5E' is what to use for the
% fms argument
%
% fileID = 1;
%
% >> displaytable(data,colheadings,wid,fms,rowheadings,fileID);
%             |number of projec |           sales |          profit
% Jimmy Slick |               3 |          0.4502 |    5.22908E-001
% Norman Noob |               1 |          0.9972 |    2.78606E-002
%
% Example 2 - Produce a latex table
%
% colheadings = {'number of projects','sales','profit'};
% rowheadings = {'Jimmy Slick', 'Norman Noob'};
% data = [3 rand(1) rand(1); 1 rand(1) rand(1)];
% wid = 16;
% fms = {'d'};
%
% colsep = ' & ';
% rowending = ' \\';
%
% fileID = 1;
%
% >> displaytable(data,colheadings,wid,fms,rowheadings,fileID,colsep,rowending);
%
%             & number of projec &            sales &           profit \\
% Jimmy Slick &                3 &    6.948286e-001 &    3.170995e-001 \\
% Norman Noob &                1 &    9.502220e-001 &    3.444608e-002 \\
%
% Example 3 - Mixed numeric and strings
%
% colheadings = {'number of projects','sales','profit'};
% rowheadings = {'Jimmy Slick', 'Norman Noob'};
%
% data = {3, rand(1), rand(1);
%         1, 'none', 0};
%
% wid = 16;
% fms = {'d','.4f','.5E'};
%
% fileID = 1;
%
% >> displaytable(data,colheadings,wid,fms,rowheadings,fileID);
%             | number of projec |            sales |           profit
% Jimmy Slick |                3 |           0.4854 |     8.00280E-001
% Norman Noob |                1 |             none |     0.00000E+000
%
% See Also: DISP, FPRINTF

% Created by Richard Crozier 2012

    if nargin < 8 || isempty(rowending)
        rowending = '';
    end

    if nargin < 7 || isempty(colsep)
        colsep = ' | ';
    end

    % do some basic checking of input
    if nargin < 6 || isempty(fid)
        % print to the command line
        fid = 1;
    end

    if nargin < 5 || isempty(rowheadings)
        % no row headings supplied, use empty cell array
        rowheadings = {};
    end

    if nargin < 4 || isempty(fms)
        % no format specifiers supplied, use 'g' for all columns
        fms = 'g';
    end

    if nargin < 3 || isempty(wid)
        % default width is 10, this will be modified if column headers are
        % supplied
        wid = 10;
    end

    if nargin < 2
        colheadings = {};
    end

    if nargin >= 6 && ~iscellstr(rowheadings)
        error ('row headings must be cell array of strings');
    end

    % get the numbers of rows and columns of data
    [nRowD,nColD] = size(data);

    % check that sensible format specifiers have been supplied
    if ~iscellstr(fms)

        if ischar(fms)
            fms = repmat({fms},1,nColD);
        else
            error('fms must be a string or cell array of strings.');
        end

    elseif isempty(fms)

        fms = repmat({'f'},1,nColD);

    elseif numel(fms) ~= nColD

        if numel(fms) == 1
            fms = repmat(fms,1,nColD);
        else
           error('fms does not have the same number of format specifiers as there are columns.');
        end

    end

    % replace empty format specifiers with 'f'
    for i = 1:numel(fms)
        if isempty(fms{i})
            fms{i} = 'f';
        end
    end

    [nRowFms,nColFms] = size(fms);
    if(nRowFms>1)
        error ('fms can not have more than one row');
    end


    if ~isempty(rowheadings)

        if ~iscellstr(rowheadings)
            error('rowheadings must be a cell array of strings');
        end

        if numel(rowheadings) ~= size(data, 1)
            error('Rowheadings must be a cell array of strings of the same size as the number of rows in data.')
        end

        rhwid = 0;

        for r = 1:numel(rowheadings)
            % find the maximum width the first column must be to accomodate
            % the row headings
            rhwid = max(rhwid, length(rowheadings{r}));
        end

    end

    if isscalar(wid)

        tempwid = zeros(1, numel(fms));

        for i = 1:numel(fms)

            % get the number of decimal places specified
            [start_idx, end_idx, extents, matches] = regexp(fms{i}, '(?<=\.)\d+');

            if isempty(start_idx)

                % no numbers were found in the format spec, just use the
                % supplied width
                tempwid(i) = wid;

            else

                % some numbers were supplied, use the larger of the width
                % or the number of desired decimal places plus two (to
                % allow for leading number plus decimal point)
                tempwid(i) = max(wid, round(str2double(matches{1})) + 2);

            end

        end

        % replace scalar width with array of width values
        wid = tempwid;

    end

    if ~isempty(colheadings)

        [nRowCH,nColCH] = size(colheadings);
        if(nRowCH>1)
            error ('column headings can not have more than one row');
        end

        if ~iscellstr(colheadings)
            error('If not empty, colheadings must be a cell array of strings');
        end

        if(nColCH ~= nColD)
            error ('data must have same number of columns as headings');
        end

%         fmt = arrayfun(@(x) ['%',num2str(wid(x)),'s |'], 1:nColD, 'UniformOutput', false);

        if ~isempty(rowheadings)
            % TODO allow extra heading for column
            fprintf(fid, ['%s',colsep], repmat(' ', 1,rhwid));
        end

        if nargin < 3

            % only data and column headings have been supplied, so
            % determine a sensible value for the width of each column

            tempwid = zeros(size(wid));
            for i = 1:numel(colheadings)

                % get a column width which is the minimum to accept the
                % column heading length or the default width specification
                tempwid(i) = max(length(colheadings{i}), wid(i));

                if tempwid < 1
                    error('Column width is less than 1, and the column header is empty.')
                end

            end

            wid = tempwid;

        end

        % Now loop through the column headings printing them out with the
        % desired column separator
        for i = 1:numel(colheadings)

            str = sprintf(['%',num2str(wid(i)),'s'],colheadings{i});

            if i == numel(colheadings)
                % If we are at the end of a row, don't print the column
                % separator
                fprintf(fid, '%s', str(1:wid(i)) );
            else
                % print the column header and column separator
                fprintf(fid, ['%s',colsep], str(1:wid(i)) );
            end

        end

        fprintf(fid, '%s\n', rowending);

    end

    fmt = arrayfun(@(x) ['%',num2str(wid(x)),fms{x}],1:nColD,'UniformOutput',false);

    for i = 1:size(data,1)

        % first print a row header if necessary
        if ~isempty(rowheadings)
            fprintf(fid, ['%',num2str(rhwid),'s',colsep], rowheadings{i});
        end

        % now loop through the data formatting and printing as appropriate
        for j = 1:size(data,2)

            if iscell(data)
                % data is a cell array, possibly containing mixed data
                % types

                if ischar(data{i,j})

                    str = sprintf(['%',num2str(wid(j)),'s'],data{i,j});

                elseif isscalar(data{i,j})

                    % write the number
                    str = sprintf(fmt{j},data{i,j});

                    if length(str) > wid(j)

                        % convert to scientific notation as it doesn't fit
                        str =  sprintf(['%',num2str(wid(j)),'g'],data{i,j});

                        if length(str) > wid(j)
                            % fill space with #s as the number stil doesn't
                            % fit in
                            str = repmat('#', 1, wid(j));
                        end

                    end

                elseif isempty(data{i,j})

                    % indicate an empty value
                    str = sprintf(['%',num2str(wid(j)),'s'],'');

                else
                    % we can only tabulate strings and scalars, so throw an
                    % error
                    error('CROZIER:displaytable:badcellcontents', ...
                          'each cell in cell array data must contain only individual strings or scalar values.')

                end

                % print the string
                fprintf(fid, '%s', str(1:wid(j)));

                if j < size(data,2)
                    % print column separator
                    fprintf(fid, colsep);
                end

            else

                % data is a numerical matrix

                str = sprintf(fmt{j},data(i,j));

                if length(str) > wid(j)

                    % convert to scientific notation as it doesn't fit
                    str =  sprintf(['%',num2str(wid(j)),'g'],data(i,j));

                    if length(str) > wid(j)
                        % fill space with #s as the number doesn't fit in
                        str = repmat('#', 1, wid(j));
                    end

                end

                if j == size(data,2)
                    % do not print the last column separator at the end of
                    % a row
                    fprintf(fid, '%s', str(1:wid(j)));
                else
                    fprintf(fid, ['%s',colsep], str(1:wid(j)));
                end
            end

        end

        % end of line so put in a new row
        fprintf(fid, '%s\n', rowending);

    end

end % of displaytable()

Contact us