flexinterpn

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

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:

• up to 4D data handling
• variable kernel (type, size, etc., such as linear, cubic, Lanczos, Gaussian)
• automatic handling of errors in kernel amplitude (e.g. sum of kernel is not exactly = 1)
• automatic handling of 4×4 quaternion matrices (for spatial transformations)
• handling of evenly-spaced coordinate ranges (so as not having to create large arrays to specify interpolation coordinates)
• being FAST (by means of pre-computed kernel windows instead of computing the interpolation function for each coordinate)

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…)

  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

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!

Coordinates can be given as

• a list (whereas the number of columns must match the number of data dimensions, if the data is a column vector, coords must be a sampling-points-by-1 vector!):

  flexinterpn_coordinate_list.m

  % create random 3D data
  data = randn(20, 20, 20);
   
  % create random coordinate to sample from
  coords = 1 + 19 .* rand(1000, 3);
   
  % sample
  sampled = flexinterpn(data, coords);

• a coordinate range specified like this:

  flexinterpn_coordinate_range.m

  % load image file
  im1 = neuroelf_file('s', 'text2');
   
  % show interesting zoom/slide effect (like Mac screen saver)
  % - loop over a counter value
  for c = 1:0.3:50
   
      % compute zoom factor
      d = 1 - c / 100;
   
      % resample data (3D, 3rd dim is the RGB color planes!)
      id = uint8(flexinterpn_method(im1, ...
          [  Inf  ,         Inf, Inf;
           0.1 * c, 150 - 3 * c,  1 ;
               d  ,     d      ,  1 ;
          [[0.1*c, 150-3*c] + d * [224, 424], 3]], 'cubic'));
   
      % show image
      image(id);
      drawnow;
      pause(0.001);
  end

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);

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.

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).

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!)

The function has a few differences to standard interpolation:

• the function does not assume 0-value samples beyond the edges of the defined data; instead it will not sample there
• the sum of samples (multiplied by the according kernel weight) is divided by the sum of weights that were used (auto-kernel-amplitude scaling); while this corrects for problems with cut-off kernels (e.g. if a smoothing kernel is cut at the edges, this does not create a problem in a way that the overall mean value of the data will decrease), it is not possible to use this function to actually multiply by a constant value!
• given that for an internal rotation support the data must be given as 3D, there is a “hidden feature” to allow rotating 2D data by first shiftdim-ming the data so that the first dimension is of size 1