Table of Contents

flexinterpn

Note: this is a MEX (compiled) function, and it is essential for the visualization and many data processing functions of NeuroElf.

Motivation

When I started with writing visualization functions, I soon discovered that the interpn function of Matlab is accurate and, in some regards, much more flexible but also, to my dismay, not fast enough to handle real-time updates for when the cursor is moved in a 3-canonical-slicing view of the brain.

After looking up several helpful links about interpolation, I decided to write a C-based function that would allow the following operations:

With respect to that last item, I myself am happy enough with the result (although there are most likely still faster interpolation algorithms out there…)

Function reference ('help flexinterpn')

  flexinterp  - flexible data interpolation (up to 4D)
 
  FORMAT:       idata = flexinterpn(data, coords [, k, ks [, hold [, qt]]])
 
  Input fields:
 
        data        up to 4D data
        coords      CxD (number of coords -by- number of dims) coordinates
                or  4xD range specification with
                    1st row: Inf
                    2nd row: from
                    3rd row: step
                    4th row: to, such that
                    a given dim's range is sample at from:step:to
        kernel      Kx1 interpolation kernel (e.g. from sinc(x)/sinc(x/a))
        ksampling   kernel sampling rate (in samples)
        hold        default sample at invalid coordinates
        qtrf        quaternion transformation applied to a given range
 
  Output fields:
 
        idata       interpolated data
 
  Example:
 
        data = randn(16, 16, 16);
        kernel = sinc(-3:0.001:3) .* sinc(-1:1/3000:1);
        kernel([1,end]) = 0;
        meshing = 0.625:0.25:16.375;
        [x, y, z] = meshgrid(meshing, meshing, meshing);
        idata = reshape( ...
            flexinterpn(data, [x(:), y(:), z(:)], kernel(:), 1000), ...
            [64, 64, 64]);
 
  Observe that the sample result is obtained by using...
 
        idata = flexinterpn(data, repmat([Inf; 0.625; 0.25; 16.375], 1, 3), ...
            kernel(:), 1000);
 
  Note: This is a compiled function.
        Also, make sure the kernel has a 0 value at the beginning and end
        e.g. by issuing k([1, end]) = 0;
        Lastly, if different kernels are to be used (only valid for
        coordinate ranges!), both the kernel and ksampling arguments can
        be given as a 1-by-ndims(data) cell array

Arguments

data

The data argument simply is the data you want to interpolate. It is assumed to be evenly spaced (that is to say constant sampling rate in all dimensions) and the first coordinate is 1 in each dimension (using Matlab's 1-based indexing; automatically taken care of by the function internally). All numeric datatypes are supported; that is to say, there is no need to convert any data to double/single precision before using flexinterpn!

coords

Coordinates can be given as

kernel

The kernel can be any (hopefully symmetric, but even that is not a must!) double array (or cell array of double arrays, if different kernels are to be used!). The number of elements must be odd (with one central elements, which usually is the maximum value of the entire kernel). The number of elements must be selected so that the following condition holds true: mod(numel(kernel) - 1, kernelsize) == 0.

If, say, the kernel has the following values (which, BTW, is a very low-resolution cubic interpolation kernel)

kernel = [0; -0.025; -0.0625; -0.07; 0; 0.227; 0.5625; 0.867; 1; 0.867; 0.5625; 0.227; 0; -0.07; -0.0625; -0.025; 0];

the kernel sampling (in this case) is 4. If now a coordinate is requested that cannot be satisfied by the kernel directly (e.g. 3.2, which is close to 3.25 but not just…), then the kernel weights themselves are linearly interpolated. That means that the kernel must be given in a “good enough” quality so that linearly interpolating this kernel will not lead to quality loss! The flexinterpn_method function (which will create the more typical kernels) internally uses a 4096-bin-per-sample resolution (e.g. the cubic interpolation kernel has 16385 values!).

If a cell array is supplied, the ksampling argument must also be a cell array (see below).

Note: the first and last element of the (or rather each) kernel MUST be set to zero to avoid symmetry problems! (this condition is met by all of the standard kernels supplied by flexinterpn_method):

kernel([1, end]) = 0;

To obtain one of these “standard” kernels, use the following syntax:

flexinterpn_getstandardkernel.m
% call flexinterpn_method with single data point
% the kernel (second output) is then a 1-by-2 cell array suitable for
% argument expansion in calls to flexinterpn
[null, kernel] = flexinterpn_method(1, 1, 'cubic');
 
% visualize kernel
plot(k{1});
 
% call flexinterpn with that kernel
r = rand(10, 10);
ri = flexinterpn(r, [Inf, Inf; 1, 1; 1, 1; 10, 10], k{:}, 0);

ksampling

Sampling rate of kernel (ratio between coordinates in data and elements in kernel).

Here is a more complex example… To spatially smooth a VTC, a Gaussian kernel is created. Next, this kernel is applied only in dims 2 through 4):

flexinterpn_3Din4Dsmoothing.m
% create gaussian kernel with size 2 (in voxels!)
% values in the kernel smaller than 1e-4 are discarded
% as they don't play any big role for the sum of weights
k = smoothkern(2, 1e-4);
 
% a VTC is loaded (datatype: uint16)
vtc = xff('*.vtc');
 
% the data is smoothed in dims 2 through 4 -> store as single!
vtc.VTCData = single(flexinterpn(vtc.VTCData, ...
    [Inf, Inf, Inf, Inf; ones(2, 4); size(vtc.VTCData)], ...
    {[0; 1; 0], k, k, k}, {1, 1, 1, 1}, 0));
 
% some patches necessary
vtc.FileVersion = 3;
vtc.DataType = 2;

Please note that this is a rather unusual and clearly not intended use, but it demonstrates the flexibility of this interpolation function.

hold

Value for samples with invalid coordinates or if all source values are invalid (inf/nan). Default is 0 (and very few applications would need to change that, but it is a standard in interpolation algorithms, so it is implemented).

qtrf

If the input data is 3D and the coordinates are given as a range (see coords section), a quaternion matrix can be given, which will then be applied to each of the coordinates in turn. The quaternion is automatically patched so that it takes the 1-base/0-base difference between Matlab and C-code into account (in other words, please use matrices such as SPM does!)

Notes

The function has a few differences to standard interpolation: