diff --git a/include/vigra/harmonics.hxx b/include/vigra/harmonics.hxx
new file mode 100644
index 000000000..2f4a52728
--- /dev/null
+++ b/include/vigra/harmonics.hxx
@@ -0,0 +1,1661 @@
+/************************************************************************/
+/*                                                                      */
+/*        Copyright 2009-2010 by Ullrich Koethe and Janis Fehr          */
+/*                                                                      */
+/*    This file is part of the VIGRA computer vision library.           */
+/*    The VIGRA Website is                                              */
+/*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
+/*    Please direct questions, bug reports, and contributions to        */
+/*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
+/*        vigra@informatik.uni-hamburg.de                               */
+/*                                                                      */
+/*    Permission is hereby granted, free of charge, to any person       */
+/*    obtaining a copy of this software and associated documentation    */
+/*    files (the "Software"), to deal in the Software without           */
+/*    restriction, including without limitation the rights to use,      */
+/*    copy, modify, merge, publish, distribute, sublicense, and/or      */
+/*    sell copies of the Software, and to permit persons to whom the    */
+/*    Software is furnished to do so, subject to the following          */
+/*    conditions:                                                       */
+/*                                                                      */
+/*    The above copyright notice and this permission notice shall be    */
+/*    included in all copies or substantial portions of the             */
+/*    Software.                                                         */
+/*                                                                      */
+/*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
+/*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
+/*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
+/*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
+/*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
+/*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
+/*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
+/*    OTHER DEALINGS IN THE SOFTWARE.                                   */
+/*                                                                      */
+/************************************************************************/
+
+#ifndef VIGRA_HARMONICS_HXX
+#define VIGRA_HARMONICS_HXX
+
+#include <complex>
+#include "config.hxx"
+#include "error.hxx"
+#include "utilities.hxx"
+#include "mathutil.hxx"
+#include "array_vector.hxx"
+#include "matrix.hxx"
+#include "tinyvector.hxx"
+#include "quaternion.hxx"
+#include "wigner-matrix.hxx"
+#include "clebsch-gordan.hxx"
+#include "multi_fft.hxx"
+#include "bessel.hxx"
+
+
+namespace vigra {
+
+namespace detail  {
+
+// computes the normalization for SH base functions
+inline double realSH(double l, double m)
+{
+    return std::sqrt((2.0*l + 1.0) / (4.0*M_PI*facLM(l,m)));
+
+}
+
+template <typename REAL>
+inline REAL fac(REAL in)
+{   
+    REAL temp = 1;
+    for (int i=2;i<=in;i++)
+    {   
+        temp *= i;
+    }     
+    return temp;
+}           
+
+
+template <typename REAL, typename T> 
+TinyVector<REAL, 3> centerOfBB(const MultiArray<3,T>  &A) 
+{
+    return TinyVector<REAL, 3>(A.shape()) /= 2.0;                        
+}
+
+
+template <typename REAL>
+void 
+eulerAngles(REAL sphereRadius_um,
+            MultiArray<3,REAL>& phi,
+            MultiArray<3,REAL>& theta,
+            MultiArray<3,REAL>& psi,
+            REAL gaussWidthAtHalfMaximum_um,
+            TinyVector<REAL,3> voxelSize=TinyVector<REAL,3>(1.0))
+{
+    REAL radiusLev = sphereRadius_um /voxelSize[0] + gaussWidthAtHalfMaximum_um*3;
+    REAL radiusRow = sphereRadius_um /voxelSize[1] + gaussWidthAtHalfMaximum_um*3;
+    REAL radiusCol = sphereRadius_um /voxelSize[2] + gaussWidthAtHalfMaximum_um*3;
+
+    int intRadiusLev = (int)std::ceil( radiusLev);
+    int intRadiusRow = (int)std::ceil( radiusRow);
+    int intRadiusCol = (int)std::ceil( radiusCol);
+
+    MultiArrayShape<3>::type NewShape( intRadiusLev*2 + 1, intRadiusRow*2 + 1, intRadiusCol*2 + 1 );
+
+    TinyVector<int,3> M;
+    M[0]= NewShape[0] / 2;
+    M[1]= NewShape[1] / 2;
+    M[2]= NewShape[2] / 2;
+
+    phi.reshape( NewShape,0 );
+    theta.reshape( NewShape,0 );
+    psi.reshape( NewShape,0 );
+
+
+    for (int z=0;z<NewShape[0];++z)
+    {
+        for (int y=0;y<NewShape[1];++y)
+        {
+            for (int x=0;x<NewShape[2];++x)
+            {
+                REAL X = (REAL)x *voxelSize[2];
+                REAL Y = (REAL)y *voxelSize[1];
+                REAL Z = (REAL)z *voxelSize[0];
+                int   iZ = z;
+
+                // calculate psi
+                phi(z,y,x ) = std::atan2( Y, X );
+
+                // calculate theta
+                REAL r = std::sqrt( X*X + Y*Y + Z*Z );
+                REAL alpha = (REAL)std::asin( Z / r );
+
+                if( iZ == 0 )
+                {
+                    theta( z,y,x ) = M_PI * 0.5;
+                }
+                else
+                {
+                    theta( z,y,x ) = M_PI * 0.5 - alpha;
+                }
+            }
+        }
+    }
+}
+
+} // namespace detail
+
+template <typename REAL>
+MultiArray<3,REAL>
+binarySphereREAL(REAL radius_um, 
+                 REAL gaussWidthAtHalfMaximum_um, 
+                 TinyVector<REAL,3> voxelSize=TinyVector<REAL,3>(1.0))
+{
+    REAL kernelRadius_um = radius_um;// + gaussWidthAtHalfMaximum_um*3;
+    REAL radiusLev = kernelRadius_um /voxelSize[0] + gaussWidthAtHalfMaximum_um*3 ;
+    REAL radiusRow = kernelRadius_um /voxelSize[1] + gaussWidthAtHalfMaximum_um*3;
+    REAL radiusCol = kernelRadius_um /voxelSize[2] + gaussWidthAtHalfMaximum_um*3;
+
+    int intRadiusLev = (int)std::ceil( radiusLev);
+    int intRadiusRow = (int)std::ceil( radiusRow);
+    int intRadiusCol = (int)std::ceil( radiusCol);
+
+    MultiArrayShape<3>::type outshape(intRadiusLev*2 + 1, intRadiusRow*2 + 1, intRadiusCol*2 + 1);
+    MultiArray<3,REAL> output( outshape);
+
+
+    for( int m = 0; m < outshape[0]; ++m)
+    {
+        REAL z_um = (m - intRadiusLev) *voxelSize[0];
+        REAL sqr_z_um = z_um * z_um;
+        for( int r = 0; r < outshape[1]; ++r)
+        {
+            REAL y_um = (r - intRadiusRow) *voxelSize[1];
+            REAL sqr_y_um = y_um * y_um;
+            for( int c = 0; c < outshape[2]; ++c)
+            {
+                REAL x_um = (c - intRadiusCol) *voxelSize[2];
+                REAL sqr_x_um = x_um * x_um;
+                REAL dist_um = sqrt( sqr_z_um + sqr_y_um + sqr_x_um);
+
+                if( fabs(dist_um - radius_um) <voxelSize[2]/2)
+                {
+                    output(m,r,c) = 1;
+                }
+                else
+                {
+                    output(m,r,c) = 0;
+                }
+            }
+        }
+    }
+    return output;
+}
+
+/** ------------------------------------------------------
+        sphereSurfHarmonic
+----------------------------------------------------------
+\brief computes a Spherical harmonic base function
+
+\pqrqm out returns 3D SH base function
+\param sphereRadius_um radius of the base function
+\param gaussWidthAtHalfMaximum_um gaussian smothing of the spherical surface
+\param l expansion band, l =[0,l_max]
+\param m expansion sub-band, m=[-l,l]
+\param full bool, if true, the volume of the sphere is filled, otherwise only surface function  
+\param voxelsize optional parameter used to compute base functions for non-equdistant volume samplings
+*/
+
+namespace detail {
+
+template <class REAL>
+inline
+void avoidNans(REAL & temp, REAL eps = 0.00000001)
+{
+    if (temp == 1.0)
+        temp -= eps;
+    if (temp == -1.0) 
+        temp += eps;
+}
+
+} // namespace detail
+
+template <typename REAL>
+void 
+sphereSurfHarmonic(MultiArray<3,FFTWComplex<REAL> >& output, 
+                   REAL sphereRadius_um, REAL gaussWidthAtHalfMaximum_um, 
+                   int l, int m, 
+                   bool full, 
+                   TinyVector<REAL,3> voxelSize=TinyVector<REAL,3>(1.0))
+{
+    if (gaussWidthAtHalfMaximum_um <=1)
+        gaussWidthAtHalfMaximum_um =1;
+
+    REAL radiusLev = sphereRadius_um /voxelSize[0];
+    REAL radiusRow = sphereRadius_um /voxelSize[1];
+    REAL radiusCol = sphereRadius_um /voxelSize[2];
+
+    radiusLev += gaussWidthAtHalfMaximum_um*3;
+    radiusRow += gaussWidthAtHalfMaximum_um*3;
+    radiusCol += gaussWidthAtHalfMaximum_um*3;
+
+    int intRadiusLev = (int)std::ceil( radiusLev);
+    int intRadiusRow = (int)std::ceil( radiusRow);
+    int intRadiusCol = (int)std::ceil( radiusCol);
+
+    MultiArrayShape<3>::type outshape(intRadiusLev*2 + 1,intRadiusRow*2 + 1, intRadiusCol*2 + 1);
+    output.reshape( outshape, 0);
+
+    REAL sigmaFactor = REAL(-2.0*std::log(0.5) / (gaussWidthAtHalfMaximum_um * gaussWidthAtHalfMaximum_um));
+    for( int s = 0; s < outshape[0]; ++s)
+    {
+        REAL z_um = (s - intRadiusLev) *voxelSize[0];
+        REAL sqr_z_um = z_um * z_um;
+        for( int r = 0; r < outshape[1]; ++r)
+        {
+            REAL y_um = (r - intRadiusRow) *voxelSize[1];
+            REAL sqr_y_um = y_um * y_um;
+            for( int c = 0; c < outshape[2]; ++c)
+            {
+                REAL x_um = (c - intRadiusCol) *voxelSize[2];
+                REAL sqr_x_um = x_um * x_um;
+                REAL dist_um = sqrt( sqr_z_um + sqr_y_um + sqr_x_um);
+                REAL gauss_x = 0;
+                if (!full||(dist_um > sphereRadius_um))
+                {
+                    gauss_x=(dist_um - sphereRadius_um);
+                }
+                else
+                {
+                    gauss_x=1;
+                }
+                if( x_um*x_um+y_um*y_um == 0)
+                {
+                    y_um += REAL(0.00001); //avoid nans
+                }
+                REAL theta;
+                REAL temp = z_um/sqrt( (REAL) (x_um*x_um+y_um*y_um+z_um*z_um));
+                detail::avoidNans(temp);
+                theta = std::acos(temp);
+
+                REAL phi;
+                if (y_um>=0)
+                {
+                    REAL temp = x_um/sqrt( (REAL)(x_um*x_um+y_um*y_um));
+                    detail::avoidNans(temp);
+                    phi = std::acos(temp);
+                }
+                else
+                {
+                    REAL temp = x_um/sqrt( (REAL)(x_um*x_um+y_um*y_um));
+                    detail::avoidNans(temp);
+                    phi = REAL(2.0*M_PI - std::acos(temp));
+                }
+
+                FFTWComplex<REAL> SHfactor;
+                SHfactor.real()= REAL(detail::realSH(l,m)*legendre(l,m,std::cos(theta)) * std::cos(m * phi));
+                SHfactor.imag()= REAL(detail::realSH(l,m)*legendre(l,m,std::cos(theta)) * std::sin(m * phi));
+                output(s,r,c) = ((REAL)std::exp( -0.5 * gauss_x * gauss_x * sigmaFactor)) * SHfactor;
+            }
+        }
+    }
+}
+
+template <typename REAL>
+void 
+sphereVecHarmonic(MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > >  & res,
+                  REAL radius, REAL gauss, 
+                  int l, int k, int m)
+{
+    //std::cerr<<"computing VH: "<<l<<" "<<k<<" "<<m<<"\n";
+    //1-m
+    MultiArray<3,FFTWComplex<REAL> > tmpSH;
+
+    InvariantViolation err("");
+    try
+    {
+        if (abs(1-m)>l) 
+            throw err;
+        FFTWComplex<REAL> cg;
+        cg.real() = clebschGordan(l+k, m, l, 1-m, 1, 1);
+        cg.imag() = 0;
+        sphereSurfHarmonic(tmpSH,radius, gauss, l, 1-m, false);
+        typename MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > >::iterator p = res.begin();
+        typename MultiArray<3,FFTWComplex<REAL> >::iterator q = tmpSH.begin();
+        for (;q!=tmpSH.end();++p,++q)
+            (*p)[0] = cg * *q;
+    }
+    catch(InvariantViolation &err) //in case clebsh gordan dilivers invalid combination
+    {
+        //std::cerr<<"no Z\n";
+        typename MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > >::iterator p = res.begin();
+        for (; p!=res.end();++p)
+            (*p)[0] = (REAL)0;
+    }
+    //-m
+    try
+    {
+        if (abs(-m)>l) 
+            throw err;
+        FFTWComplex<REAL> cg;
+        cg.real() = clebschGordan(l+k, m, l, -m, 1, 0);
+        sphereSurfHarmonic(tmpSH,radius, gauss, l, 1-m, false);
+        typename MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > >::iterator p = res.begin();
+        typename MultiArray<3,FFTWComplex<REAL> >::iterator q = tmpSH.begin();
+        for (;q!=tmpSH.end();++p,++q)
+            (*p)[1] = cg * *q;
+    }
+    catch(InvariantViolation &err) //in case clebsh gordan dilivers invalid combination
+    {
+        //std::cerr<<"no Y\n";
+        typename MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > >::iterator p = res.begin();
+        for (; p!=res.end();++p)
+            (*p)[1] = (REAL)0;
+
+    }
+    //-(1+m)
+    try
+    {
+        if (abs(-(1+m))>l) 
+            throw err;
+        FFTWComplex<REAL> cg;
+        cg.real() = clebschGordan(l+k, m, l, -(1+m), 1, -1);
+        cg.imag() = 0;
+        sphereSurfHarmonic(tmpSH,radius, gauss, l, -(m+1), false);
+        typename MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > >::iterator p = res.begin();
+        typename MultiArray<3,FFTWComplex<REAL> >::iterator q = tmpSH.begin();
+        for (;q!=tmpSH.end();++p,++q)
+            (*p)[0] = cg * *q;
+    }
+    catch(InvariantViolation &err) //in case clebsh gordan dilivers invalid combination
+    {
+        //std::cerr<<"no X\n";
+        typename MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > >::iterator p = res.begin();
+        for (; p!=res.end();++p)
+            (*p)[2] = (REAL)0;
+    }
+}
+
+template <typename REAL>
+inline REAL bessel_zero_Jnu(unsigned int l, unsigned int n)
+{
+    vigra_invariant(l <= 10 && n <= 10,
+        "bessel_zero_Jnu(): max implemented band is 10.");
+
+    //cache of the first 10 bessel zero points (max l=10)
+    static REAL cache[110] = {
+        0,0,0,0,0,0,0,0,0,0, 
+        REAL(2.4048255576957729), REAL(3.8317059702075125), 
+        REAL(5.1356223018406828), REAL(6.3801618959239841), 
+        REAL(7.5883424345038035), REAL(8.7714838159599537), 
+        REAL(9.9361095242176845), REAL(11.086370019245084), 
+        REAL(12.225092264004656), REAL(13.354300477435331), 
+        REAL(5.5200781102863106), REAL(7.0155866698156188), 
+        REAL(8.4172441403998643), REAL(9.7610231299816697), 
+        REAL(11.064709488501185), REAL(12.338604197466944), 
+        REAL(13.589290170541217), REAL(14.821268727013171), 
+        REAL(16.03777419088771), REAL(17.241220382489129), 
+        REAL(8.6537279129110125), REAL(10.173468135062722), 
+        REAL(11.61984117214906), REAL(13.015200721698434), 
+        REAL(14.37253667161759), REAL(15.700174079711671), 
+        REAL(17.003819667816014), REAL(18.287582832481728), 
+        REAL(19.554536430997054), REAL(20.807047789264107), 
+        REAL(11.791534439014281), REAL(13.323691936314223), 
+        REAL(14.795951782351262), REAL(16.223466160318768), 
+        REAL(17.615966049804832), REAL(18.98013387517992), 
+        REAL(20.320789213566506), REAL(21.6415410198484), 
+        REAL(22.945173131874618), REAL(24.233885257750551), 
+        REAL(14.930917708487787), REAL(16.470630050877634), 
+        REAL(17.959819494987826), REAL(19.409415226435012), 
+        REAL(20.826932956962388), REAL(22.217799896561267), 
+        REAL(23.586084435581391), REAL(24.934927887673023), 
+        REAL(26.266814641176644), REAL(27.583748963573008), 
+        REAL(18.071063967910924), REAL(19.615858510468243), 
+        REAL(21.116997053021844), REAL(22.582729593104443), 
+        REAL(24.01901952477111), REAL(25.430341154222702), 
+        REAL(26.820151983411403), REAL(28.1911884594832), 
+        REAL(29.54565967099855), REAL(30.885378967696674), 
+        REAL(21.211636629879258), REAL(22.760084380592772), 
+        REAL(24.270112313573105), REAL(25.748166699294977), 
+        REAL(27.19908776598125), REAL(28.626618307291139), 
+        REAL(30.033722386570467), REAL(31.422794192265581), 
+        REAL(32.795800037341465), REAL(34.154377923855094), 
+        REAL(24.352471530749302), REAL(25.903672087618382), 
+        REAL(27.420573549984557), REAL(28.908350780921758), 
+        REAL(30.371007667117247), REAL(31.811716724047763), 
+        REAL(33.233041762847122), REAL(34.637089352069324), 
+        REAL(36.025615063869573), REAL(37.400099977156586), 
+        REAL(27.493479132040257), REAL(29.046828534916855), 
+        REAL(30.569204495516395), REAL(32.06485240709771), 
+        REAL(33.53713771181922), REAL(34.988781294559296), 
+        REAL(36.422019668258457), REAL(37.838717382853609), 
+        REAL(39.240447995178137), REAL(40.628553718964525), 
+        REAL(30.634606468431976), REAL(32.189679910974405), 
+        REAL(33.716519509222699), REAL(35.218670738610115), 
+        REAL(36.699001128744648), REAL(38.15986856196713), 
+        REAL(39.603239416075404), REAL(41.030773691585537), 
+        REAL(42.443887743273557), REAL(43.84380142033735) 
+    };
+    return cache[n*10+l];
+}
+
+template <typename REAL>
+void 
+sphereFullHarmonic(MultiArray<3,FFTWComplex<REAL> >& output, 
+                   REAL sphereRadius_um, 
+                   int n, int l, int m, 
+                   TinyVector<REAL,3> voxelSize=TinyVector<REAL,3>(1.0))
+{
+    REAL radiusLev = sphereRadius_um / voxelSize[0] +3;
+    REAL radiusRow = sphereRadius_um / voxelSize[1] +3;
+    REAL radiusCol = sphereRadius_um / voxelSize[2] +3;
+    int intRadiusLev = (int)std::ceil( radiusLev);
+    int intRadiusRow = (int)std::ceil( radiusRow);
+    int intRadiusCol = (int)std::ceil( radiusCol);
+
+
+    //precompute SH parts for l and m
+    MultiArray<3,FFTWComplex<REAL> > SH;
+    sphereSurfHarmonic<REAL>(SH,sphereRadius_um, 1, l, m, true);
+
+    MultiArrayShape<3>::type outshape(intRadiusLev*2 + 1,intRadiusRow*2 + 1, intRadiusCol*2 + 1);
+    output.reshape( outshape, 0);
+
+    REAL xnl=sphereRadius_um;
+    if(n>0)
+    {
+        xnl=bessel_zero_Jnu<REAL>(l,n);
+    }
+    REAL k=xnl/sphereRadius_um;
+    REAL J2=(REAL)std::pow(besselJ(l+1,xnl),2.0);
+    REAL N=(sphereRadius_um*sphereRadius_um*sphereRadius_um)/2*J2;
+
+    REAL sigmaFactor = REAL(-2.0*std::log(0.5) / 4.0);
+    FFTWComplex<REAL> I(0,1);
+    for( int z = 0; z < outshape[0]; ++z)
+    {
+        REAL z_um = (z - intRadiusLev) * voxelSize[0];
+        REAL sqr_z_um = z_um * z_um;
+        for( int y = 0; y < outshape[1]; ++y)
+        {
+            REAL y_um = (y - intRadiusRow) * voxelSize[1];
+            REAL sqr_y_um = y_um * y_um;
+            for( int x = 0; x < outshape[2]; ++x)
+            {
+                REAL x_um = (x - intRadiusCol) * voxelSize[2];
+                REAL sqr_x_um = x_um * x_um;
+                REAL r = sqrt( sqr_z_um+sqr_y_um + sqr_x_um);
+                REAL gauss_x = (r - sphereRadius_um);
+
+                FFTWComplex<REAL> Phi = SH(z,y,x);
+                REAL J1=(REAL)besselJ(l,k*r);
+                FFTWComplex<REAL> R = 1/sqrt(N)*J1;
+                FFTWComplex<REAL> Psi = R*Phi;
+
+                if(r<=sphereRadius_um)
+                    output(z,y,x) = Psi;
+                else
+                    output(z,y,x) = (REAL)exp( -0.5 * gauss_x * gauss_x * sigmaFactor) * Psi;
+            }
+        }
+    }
+}
+
+template <typename REAL>
+void 
+sphereFullVecHarmonic(MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > >& output, 
+                      REAL sphereRadius_um, 
+                      int n, int l, int k, int m, 
+                      TinyVector<REAL,3> voxelSize=TinyVector<REAL,3>(1.0))
+{
+    REAL radiusLev = sphereRadius_um / voxelSize[0] +3;
+    REAL radiusRow = sphereRadius_um / voxelSize[1] +3;
+    REAL radiusCol = sphereRadius_um / voxelSize[2] +3;
+    int intRadiusLev = (int)std::ceil( radiusLev);
+    int intRadiusRow = (int)std::ceil( radiusRow);
+    int intRadiusCol = (int)std::ceil( radiusCol);
+
+
+    //precompute VH parts for l and m
+    MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > > VH;
+    sphereVecHarmonic(VH, sphereRadius_um, 1, l, k, m);
+
+    MultiArrayShape<3>::type outshape(intRadiusLev*2 + 1,intRadiusRow*2 + 1, intRadiusCol*2 + 1);
+    TinyVector<FFTWComplex<REAL>,3 > zero(0,0,0);
+    output.reshape( outshape, zero);
+
+    REAL xnl=sphereRadius_um;
+    if(n>0)
+    {
+            xnl=0;//gsl_sf_bessel_zero_Jnu(l,n);
+            //std::cerr<<n<<" "<<m<<" "<<xnm<<"\n";
+    }
+    REAL K=xnl/sphereRadius_um;
+    REAL J2=std::pow(besselJ(l+1,xnl),2.0);
+    REAL N=(sphereRadius_um*sphereRadius_um*sphereRadius_um)/2*J2;
+
+    REAL sigmaFactor = -2*log(0.5) / (4);
+    FFTWComplex<REAL> I(0,1);
+    for( int z = 0; z < outshape[0]; ++z)
+    {
+        REAL z_um = (z - intRadiusLev) * voxelSize[0];
+        REAL sqr_z_um = z_um * z_um;
+        for( int y = 0; y < outshape[1]; ++y)
+        {
+            REAL y_um = (y - intRadiusRow) * voxelSize[1];
+            REAL sqr_y_um = y_um * y_um;
+            for( int x = 0; x < outshape[2]; ++x)
+            {
+                REAL x_um = (x - intRadiusCol) * voxelSize[2];
+                REAL sqr_x_um = x_um * x_um;
+                REAL r = sqrt( sqr_z_um+sqr_y_um + sqr_x_um);
+                REAL gauss_x = (r - sphereRadius_um);
+                TinyVector<FFTWComplex<REAL>,3> Phi = VH(z,y,x);
+                REAL J1=besselJ(l,K*r);
+                FFTWComplex<REAL> R = 1/sqrt(N)*J1;
+                TinyVector<FFTWComplex<REAL>,3> Psi;
+                Psi[0]= R*Phi[0];
+                Psi[1]= R*Phi[1];
+                Psi[2]= R*Phi[2];
+
+                if(r<=sphereRadius_um)
+                {
+                    output(z,y,x) = Psi;
+                }
+                else
+                {
+                    Psi[0]*=((FFTWComplex<REAL>)exp( -0.5 * gauss_x * gauss_x * sigmaFactor));
+                    Psi[1]*=((FFTWComplex<REAL>)exp( -0.5 * gauss_x * gauss_x * sigmaFactor));
+                    Psi[2]*=((FFTWComplex<REAL>)exp( -0.5 * gauss_x * gauss_x * sigmaFactor));
+                    output(z,y,x) =  Psi;
+                }
+            }
+        }
+    }
+}
+
+/** ------------------------------------------------------
+        computeSHbaseF
+----------------------------------------------------------
+\brief pre-computes Spherical harmonic base functions
+
+\param radius radius of the spherical expansion
+\param gauss smothing of the spherical surface
+\param band maximum expansion band
+\param SHbaseF holds precomputed SH base functions  
+*/
+template <typename REAL>
+void 
+computeSHbaseF(REAL radius, REAL gauss, unsigned int band, 
+               std::vector<std::vector<MultiArray<3,FFTWComplex<REAL> > > > &SHbaseF)
+{
+    SHbaseF.resize(band + 1);
+
+    for (int l = 0; l <= (int)band; l++)
+    {
+        SHbaseF[l].resize(2*l + 1);
+        for (int m = -l; m <= l; m++)
+        {
+            MultiArray<3, FFTWComplex<REAL> > Coeff;
+            sphereSurfHarmonic(Coeff,radius, gauss, l, m, false);
+            SHbaseF[l][m+l].reshape(Coeff.shape());
+            SHbaseF[l][m+l]=Coeff;
+        }
+    }
+}
+
+template <typename REAL>
+void 
+computePHbaseF(REAL radius, unsigned int band, 
+               std::vector<std::vector<std::vector<MultiArray<3,FFTWComplex<REAL> > > > > &PHbaseF, 
+               bool realdata=false)
+{
+    PHbaseF.resize(band + 1);
+    //n=0 undefined
+    int mMin;
+    int mOff;
+    //for real data one only needs posstive coeefs (symmetry)
+    for (int n = 1; n <= (int)band; n++)
+    {
+        PHbaseF[n].resize(band + 1);
+        for (int l = 0; l <= (int)band; l++)
+        {
+            if (realdata)
+            {
+                mMin = 0;
+                mOff = 0;
+                PHbaseF[n][l].resize(l + 1);
+            }
+            else
+            {
+                mMin = -l;
+                mOff = l;
+                PHbaseF[n][l].resize(2*l + 1);
+            }
+            for (int m = mMin; m <= l; m++)
+            {
+                //std::cerr<<n<<" "<<l<<" "<<m+mOff<<"\n";
+                MultiArray< 3,FFTWComplex<REAL> > Coeff;
+                sphereFullHarmonic(Coeff, radius, n,l,m);
+                PHbaseF[n][l][m+mOff].reshape(Coeff.shape());
+                PHbaseF[n][l][m+mOff]=Coeff;
+            }
+        }
+    }
+}
+
+template <typename REAL>
+void 
+computeVHbaseF(REAL radius, REAL gauss, unsigned int band, 
+               std::vector<std::vector<std::vector<MultiArray<3, TinyVector<FFTWComplex<REAL>,3> > > > > &VHbaseF)
+{
+    VHbaseF.resize(band + 1);
+
+    for (int l = 0; l <= band; l++)
+    {
+        VHbaseF[l].resize(3);
+        for (int k=-1;k<=1;k++)
+        {
+            VHbaseF[l][k+1].resize(2*band + 1);
+            for (int m = -l; m <= l; m++)
+            {
+                MultiArray<3, TinyVector< FFTWComplex<REAL> ,3 > > Coeff;
+
+                sphereVecHarmonic(Coeff,radius, gauss, l, k, m);
+                VHbaseF[l][k+1][m+l].reshape(Coeff.shape());
+                VHbaseF[l][k+1][m+l]=Coeff;
+            }
+        }
+    }
+}
+
+template <typename REAL>
+void 
+computeVPHbaseF(REAL radius, unsigned int band, 
+                std::vector<std::vector<std::vector<std::vector<MultiArray<3, TinyVector< FFTWComplex<REAL>,3 > > > > > >&VHbaseF)
+{
+    VHbaseF.resize(band + 1);
+    for(int n=0; n>= band; n++)
+    {
+        VHbaseF[n].resize(band + 1);
+
+        for (int l = 0; l <= band; l++)
+        {
+            VHbaseF[n][l].resize(3);
+            for (int k=-1;k<=1;k++)
+            {
+                VHbaseF[l][k+1].resize(2*band + 1);
+                for (int m = -l; m <= l; m++)
+                {
+                    MultiArray<3, TinyVector< FFTWComplex<REAL> ,3 > > Coeff;
+
+                    sphereFullVecHarmonic(Coeff,radius, n, l, k, m);
+                    VHbaseF[n][l][k+1][m+l].reshape(Coeff.shape());
+                    VHbaseF[n][l][k+1][m+l]=Coeff;
+                }
+            }
+        }
+    }
+}
+
+
+template <typename REAL>
+MultiArray<3,REAL> 
+sphereSurfGauss(REAL sphereRadius_um, REAL gaussWidthAtHalfMaximum_um, 
+                TinyVector<REAL,3> voxelSize=TinyVector<REAL,3>(1.0))
+{
+    REAL kernelRadius_um = sphereRadius_um;;
+    REAL radiusLev = kernelRadius_um /voxelSize[0] + gaussWidthAtHalfMaximum_um*3;
+    REAL radiusRow = kernelRadius_um /voxelSize[1] + gaussWidthAtHalfMaximum_um*3;
+    REAL radiusCol = kernelRadius_um /voxelSize[2] + gaussWidthAtHalfMaximum_um*3;
+
+    int intRadiusLev = (int)std::ceil( radiusLev);
+    int intRadiusRow = (int)std::ceil( radiusRow);
+    int intRadiusCol = (int)std::ceil( radiusCol);
+
+    MultiArrayShape<3>::type outShape(intRadiusLev*2 + 1, intRadiusRow*2 + 1, intRadiusCol*2 + 1);
+    MultiArray<3,REAL> output( outShape );
+
+    REAL sigmaFactor = -2*log(0.5) / (gaussWidthAtHalfMaximum_um * gaussWidthAtHalfMaximum_um);
+
+    for( int m = 0; m < outShape[0]; ++m)
+    {
+        REAL z_um = (m - intRadiusLev) *voxelSize[0];
+        REAL sqr_z_um = z_um * z_um;
+        for( int r = 0; r < outShape[1]; ++r)
+        {
+            REAL y_um = (r - intRadiusRow) *voxelSize[1];
+            REAL sqr_y_um = y_um * y_um;
+            for( int c = 0; c < outShape[2]; ++c)
+            {
+                REAL x_um = (c - intRadiusCol) *voxelSize[2];
+                REAL sqr_x_um = x_um * x_um;
+                REAL dist_um = sqrt( sqr_z_um + sqr_y_um + sqr_x_um);
+                REAL gauss_x = (dist_um - sphereRadius_um);
+                output(m,r,c) = exp( -0.5 * gauss_x * gauss_x * sigmaFactor);
+            }
+        }
+    }
+
+    REAL kernelSum = 0;
+    typename MultiArray<3,REAL>::iterator p=output.begin();
+    for (;p!=output.end();++p)
+        kernelSum+=*p;
+    output *= 1.0/kernelSum;
+
+    return output;
+}
+
+template <typename REAL>
+void 
+reconstSH(REAL radius, REAL gauss, unsigned int band, 
+          MultiArray<3,REAL >& reconstruct, 
+          const std::vector<std::vector<FFTWComplex<REAL> >  >& SH_A, 
+          const std::vector<std::vector<MultiArray<3,FFTWComplex<REAL> > > >& SHbaseF)
+{
+    //FIXME reconstSH currently only works for real data
+    reconstruct.reshape(SHbaseF[0][0].shape());
+
+    MultiArray<3,FFTWComplex<REAL> > T;
+    std::string name;
+    for (int l=0;l<=(int)band;l++)
+    {
+        for (int m=-l;m<=l;m++)
+        {
+            typename MultiArray<3,FFTWComplex<REAL> >::iterator p=SHbaseF[l][l+m].begin();
+            typename MultiArray<3,REAL >::iterator q=reconstruct.begin();
+            for (;q!=reconstruct.end();++p,++q)
+                *q += real(*p * conj(SH_A[l][l+m]));
+        }
+    }
+}
+
+template <typename REAL>
+void 
+reconstPH(REAL radius, unsigned int band, 
+          MultiArray<3, REAL >& reconstruct, 
+          const std::vector<std::vector<std::vector<FFTWComplex<REAL> >  > >& PH_A, 
+          std::vector<std::vector<std::vector<MultiArray<3,FFTWComplex<REAL> > > > > &PHbaseF)
+{
+    //FIXME reconstPH currently only works for real data
+    reconstruct.reshape(PHbaseF[1][0][0].shape(), 0);
+
+    for (int n=1;n<(int)PH_A.size();n++)
+    {
+        for (int l=0;l<(int)PH_A[n].size();l++)
+        {
+            for (int m=0;m<(int)PH_A[n][l].size();m++)
+            {
+                typename MultiArray<3,FFTWComplex<REAL> >::iterator p = PHbaseF[n][l][m].begin();
+                typename MultiArray<3, REAL >::iterator q=reconstruct.begin();
+                for (;q!=reconstruct.end();++q,++p)
+                    *q += real(*p * conj(PH_A[n][l][m]));
+            }
+        }
+    }
+}
+
+template <typename REAL>
+void 
+reconstVPH(REAL radius, unsigned int band, 
+           MultiArray<3, TinyVector<REAL,3> >& reconstruct, 
+           const std::vector<std::vector<std::vector<std::vector<FFTWComplex<REAL> > > > >& VPH_A)
+{
+    MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > > base_tmp;
+    sphereVecHarmonic(base_tmp, radius, 1, 0, 0, 0);
+    reconstruct.reshape(base_tmp.shape());
+    MultiArray<3, TinyVector<FFTWComplex<REAL>,3> > tmp(reconstruct.shape(), 0);
+    FFTWComplex<REAL> zero(0,0);
+    TinyVector<FFTWComplex<REAL>,3> Zero(zero,zero,zero);
+
+    for (int n=1;n<=(int)band;n++)
+    {
+        for (int l=0;l<=(int)band;l++)
+        {
+            for (int k=-1;k<=1;++k)
+            {
+                for (int m=-(l+k);m<=(l+k);m++)
+                {
+                    sphereVecHarmonic(base_tmp, radius, n, l, k, m);
+                    typename MultiArray<3, TinyVector<FFTWComplex<REAL>,3> >::iterator p=tmp.begin();
+                    typename MultiArray< 3, TinyVector<FFTWComplex<REAL>,3 > >::iterator q=base_tmp.begin();
+                    for (;q!=base_tmp.end();++p,++q)
+                    {
+                        (*p)[0] += conj(VPH_A[n][l][k+1][(l+k)+m]) * (*q)[0];
+                        (*p)[1] += conj(VPH_A[n][l][k+1][(l+k)+m]) * (*q)[1];
+                        (*p)[2] += conj(VPH_A[n][l][k+1][(l+k)+m]) * (*q)[2];
+                    }
+                }
+            }
+        }
+    }
+    //reconst real vec directions
+    typename MultiArray<3, TinyVector<REAL,3> >::iterator p=reconstruct.begin();
+    typename MultiArray<3, TinyVector<FFTWComplex<REAL>,3> >::iterator q=tmp.begin();
+    for (;q!=tmp.end();++p,++q)
+    {
+        (*p)[0] = real((*q)[1]);
+        (*p)[1] = -sqrt(2)*real((*q)[0]);
+        (*p)[2] = sqrt(2)*imag((*q)[0]);
+    }
+}
+
+template <typename REAL>
+void 
+reconstVH(REAL radius, REAL gauss, unsigned int band, 
+          MultiArray<3, TinyVector<REAL,3> >& reconstruct,
+          const std::vector<std::vector<std::vector<FFTWComplex<REAL> > > >& VH_A)
+{
+
+    MultiArray< 3, TinyVector<FFTWComplex<REAL>,3 > > base_tmp;
+    sphereVecHarmonic(base_tmp, radius, gauss, 0, 0, 0);
+    reconstruct.reshape(base_tmp.shape());
+    FFTWComplex<REAL> zero(0,0);
+    MultiArray<3, TinyVector<FFTWComplex<REAL>,3> > tmp(reconstruct.shape(), 0);
+    TinyVector<FFTWComplex<REAL>,3> Zero(zero,zero,zero);
+
+    for (int l=0;l<=band;l++)
+    {
+        for (int k=-1;k<=1;++k)
+        {
+            for (int m=-(l+k);m<=(l+k);m++)
+            {
+                sphereVecHarmonic(base_tmp, radius, gauss, l, k, m);
+                typename MultiArray<3, TinyVector<FFTWComplex<REAL>,3> >::iterator p=tmp.begin();
+                typename MultiArray< 3, TinyVector<FFTWComplex<REAL>,3 > >::iterator q=base_tmp.begin();
+                for (;q!=base_tmp.end();++p,++q)
+                {
+                    (*p)[0] += conj(VH_A[l][k+1][(l+k)+m]) * (*q)[0];
+                    (*p)[1] += conj(VH_A[l][k+1][(l+k)+m]) * (*q)[1];
+                    (*p)[2] += conj(VH_A[l][k+1][(l+k)+m]) * (*q)[2];
+                }
+            }
+        }
+    }
+
+    //reconst real vec directions
+    typename MultiArray<3, TinyVector<REAL,3> >::iterator p=reconstruct.begin();
+    typename MultiArray<3, TinyVector<FFTWComplex<REAL>,3> >::iterator q=tmp.begin();
+    for (;q!=tmp.end();++p,++q)
+    {
+        (*p)[0] = real((*q)[1]);
+        (*p)[1] = -sqrt(2)*real((*q)[0]);
+        (*p)[2] = sqrt(2)*imag((*q)[0]);
+    }
+}
+
+/** ------------------------------------------------------
+    SHpos
+----------------------------------------------------------
+\brief computes singe local Spherical harmonic expansion at given 3D position
+
+\param SH_A holds the returned SH coefficients  
+\param radius radius of the spherical expansion
+\param band maximum expansion band
+\param A input volume data
+\param pos 3D position of the SH expansion
+\param SHbaseF precomputed SH base functions (-> see computeSHbaseF) 
+*/
+template <typename REAL>
+void 
+SHpos(std::vector<std::vector<FFTWComplex<REAL> >  > &SH_A, 
+      REAL radius, unsigned int band,
+      const MultiArray<3,REAL> &A, 
+      const TinyVector<REAL, 3> &pos, 
+      std::vector<std::vector<MultiArray<3,FFTWComplex<REAL> > > > const &SHbaseF)
+{
+    SH_A.resize(SHbaseF.size());
+
+    // SH coeffs
+    for (int l = 0; l <= (int)band; l++)
+    {
+        SH_A[l].resize(SHbaseF[l].size());
+        for (int m = 0; m < (int)SHbaseF[l].size(); m++)
+        {
+            MultiArrayShape<3>::type coffShape(SHbaseF[l][m].shape());
+
+            int xa = (int) floor(pos[2] - coffShape[2] / 2);
+            int xe = xa + coffShape[2] - 1;
+            int ya = (int) floor(pos[1] - coffShape[1] / 2);
+            int ye = ya + coffShape[1] - 1;
+            int za = (int) floor(pos[0] - coffShape[0] / 2);
+            int ze = za + coffShape[0] - 1;
+
+            SH_A[l][m] = (REAL)0;
+            int sz=0;
+            for (int z=ze;z>=za;z--,sz++)
+            {
+                int sy=0;
+                for (int y=ye;y>=ya;y--,sy++)
+                {
+                    int sx=0;
+                    for (int x=xe;x>=xa;x--,sx++)
+                    SH_A[l][m]+=A(z,y,x)*(SHbaseF[l][m](sz,sy,sx));
+                }
+            }
+        }
+    }
+}
+
+/** ------------------------------------------------------
+        SHcenter
+----------------------------------------------------------
+\brief computes singe local Spherical harmonic expansion at the center of te given 3D volume
+
+\param SH_A holds the returned SH coefficients  
+\param radius radius of the spherical expansion
+\param band maximum expansion band
+\param A input volume data
+\param SHbaseF precomputed SH base functions (-> see computeSHbaseF) 
+*/
+
+template <typename REAL>
+void 
+SHcenter(std::vector<std::vector<FFTWComplex<REAL> >  > &SH_A, 
+         REAL radius, unsigned int band,
+         const MultiArray<3,REAL> &A, 
+         std::vector<std::vector<MultiArray<3,FFTWComplex<REAL> > > > const &SHbaseF)
+{
+    SHpos(SH_A, radius, band, A,  detail::centerOfBB<REAL>(A), SHbaseF);
+}
+
+/** ------------------------------------------------------
+        SH_Iterator
+---------------------------------------------------------- 
+\brief iterator class to access the cascaded vector representations of SH base functions and SH coefficients
+*/
+template <typename T>
+class SH_Iterator
+{
+  public:
+    typedef MultiArray<3,T >* pointer;
+    typedef MultiArray<3,T >& reference;
+    typedef MultiArray<3,T > value_type;
+    typedef int difference_type;
+    typedef std::random_access_iterator_tag iterator_category;
+
+    SH_Iterator(std::vector<std::vector<MultiArray<3,T > > >& data, 
+                 int l, int m, std::string name="name")
+    :   _data(data),
+        _l(l),
+        _m(m),
+        _name(name)
+    {}
+    
+    void operator++()
+    {
+        if (_m< (int)_data[_l].size()-1)
+        {
+            //std::cerr<<_name<<"++ "<<_m<<" "<<_data[_l].size()<<"\n"<<std::flush;
+            _m++;
+        }
+        else
+        {
+            if (_l< (int)_data.size())
+            {
+                _m=0;
+                _l++;
+            }
+        }
+    }
+    
+    MultiArray<3,T >& operator*() const
+    {
+        return _data[_l][_m];
+    }
+
+    const MultiArray<3,T >* operator->() const
+    {
+        //std::cerr<<_name<<" "<<_l<<" "<<_m<<" "<<_data[_l].size()<<" "<<_data[_l][_m].shape()<<"\n";
+        return &_data[_l][_m];
+    }
+
+    int getL()
+    {
+        return _l;
+    }
+    
+    int getM()
+    {
+        return _m;
+    }
+    
+    bool operator==(SH_Iterator& A)
+    {
+        return ((this->_l==A.getL())&&(this->_m==A.getM()));
+    }
+    
+    bool operator!=(SH_Iterator& A)
+    {
+        return !((this->_l==A.getL())&&(this->_m==A.getM()));
+    }
+
+  private:
+    std::vector<std::vector<MultiArray<3,T > > > &_data;
+    int _l;
+    int _m;
+    std::string _name;
+};
+
+template <typename REAL>
+class PH_Iterator
+{
+  public:
+    typedef MultiArray<3,FFTWComplex<REAL> >* pointer;
+    typedef MultiArray<3,FFTWComplex<REAL> >& reference;
+    typedef MultiArray<3,FFTWComplex<REAL> > value_type;
+    typedef int difference_type;
+    typedef std::random_access_iterator_tag iterator_category;
+
+    PH_Iterator(std::vector<std::vector<std::vector<MultiArray<3,FFTWComplex<REAL> > > > >& data, 
+                int k, int l, int m)
+    : _data(data),
+      _k(k),
+      _l(l),
+      _m(m)
+    {}
+    
+    void operator++()
+    {
+        if (_m<_data[_k][_l].size()-1)
+        {
+            _m++;
+        }
+        else
+        {
+            if (_l<_data[_k].size()-1)
+            {
+                _m=0;
+                _l++;
+            }
+            else
+            {
+                if (_k<_data.size())
+                {
+                    _l=0;
+                    _k++;
+                }
+            }
+        }
+    }
+    
+    MultiArray<3,FFTWComplex<REAL> >& operator*() const
+    {
+        return _data[_k][_l][_m];
+    }
+
+    const MultiArray<3,FFTWComplex<REAL> >* operator->() const
+    {
+        return &_data[_k][_l][_m];
+    }
+    
+    int getK()
+    {
+        return _k;
+    }
+    
+    int getL()
+    {
+        return _l;
+    }
+    
+    int getM()
+    {
+        return _m;
+    }
+    
+    bool operator==(PH_Iterator& A)
+    {
+        return ((this->_k==A.getL())&&(this->_l==A.getL())&&(this->_m==A.getM()));
+    }
+    
+    bool operator!=(PH_Iterator& A)
+    {
+        return !((this->_k==A.getK())&&(this->_l==A.getL())&&(this->_m==A.getM()));
+    }
+
+  private:
+    std::vector<std::vector<std::vector<MultiArray<3,FFTWComplex<REAL> > > > > &_data;
+    int _k;
+    int _l;
+    int _m;
+};
+
+/** ------------------------------------------------------
+        Array2SH
+----------------------------------------------------------
+\brief computes  local Spherical harmonic expansion at all positions of the given 3D volume
+
+\param SH_A holds the returned SH coefficients  
+\param radius radius of the spherical expansion
+\param band maximum expansion band
+\param A input volume data
+\param SHbaseF precomputed SH base functions (-> see computeSHbaseF) 
+*/
+template <typename REAL>
+void 
+Array2SH(std::vector<std::vector<MultiArray<3,FFTWComplex<REAL> >  > > &SH_A, 
+         unsigned int band, REAL radius, 
+         std::vector<std::vector<MultiArray<3,FFTWComplex<REAL> > > > &SHbaseF,
+         const MultiArray<3,REAL> &A)
+//FIXME add const iterator and make SHbaseF const
+{
+
+    SH_A.resize(SHbaseF.size());
+    for (int l=0;l<(int)SHbaseF.size();l++)
+    {
+        SH_A[l].resize(SHbaseF[l].size());
+        for (int m=0;m<(int)SHbaseF[l].size();m++)
+        {
+            SH_A[l][m].reshape(A.shape(),0);
+        }
+    }
+
+    SH_Iterator< FFTWComplex<REAL> > SHbaseF_Iter(SHbaseF,0,0,"BaseF");
+    SH_Iterator< FFTWComplex<REAL> > SHbaseF_Iter_end(SHbaseF,SHbaseF.size(),0);
+    SH_Iterator< FFTWComplex<REAL> > SH_A_Iter(SH_A,0,0,"SH_A");
+    convolveFFTComplexMany(A, SHbaseF_Iter, SHbaseF_Iter_end, SH_A_Iter,false);
+}
+
+template <typename REAL>
+void 
+PHpos(std::vector<std::vector<std::vector<FFTWComplex<REAL> >  > > &PH_A, 
+      REAL radius, unsigned int band,
+      const MultiArray<3,FFTWComplex<REAL> > &A, 
+      std::vector<std::vector<std::vector<MultiArray<3,FFTWComplex<REAL> > > > > &PHbaseF,
+      const TinyVector<REAL, 3> &pos)
+{
+    PH_A.resize(PHbaseF.size());
+
+    for (int n=1;n<(int)PHbaseF.size();n++)
+    {
+        PH_A[n].resize(PHbaseF[n].size());
+        for (int l = 0; l < (int)PHbaseF[n].size(); l++)
+        {
+            PH_A[n][l].resize(PHbaseF[n][l].size());
+            for (int m = 0; m < (int)PHbaseF[n][l].size(); m++)
+            {
+                MultiArrayShape<3>::type coffShape(PHbaseF[n][l][m].shape());
+                int xa = (int) floor(pos[2] - coffShape[2] / 2);
+                int xe = xa + coffShape[2] - 1;
+                int ya = (int) floor(pos[1] - coffShape[1] / 2);
+                int ye = ya + coffShape[1] - 1;
+                int za = (int) floor(pos[0] - coffShape[0] / 2);
+                int ze = za + coffShape[0] - 1;
+                PH_A[n][l][m] = (REAL)0;
+                int sz=0;
+                int sy=0;
+                int sx=0;
+                for (int z=za;z<ze;z++,sz++)
+                {
+                    sy=0;
+                    for (int y=ya;y<ye;y++,sy++)
+                    {
+                        sx=0;
+                        for (int x=xa;x<xe;x++,sx++)
+                            PH_A[n][l][m]+=A(z,y,x)*PHbaseF[n][l][m](sz,sy,sx);
+                    }
+                }
+            }
+        }
+    }
+}
+
+template <typename REAL>
+void 
+PHcenter(std::vector<std::vector<std::vector<FFTWComplex<REAL> >  > > &PH_A, 
+         REAL radius, unsigned int band,
+         const MultiArray<3,FFTWComplex<REAL> > &A, 
+         std::vector<std::vector<std::vector<MultiArray<3,FFTWComplex<REAL> > > > > &PHbaseF)
+{
+    PHpos(PH_A, radius, band, A, PHbaseF, detail::centerOfBB<REAL>(A));
+}
+
+template <typename REAL>
+void 
+Array2PH(std::vector<std::vector<std::vector<MultiArray<3,FFTWComplex<REAL> >  > > > &PH_A, 
+         unsigned int band, REAL radius, 
+         bool realData, 
+         const MultiArray<3,FFTWComplex<REAL> > &A, 
+         fftwf_plan forward_plan, fftwf_plan backward_plan)
+{
+    std::vector<std::vector<std::vector<MultiArray<3,FFTWComplex<REAL> > > > > PHbaseF;
+    computePHbaseF(radius, band, PHbaseF, realData);
+
+    PH_A.resize(band+1);
+
+    for (int n=1;n<=band;n++)
+    {
+        PH_A[n].resize(band+1);
+        for (int l=0;l<=band;l++)
+            PH_A[n][l].resize(2*band+1);
+    }
+
+    PH_Iterator<REAL> PHbaseF_Iter(PHbaseF,0,0,0);
+    PH_Iterator<REAL> PHbaseF_Iter_end(PHbaseF,3,PHbaseF.size(),0);
+    PH_Iterator<REAL> PH_A_Iter(PH_A,0,0,0);
+    convolveFFTComplexMany(A, PHbaseF_Iter, PHbaseF_Iter_end, PH_A_Iter,false);
+} 
+
+/*
+template <typename REAL>
+void Array2VH(std::vector<std::vector<std::vector<MultiArray<3, FFTWComplex<REAL> > > > > &VH_A, unsigned int band, REAL gauss, REAL radius, const MultiArray<3, TinyVector<REAL,3> > &A)
+{
+    //transform input to C^(2j+1)
+    MultiArray< 3,FFTWComplex<REAL> >  inputZ(A.shape());
+    MultiArray< 3,FFTWComplex<REAL> >  inputY(A.shape());
+    MultiArray< 3,FFTWComplex<REAL> >  inputX(A.shape());
+    MultiArray< 3,FFTWComplex<REAL> >::iterator z=inputZ.begin();
+    MultiArray< 3,FFTWComplex<REAL> >::iterator y=inputY.begin();
+    MultiArray< 3,FFTWComplex<REAL> >::iterator x=inputX.begin();
+
+    for (MultiArray<3,TinyVector<REAL,3> >::const_iterator p=A.begin();p!=A.end();++p,++z,++y,++x)
+    {
+        z->real() = -(*p)[1];
+        z->imag() = -(*p)[2];
+        *z *= 1/sqrt(2);
+        y->real() =  (*p)[0];
+        y->imag() =  0;
+        x->real() = (*p)[1];
+        x->imag() = -(*p)[2];
+        *x *= 1/sqrt(2);
+    }
+   VH_A.resize(band + 1);
+   REAL norm = M_PI/2*pow((REAL)radius,2.0);
+
+
+    for(int k=-1;k<=1;++k)
+    {
+        VH_A[0][1+k].resize(3);
+        for (int m=-1;m<=1;++m)
+        {
+            MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > > vh;
+            sphereVecHarmonic(vh, radius, gauss, 0, k, m);
+            MultiArray<3, FFTWComplex<REAL> > vhz(vh.shape());
+            MultiArray<3, FFTWComplex<REAL> > vhy(vh.shape());
+            MultiArray<3, FFTWComplex<REAL> > vhx(vh.shape());
+            MultiArray<3, FFTWComplex<REAL> >::iterator z = vhz.begin();
+            MultiArray<3, FFTWComplex<REAL> >::iterator y = vhy.begin();
+            MultiArray<3, FFTWComplex<REAL> >::iterator x = vhx.begin();
+            for(MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > >::iterator p= vh.begin();p!=vh.end();p++,z++,y++,x++)
+            {
+                *z = (*p)[0];
+                *y = (*p)[1];
+                *x = (*p)[2];
+            }
+            VH_A[0][1+k][1+m].reshape(A.shape());
+            VH_A[0][1+k][1+m] = convolveComplex(inputZ, vhz);
+            VH_A[0][1+k][1+m] += convolveComplex(inputY, vhy);
+            VH_A[0][1+k][1+m] += convolveComplex(inputX, vhx);
+
+            for (MultiArray<3, FFTWComplex<REAL> >::iterator p=VH_A[0][1+k][1+m].begin();p!=VH_A[0][1+k][1+m].end();++p)
+                *p*= (complex<REAL>)pow(-1.0,(REAL) 0)/(3*norm);
+        }
+    }
+
+    for (int l=1;l<=band;l++)
+    {
+        VH_A[l].resize(3);
+        for(int k=-1;k<=1;++k)
+        {
+            VH_A[l][1+k].resize(2*(l+k)+1);
+            for (int m=-(l+k);m<=(l+k);++m)
+            {
+                MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > > vh;
+                sphereVecHarmonic(vh, radius, gauss, 0, k, m);
+                MultiArray<3, FFTWComplex<REAL> > vhz(vh.shape());
+                MultiArray<3, FFTWComplex<REAL> > vhy(vh.shape());
+                MultiArray<3, FFTWComplex<REAL> > vhx(vh.shape());
+                MultiArray<3, FFTWComplex<REAL> >::iterator z = vhz.begin();
+                MultiArray<3, FFTWComplex<REAL> >::iterator y = vhy.begin();
+                MultiArray<3, FFTWComplex<REAL> >::iterator x = vhx.begin();
+                for(MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > >::iterator p= vh.begin();p!=vh.end();p++,z++,y++,x++)
+                {
+                    *z = (*p)[0];
+                    *y = (*p)[1];
+                    *x = (*p)[2];
+                }
+
+                VH_A[l][1+k][1+m].reshape(A.shape());
+
+
+                VH_A[l][1+k][1+m] = convolveComplex(inputZ, vhz);
+                VH_A[l][1+k][1+m] += convolveComplex(inputY, vhy);
+                VH_A[l][1+k][1+m] += convolveComplex(inputX, vhx);
+
+                for (MultiArray<3, FFTWComplex<REAL> >::iterator p=VH_A[l][1+k][1+m].begin();p!=VH_A[l][1+k][1+m].end();++p)
+                    *p*= (complex<REAL>)pow(-1.0,(REAL) 0)/(3*norm);
+
+            }
+
+        }
+    }
+}
+
+template <typename REAL>
+void Array2VPH(std::vector<std::vector<std::vector<std::vector<MultiArray<3, FFTWComplex<REAL> > > > > >&VH_A, unsigned int band, REAL radius, const MultiArray<3, TinyVector<REAL,3> > &A)
+{
+    //transform input to C^(2j+1)
+    MultiArray< 3,FFTWComplex<REAL> >  inputZ(A.shape());
+    MultiArray< 3,FFTWComplex<REAL> >  inputY(A.shape());
+    MultiArray< 3,FFTWComplex<REAL> >  inputX(A.shape());
+    MultiArray< 3,FFTWComplex<REAL> >::iterator z=inputZ.begin();
+    MultiArray< 3,FFTWComplex<REAL> >::iterator y=inputY.begin();
+    MultiArray< 3,FFTWComplex<REAL> >::iterator x=inputX.begin();
+
+    for (MultiArray<3,TinyVector<REAL,3> >::const_iterator p=A.begin();p!=A.end();++p,++z,++y,++x)
+    {
+        real(*z) = -(*p)[1];
+        imag(*z) = -(*p)[2];
+        *z *= 1/sqrt(2);
+        real(*y) =  (*p)[0];
+        imag(*y) =  0;
+        real(*x) = (*p)[1];
+        imag(*x) = -(*p)[2];
+        *x *= 1/sqrt(2);
+    }
+
+    for(int n=1;n<=band;n++)
+    {
+        VH_A.resize(band+1);
+        //0th coeff
+        VH_A[n][0].resize(3);
+
+        REAL norm = M_PI/2*pow((REAL)radius,2.0);
+
+        for(int k=-1;k<=1;++k)
+        {
+            VH_A[n][0][1+k].resize(3);
+            for (int m=-1;m<=1;++m)
+            {
+                MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > > vh;
+                sphereFullVecHarmonic(vh, radius,  n, 0, k, m);
+                MultiArray<3, FFTWComplex<REAL> > vhz(vh.shape());
+                MultiArray<3, FFTWComplex<REAL> > vhy(vh.shape());
+                MultiArray<3, FFTWComplex<REAL> > vhx(vh.shape());
+                MultiArray<3, FFTWComplex<REAL> >::iterator z = vhz.begin();
+                MultiArray<3, FFTWComplex<REAL> >::iterator y = vhy.begin();
+                MultiArray<3, FFTWComplex<REAL> >::iterator x = vhx.begin();
+                for(MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > >::iterator p= vh.begin();p!=vh.end();p++,z++,y++,x++)
+                {
+                    *z = (*p)[0];
+                    *y = (*p)[1];
+                    *x = (*p)[2];
+                }
+                VH_A[n][0][1+k][1+m].reshape(A.shape());
+                VH_A[n][0][1+k][1+m] = convolveComplex(inputZ, vhz);
+                VH_A[n][0][1+k][1+m] += convolveComplex(inputY, vhy);
+                VH_A[n][0][1+k][1+m] += convolveComplex(inputX, vhx);
+
+                for (MultiArray<3, FFTWComplex<REAL> >::iterator p=VH_A[n][0][1+k][1+m].begin();p!=VH_A[n][0][1+k][1+m].end();++p)
+                    *p*= (complex<REAL>)pow(-1.0,(REAL) 0)/(3*norm);
+
+            }
+        }
+
+        for (int l=1;l<=band;l++)
+        {
+            VH_A[n][l].resize(3);
+            for(int k=-1;k<=1;++k)
+            {
+                VH_A[n][l][1+k].resize(2*(l+k)+1);
+                for (int m=-(l+k);m<=(l+k);++m)
+                {
+                    MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > > vh;
+                    sphereFullVecHarmonic(vh, radius,  n, 0, k, m);
+                    MultiArray<3, FFTWComplex<REAL> > vhz(vh.shape());
+                    MultiArray<3, FFTWComplex<REAL> > vhy(vh.shape());
+                    MultiArray<3, FFTWComplex<REAL> > vhx(vh.shape());
+                    typename MultiArray<3, FFTWComplex<REAL> >::iterator z = vhz.begin();
+                    typename MultiArray<3, FFTWComplex<REAL> >::iterator y = vhy.begin();
+                    typename MultiArray<3, FFTWComplex<REAL> >::iterator x = vhx.begin();
+                    typename MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > >::iterator p= vh.begin();
+                    for(;p!=vh.end();p++,z++,y++,x++)
+                    {
+                        *z = (*p)[0];
+                        *y = (*p)[1];
+                        *x = (*p)[2];
+                    }
+
+                    VH_A[n][l][1+k][1+m].reshape(A.shape());
+
+
+                    VH_A[n][l][1+k][1+m] = convolveComplex(inputZ, vhz);
+                    VH_A[n][l][1+k][1+m] += convolveComplex(inputY, vhy);
+                    VH_A[n][l][1+k][1+m] += convolveComplex(inputX, vhx);
+
+                    typename MultiArray<3, FFTWComplex<REAL> >::iterator p=VH_A[n][l][1+k][1+m].begin();
+                    for (;p!=VH_A[n][l][1+k][1+m].end();++p)
+                        *p*= (complex<REAL>)pow(-1.0,(REAL) 0)/(3*norm);
+
+                }
+            }
+        }
+    }
+}
+*/
+
+template <typename REAL>
+void 
+VHpos(std::vector<std::vector<std::vector<FFTWComplex<REAL> > > > &VH_A, 
+      unsigned int band, REAL gauss, REAL radius, 
+      const MultiArray<3,TinyVector<REAL,3> > &A, 
+      const TinyVector<REAL, 3> &pos)
+{
+    //transform input to C^(2j+1)
+    MultiArray< 3,FFTWComplex<REAL> >  inputZ(A.shape());
+    MultiArray< 3,FFTWComplex<REAL> >  inputY(A.shape());
+    MultiArray< 3,FFTWComplex<REAL> >  inputX(A.shape());
+    typename MultiArray< 3,FFTWComplex<REAL> >::iterator z=inputZ.begin();
+    typename MultiArray< 3,FFTWComplex<REAL> >::iterator y=inputY.begin();
+    typename MultiArray< 3,FFTWComplex<REAL> >::iterator x=inputX.begin();
+
+    for (typename MultiArray<3,TinyVector<REAL,3> >::const_iterator p=A.begin();p!=A.end();++p,++z,++y,++x)
+    {
+        z->real() = -(*p)[1];
+        z->imag() = -(*p)[2];
+        *z *= 1/sqrt(2);
+        y->real() =  (*p)[0];
+        y->imag() =  0;
+        x->real() = (*p)[1];
+        x->imag() = -(*p)[2];
+        *x *= 1/sqrt(2);
+    }
+    VH_A.resize(band + 1);
+    REAL norm = M_PI/2*std::pow((REAL)radius,2.0);
+
+    // SH coeffs
+    int xa=0,xe=0,ya=0,ye=0,za=0,ze=0;
+
+    VH_A[0].resize(3);
+    for(int k=-1;k<=1;++k)
+    {
+        VH_A[0][1+k].resize(3);
+        for (int m=-1;m<=1;++m)
+        {
+            MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > > vh;
+            sphereVecHarmonic(vh, radius, gauss, 0, k, m);
+
+            xa = (int) floor(pos[2] - vh.shape()[2] / 2);
+            xe = xa + vh.shape()[2] - 1;
+            ya = (int) floor(pos[1] - vh.shape()[1] / 2);
+            ye = ya + vh.shape()[1] - 1;
+            za = (int) floor(pos[0] - vh.shape()[0] / 2);
+            ze = za + vh.shape()[0] - 1;
+
+            int zs=0;
+            int ys=0;
+            int xs=0;
+            for (int z=za;z!=ze;z++,zs++)
+            {
+                for (int y=ya;y<=ye;y++,ys++)
+                {
+                    for (int x=xa;x<=xe;x++,xs++)
+                    {
+                        VH_A[0][1+k][1+m] += inputZ(z,y,x) * vh(zs,ys,xs)[0]/norm + inputZ(z,y,x) * vh(zs,ys,xs)[1]/norm + inputZ(z,y,x) * vh(zs,ys,xs)[2]/norm;
+                    }
+                }
+            }
+        }
+    }
+    for (int l = 1; l <= band; l++)
+    {
+        VH_A[l].resize(3);
+        for (int k=-1;k<=1;++k)
+        {
+            VH_A[l][k+1].resize(2 * (l+k) + 1);
+            for (int m = -(l+k); m <= (l+k); m++)
+            {
+                MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > > vh;
+                sphereVecHarmonic(vh, radius, gauss, 0, k, m);
+
+                xa = (int) floor(pos[2] - vh.shape()[2] / 2);
+                xe = xa + vh.shape()[2] - 1;
+                ya = (int) floor(pos[1] - vh.shape()[1] / 2);
+                ye = ya + vh.shape()[1] - 1;
+                za = (int) floor(pos[0] - vh.shape()[0] / 2);
+                ze = za + vh.shape()[0] - 1;
+
+                int zs=0;
+                for (int z=za;z!=ze;z++,zs++)
+                {
+                    int ys=0;
+                    for (int y=ya;y<=ye;y++,ys++)
+                    {
+                        int xs=0;
+                        for (int x=xa;x<=xe;x++,xs++)
+                        {
+                            VH_A[l][1+k][1+m] += inputZ(z,y,x) * vh(zs,ys,xs)[0]/norm + inputZ(z,y,x) * vh(zs,ys,xs)[1]/norm + inputZ(z,y,x) * vh(zs,ys,xs)[2]/norm;
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+template <typename REAL>
+void 
+VHcenter(std::vector<std::vector<std::vector< FFTWComplex<REAL> > > >&VH_A, 
+         unsigned int band, REAL gauss, REAL radius, 
+         const MultiArray<3, TinyVector<REAL,3> > &A)
+{
+    VHpos(VH_A, band, gauss, radius, A, detail::centerOfBB(A));
+}
+
+template <typename REAL>
+void 
+VPHpos(std::vector<std::vector<std::vector<std::vector<FFTWComplex<REAL> > > > >&VH_A, 
+       unsigned int band, REAL radius, 
+       const MultiArray<3, TinyVector<REAL,3> > &A, const TinyVector<REAL, 3> pos)
+{
+    REAL gauss = 1;
+    //transform input to C^(2j+1)
+    MultiArray< 3,FFTWComplex<REAL> >  inputZ(A.shape());
+    MultiArray< 3,FFTWComplex<REAL> >  inputY(A.shape());
+    MultiArray< 3,FFTWComplex<REAL> >  inputX(A.shape());
+    typename MultiArray< 3,FFTWComplex<REAL> >::iterator z=inputZ.begin();
+    typename MultiArray< 3,FFTWComplex<REAL> >::iterator y=inputY.begin();
+    typename MultiArray< 3,FFTWComplex<REAL> >::iterator x=inputX.begin();
+
+    typename MultiArray<3,TinyVector<REAL,3> >::const_iterator p=A.begin();
+    for (;p!=A.end();++p,++z,++y,++x)
+    {
+        z->real() = -(*p)[1];
+        z->imag() = -(*p)[2];
+        *z *= 1/sqrt(2);
+        y->real() =  (*p)[0];
+        y->imag() =  0;
+        x->real() = (*p)[1];
+        x->imag() = -(*p)[2];
+        *x *= 1/sqrt(2);
+    }
+
+    VH_A.resize(band + 1);
+    for(int n=1;n<=band;n++)
+    {
+        VH_A[n].resize(band + 1);
+        REAL norm = M_PI/2*std::pow((REAL)radius,2.0);
+
+        // SH coeffs
+        int xa=0,xe=0,ya=0,ye=0,za=0,ze=0;
+
+        VH_A[n][0].resize(3);
+        for(int k=-1;k<=1;++k)
+        {
+            VH_A[n][0][1+k].resize(3);
+            for (int m=-1;m<=1;++m)
+            {
+                MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > > vh;
+                sphereVecHarmonic(vh, radius, gauss, 0, k, m);
+
+                xa = (int) floor(pos[2] - vh.shape()[2] / 2);
+                xe = xa + vh.shape()[2] - 1;
+                ya = (int) floor(pos[1] - vh.shape()[1] / 2);
+                ye = ya + vh.shape()[1] - 1;
+                za = (int) floor(pos[0] - vh.shape()[0] / 2);
+                ze = za + vh.shape()[0] - 1;
+
+                int zs=0;
+                for (int z=za;z!=ze;z++,zs++)
+                {
+                    int ys=0;
+                    for (int y=ya;y<=ye;y++,ys++)
+                    {
+                        int xs=0;
+                        for (int x=xa;x<=xe;x++,xs++)
+                        {
+                            VH_A[n][0][1+k][1+m] += inputZ(z,y,x) * vh(zs,ys,xs)[0]/norm + inputZ(z,y,x) * vh(zs,ys,xs)[1]/norm + inputZ(z,y,x) * vh(zs,ys,xs)[2]/norm;
+                        }
+                    }
+                }
+            }
+        }
+        for (int l = 1; l <= band; l++)
+        {
+            VH_A[n][l].resize(3);
+            for (int k=-1;k<=1;++k)
+            {
+                VH_A[n][l][k+1].resize(2 * (l+k) + 1);
+                for (int m = -(l+k); m <= (l+k); m++)
+                {
+                    MultiArray<3, TinyVector<FFTWComplex<REAL>,3 > > vh;
+                    sphereVecHarmonic(vh, radius, gauss, 0, k, m);
+
+                    xa = (int) floor(pos[2] - vh.shape()[2] / 2);
+                    xe = xa + vh.shape()[2] - 1;
+                    ya = (int) floor(pos[1] - vh.shape()[1] / 2);
+                    ye = ya + vh.shape()[1] - 1;
+                    za = (int) floor(pos[0] - vh.shape()[0] / 2);
+                    ze = za + vh.shape()[0] - 1;
+
+                    int zs=0;
+                    int ys=0;
+                    int xs=0;
+                    for (int z=za;z!=ze;z++,zs++)
+                    {
+                        for (int y=ya;y<=ye;y++,ys++)
+                        {
+                            for (int x=xa;x<=xe;x++,xs++)
+                            {
+                                VH_A[n][l][1+k][1+m] += inputZ(z,y,x) * vh(zs,ys,xs)[0]/norm + inputZ(z,y,x) * vh(zs,ys,xs)[1]/norm + inputZ(z,y,x) * vh(zs,ys,xs)[2]/norm;
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+template <typename REAL>
+void 
+VPHcenter(std::vector<std::vector<std::vector<std::vector< FFTWComplex<REAL> > > > >&VH_A, 
+          unsigned int band, REAL radius, 
+          const MultiArray<3, TinyVector<REAL,3> > &A)
+{
+    VPHpos(VH_A, band, radius, A, detail::centerOfBB(A));
+}
+
+
+} // namespace vigra 
+
+#endif // VIGRA_INVARIANT_FEATURES3D_HXX
diff --git a/include/vigra/invariant_features3D.hxx b/include/vigra/invariant_features3D.hxx
index 5c15d974a..2bf0dca55 100644
--- a/include/vigra/invariant_features3D.hxx
+++ b/include/vigra/invariant_features3D.hxx
@@ -45,26 +45,385 @@
 #include "matrix.hxx"
 #include "tinyvector.hxx"
 #include "quaternion.hxx"
+#include "harmonics.hxx"
+#include "clebsch-gordan.hxx"
+#include "multi_fft.hxx"
+#include "wigner-matrix.hxx"
+#include "functortraits.hxx"
+#include <vigra/multi_math.hxx>
 
 namespace vigra {
 
-namespace detail  {
+/** \addtogroup InvariantFeatures local invariant features
+*/
+//@{
 
-// computes the normalization for SH base functions
-inline double realSH(double l, double m)
+/********************************************************/
+/*                                                      */
+/*           Spherical Harmonic cross correlation       */
+/*                                                      */
+/********************************************************/
+
+/** \brief Computes the complete cross correlation over all possible 3D
+    rotations of two Spherical Harmonic (SH) representations. See harmonics.hxx
+    on how to compute SH representations of 3D volume data.
+
+    Returns (normalized) correlation value as well as the roation offset at maximum correlation
+    of the two signals.
+
+    Uses the fast Fourier space corrleation algorithm presented in:
+
+    <i>Fehr, J., Reisert, M. & Burkhardt, H. "Fast and Accurate Rotation Estimation on the 2-Sphere without Correspondences", Computer Vision -- ECCV 2008, LNCS 5303.</i>
+
+    More details and theoretic analysis in:
+
+    <i>2009: Fehr, J. "Local invariant Features for 3D Image Analysis", Dissertation, 2009, pages:47-53</i> 
+
+    NOTE: this is actually not an invariant feature! Used for one-on-one euqivariant comparison. 
+    Use SHautocorr if a feature is needed. 
+
+    The Euler angle representation is given in ZYZ' notation. 
+
+    \param SH_A	first SH expansion
+    \param SH_B	second SH expansion
+    \param padsize size of the padding for FFT sinc interpolation. Larger padding increases the
+    accuracy of the rotation estimation at increased cost. Theoretic limit of the angular resolution is: Err >= 180°/(2*b_max+padzize), where b_max is the maximum band of the SH expansion. Reasonable
+    values for padsize are [0,32,64,...,512] 
+    \param normalize bool:compute normalized cross corelation with perfect match at max==1 or unnormalized
+    \param max returns maximum correlation value
+    \param phi returns rotation offset in first Euler angle [0...2PI]
+    \param theta returns rotation offset in second Euler angle [0...PI]
+    \param psi returns rotation offset in third Euler angle [0...2PI]
+    \param corMat returns the 3D (theta,phi,psi) correlation space
+    \WARNING: there are natural ambiguities of the rotation representation in Euler angles, especially for theta=0  
+
+*/
+template <typename REAL>
+void 
+SHxcorr(std::vector<std::vector<FFTWComplex<REAL> >  > const &SH_A, 
+        std::vector<std::vector<FFTWComplex<REAL> >  > const &SH_B,
+        unsigned int padsize, 
+        bool normalize, 
+        REAL &max, REAL &phi, REAL &theta, REAL &psi, 
+        MultiArray<3,REAL > &corMat)
+{
+    int band=SH_A.size()-1;
+    TinyVector<int,3> baseshape(2*band+padsize,2*band+padsize,2*band+padsize);
+    TinyVector<int,3> padshape(fftwBestPaddedShape(baseshape));
+    MultiArray<3,FFTWComplex<REAL> > C_fft;
+    C_fft.reshape(padshape,0);
+
+    int size = padshape[0];
+
+    WignerMatrix<REAL> W(band);
+    W.compute_D(band);
+    FFTWComplex<REAL> Im(0, 1);
+
+    for (int l=0;l<=band;l++)
+    {
+        for (int m1=-l;m1<=l;m1++)
+        {
+            for (int m2=-l;m2<=l;m2++)
+            {
+                for (int m3=-l;m3<=l;m3++)
+                {
+                    C_fft((m2+size+1)%(size),(m1+size+1)%(size),(m3+size+1)%(size)) +=  
+                    (W.get_d(l,m1,m2) * W.get_d(l,m2,m3) * (SH_A[l][l+m1]) * conj(SH_B[l][l+m3])) *
+                       pow(Im, (m1 + 2*m2 + m3));
+                }
+            }
+        }
+    }
+
+    MultiArray<3,FFTWComplex<REAL> > res(C_fft.shape());
+
+    //FFT^-1
+    fourierTransformInverse(C_fft,res);
+    
+    REAL variance = 1;
+    if(normalize)
+    {
+        REAL varianceA=0, varianceB=0;
+        for (int l=0;l<=band;l++)
+        {
+            REAL tmpA=0, tmpB=0;
+            for (int m=-l;m<=l;m++)
+            {
+                tmpA += (REAL)(std::pow((double)real(SH_A[l][m+l]),2.0) + 
+                               std::pow((double)imag(SH_A[l][m+l]),2.0));
+                tmpB += (REAL)(std::pow((double)real(SH_B[l][m+l]),2.0) + 
+                               std::pow((double)imag(SH_B[l][m+l]),2.0));
+            }
+            varianceA += tmpA;
+            varianceB += tmpB;
+        }
+
+        variance = std::sqrt(varianceA*varianceB)/(size*size*size);
+        if (variance == 0) 
+        {
+            variance = FLT_MIN;
+        }
+    }
+
+    corMat.reshape(res.shape(),0);
+    
+    TinyVector<int, 3> hotpos1;
+    max = 0;
+    for (int z=0;z<corMat.shape()[0];z++)
+    {
+        for (int y=0;y<corMat.shape()[1];y++)
+        {
+            for (int x=0;x<corMat.shape()[2];x++)
+            {
+                REAL tmp = abs(res(z,y,x)/variance);
+                if (tmp>max)
+                {
+                    max=tmp;
+                    hotpos1[0]=z;
+                    hotpos1[1]=y;
+                    hotpos1[2]=x;
+                }
+                corMat(z,y,x) = tmp;
+            }
+        }
+    }
+    REAL grid = REAL(2.0*M_PI/(size));
+
+    //phi
+    if(hotpos1[1]*grid >M_PI)
+        phi = REAL(M_PI + (2*M_PI -hotpos1[1]*grid ));
+    else
+        phi = REAL(M_PI - hotpos1[1]*grid);
+
+    //theta     
+    if(hotpos1[0]*grid >M_PI)
+        theta = REAL(2*M_PI - (hotpos1[0])*grid);
+    else
+        theta = REAL(hotpos1[0]*grid);
+
+    //psi
+    if(hotpos1[2]*grid >M_PI)
+        psi = REAL(M_PI + (2*M_PI -hotpos1[2]*grid ));
+    else
+        psi = REAL(M_PI - hotpos1[2]*grid);
+}
+
+//@{
+
+/********************************************************/
+/*                                                      */
+/*           Spherical Harmonic auto correlation       */
+/*                                                      */
+/********************************************************/
+
+/** \brief Computes integrates over  complete auto correlation over all possible 3D
+    rotations of a Spherical Harmonic (SH) representation. See harmonics.hxx
+    on how to compute SH representations of 3D volume data.
+
+    Returns scalar invariant feature.
+
+    More details and theoretic analysis in:
+
+    <i>2009: Fehr, J. "Local invariant Features for 3D Image Analysis", Dissertation, 2009, pages:73-74</i> 
+
+    \param SH_A  SH expansion
+    \param padsize size of the padding for FFT sinc interpolation (see SHxcorr for details). 
+    \param f functor with scalar non-linear mapping of the correlation values before integration. 
+    Reasonable functors could be f(x)=x^2, sqrt(x), 1/x, ...
+    \param res return feature value
+*/
+template <typename REAL, typename FUNCTOR>
+void 
+SHautocorr(std::vector<std::vector<FFTWComplex<REAL> >  > const &SH_A, 
+           unsigned int padsize, FUNCTOR &f, REAL &res)
 {
-    return std::sqrt((2.0*l + 1.0) / (4.0*M_PI*facLM(l,m)));
+    REAL cmax,phi,theta,psi=0;
+    MultiArray<3,float > corMat;
+    SHxcorr(SH_A,SH_A,padsize,true,cmax,phi,theta,psi,corMat);
+    res=0;
+    for ( MultiArray<3,float >::iterator p=corMat.begin();p!=corMat.end();p++)
+        res *= f(*p);
+}
+
+/********************************************************/
+/*                                                      */
+/*           Spherical Harmonic power-specturm          */
+/*                                                      */
+/********************************************************/
+
+/** \brief Computes the power spectrum of the SH representation. 
+    See harmonics.hxx on how to compute SH representations of 3D volume data.
+    
+    Returns vectorial invariant feature.
+
+    Well known Idea, one of the first papers using the SH power-spectrum as a feature was:
+
+    <i>M. Kazhdan. Rotation invariant spherical harmonic representation of 3d shape descriptors. Symp. on Geom. Process., 2003.</i>
+ 
+    More details and theoretic analysis in:
 
+    <i>2009: Fehr, J. "Local invariant Features for 3D Image Analysis", Dissertation, 2009, pages:68-69</i> 
+
+    \param SH_A  SH expansion
+    \param res return vector of feature values
+*/
+template <typename REAL>
+void 
+SHabs(std::vector<std::vector<FFTWComplex<REAL> >  > const &SH_A, 
+      std::vector<REAL> &res)
+{
+    res.resize(SH_A.size());
+    for (int l=0;l<(int)SH_A.size();l++)
+    {
+        res[l]=0;
+        for (int m=-l;m<=l;m++)
+            res[l] += (REAL)(std::pow((double)real(SH_A[l][m+l]),2.0)+ 
+                             std::pow((double)imag(SH_A[l][m+l]),2.0));
+    }   
 }
 
-template<int N, class T, class C>
-TinyVector<float, N> centerOfBB(MultiArrayView<N, T, C> const & A)
+//this version computes feature at all positions, use Array2SH to get SH_A 
+template <typename REAL>
+void 
+SHabs(std::vector<std::vector<MultiArray<3,FFTWComplex<REAL> > > >const &SH_A, 
+      std::vector<MultiArray<3,REAL> >&res)
 {
-    return TinyVector<float, N>(A.shape()) /= 2.0;                        
+    using namespace multi_math;
+    res.resize(SH_A.size());
+    for (int l=0;l<(int)SH_A.size();l++)
+    {
+        res[l].reshape(SH_A[l][0].shape(),0);
+        for (int m=-l;m<=l;m++)
+        {
+            typename MultiArray<3,REAL >::iterator q=res[l].begin();
+            typename MultiArray<3,FFTWComplex<REAL> >::iterator p=SH_A[l][m+l].begin();
+            for (;p<SH_A[l][m+l].end();p++,q++)
+            {
+                *q += REAL(std::pow((double)real(*p),2.0)+ std::pow((double)imag(*p),2.0));
+            }
+        }
+    }
 }
 
-} // namespace detail
 
+/********************************************************/
+/*                                                      */
+/*           Spherical Harmonic bi-specturm             */
+/*                                                      */
+/********************************************************/
+
+/** \brief Computes the bi-spectrum of the SH representation. 
+    See harmonics.hxx on how to compute SH representations of 3D volume data.
+    
+    Returns vectorial invariant feature.
+
+    This is a complete feature! Turns out to be less robust for real world data, hence SH expansion
+    should be limited to low bands only. Original paper:
+
+    <i>S. Wenndt and S. Shamsunder. Bispectrum features for robust speaker identification. Acoustics, Speech, and Signal Processing, IEEE International Conference on, 2:1095, 1997.</i>
+ 
+    More details and theoretic analysis in:
+
+    <i>2009: Fehr, J. "Local invariant Features for 3D Image Analysis", Dissertation, 2009, pages:75-76</i> 
+
+    \param SH_A  SH expansion
+    \param res return vector of feature values
+*/
+
+template <typename REAL>
+void 
+SHbispec(std::vector<std::vector<FFTWComplex<REAL> >  > const &SH_A, 
+         std::vector<FFTWComplex<REAL> > &res)
+{
+    res.clear();
+    for (int l=0;l<(int)SH_A.size();l++)
+    {
+        for (int l1=0;l1<=l;l1++)
+        {
+            for (int l2=0;l2<=l;l2++)
+            {
+                // some l,l1,l2 combinations are not allowed (triangular inequallity) - rejected by Clebschgordan
+                try
+                {
+                    FFTWComplex<REAL> f;
+                    f = (REAL)0.0;
+
+                    for(int m=-l;m<=l;m++)
+                    {
+                        FFTWComplex<REAL> tmp_feature;
+                        tmp_feature = (REAL)0.0;
+                        for(int m1=std::max(-l1,m-l2);m1<=std::min(l1,m+l2);m1++)
+                        {
+                            FFTWComplex<REAL> s;
+                            s = (REAL)0.0;
+                            s.real() = (REAL)clebschGordan(l1, m1, l2, m-m1, l, m);
+                            tmp_feature += s * conj(SH_A[l1][m1+l1]) * conj(SH_A[l2][m-m1+l2]);
+                        }
+                        f += SH_A[l][m+l] * tmp_feature;
+
+                    }
+
+                    res.push_back(f);
+                }
+                catch (ContractViolation &){}
+            }
+        }
+    }
+}
+
+//this version computes feature at all positions, use Array2SH to get SH_A 
+template <typename REAL>
+void 
+SHbispec(std::vector<std::vector<MultiArray<3,FFTWComplex<REAL> > > >const &SH_A, 
+         std::vector<MultiArray<3,FFTWComplex<REAL> > > &res)
+{
+    res.clear();
+    for (int l=0;l<(int)SH_A.size();l++)
+    {
+        for (int l1=0;l1<=l;l1++)
+        {
+            for (int l2=0;l2<=l;l2++)
+            {
+                // some l,l1,l2 combinations are not allowed (triangular inequallity) - rejected by Clebschgordan
+                try
+                {
+                    MultiArray<3,FFTWComplex<REAL> > f;
+                    f.reshape(SH_A[0][0].shape(),0);
+
+                    for(int m=-l;m<=l;m++)
+                    {
+                        MultiArray<3,FFTWComplex<REAL> > tmp_feature;
+                        tmp_feature.reshape(SH_A[0][0].shape(),0);;
+                        for(int m1=std::max(-l1,m-l2);m1<=std::min(l1,m+l2);m1++)
+                        {
+                            FFTWComplex<REAL> s;
+                            s = (REAL)0.0;
+                            s.real() = (REAL)clebschGordan(l1, m1, l2, m-m1, l, m);
+                            typename MultiArray<3,FFTWComplex<REAL> >::iterator h = 
+                                SH_A[l2][m-m1+l2].begin();
+                            typename MultiArray<3,FFTWComplex<REAL> >::iterator q =
+                                SH_A[l1][m1+l1].begin();
+                            typename MultiArray<3,FFTWComplex<REAL> >::iterator p = 
+                                tmp_feature.begin();
+                            for (;p<tmp_feature.end();p++,q++,h++)
+                                *p += s * conj(*q) * conj(*h);
+                        }
+                        typename MultiArray<3,FFTWComplex<REAL> >::iterator h=SH_A[l][m+l].begin();
+                        typename MultiArray<3,FFTWComplex<REAL> >::iterator q=f.begin();
+                        typename MultiArray<3,FFTWComplex<REAL> >::iterator p = tmp_feature.begin();
+                        for (;p<tmp_feature.end();p++,q++,h++)
+                            *q += *h * *p;
+
+                    }
+
+                    res.push_back(f);
+
+                }
+                catch (ContractViolation &){}
+            }
+        }
+    }
+}
 
 } // namespace vigra 
 
diff --git a/include/vigra/mathutil.hxx b/include/vigra/mathutil.hxx
index f87d6a3b7..e5ebe1e75 100644
--- a/include/vigra/mathutil.hxx
+++ b/include/vigra/mathutil.hxx
@@ -1020,14 +1020,14 @@ REAL legendre(unsigned int l, int m, REAL x)
     {
         m = -m;
         REAL s = odd(m)
-                   ? -1.0
-                   :  1.0;
-        return legendre(l,m,x) * s / detail::facLM<REAL>(l,m);
+                   ? REAL(-1.0)
+                   : REAL( 1.0);
+        return legendre(l,m,x) * s / (REAL)detail::facLM<double>(l,m);
     }
     REAL result = 1.0;
     if (m > 0)
     {
-        REAL r = std::sqrt( (1.0-x) * (1.0+x) );
+        REAL r = (REAL)std::sqrt( (1.0-x) * (1.0+x) );
         REAL f = 1.0;
         for (int i=1; i<=m; i++)
         {
@@ -1038,13 +1038,13 @@ REAL legendre(unsigned int l, int m, REAL x)
     if((int)l == m) 
         return result;
 
-    REAL result_1 = x * (2.0 * m + 1.0) * result;
+    REAL result_1 = REAL(x * (2.0 * m + 1.0) * result);
     if((int)l == m+1) 
         return result_1;
     REAL other = 0.0;
     for(unsigned int i = m+2; i <= l; ++i)
     {
-        other = ( (2.0*i-1.0) * x * result_1 - (i+m-1.0)*result) / (i-m);
+        other = REAL(( (2.0*i-1.0) * x * result_1 - (i+m-1.0)*result) / (i-m));
         result = result_1;
         result_1 = other;
     }
diff --git a/include/vigra/multi_fft.hxx b/include/vigra/multi_fft.hxx
index d9576ed27..e67a0552f 100644
--- a/include/vigra/multi_fft.hxx
+++ b/include/vigra/multi_fft.hxx
@@ -815,10 +815,10 @@ fftEmbedKernel(MultiArrayView<N, Real, C1> kernel,
     FFTEmbedKernel<(int)N-1>::exec(out, kernel.shape(), srcPoint, destPoint, false);
 }
 
-template <unsigned int N, class Real, class C1, class C2>
+template <unsigned int N, class T1, class C1, class T2, class C2>
 void 
-fftEmbedArray(MultiArrayView<N, Real, C1> in,
-              MultiArrayView<N, Real, C2> out)
+fftEmbedArray(MultiArrayView<N, T1, C1> in,
+              MultiArrayView<N, T2, C2> out)
 {
     typedef typename MultiArrayShape<N>::type Shape;
     
@@ -829,7 +829,7 @@ fftEmbedArray(MultiArrayView<N, Real, C1> in,
     
     out.subarray(leftDiff, right) = in;
     
-    typedef typename MultiArrayView<N, Real, C2>::traverser Traverser;
+    typedef typename MultiArrayView<N, T2, C2>::traverser Traverser;
     typedef MultiArrayNavigator<Traverser, N> Navigator;
     typedef typename Navigator::iterator Iterator;
     
@@ -1153,6 +1153,16 @@ class FFTWConvolvePlan
         init(in, kernel, out, planner_flags);
     }
    
+    template <class C1, class C2, class C3>
+    FFTWConvolvePlan(MultiArrayView<N, Real, C1> in,
+                     MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
+                     MultiArrayView<N, FFTWComplex<Real>, C3> out, 
+                     bool fourierDomainKernel,
+                     unsigned int planner_flags = FFTW_ESTIMATE)
+    {
+        init(in, kernel, out, fourierDomainKernel, planner_flags);
+    }
+
     template <class C1, class C2, class C3>
     FFTWConvolvePlan(MultiArrayView<N, FFTWComplex<Real>, C1> in,
                      MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
@@ -1197,6 +1207,19 @@ class FFTWConvolvePlan
         initFourierKernel(in.shape(), kernel.shape(), planner_flags);
     }
     
+    template <class C1, class C2, class C3>
+    void init(MultiArrayView<N, Real, C1> in, 
+              MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
+              MultiArrayView<N, FFTWComplex<Real>, C3> out, 
+              bool fourierDomainKernel,
+              unsigned int planner_flags = FFTW_ESTIMATE)
+    {
+        vigra_precondition(in.shape() == out.shape(),
+            "FFTWConvolvePlan::init(): input and output must have the same shape.");
+        useFourierKernel = fourierDomainKernel;
+        initComplex(in.shape(), kernel.shape(), planner_flags);
+    }
+    
     template <class C1, class C2, class C3>
     void init(MultiArrayView<N, FFTWComplex<Real>, C1> in, 
               MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
@@ -1242,37 +1265,32 @@ class FFTWConvolvePlan
     }
      
     template <class C1, class KernelIterator, class OutIterator>
-    void initMany(MultiArrayView<N, FFTWComplex<Real>, C1> in, 
-                  KernelIterator kernels, KernelIterator kernelsEnd,
-                  OutIterator outs,
-                  bool fourierDomainKernels,
-                  unsigned int planner_flags = FFTW_ESTIMATE)
+    void initManyComplex(MultiArrayView<N, FFTWComplex<Real>, C1> in, 
+                         KernelIterator kernels, KernelIterator kernelsEnd,
+                         OutIterator outs,
+                         bool fourierDomainKernels,
+                         unsigned int planner_flags = FFTW_ESTIMATE)
     {
-        typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
-        typedef typename KernelArray::value_type KernelValue;
-        typedef typename std::iterator_traits<OutIterator>::value_type OutArray;
-        typedef typename OutArray::value_type OutValue;
-
-        vigra_precondition((IsSameType<KernelValue, Complex>::value),
-             "FFTWConvolvePlan::initMany(): kernels have unsuitable value_type.");
-        vigra_precondition((IsSameType<OutValue, Complex>::value),
-             "FFTWConvolvePlan::initMany(): outputs have unsuitable value_type.");
-
-        useFourierKernel = fourierDomainKernels;
-        
-        Shape paddedShape = checkShapesComplex(in.shape(), kernels, kernelsEnd, outs);
-    
-        CArray newFourierArray(paddedShape), newFourierKernel(paddedShape);
-    
-        FFTWPlan<N, Real> fplan(newFourierArray, newFourierArray, FFTW_FORWARD, planner_flags);
-        FFTWPlan<N, Real> bplan(newFourierArray, newFourierArray, FFTW_BACKWARD, planner_flags);
+        initManyComplexImpl(in, kernels, kernelsEnd, outs, fourierDomainKernels, planner_flags);
+    }
     
-        forward_plan = fplan;
-        backward_plan = bplan;
-        fourierArray.swap(newFourierArray);
-        fourierKernel.swap(newFourierKernel);
+    template <class C1, class KernelIterator, class OutIterator>
+    void initManyComplex(MultiArrayView<N, Real, C1> in, 
+                         KernelIterator kernels, KernelIterator kernelsEnd,
+                         OutIterator outs,
+                         bool fourierDomainKernels,
+                         unsigned int planner_flags = FFTW_ESTIMATE)
+    {
+        initManyComplexImpl(in, kernels, kernelsEnd, outs, fourierDomainKernels, planner_flags);
     }
     
+    template <class T, class C1, class KernelIterator, class OutIterator>
+    void initManyComplexImpl(MultiArrayView<N, T, C1> in, 
+                             KernelIterator kernels, KernelIterator kernelsEnd,
+                             OutIterator outs,
+                             bool fourierDomainKernels,
+                             unsigned int planner_flags = FFTW_ESTIMATE);
+    
     void init(Shape inOut, Shape kernel,
               unsigned int planner_flags = FFTW_ESTIMATE);
     
@@ -1305,16 +1323,47 @@ class FFTWConvolvePlan
                  MultiArrayView<N, Real, C3> out);
 
     template <class C1, class C2, class C3>
-    void execute(MultiArrayView<N, FFTWComplex<Real>, C1> in,
-                 MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
-                 MultiArrayView<N, FFTWComplex<Real>, C3> out);
+    void executeComplex(MultiArrayView<N, Real, C1> in,
+                        MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
+                        MultiArrayView<N, FFTWComplex<Real>, C3> out)
+    {
+        executeComplexImpl(in, kernel, out);
+    }
 
+    template <class C1, class C2, class C3>
+    void executeComplex(MultiArrayView<N, FFTWComplex<Real>, C1> in,
+                        MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
+                        MultiArrayView<N, FFTWComplex<Real>, C3> out)
+    {
+        executeComplexImpl(in, kernel, out);
+    }
 
-    template <class C1, class KernelIterator, class OutIterator>
-    void executeMany(MultiArrayView<N, FFTWComplex<Real>, C1> in, 
-                     KernelIterator kernels, KernelIterator kernelsEnd,
-                     OutIterator outs);
+    template <class T, class C1, class C2, class C3>
+    void executeComplexImpl(MultiArrayView<N, T, C1> in,
+                            MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
+                            MultiArrayView<N, FFTWComplex<Real>, C3> out);
+
+    template <class T, class C1, class KernelIterator, class OutIterator>
+    void executeManyComplexImpl(MultiArrayView<N, T, C1> in, 
+                                KernelIterator kernels, KernelIterator kernelsEnd,
+                                OutIterator outs);
                      
+    template <class C1, class KernelIterator, class OutIterator>
+    void executeManyComplex(MultiArrayView<N, FFTWComplex<Real>, C1> in, 
+                            KernelIterator kernels, KernelIterator kernelsEnd,
+                            OutIterator outs)
+    {
+        executeManyComplexImpl(in, kernels, kernelsEnd, outs);
+    }
+    
+    template <class C1, class KernelIterator, class OutIterator>
+    void executeManyComplex(MultiArrayView<N, Real, C1> in, 
+                            KernelIterator kernels, KernelIterator kernelsEnd,
+                            OutIterator outs)
+    {
+        executeManyComplexImpl(in, kernels, kernelsEnd, outs);
+    }
+    
     template <class C1, class KernelIterator, class OutIterator>
     void executeMany(MultiArrayView<N, Real, C1> in, 
                      KernelIterator kernels, KernelIterator kernelsEnd,
@@ -1457,6 +1506,40 @@ FFTWConvolvePlan<N, Real>::initComplex(Shape in, Shape kernel,
     fourierKernel.swap(newFourierKernel);
 }
 
+template <unsigned int N, class Real>
+template <class T, class C1, class KernelIterator, class OutIterator>
+void 
+FFTWConvolvePlan<N, Real>::initManyComplexImpl(MultiArrayView<N, T, C1> in, 
+                                               KernelIterator kernels, KernelIterator kernelsEnd,
+                                               OutIterator outs,
+                                               bool fourierDomainKernels,
+                                               unsigned int planner_flags)
+{
+    typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
+    typedef typename KernelArray::value_type KernelValue;
+    typedef typename std::iterator_traits<OutIterator>::value_type OutArray;
+    typedef typename OutArray::value_type OutValue;
+
+    vigra_precondition((IsSameType<KernelValue, Complex>::value),
+         "FFTWConvolvePlan::initMany(): kernels have unsuitable value_type.");
+    vigra_precondition((IsSameType<OutValue, Complex>::value),
+         "FFTWConvolvePlan::initMany(): outputs have unsuitable value_type.");
+
+    useFourierKernel = fourierDomainKernels;
+    
+    Shape paddedShape = checkShapesComplex(in.shape(), kernels, kernelsEnd, outs);
+
+    CArray newFourierArray(paddedShape), newFourierKernel(paddedShape);
+
+    FFTWPlan<N, Real> fplan(newFourierArray, newFourierArray, FFTW_FORWARD, planner_flags);
+    FFTWPlan<N, Real> bplan(newFourierArray, newFourierArray, FFTW_BACKWARD, planner_flags);
+
+    forward_plan = fplan;
+    backward_plan = bplan;
+    fourierArray.swap(newFourierArray);
+    fourierKernel.swap(newFourierKernel);
+}
+
 template <unsigned int N, class Real>
 template <class C1, class C2, class C3>
 void 
@@ -1529,14 +1612,14 @@ FFTWConvolvePlan<N, Real>::execute(MultiArrayView<N, Real, C1> in,
 }
 
 template <unsigned int N, class Real>
-template <class C1, class C2, class C3>
+template <class T, class C1, class C2, class C3>
 void 
-FFTWConvolvePlan<N, Real>::execute(MultiArrayView<N, FFTWComplex<Real>, C1> in, 
-                                    MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
-                                    MultiArrayView<N, FFTWComplex<Real>, C3> out)
+FFTWConvolvePlan<N, Real>::executeComplexImpl(MultiArrayView<N, T, C1> in, 
+                                              MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
+                                              MultiArrayView<N, FFTWComplex<Real>, C3> out)
 {
     vigra_precondition(in.shape() == out.shape(),
-        "FFTWConvolvePlan::execute(): input and output must have the same shape.");
+        "FFTWConvolvePlan::executeComplex(): input and output must have the same shape.");
     
     Shape paddedShape = fourierArray.shape(),
           diff = paddedShape - in.shape(), 
@@ -1640,11 +1723,11 @@ FFTWConvolvePlan<N, Real>::executeManyImpl(MultiArrayView<N, Real, C1> in,
 }
 
 template <unsigned int N, class Real>
-template <class C1, class KernelIterator, class OutIterator>
+template <class T, class C1, class KernelIterator, class OutIterator>
 void 
-FFTWConvolvePlan<N, Real>::executeMany(MultiArrayView<N, FFTWComplex<Real>, C1> in, 
-                                       KernelIterator kernels, KernelIterator kernelsEnd,
-                                       OutIterator outs)
+FFTWConvolvePlan<N, Real>::executeManyComplexImpl(MultiArrayView<N, T, C1> in, 
+                                                  KernelIterator kernels, KernelIterator kernelsEnd,
+                                                  OutIterator outs)
 {
     typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
     typedef typename KernelArray::value_type KernelValue;
@@ -1903,6 +1986,16 @@ convolveFFT(MultiArrayView<N, Real, C1> in,
     FFTWConvolvePlan<N, Real>(in, kernel, out).execute(in, kernel, out);
 }
 
+template <unsigned int N, class Real, class C1, class C2, class C3>
+void
+convolveFFTComplex(MultiArrayView<N, Real, C1> in,
+            MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
+            MultiArrayView<N, FFTWComplex<Real>, C3> out,
+            bool fourierDomainKernel)
+{
+    FFTWConvolvePlan<N, Real>(in, kernel, out, fourierDomainKernel).executeComplex(in, kernel, out);
+}
+
 template <unsigned int N, class Real, class C1, class C2, class C3>
 void
 convolveFFTComplex(MultiArrayView<N, FFTWComplex<Real>, C1> in,
@@ -1910,7 +2003,7 @@ convolveFFTComplex(MultiArrayView<N, FFTWComplex<Real>, C1> in,
             MultiArrayView<N, FFTWComplex<Real>, C3> out,
             bool fourierDomainKernel)
 {
-    FFTWConvolvePlan<N, Real>(in, kernel, out, fourierDomainKernel).execute(in, kernel, out);
+    FFTWConvolvePlan<N, Real>(in, kernel, out, fourierDomainKernel).executeComplex(in, kernel, out);
 }
 
 template <unsigned int N, class Real, class C1, 
@@ -1925,6 +2018,19 @@ convolveFFTMany(MultiArrayView<N, Real, C1> in,
     plan.executeMany(in, kernels, kernelsEnd, outs);
 }
 
+template <unsigned int N, class Real, class C1, 
+          class KernelIterator, class OutIterator>
+void 
+convolveFFTComplexMany(MultiArrayView<N, Real, C1> in, 
+                KernelIterator kernels, KernelIterator kernelsEnd,
+                OutIterator outs,
+                bool fourierDomainKernel)
+{
+    FFTWConvolvePlan<N, Real> plan;
+    plan.initManyComplex(in, kernels, kernelsEnd, outs, fourierDomainKernel);
+    plan.executeManyComplex(in, kernels, kernelsEnd, outs);
+}
+
 template <unsigned int N, class Real, class C1, 
           class KernelIterator, class OutIterator>
 void 
@@ -1934,8 +2040,8 @@ convolveFFTComplexMany(MultiArrayView<N, FFTWComplex<Real>, C1> in,
                 bool fourierDomainKernel)
 {
     FFTWConvolvePlan<N, Real> plan;
-    plan.initMany(in, kernels, kernelsEnd, outs, fourierDomainKernel);
-    plan.executeMany(in, kernels, kernelsEnd, outs);
+    plan.initManyComplex(in, kernels, kernelsEnd, outs, fourierDomainKernel);
+    plan.executeManyComplex(in, kernels, kernelsEnd, outs);
 }
 
 } // namespace vigra
diff --git a/include/vigra/wigner-matrix.hxx b/include/vigra/wigner-matrix.hxx
index 4f7befcb5..2bc48fa04 100644
--- a/include/vigra/wigner-matrix.hxx
+++ b/include/vigra/wigner-matrix.hxx
@@ -46,9 +46,12 @@
 #include "tinyvector.hxx"
 #include "quaternion.hxx"
 #include "clebsch-gordan.hxx"
+#include "multi_fft.hxx"
+
 
 namespace vigra {
 
+
     /*!
      *  \class WignerMatrix 
      *  \brief computation of Wigner D matrix + rotation functions 
@@ -64,8 +67,7 @@ class WignerMatrix
 {   
   public:
   
-    // FIXME: should we rather use FFTWComplex?
-    typedef std::complex<Real> Complex;
+    typedef FFTWComplex<Real> Complex;
     typedef ArrayVector<ArrayVector<ArrayVector<Complex> > > NestedArray;
     
          /*! \brief constructor
@@ -80,7 +82,6 @@ class WignerMatrix
           *
           * \param band	expansion band
           *
-             FIXME: compute_D(l, 0.0, M_PI / 2.0, 0.0) creates the transposed matrix!
           */
     void compute_D(int band);
 
@@ -100,17 +101,17 @@ class WignerMatrix
           * \param n      	
           * \param m    
           */
-    Complex get_D(int l, int n, int m) const
+    Complex get_d(int l, int h, int m) const
     {
         if (l>0)
         {
             std::string message = std::string("WignerMatrix::get_D(): index out of bounds: l=");
-            message << l << " l_max=" << D.size() << " m=" << m << " n=" << n << "\n";
+            message << l << " l_max=" << D.size() << " m=" << m << " n=" << h << "\n";
                                   
-            vigra_precondition(l < D.size() && m+l <= 2*l+1 &&
-                               n+l <= 2*l+1 && m+l >= 0 && n+l >= 0,
+            vigra_precondition(l < (int)D.size() && m+l <= 2*l+1 &&
+                               h+l <= 2*l+1 && m+l >= 0 && h+l >= 0,
                                message.c_str());
-            return D[l](n+l, m+l);
+            return D[l](h+l, m+l);
         }
         else 
         {
@@ -140,21 +141,13 @@ class WignerMatrix
           *
           * \retval PHresult PH expansion   
           */
-    void rotatePH(NestedArray const & PH, Real phi, Real theta, Real psi,
-                  NestedArray & PHresult);
+    void rotatePH(std::vector<std::vector<std::vector<vigra::FFTWComplex<Real> > > > const & PH, Real phi, Real theta, Real psi, std::vector<std::vector<std::vector<vigra::FFTWComplex<Real> > > >  & PHresult);
+
+    void rotateSH(std::vector<std::vector<vigra::FFTWComplex<Real> > >  const & SH, Real phi, Real theta, Real psi, std::vector<std::vector<vigra::FFTWComplex<Real> > > & SHresult);
 
 
-  private:
-    // FIXME: is function is not called (and cannot be called from outside the class)
-    TinyVector<double,3> 
-    rot(TinyVector<double,3> const & vec, TinyVector<double,3> const & axis, double angle)
-    {
-        typedef Quaternion<double> Q;
-        Q qr = Q::createRotation(angle, axis),
-          qv(0.0, vec);
-        return (qr * qv * conj(qr)).v();
-    }
 
+  private:
     int l_max;
     int cg1cnt;
     ArrayVector<double> CGcoeff;
@@ -194,27 +187,27 @@ template <class Real>
 void
 WignerMatrix<Real>::compute_D(int l, Real phi, Real theta, Real psi)
 {
-    double s = std::sin(theta);
-    double c = std::cos(theta);
+    Real s = std::sin(theta);
+    Real c = std::cos(theta);
     
-    Complex i(0.0, 1.0);
-    Complex eiphi = std::exp(i*phi);
-    Complex emiphi = std::exp(-i*phi);
-    Complex eipsi = std::exp(i*psi);
-    Complex emipsi = std::exp(-i*psi);
+    FFTWComplex<Real> i(0.0, 1.0);
+    FFTWComplex<Real> eiphi = exp(i*phi);
+    FFTWComplex<Real> emiphi = exp(-i*phi);
+    FFTWComplex<Real> eipsi = exp(i*psi);
+    FFTWComplex<Real> emipsi = exp(-i*psi);
     
     if (D.size() < (std::size_t)(l+1)) 
         D.resize(l+1);
     D[1].reshape(MultiArrayShape<2>::type(3,3));
     
     D[1](0,0) = emipsi * Complex(Real(0.5*(1.0+c))) * emiphi;
-    D[1](0,1) = Complex(Real(-s/M_SQRT2)) * emiphi;
-    D[1](0,2) = eipsi * Complex(Real(0.5*(1.0-c))) * emiphi;
-    D[1](1,0) = emipsi * Complex(Real(s/M_SQRT2));
+    D[1](1,0) = Complex(Real(-s/M_SQRT2)) * emiphi;
+    D[1](2,0) = eipsi * Complex(Real(0.5*(1.0-c))) * emiphi;
+    D[1](0,1) = emipsi * Complex(Real(s/M_SQRT2));
     D[1](1,1) = Complex(Real(c));
-    D[1](1,2) = eipsi * Complex(Real(-s/M_SQRT2));
-    D[1](2,0) = emipsi * Complex(Real(0.5*(1.0-c))) * eiphi; 
-    D[1](2,1) = Complex(Real(s/M_SQRT2)) * eiphi;
+    D[1](2,1) = eipsi * Complex(Real(-s/M_SQRT2));
+    D[1](0,2) = emipsi * Complex(Real(0.5*(1.0-c))) * eiphi; 
+    D[1](1,2) = Complex(Real(s/M_SQRT2)) * eiphi;
     D[1](2,2) = eipsi * Complex(Real(0.5*(1.0+c))) * eiphi;
 
     l_max = 1;
@@ -237,8 +230,6 @@ WignerMatrix<Real>::compute_D(int l)
     if (l==1)
     {
         //precompute D0 =1 and D1 = (90 degree rot)
-        // FIXME: signs are inconsistent with above explicit formula for 
-        //        theta = pi/2, phi=0, psi=0 (sine terms should be negated)
         D[1].reshape(MultiArrayShape<2>::type(3,3));
         D[1](0,0) = Real(0.5);
         D[1](0,1) = Real(0.5*M_SQRT2);
@@ -286,28 +277,57 @@ WignerMatrix<Real>::compute_D(int l)
     }
 }
 
+template <class Real>
+void
+WignerMatrix<Real>::rotateSH(std::vector<std::vector< FFTWComplex<Real> > > const  & SH, 
+                             Real phi, Real theta, Real psi, 
+                             std::vector<std::vector< FFTWComplex<Real> > > & SHresult)
+{
+    int band = SH.size()-1;
+    compute_D(band, phi, theta, psi);
+
+    SHresult.resize(SH.size());
+
+    for (int l=0; l<(int)SH.size(); l++)
+    {
+        SHresult[l].resize(SH[l].size());
+        for(int m=-l; m<=l; m++)
+        {
+            Complex tmp = 0;
+            for (int h=-l; h<=l; h++)
+            {
+                tmp += get_d(l,h,m) * SH[l][h+l];
+            }
+
+            SHresult[l][m+l] = tmp;
+        }
+    }
+}
+
+
 template <class Real>
 void 
-WignerMatrix<Real>::rotatePH(NestedArray const & PH, Real phi, Real theta, Real psi,
-                             NestedArray & PHresult)
+WignerMatrix<Real>::rotatePH(std::vector<std::vector<std::vector< FFTWComplex<Real> > > > const  & PH, 
+                             Real phi, Real theta, Real psi,
+                             std::vector<std::vector<std::vector< FFTWComplex<Real> > > > & PHresult)
 {
-    int band = PH[1].size()-1;
+    int band = PH.size()-1;
     compute_D(band, phi, theta, psi);
 
     PHresult.resize(PH.size());
 
-    for(int n=1; n<=band; n++)
+    for(int n=1; n<(int)PH.size(); n++)
     {
-        PHresult[n].resize(band+1);
-        for (int l=0; l<=band; l++)
+        PHresult[n].resize(PH[n].size());
+        for (int l=0; l<(int)PH[n].size(); l++)
         {
-            PHresult[n][l].resize(2*band+1);
+            PHresult[n][l].resize(PH[n][l].size());
             for(int m=-l; m<=l; m++)
             {
                 Complex tmp = 0;
                 for (int h=-l; h<=l; h++)
                 {
-                    tmp += get_D(l,h,m) * PH[n][l][h+l];
+                    tmp += get_d(l,h,m) * PH[n][l][h+l];
                 }
 
                 PHresult[n][l][m+l] = tmp;
diff --git a/test/features/CMakeLists.txt b/test/features/CMakeLists.txt
index 6f2579000..870ce4b92 100644
--- a/test/features/CMakeLists.txt
+++ b/test/features/CMakeLists.txt
@@ -1,7 +1,12 @@
 if(FFTW3_FOUND)
     INCLUDE_DIRECTORIES(${FFTW3_INCLUDE_DIR})
+    
+    IF(HDF5_FOUND)
+        ADD_DEFINITIONS(-DHasHDF5 ${HDF5_CPPFLAGS})
+        INCLUDE_DIRECTORIES(${HDF5_INCLUDE_DIR})
+    ENDIF(HDF5_FOUND)
 
-    VIGRA_ADD_TEST(test_features test.cxx LIBRARIES vigraimpex ${FFTW3_LIBRARIES})
+    VIGRA_ADD_TEST(test_features test.cxx LIBRARIES vigraimpex ${FFTW3_LIBRARIES} ${FFTW3F_LIBRARIES})
 
     # VIGRA_COPY_TEST_DATA(ghouse.gif)
 else()
diff --git a/test/features/test.cxx b/test/features/test.cxx
index bd0163a5a..3628af58c 100644
--- a/test/features/test.cxx
+++ b/test/features/test.cxx
@@ -42,13 +42,16 @@
 #include <vigra/wigner-matrix.hxx>
 #include <vigra/multi_pointoperators.hxx>
 #include "wigner-matrix-reference.hxx"
+#include <vigra/harmonics.hxx>
+#include <vigra/invariant_features3D.hxx>
+#include <vigra/hdf5impex.hxx>
 
 using namespace vigra;
 
 template <class Real>
 struct ComplexRealFunctor
 {
-    Real operator()(std::complex<Real> const & c) const
+    Real operator()(FFTWComplex<Real> const & c) const
     {
         return c.real();
     }
@@ -57,7 +60,7 @@ struct ComplexRealFunctor
 template <class Real>
 struct ComplexImagFunctor
 {
-    Real operator()(std::complex<Real> const & c) const
+    Real operator()(FFTWComplex<Real> const & c) const
     {
         return c.imag();
     }
@@ -66,8 +69,396 @@ struct ComplexImagFunctor
 struct InvariantFeaturesTest
 {
 
-    InvariantFeaturesTest()
+    void SHbaseTest()
     {
+        std::cerr<<"SHbaseTest\n";
+        std::vector<std::vector<vigra::MultiArray<3,FFTWComplex<float> > > > SHbaseF;
+        vigra::computeSHbaseF<float>(10,2,5,SHbaseF);
+
+	//FIXME: this is a qualitative test, find objective criteria 
+        //HDF5File file("SHbaseTest.h5",HDF5File::New);
+	for (int l=0;l<SHbaseF.size();l++)
+	    for (int m=0;m<SHbaseF[l].size();m++)
+	    {
+		std::string name = "r_"+vigra::asString(l)+"_"+vigra::asString(-l+m);
+		vigra::MultiArray<3, float > tmp;
+		tmp.reshape(SHbaseF[l][m].shape(),0);
+		for (int z=0;z<SHbaseF[l][m].shape()[0];++z)
+		    for (int y=0;y<SHbaseF[l][m].shape()[1];++y)
+			for (int x=0;x<SHbaseF[l][m].shape()[2];++x)
+			    tmp(z,y,x)=real(SHbaseF[l][m](z,y,x));
+		//file.write(name, tmp);
+		name = "c_"+vigra::asString(l)+"_"+vigra::asString(-l+m);
+		for (int z=0;z<SHbaseF[l][m].shape()[0];++z)
+		    for (int y=0;y<SHbaseF[l][m].shape()[1];++y)
+			for (int x=0;x<SHbaseF[l][m].shape()[2];++x)
+			    tmp(z,y,x)=imag(SHbaseF[l][m](z,y,x));
+		//file.write(name, tmp);
+
+	    }
+    }
+
+    void SHcenterTest()
+    {
+        std::cerr<<"SHCenterTest\n";
+        typedef MultiArray<3, float>::difference_type Shape;
+        vigra::MultiArrayShape<3>::type testShape(51,51,51);
+        vigra::MultiArray<3,float> testVolume;
+        vigra::MultiArray<3,float> testVolumeRot;
+        testVolume.reshape(testShape,0);
+        testVolumeRot.reshape(testShape,0);
+        vigra::MultiArrayView<3,float> testView = testVolume.subarray(Shape(0,0,0), Shape(25, 25, 25));
+        vigra::MultiArrayView<3,float> testViewRot = testVolumeRot.subarray(Shape(26,0,0), Shape(51, 25, 25));
+        testView+=1;
+        testViewRot+=1;
+
+	std::vector<std::vector<vigra::MultiArray<3,FFTWComplex<float> > > > SHbaseF;
+        vigra::computeSHbaseF<float>(10,2,5,SHbaseF);
+
+        std::vector<std::vector<FFTWComplex<float> >  >  SH_A;
+        std::vector<std::vector<FFTWComplex<float> >  >  SH_A_rot;
+        std::vector<std::vector<FFTWComplex<float> >  >  SH_A_rerot;
+
+        vigra::SHcenter<float>(SH_A, 10, 5, testVolume, SHbaseF);
+
+        //check if all coefficients are real for m=0
+        float cm0=0;
+        for (int l=0;l<SHbaseF.size();l++)
+	    cm0+=imag(SH_A[l][l]);
+        should(abs(cm0) < 1e-4f);
+
+        vigra::MultiArray<3,float> reconstVolume;
+        vigra::MultiArray<3,float> reconstVolumeRot;
+        vigra::MultiArray<3,float> reconstVolumeReRot;
+        reconstVolume.reshape(testShape,0);
+        reconstVolumeRot.reshape(testShape,0);
+        reconstVolumeReRot.reshape(testShape,0);
+
+        vigra::WignerMatrix<float> W(5);
+
+        //check if zero rotation changes coeffs
+        W.rotateSH(SH_A, (float) 0.0, (float) 0.0, (float) 0.0, SH_A_rot);
+        float diff = 0;
+	for (int l=0;l<SHbaseF.size();l++)
+	    for (int m=0;m<SHbaseF[l].size();m++)
+		diff+=abs(SH_A[l][m]-SH_A_rot[l][m]);
+        should(abs(diff) < 1e-4f);
+
+        //check if rotation and unrotation  changes coeffs
+        W.rotateSH(SH_A, (float) M_PI/4, (float) M_PI/8, (float) M_PI/2, SH_A_rot);
+	W.rotateSH(SH_A_rot, (float) -M_PI/2, (float) -M_PI/8, (float) -M_PI/4, SH_A_rerot);
+        diff = 0;
+	for (int l=0;l<SHbaseF.size();l++)
+	    for (int m=0;m<SHbaseF[l].size();m++)
+	    {
+		//std::cerr<<l<<" "<<m-l<<" "<<SH_A[l][m]<<" "<<SH_A_rerot[l][m]<<"\n ";
+		diff+=abs(SH_A[l][m]-SH_A_rerot[l][m]);
+	    }
+        should(abs(diff) < 1e-2f);
+
+	
+        //check for coeff diff after rot in spacial domain and unrot in harmonics
+        vigra::SHcenter<float>(SH_A_rot, 10, 5, testVolumeRot, SHbaseF);
+        W.rotateSH(SH_A, (float) M_PI/2, (float) 0.0, (float) 0, SH_A_rerot);
+        diff = 0;
+	for (int l=0;l<SHbaseF.size();l++)
+	    for (int m=0;m<SHbaseF[l].size();m++)
+	    {
+		//std::cerr<<l<<" "<<m-l<<" "<<SH_A_rot[l][m]<<" "<<SH_A_rerot[l][m]<<"\n ";
+		diff+=abs(SH_A_rot[l][m]-SH_A_rerot[l][m]);
+	    }
+	std::cerr<<"diff:"<<abs(diff)<<"\n";
+        //should(abs(diff) < 1);//FIXME
+
+        vigra::reconstSH<float>(10,2,5,reconstVolume,SH_A,SHbaseF);
+        vigra::reconstSH<float>(10,2,5,reconstVolumeRot,SH_A_rot,SHbaseF);
+        vigra::reconstSH<float>(10,2,5,reconstVolumeReRot,SH_A_rerot,SHbaseF);
+
+        //HDF5File file("SHreconstTest.h5",HDF5File::New);
+        //file.write("reconst", reconstVolume);
+        //file.write("reconstRot", reconstVolumeRot);
+        //file.write("reconstReRot", reconstVolumeReRot);
+        //file.write("orig",testVolume);
+    }
+
+
+    struct pow2f
+    {
+	float operator()(float t) const
+	{
+	    return pow( (double)t,2.0);
+	}
+    };
+
+    void SHfeatureTest()
+    {
+        std::cerr<<"SHCenterTest\n";
+        typedef MultiArray<3, float>::difference_type Shape;
+        vigra::MultiArrayShape<3>::type testShape(51,51,51);
+        vigra::MultiArray<3,float> testVolume;
+        vigra::MultiArray<3,float> testVolumeRot;
+        testVolume.reshape(testShape,0);
+        testVolumeRot.reshape(testShape,0);
+        vigra::MultiArrayView<3,float> testView = testVolume.subarray(Shape(0,0,0), Shape(25, 25, 25));
+        vigra::MultiArrayView<3,float> testViewRot = testVolumeRot.subarray(Shape(26,0,0), Shape(51, 25, 25));
+        testView+=1;
+        testViewRot+=1;
+
+        std::vector<std::vector<vigra::MultiArray<3,FFTWComplex<float> > > > SHbaseF;
+        vigra::computeSHbaseF<float>(10,2,5,SHbaseF);
+
+        std::vector<std::vector<FFTWComplex<float> >  >  SH_A;
+	std::vector<std::vector<FFTWComplex<float> >  >  SH_A_rot;
+        std::vector<std::vector<FFTWComplex<float> >  >  SH_A_real_rot;
+
+        vigra::SHcenter<float>(SH_A, 10, 5, testVolume, SHbaseF);
+	vigra::WignerMatrix<float> W(5);
+	W.rotateSH(SH_A, (float) M_PI/4, (float) M_PI/8, (float) M_PI/2, SH_A_rot);
+	vigra::SHcenter<float>(SH_A_real_rot, 10, 5, testVolumeRot, SHbaseF);
+
+	//test if SHabs is the same 
+	std::vector<float> res_A, res_A_rot, res_A_real_rot;
+	vigra::SHabs(SH_A, res_A);
+	vigra::SHabs(SH_A_rot, res_A_rot);
+	vigra::SHabs(SH_A_real_rot, res_A_real_rot);
+	float diff1=0,diff2=0;
+	for (int l=0;l<res_A.size();l++)
+	{
+	    //std::cerr<<res_A[l]<<" "<<res_A_rot[l]<<" "<<res_A_real_rot[l]<<"\n";
+	    diff1+=abs(res_A[l]-res_A_rot[l]);
+	    diff2+=abs(res_A[l]-res_A_real_rot[l]);
+	}
+	should(abs(diff1) < 1e-1f);//should be zero, but numerical problems
+	should(abs(diff2) < 1);//real rotations not as accurate (due to band limitation)
+
+
+        //test if SHbispec is the same 
+        std::vector<FFTWComplex<float> > resC_A, resC_A_rot, resC_A_real_rot;
+	vigra::SHbispec(SH_A, resC_A);
+        vigra::SHbispec(SH_A_rot, resC_A_rot);
+        vigra::SHbispec(SH_A_real_rot, resC_A_real_rot);
+        diff1=0;
+	diff2=0;
+        for (int l=0;l<resC_A.size();l++)
+        {
+            //std::cerr<<resC_A[l]<<" "<<resC_A_rot[l]<<" "<<resC_A_real_rot[l]<<"\n";
+            diff1+=abs(resC_A[l]-resC_A_rot[l]);
+            diff2+=abs(resC_A[l]-resC_A_real_rot[l]);
+        }
+        //std::cerr<<abs(diff1)/resC_A.size()<<" "<<abs(diff2)/resC_A.size()<<"\n";
+	should(abs(diff1)/resC_A.size() < 1);//should be zero, but numerical problems
+        should(abs(diff2)/resC_A.size() < 100000);//feature very sensitive: large values for real data 
+
+	//test SHcorr
+        //check if normalzed cross correlation has max=1 at zero rotation offset
+        float cmax,phi,theta,psi=0;
+        vigra::MultiArray<3,float > corMat;
+        SHxcorr(SH_A,SH_A,64,true,cmax,phi,theta,psi, corMat);
+        should( (abs(phi+psi)<0.1 || abs(phi+psi-2*M_PI) < 0.1) && (theta<0.1) && (abs(cmax-1.0)< 1e-4f) ); //note: at theta=0 abitrary combinations of phi+psi%2PI are possible
+        //std::cerr<<"xcorr: "<<cmax<<" "<<phi<<" "<<theta<<" "<<psi<<" : "<<phi+psi<<"\n";
+
+	//check if normalzed cross correlation predicts correct correlation offset
+	W.rotateSH(SH_A, (float) M_PI/2, (float) M_PI/2, (float) M_PI/2, SH_A_rot);
+	SHxcorr(SH_A,SH_A_rot,64,true,cmax,phi,theta,psi, corMat);
+	//std::cerr<<"xcorr: "<<cmax<<" "<<phi<<" "<<theta<<" "<<psi<<"\n";
+	should( (abs(phi-M_PI/2)<1e-1f) && (abs(theta-M_PI/2)<1e-1f) && (abs(psi-M_PI/2)<1e-1f) && (abs(cmax-1)<1e-1f) ); 
+
+	//HDF5File corrfile("SHCorTest.h5",HDF5File::New);
+        //corrfile.write("cormat", corMat);
+	
+	//test invariance of SHautocorr
+	float res1,res2;
+	pow2f myfunctor;
+	SHautocorr(SH_A, 32, myfunctor, res1);
+	SHautocorr(SH_A_rot, 32, myfunctor, res2);
+	should( abs(res1-res2)<1e-3f);
+    }
+
+    void SHallTest()
+    {
+        std::cerr<<"SHallTest\n";
+        typedef MultiArray<3, float>::difference_type Shape;
+        vigra::MultiArrayShape<3>::type testShape(51,51,51);
+        vigra::MultiArray<3,float> testVolume;
+        testVolume.reshape(testShape,0);
+        vigra::MultiArrayView<3,float> testView = testVolume.subarray(Shape(0,0,0), Shape(25, 25, 25));
+        testView+=1;
+
+        std::vector<std::vector<vigra::MultiArray<3,FFTWComplex<float> > > > SHbaseF;
+        vigra::computeSHbaseF<float>(10,2,5,SHbaseF);
+
+        std::vector<std::vector<FFTWComplex<float> >  >  SH_A;
+        vigra::SHcenter<float>(SH_A, 10, 5, testVolume, SHbaseF);
+
+	std::vector<std::vector<vigra::MultiArray<3,FFTWComplex<float> >  > > SH_All;
+	vigra::Array2SH<float>(SH_All,5, 10, SHbaseF, testVolume);
+	
+	//test if computation in fourier space equals computation in spacial domain
+	float diff = 0;
+        for (int l=0;l<SHbaseF.size();l++)
+            for (int m=0;m<SHbaseF[l].size();m++)
+            {
+                //std::cerr<<l<<" "<<m-l<<" "<<SH_A[l][m]<<" "<<SH_All[l][m](25,25,25)<<"\n ";
+                diff+=abs(SH_A[l][m]-SH_All[l][m](25,25,25));
+            }
+        should(abs(diff) < 1e-2f);
+
+	//test SHabs all
+	std::vector<vigra::MultiArray<3,float> > res_All;
+	SHabs(SH_All,res_All);
+	std::vector<float> res_A;
+	vigra::SHabs(SH_A, res_A);
+	diff=0;
+	for (int i=0;i<SH_A.size();i++)
+	    diff+=abs(res_A[i]-res_All[i](25,25,25));
+	should(abs(diff) < 1);
+
+	//test SHbispec all;
+	std::vector<vigra::MultiArray<3,FFTWComplex<float> > > res_All2;
+	SHbispec(SH_All,res_All2);
+	std::vector<FFTWComplex<float> > res_A2;
+	SHbispec(SH_A,res_A2);
+	diff=0;
+	for (int i=0;i<SH_A.size();i++)
+	    diff+=abs(res_A[i]-res_All[i](25,25,25));
+	should(abs(diff) < 1);
+
+
+    }
+
+    void PHbaseTest()
+    {
+	std::cerr<<"PHbaseTest\n";
+	std::vector<std::vector<std::vector<vigra::MultiArray<3,FFTWComplex<float> > > > > PHbaseF;
+	vigra::computePHbaseF<float>(10,5,PHbaseF,false);
+	
+	//FIXME: this is a qualitative test, find objective criteria
+	//HDF5File file("PHbaseTest.h5",HDF5File::New);
+	for (int n=0;n<PHbaseF.size();n++)
+	    for (int l=0;l<PHbaseF[n].size();l++)
+		for (int m=0;m<PHbaseF[n][l].size();m++)
+		{
+		    std::string name = "r_"+vigra::asString(n)+"_"+vigra::asString(l)+"_"+vigra::asString(-l+m);
+		    //std::cerr<<name<<" "<<PHbaseF[n][l][m].shape()<<" "<<k<<"\n";
+		    vigra::MultiArray<3, float > tmp;
+		    tmp.reshape(PHbaseF[n][l][m].shape(),0);
+		    for (int z=0;z<PHbaseF[n][l][m].shape()[0];++z)
+			for (int y=0;y<PHbaseF[n][l][m].shape()[1];++y)
+			    for (int x=0;x<PHbaseF[n][l][m].shape()[2];++x)
+				tmp(z,y,x)=real(PHbaseF[n][l][m](z,y,x));
+		    //file.write(name, tmp);
+		    name = "c_"+vigra::asString(n)+"_"+vigra::asString(l)+"_"+vigra::asString(-l+m);
+		    for (int z=0;z<PHbaseF[n][l][m].shape()[0];++z)
+		        for (int y=0;y<PHbaseF[n][l][m].shape()[1];++y)
+		            for (int x=0;x<PHbaseF[n][l][m].shape()[2];++x)
+		                tmp(z,y,x)=imag(PHbaseF[n][l][m](z,y,x));
+		    //file.write(name, tmp);
+
+		}
+    }
+
+    void PHcenterTest()
+    {
+	std::cerr<<"PHCenterTest\n";
+	typedef MultiArray<3, float>::difference_type Shape;
+	vigra::MultiArrayShape<3>::type testShape(51,51,51);
+	vigra::MultiArray<3,float> testVolume;
+	vigra::MultiArray<3,float> testVolumeRot;
+	testVolume.reshape(testShape,0);
+	testVolumeRot.reshape(testShape,0);
+	vigra::MultiArrayView<3,float> testView = testVolume.subarray(Shape(0,0,0), Shape(25, 25, 25));
+	vigra::MultiArrayView<3,float> testViewRot = testVolumeRot.subarray(Shape(26,0,0), Shape(51, 25, 25));
+	testView+=1;
+	testViewRot+=1;
+
+	std::vector<std::vector<std::vector<vigra::MultiArray<3,FFTWComplex<float> > > > > PHbaseF;
+	vigra::computePHbaseF<float>(10,5,PHbaseF);
+
+	std::vector<std::vector<std::vector<FFTWComplex<float> >  > > PH_A;
+	std::vector<std::vector<std::vector<FFTWComplex<float> >  > > PH_A_rot;
+        std::vector<std::vector<std::vector<FFTWComplex<float> >  > > PH_A_rerot;
+
+	vigra::PHcenter<float>(PH_A, 10, 5, testVolume, PHbaseF);
+	
+	//check if all coefficients are real for m=0
+	float cm0=0;
+	for (int n=1;n<PHbaseF.size();n++)
+            for (int l=0;l<PHbaseF[n].size();l++)
+		    cm0+=imag(PH_A[n][l][l]);
+	should(abs(cm0) < 1e-4f);
+
+	vigra::MultiArray<3,float> reconstVolume;
+	vigra::MultiArray<3,float> reconstVolumeRot;
+	vigra::MultiArray<3,float> reconstVolumeReRot;
+	reconstVolume.reshape(testShape,0);
+	reconstVolumeRot.reshape(testShape,0);
+	reconstVolumeReRot.reshape(testShape,0);
+	
+	vigra::WignerMatrix<float> W(5);
+	
+	//check if zero rotation changes coeffs
+	W.rotatePH(PH_A, (float) 0.0, (float) 0.0, (float) 0.0, PH_A_rot);
+        float diff = 0;
+        for (int n=1;n<PHbaseF.size();n++)
+            for (int l=0;l<PHbaseF[n].size();l++)
+                for (int m=0;m<PHbaseF[n][l].size();m++)
+                    diff+=abs(PH_A[n][l][m]-PH_A_rot[n][l][m]);
+	should(abs(diff) < 1e-4f);
+
+        //check if rotation keeps all coefficients real for l=0, m=0
+        W.rotatePH(PH_A, (float) M_PI/4, (float) M_PI/8, (float) M_PI/2, PH_A_rot);
+        cm0=0;
+        for (int n=1;n<PHbaseF.size();n++)
+	{
+	    //std::cerr<<n<<" "<<PH_A_rot[n][0][0]<<"\n ";
+            cm0+=imag(PH_A_rot[n][0][0]);
+	}
+        should(abs(cm0) < 1e-4f);
+
+
+	//check if rotation and unrotation  changes coeffs
+	W.rotatePH(PH_A_rot, (float) -M_PI/2, (float) -M_PI/8, (float) -M_PI/4, PH_A_rerot);
+        diff = 0;
+        for (int n=1;n<PHbaseF.size();n++)
+            for (int l=0;l<PHbaseF[n].size();l++)
+                for (int m=0;m<PHbaseF[n][l].size();m++)
+                {
+                    //std::cerr<<n<<" "<<l<<" "<<m-l<<" "<<PH_A_rot[n][l][m]<<" "<<PH_A_rerot[n][l][m]<<"\n ";
+                    diff+=abs(PH_A[n][l][m]-PH_A_rerot[n][l][m]);
+                }
+        should(abs(diff) < 1e-4f);
+
+	//check for coeff diff after rot in spacial domain and unrot un harmonics
+	vigra::PHcenter<float>(PH_A_rot, 10, 5, testVolumeRot, PHbaseF);
+	W.rotatePH(PH_A_rot, (float) 0.0, (float) M_PI/2, (float) 0.0, PH_A_rerot);
+	diff = 0;
+        for (int n=1;n<PHbaseF.size();n++)
+            for (int l=0;l<PHbaseF[n].size();l++)
+                for (int m=0;m<PHbaseF[n][l].size();m++)
+                {
+                    //std::cerr<<n<<" "<<l<<" "<<m-l<<" "<<PH_A[n][l][m]<<" "<<PH_A_rerot[n][l][m]<<"\n ";
+                    diff+=abs(PH_A[n][l][m]-PH_A_rerot[n][l][m]);
+                }
+        //should(abs(diff) < 1);
+	//FIXME
+
+	vigra::reconstPH<float>(10,5,reconstVolume,PH_A,PHbaseF);
+	vigra::reconstPH<float>(10,5,reconstVolumeRot,PH_A_rot,PHbaseF);
+	vigra::reconstPH<float>(10,5,reconstVolumeReRot,PH_A_rerot,PHbaseF);
+	
+	//HDF5File file("PHreconstTest.h5",HDF5File::New);
+	//file.write("reconst", reconstVolume);
+	//file.write("reconstRot", reconstVolumeRot);
+	//file.write("reconstReRot", reconstVolumeReRot);
+	//file.write("orig",testVolume);
+	
+	
+    }
+
+
+    void PHallTest()
+    {
+
     }
 
     void wignerMatrixTest()
@@ -104,15 +495,13 @@ struct InvariantFeaturesTest
             M diff(2*l+1, 2*l+1);
             FindMinMax<float> minmax;
 
-            transformMultiArray(srcMultiArrayRange(wigner.get_D(l)),
-                                destMultiArray(diff), ComplexImagFunctor<float>());
+	    transformMultiArray(srcMultiArrayRange(wigner.get_D(l)), destMultiArray(diff), ComplexImagFunctor<float>());
             inspectMultiArray(srcMultiArrayRange(diff), minmax);
             shouldEqual(minmax.min, 0.0f);
             shouldEqual(minmax.max, 0.0f);
 
-            transformMultiArray(srcMultiArrayRange(wigner.get_D(l)),
-                                destMultiArray(diff), ComplexRealFunctor<float>());
-            diff -= ref[l];
+	    transformMultiArray(srcMultiArrayRange(wigner.get_D(l)), destMultiArray(diff), ComplexRealFunctor<float>());
+	    diff -= ref[l];
             inspectMultiArray(srcMultiArrayRange(diff), minmax);
             should(minmax.min > -1e-4f);
             should(minmax.max <  1e-4f);
@@ -128,22 +517,20 @@ struct InvariantFeaturesTest
             M diff(2*l+1, 2*l+1);
             FindMinMax<float> minmax;
 
-            transformMultiArray(srcMultiArrayRange(wigner.get_D(l)),
-                                destMultiArray(diff), ComplexImagFunctor<float>());
+			transformMultiArray(srcMultiArrayRange(wigner.get_D(l)),
+				                destMultiArray(diff), ComplexImagFunctor<float>());
             inspectMultiArray(srcMultiArrayRange(diff), minmax);
             shouldEqual(minmax.min, 0.0f);
             shouldEqual(minmax.max, 0.0f);
             
-            // FIXME: transpose() shouldn't be necessary below
-            transformMultiArray(srcMultiArrayRange(transpose(wigner2.get_D(l))),
-                                destMultiArray(diff), ComplexRealFunctor<float>());
-            diff -= ref[l];
-            inspectMultiArray(srcMultiArrayRange(diff), minmax);
-            should(minmax.min > -1e-4f);
+	    transformMultiArray(srcMultiArrayRange((wigner2.get_D(l))), destMultiArray(diff), ComplexRealFunctor<float>());
+	    diff -= ref[l];
+	    inspectMultiArray(srcMultiArrayRange(diff), minmax);
+	    should(minmax.min > -1e-4f);
             should(minmax.max <  1e-4f);
         }
         
-        // FIXME: compute_D() with arbitrary angles, rot(), and rotatePH() are untested!
+        // FIXME: compute_D() with arbitrary angles, rot()
     }
 
 };
@@ -155,6 +542,12 @@ struct FeaturesTestSuite
     : vigra::test_suite("FeaturesTestSuite")
     {
         add( testCase( &InvariantFeaturesTest::wignerMatrixTest));
+	add( testCase( &InvariantFeaturesTest::SHbaseTest));
+	add( testCase( &InvariantFeaturesTest::SHcenterTest));
+	add( testCase( &InvariantFeaturesTest::SHallTest));
+	add( testCase( &InvariantFeaturesTest::SHfeatureTest));
+	add( testCase( &InvariantFeaturesTest::PHbaseTest));
+	add( testCase( &InvariantFeaturesTest::PHcenterTest));
     }
 };