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

flexinterpn.txt · Last modified: 2011/05/29 04:44 by jochen