Wiki

Clone wiki

sta-toolbox / Home

Welcome

This is the WIKI of the STA-toolbox.
The STA-toolbox is free software. License: GPL 3

A 64bit Ubuntu Linux is the recommended platform. The toolbox has been developed on Linux, but has been successfully compiled on Windows as well.

Feel free to check out the code or download the ZIP file containing the most current version.

welcome

Table of Contents



#Installation

Please consider the README.txt for details.


#Documentation ## C and C++ API The STA-toolbox provides a C api (stensor.h) and a C++ api (stafield.h). The first one is the counterpart of the matlab stafieldStruct interface, and the latter one the counterpart of the matlab stafield class (see matlab examples for details). You can find further details regarding the C and C++ here: stafield class (C / C++ api)

Examples (Matlab)

Preliminaries

The toolbox uses a data container to represent spherical tensor fields. We call this container a stafield. A stafield is a multidimensional array with attributes. The attributes are representing the properties of the spherical tensor field. The tensor field data is stored in a five dimensional array of size 2 x N x X x Y x Z. The dimensions are consistent with the C and C++ implementations. The first dimension represents the real and imaginary part of the field. The second component the tensor components. The remaining three dimensions are the image dimensions. In matlab (and cotave), the stafield is represented by the stafieldStruct structure. Note that there exists a stafield class in matlab as well. Since classes are not (yet) supported in octave, we will focus on the structure here. The stafield class implements all basic tensor operators as member functions, while for the stafieldStruct, the tensor operations are external function calls (like C++ vs C).

Given an 3D image, you can initialize a stafield with the stafieldStruct:

#!matlab

shape=[64,64,64];
img=randn(shape,'single');

ofield=stafieldStruct(img);
IMPORTANT: if your initial image is in double precision, then the toolbox uses double precision. If in single precision, then it uses single precision. If there is no specific reason, then single precision is highly recommend (in order to safe memory). So please don't forget to cast your image into single precision.

As you can see, a stafield has several attributes:

#!matlab
 ofield

ofield = 

         storage: 'STA_FIELD_STORAGE_R'
               L: 0
            type: 'STA_OFIELD_SINGLE'
            data: [5-D single]
           shape: [64 64 64]
    element_size: [1 1 1]
* element_size defines the extents of a voxel * L defines the tensor rank * type : we distinguish between four different field types. (1) 'STA_OFIELD_SINGLE' is the default and defines a single tensor field of order L. (2) 'STA_OFIELD_FULL' is a collection of tensor fields of order 0 up to order L. You can consider it as a container for a bandlimited expansion of an orientation field. (3,4) 'STA_OFIELD_ODD' and 'STA_OFIELD_EVEN' are basically full fields, however, the even or odd orders are missing. These fields are used to represent e.g. symmetric orientation functions. * storage : For many of our applications, the spherical tensors have a symmetry. For such tensors, it is sufficient to store only L+1 coefficients instead of all 2L+1. 'STA_FIELD_STORAGE_C' denotes the "full" tensor, 'STA_FIELD_STORAGE_R' the tensor with intra-symmetries. 'STA_FIELD_STORAGE_RF' is the 'STA_FIELD_STORAGE_R' version in Fourier space.

Cartesian 2 spherical tensor field conversion

The toolbox provides functions to convert spherical to Cartesian tensor fields and vice versa. One a field has been converted into a spherical tensor field, all our tools can be used for further analysis and visualization. The following example shows how to transform vector and matrix valued tensor fields to their spherical counterparts:

#!matlab
  %% drawing a smoothed cross
    img = zeros([65,65,65]); 
    img2D = zeros([65,65]); 
    img2D (33,:)=1; 
    img2D (:,33)=1; 
    img2D =imrotate(img2D,45,'bicubic','crop');
    img(:,:,33)=img2D;
    img=img(26:end-25,26:end-25,26:end-25);
    img=exp(-bwdist(img>0.99).^2/4);
    img=permute(img,[3,1,2]);

    %% showing cross
    figure(1);
    imagesc(squeeze(img(ceil(end/2),:,:)))
    axis equal
    colormap gray

    %%computing the image gradient (Cartesian)
    G=sta_Grad(img);
    size(G)
    %and the Hessian matrix (only six elements)
    H=sta_Hessian(img);
    size(H)
    %Gradient as asymetric second order tensor
    GM=sta_Grad(img,2);
    size(GM)
    %full Hessian matrix (9 elelemts)
    HF=sta_Hessian(img,false);
    size(HF)
    %creating a full 9 DOF tensor based on Hessian and gradient
    T=HF+GM;

    %plotting the spherical gradient
    fieldG=sta_c2s(G);
    fieldG
    figure(2);
    sta_glyph(fieldG,'gamma',-0.8,'background',img);
    axis equal
    colormap jet

    %plotting the spherical counterparts of the Hessian
    fieldH=sta_c2s(H);
    fieldH
    figure(3);
    sta_glyph(fieldH,'gamma',-0.8,'background',img);
    axis equal
    colormap jet

    %plotting the spherical counterparts of the Hessian (traceless part)
    figure(4);
    sta_glyph(stafieldStruct(fieldH,2),'gamma',-0.8,'background',img);
    axis equal
    colormap jet

    %plotting the combination of Hessian and gradient
    fieldT=sta_c2s(T);
    fieldT
    figure(5);
    sta_glyph(fieldT,'gamma',-0.8,'background',img);
    axis equal
    colormap jet

    %extracting and plotting the gradient from the tensor field
    figure(6);
    sta_glyph(stafield(fieldT,1),'gamma',-0.8,'background',img);
    axis equal
    colormap jet
s2c

Basic Tensor Operations

The toolbox implements five basic tensor operations:

operation matlab function
spherical products sta_prod
spherical derivatives sta_deriv
Laplace sta_lap
tensor FFT sta_fft, sta_ifft
multiplication sta_mult

Orientation Fields

In many of our applications, the tensor fields are expansion coefficients of local orientation functions. These orientation functions are representing local image features in an angular dependent manner. In these cases, the spherical tensors are spherical harmonic expansion coefficients of the orientation function. Features can be quickly computed using a specific set of basis functions. For instance, spherical Gaussian derivatives up to order 3 can be computed with:

#!matlab 

 % img to tensor field
 ifield=stafieldStruct(img) 
 % convolution kernel
 kernel=stafieldStruct('gauss',ifield.shape,5); 
 L=3;
 % initializing orientation field
 ofield=stafieldStruct(ifield.shape,L,'STA_OFIELD_FULL',...
   'STA_FIELD_STORAGE_R','single');
 % convolution
 ifield=sta_ifft(sta_prod(sta_fft(ifield),sta_fft(kernel),0,... 
     {'alpha',1/prod(ifield.shape)})); 
 % defining the 0 order field of the full field 
 ofield=stafieldStruct(ofield,ifield);
 % iteratively computing the higher order components
 for l=1:L, ...
  ofield=stafieldStruct(ofield,sta_deriv(stafieldStruct(ofield,l-1),1)); 
 end;
Since the projection of images to local orientation fields is a very important and frequently used step, the toolbox provides a memory efficient version. The following does the same as above, but uses a pure C version of the code:

#!matlab
  ifield=stafieldStruct(img);
  L=3;
  ofield=sta_steerFilt(ifield,'kernel','gauss','kparam',5,'L',L,'type','STA_OFIELD_FULL');

The toolbox supports the fast Gaussian Derivatives [1], SHOG [2], Gauss-Laguerre (the spherical counterpart of 3D Gaussian Hermite Moments) [3] and Gauss-Bessel functions [3] (A spherical Gabor transform).

High-Level Functions

The toolbox provides a high-level interface for several applications. These functions are based on the basic operations above. The most important ones are

description | matlab function | | ------------- | -------------| Fast computation of rotation invariant Gabor descriptor images | sta_localfeat_SGD_cfast | [3]| Fast computation of rotation invariant Gauss-Laguerre descriptor images | sta_localfeat_SLD_cfast |[3] | Expanding orientation fields in terms of tensorial harmonics | sta_th_exp | [7] | Computes the power- and bispectra from orientation fields | sta_invrts | | Computes the spherical expansion coefficients of local gradient orientation histograms (spherical HOG) |sta_shog| [2] | starting from the input rank, it quickly computes a full orientation field based on all up and down derivatives |sta_wavelet|| sta_wavelet with a contrast normalization |sta_wavelet_inorm|| projects a scalar or vector image onto a spherical harmonic kernel|sta_steerFilt|| renders back an orientation field into a scalar image|sta_voteFilt|| diffusion time evolution |sta_propagateFC|[8]| more options than sta_propagateFC |sta_propagateAFC|[8]| spatially regularized spherical deconvolution |sta_spdeconv|| anisotropic smoothing |sta_steerdeconv| [6] | training a filter for detecting objects |sta_gfilter_train|| training a filter for detecting objects |sta_gfilter_coordinates_train|| applying a trained filter to unseen objects |sta_gfilter_apply|| ...|sta_filter_train|| ...|sta_filter_apply|| ...|||

Supplementary Functions

The toolbox also provides a set of supplementary functions for data conversion and visualization. The most important once are

description | matlab function | | ------------- | -------------| A 3D stack viewer |stackview|| An ortho-viewer |sta_LabelMe|| A glyph viewer for visualizing orientation fields |sta_glyph|| converting hardy data to orientation fields (stafields) |sta_hardi2stafield|| converting stafield to a discrete orientation field |sta_stafield2ad|| converting a discrete orientation field to a stafield |sta_ad2stafield|| computes spherical harmonics |sta_sphericalHarmonic|| computes Clebsch-Gordan coefficients |sta_clebschGordan|| extracts local maxima from an orientation field|getLocalMax|| detects local maxima in images |sta_loacalmax|| Cartesian to spherical tensor field conversion |sta_c2s|| spherical tensor to Cartesian field conversion |sta_s2c|| ...|||


#Applications

Trainable filters

Suppose you want to detect some specific landmarks within images. For this task, you can use spherical tensors to create trainable voting filters. The toolbox provides a high level interface to such filters. For instance, a harmonic filter [1] can be created with

Harmonic Filter [1]

#!matlab
model=sta_gfilter_train(...
    {img,img2},...                           %two training images
    {label_img1,label_img2},...                  %two label images  
    5,...                            %L=5
    {'gauss',2},...                      %use wavelet features 
    {'gauss',3},...                  %use a voting functions with radius 3
    'options_combo',{'o2_options_power',[1,1,5,0,5],...      %second order products
    'featurefunc',@sta_wavelet_inorm);              %contrast normalized wavelets
where img1 and img2 are two training images, and label_img1, label_img2 are two binary label images. Once you have a model, you can apply the filter to new images with

#!matlab
    H=sta_gfilter_apply(img,model);

You can easily design new filters. For instance, a SHOG filter [2], which uses spherical HOG feature, can be initialized with

SHOG Filter [2]

#!matlab
model=sta_gfilter_train(...
    {img,img2},...                           %two training images
    {label_img1,label_img2},...                  %two label images  
    5,...                            %L=5
    {'sh',[0,2],[2,2]},...                       %use two Gaussian smoothed spheres features 
    {'gauss',3},...                  %use a voting functions with radius 3
    'options_combo',{'o2_options_power',[1,1,5,0,5],...      %second order products
    'featurefunc',@sta_shog);               %SHIG features

Rotation Invariant Descriptor Fields

The angular power or bi-spectrum of orientation function is rotation invariant. Once you have rotation invariant descriptor fields, you can use any classifier to classify voxels or detect structures in images in a rotation invariant way.

The power spectrum for instance, is the squared magnitude of the spherical harmonic expansion coefficients. In order to compute the bi-spectrum, you need to combine three expansion coefficients in order to form a new, scalar valued, tensor with the spherical product. Particularly for the latter exist many symmetries. The toolbox function sta_invrts is taking care of this issue. Given a set (a cell array) of orientation fields (or a single orientation field), the function computes a four dimensional descriptor image [D, X ,Y ,Z], where the first dimension contains the (locally) rotation invariant spectra.

For instance,

#!matlab
features=sta_invtrs(ofield,'power2',true);
computes the power spectrum from the stafield ofield. For the power spectrum of the Gauss-Laguerre and Gauss-Bessel coefficients, there exist the even more memory and computation efficient functions [3]

#!matlab
features=sta_localfeat_SGD_cfast(img,params,BW,Lap);

features=sta_localfeat_SLD_cfast(img,scales,BW,Lap);
With these functions, the intermediate computation of the whole orientation field is avoided. However, note that the application of these functions is restricted to the power spectrum of those two expansions only.

Rotation Invariant Features form High Angular Resolution Images (HARDI)

The HARDI signal is itself already an orientation field. With the sta_hardi2stafield function, you can easily turn a discrete orientation field representation into a stafield. Rotation invariant spherical harmonic descriptors from the HARDI signal have been used e.g. for tissue classification [4,5]. With the sta_invrts function, you can easily create an invariant descriptor field form the HARDI data in two steps: If you extent the features with invariant descriptors from the spherical derivatives of the ofield, you can improve the results [5].

#!matlab
% converting HARDI signal to an even orientation field (default L=4)
ofield=sta_hardi2stafield(data,b_dir);
% computing the local image descriptor image using power spectrum features
d=sta_invrts(ofield, ’power2’,true);

Steerable Deconvolution

[6,8]

#!matlab
 %% line enhancement example
      image=zeros(32,32,32);
      image(16,16,[1:5 8:32])=1;
      image(16,:,16) = 1;
      kernel =  sta_fft(sta_fspecial('gauss',{'kname','gauss','kparams',1,'shape',size(image),'precision','single'}));
      data = sta_fft(stafieldStruct(single(image)));
      image = sta_ifft(sta_prod(data,kernel,0));
      image.data(1,1,:,:,:) = image.data(1,1,:,:,:) + 0.2*randn([1 1 image.shape])*max(image.data(:));
      res = sta_steerdeconv(image);

Spherical Deconvolution

[8]

#!matlab
   ds = createTestCrossing(pi/8,pi/4+pi/8,1,0.01)
           res = sta_spdeconv(ds.data,ds.dirs,ds.mask,'verbose',true,'lambda',0.01)

Further Helpful Functions

###Data Visualization

#!matlab

% rendering an orientation functions as glyphs (deformed spheres)
sta_glyph(ofield);
% start in a specific slice / specific dimension
sta_glyph(ofield,'gamma',-0.5,'init_view',3,'init_slice',14);

% with sta_glyph, you can visualize the orientation fields. For example,
% the next lines of code create the image shown below
img=zeros([32,32,16]);img(:,16,8)=1;img(8,16:end,8)=1;img(24,:,8)=1;img=exp(-bwdist(img)/2);
ofield=sta_steerFilt(stafieldStruct(img),'kernel','gaussBessel','kparam',[3,2*pi,1],'L',5,'type','STA_OFIELD_FULL');
sta_glyph(ofield,'init_view',3,'gamma',-0.5,'init_slice',8,'background',img);
colormap hot

tools

These functions are not restricted to stafields:

#!matlab

% a simple 3D image-stack viewer
stackview(img,@imagesc)

% a simple 3D slice viewer
sta_LabelMe('Image',randn([32,32,32]))

tools

###Data Conversion

#!matlab

% a function to convert HARDI data to stafields
ofield = sta_hardi2stafield(data,b_dir);


% more general : coinverting discrete orientation field 2 stafield
size(dirs) % array of directions
ans =
     3   128
size(data) % orientation field with 128 discrete directions
ans =
    32    32    16   128
ofield=sta_ad2stafield(data,dirs,'type','STA_OFIELD_EVEN','L',5);

% and backwards
data=sta_stafield2ad(ofield ,dirs);

#Literature

  1. M. Reisert and H. Burkhardt
    Harmonic filters for generic feature detection in 3D
    In Proc. of the DAGM, LNCS, Springer, 2009.

  2. H. Skibbe, M. Reisert, H. Burkhardt
    SHOG - Spherical HOG Descriptors for Rotation Invariant 3D Object Detection
    In Proc. of the DAGM, LNCS, Springer, 2011.
    SHOG filter demo

  3. H. Skibbe, M. Reisert, T. Schmidt, T. Brox, O. Ronneberger, and H. Burkhardt
    Fast rotation invariant 3d feature computation utilizing efficient local neighborhood operators
    in IEEE Trans. on PAMI, 2012.

  4. S. Schnell, D. Saur, B.W. Kreher, J. Hennig, H. Burkhardt, and V.G. Kiselev.
    Fully automated classification of HARDI in vivo data using a support vector machine
    Neuroimage, 2009.

  5. H. Skibbe, M. Reisert, H. Burkhardt
    Gaussian Neighborhood Descriptors for Brain Segmentation
    In Proc. of the MVA, 2011.

  6. M.Reisert and H. Skibbe
    Steerable deconvolution feature detection as an inverse problem
    In Proc. of the DAGM, 2011.

  7. H. Skibbe, M. Reisert, O. Ronneberger, H. Burkhardt
    Increasing the Dimension of Creativity in Rotation Invariant Feature Design Using 3D Tensorial Harmonics
    In Proc. of the DAGM, 2008

  8. M. Reisert, H. Skibbe
    Left-Invariant Diffusion on the Motion Group in terms of the Irreducible Representations of SO(3)
    In arXiv:submit/0423757 [math.AP], 2012.

#Related Projects

Updated