ProSHADE_spharmonickit.cpp

Go to the documentation of this file.

//============================================ FFTW3 + SOFT
#ifdef __cplusplus
extern "C" {
#endif
#include <fftw3.h>
#include <wrap_fftw.h>
#include <makeweights.h>
#include <s2_primitive.h>
#include <s2_cospmls.h>
#include <s2_legendreTransforms.h>
#include <s2_semi_fly.h>
#include <rotate_so3_utils.h>
#include <utils_so3.h>
#include <soft_fftw.h>
#include <rotate_so3_fftw.h>
#ifdef __cplusplus
}
#endif

//============================================ RVAPI
#include <rvapi_interface.h>
    
//============================================ ProSHADE
#include "ProSHADE.h"
#include "ProSHADE_internal.h"

void ProSHADE_internal::ProSHADE_data::mapPhaselessToSphere ( ProSHADE::ProSHADE_settings* settings,
                                                              double theta,
                                                              double phi,
                                                              double shellSz,
                                                              unsigned int manualShells,
                                                              bool keepInMemory,
                                                              bool rotDefaults )
{
    //======================================== Sanity check
    if ( !this->_phaseRemoved )
    {
        std::cerr << "!!! ProSHADE ERROR !!! Error in file " << this->_inputFileName << " !!! Attempted to map data onto shell before computing the phaseless data. Call the removePhaseFromMap function BEFORE the mapPhaselessToSphere function." << std::endl;
        
        if ( settings->htmlReport )
        {
            std::stringstream hlpSS;
            hlpSS << "<font color=\"red\">" << "Attempted to map data onto shell before computing the phaseless data. This looks like an internal bug, please report this case." << "</font>";
            rvapi_set_text                    ( hlpSS.str().c_str(),
                                               "ProgressSection",
                                               settings->htmlReportLineProgress,
                                               1,
                                               1,
                                               1 );
            settings->htmlReportLineProgress += 1;
            rvapi_flush                       ( );
        }
        
        exit ( -1 );
    }
    
    //======================================== Initialise internal values
    this->_thetaAngle                         = theta;
    this->_phiAngle                           = phi;
    this->_xSamplingRate                      = this->_xRange / static_cast<double> ( this->_maxMapU );
    this->_ySamplingRate                      = this->_yRange / static_cast<double> ( this->_maxMapV );
    this->_zSamplingRate                      = this->_zRange / static_cast<double> ( this->_maxMapW );
    
    //======================================== Initialise local variables
    double x, y, z;
    int xBottom, yBottom, zBottom;
    int xTop, yTop, zTop;
    
    std::array<double,4> c000, c001, c010, c011, c100, c101, c110, c111;
    std::array<double,4> c00, c01, c10, c11;
    std::array<double,4> c0, c1;
    double xd, yd, zd;
    
    //======================================== How many shells can the map fill?
    if ( this->_noShellsWithData == 0 )
    {
        //==================================== Automatically determine maximum number of shells
        this->_shellSpacing                   = shellSz;
        this->_shellPlacement                 = std::vector<double> ( std::floor ( this->_maxMapRange / this->_shellSpacing ) );
        
        for ( unsigned int shellIter = 0; shellIter < static_cast<unsigned int> ( this->_shellPlacement.size() ); shellIter++ )
        {
            this->_shellPlacement.at(shellIter)   = (shellIter+1) * this->_shellSpacing;
        }
    }
    
    if ( !rotDefaults )
    {
        unsigned int shellNoHlp = 0;
        for ( unsigned int sh = 0; sh < static_cast<unsigned int> ( this->_shellPlacement.size() ); sh++ )
        {
            if ( ( ( (((this->_maxMapU+1)/2)-1) * this->_xSamplingRate ) > this->_shellPlacement.at(sh) ) ||
                 ( ( (((this->_maxMapV+1)/2)-1) * this->_ySamplingRate ) > this->_shellPlacement.at(sh) ) ||
                 ( ( (((this->_maxMapW+1)/2)-1) * this->_zSamplingRate ) > this->_shellPlacement.at(sh) ) )
            { shellNoHlp = sh; }
            else
            { break; }
        }
        if ( shellNoHlp > 0 ) { shellNoHlp -= 1; }
        
        if ( manualShells == 0 ) { this->_noShellsWithData = shellNoHlp; }
        else
        {
            this->_noShellsWithData = manualShells;
            if ( manualShells > shellNoHlp )
            {
                std::cerr << "!!! ProSHADE ERROR !!! The required number of shells is larger than the data allows. Either use automatic shell number determination (i.e. manualShells = 0) or decrease the number of required shells." << std::endl;
                
                if ( settings->htmlReport )
                {
                    std::stringstream hlpSS;
                    hlpSS << "<font color=\"red\">" << "The required number of shells is larger than the data allows. Either use automatic shell number determination (i.e. manualShells = 0) or decrease the number of required shells." << "</font>";
                    rvapi_set_text                    ( hlpSS.str().c_str(),
                                                       "ProgressSection",
                                                       settings->htmlReportLineProgress,
                                                       1,
                                                       1,
                                                       1 );
                    settings->htmlReportLineProgress += 1;
                    rvapi_flush                       ( );
                }
                
                exit ( -1 );
            }
        }
    }
    else
    {
        //==================================== Make sure shells cover the corners
        unsigned int shellNoHlp               = 0;
        double maxDist                        = std::max ( sqrt( pow( ( (((this->_maxMapU+1)/2)-1) * this->_xSamplingRate ), 2.0 ) + pow ( ( (((this->_maxMapV+1)/2)-1) * this->_ySamplingRate ), 2.0 ) ),
                                                           std::max ( sqrt( pow( ( (((this->_maxMapU+1)/2)-1) * this->_xSamplingRate ), 2.0 ) + pow ( ( (((this->_maxMapW+1)/2)-1) * this->_zSamplingRate ), 2.0 ) ),
                                                                      sqrt( pow( ( (((this->_maxMapW+1)/2)-1) * this->_zSamplingRate ), 2.0 ) + pow ( ( (((this->_maxMapV+1)/2)-1) * this->_ySamplingRate ), 2.0 ) ) ) );
        for ( unsigned int sh = 0; sh < static_cast<unsigned int> ( this->_shellPlacement.size() ); sh++ )
        {
            if ( maxDist > this->_shellPlacement.at(sh) ) { shellNoHlp = sh; }
            else { break; }
        }
        shellNoHlp                           -= 1;
        
        if ( manualShells == 0 ) { this->_noShellsWithData = shellNoHlp; }
        else
        {
            this->_noShellsWithData = manualShells;
            if ( manualShells > shellNoHlp )
            {
                std::cerr << "!!! ProSHADE ERROR !!! The required number of shells is larger than the data allows. Either use automatic shell number determination (i.e. manualShells = 0) or decrease the number of required shells." << std::endl;
                
                if ( settings->htmlReport )
                {
                    std::stringstream hlpSS;
                    hlpSS << "<font color=\"red\">" << "The required number of shells is larger than the data allows. Either use automatic shell number determination (i.e. manualShells = 0) or decrease the number of required shells." << "</font>";
                    rvapi_set_text                    ( hlpSS.str().c_str(),
                                                       "ProgressSection",
                                                       settings->htmlReportLineProgress,
                                                       1,
                                                       1,
                                                       1 );
                    settings->htmlReportLineProgress += 1;
                    rvapi_flush                       ( );
                }
                
                exit ( -1 );
            }
        }
    }
    
    if ( this->_noShellsWithData < 2 )
    {
        std::cerr << "!!! ProSHADE ERROR !!! Did not find enough shells to which the data could be mapped. Either the structure is rather small and the default values do not work - use the \'sphDistance\' option to adjust, or this is a bug. Terminating..." << std::endl;
        
        if ( settings->htmlReport )
        {
            std::stringstream hlpSS;
            hlpSS << "<font color=\"red\">" << "Did not find enough shells to which the data could be mapped. Either the structure is rather small and the default values do not work - use the \'sphDistance\' option to adjust." << "</font>";
            rvapi_set_text                    ( hlpSS.str().c_str(),
                                               "ProgressSection",
                                               settings->htmlReportLineProgress,
                                               1,
                                               1,
                                               1 );
            settings->htmlReportLineProgress += 1;
            rvapi_flush                       ( );
        }
        
        exit ( -1 );
    }
    
    //======================================== Initialise internal values ( and allocate memory )
    this->_shellMappedData                    = new double*[this->_noShellsWithData];
    for ( unsigned int sh = 0; sh < this->_noShellsWithData; sh++ )
    {
        this->_shellMappedData[sh]            = new double[static_cast<unsigned int> ( this->_thetaAngle * this->_phiAngle )];
    }
    
    //======================================== Find pixelisation cutOffs
    std::vector<double> lonCO ( this->_thetaAngle + 1 );
    for ( unsigned int iter = 0; iter <= this->_thetaAngle; iter++ ) { lonCO.at(iter) = static_cast<double> ( iter ) * ( ( static_cast<double> ( M_PI ) * 2.0 ) / static_cast<double> ( this->_thetaAngle ) ) - ( static_cast<double> ( M_PI ) ); }
    std::vector<double> latCO ( this->_phiAngle + 1 );
    for ( unsigned int iter = 0; iter <= this->_phiAngle; iter++ ) { latCO.at(iter) = ( static_cast<double> ( iter ) * ( static_cast<double> ( M_PI ) / static_cast<double> ( this->_phiAngle ) ) - ( static_cast<double> ( M_PI ) / 2.0 ) ); }
    
    //======================================== Tri-linear interpolation
    unsigned int posIter;
    
    for ( unsigned int sh = 0; sh < this->_noShellsWithData; sh++ )
    {
        //==================================== Find pixel centres
        for ( unsigned int thIt = 0; thIt < this->_thetaAngle; thIt++ )
        {
            for ( unsigned int phIt = 0; phIt < this->_phiAngle; phIt++ )
            {
                //============================ Get grid point x, y and z
                x                             = this->_shellPlacement.at(sh) * cos( ( static_cast<double> ( lonCO.at(thIt) ) + static_cast<double> ( lonCO.at(thIt+1) ) ) / 2.0 ) * cos( ( static_cast<double> ( latCO.at(phIt) ) + static_cast<double> ( latCO.at (phIt+1) ) ) / 2.0 );
                y                             = this->_shellPlacement.at(sh) * sin( ( static_cast<double> ( lonCO.at(thIt) ) + static_cast<double> ( lonCO.at(thIt+1) ) ) / 2.0 ) * cos( ( static_cast<double> ( latCO.at(phIt) ) + static_cast<double> ( latCO.at (phIt+1) ) ) / 2.0 );
                z                             = this->_shellPlacement.at(sh) * sin( ( static_cast<double> ( latCO.at(phIt) ) + static_cast<double> ( latCO.at(phIt+1) ) ) / 2.0 );
                
                //============================ Find 8 closest point around the grid point
                xBottom                       = floor ( (x / this->_xSamplingRate) ) + ((this->_maxMapU+1)/2);
                yBottom                       = floor ( (y / this->_ySamplingRate) ) + ((this->_maxMapV+1)/2);
                zBottom                       = floor ( (z / this->_zSamplingRate) ) + ((this->_maxMapW+1)/2);
                
                xTop                          = xBottom + 1;
                yTop                          = yBottom + 1;
                zTop                          = zBottom + 1;
                
                posIter                       = zBottom  + (this->_maxMapW + 1) * ( yBottom  + (this->_maxMapV + 1) * xBottom );
                if ( posIter > ( (this->_maxMapU+1) * (this->_maxMapV+1) * (this->_maxMapW+1) ) ) { this->_shellMappedData[sh][static_cast<unsigned int> ( phIt * this->_thetaAngle + thIt)] = 0.0; continue; }
                c000[0]                       = (this->_densityMapCorCoords[posIter])[0] * this->_xSamplingRate;
                c000[1]                       = (this->_densityMapCorCoords[posIter])[1] * this->_ySamplingRate;
                c000[2]                       = (this->_densityMapCorCoords[posIter])[2] * this->_zSamplingRate;
                c000[3]                       = this->_densityMapCor[posIter];
                
                posIter                       = zTop  + (this->_maxMapW + 1) * ( yBottom  + (this->_maxMapV + 1) * xBottom );
                if ( posIter > ( (this->_maxMapU+1) * (this->_maxMapV+1) * (this->_maxMapW+1) ) ) { this->_shellMappedData[sh][static_cast<unsigned int> ( phIt * this->_thetaAngle + thIt)] = 0.0; continue; }
                c001[0]                       = (this->_densityMapCorCoords[posIter])[0] * this->_xSamplingRate;
                c001[1]                       = (this->_densityMapCorCoords[posIter])[1] * this->_ySamplingRate;
                c001[2]                       = (this->_densityMapCorCoords[posIter])[2] * this->_zSamplingRate;
                c001[3]                       = this->_densityMapCor[posIter];
                
                posIter                       = zBottom  + (this->_maxMapW + 1) * ( yTop  + (this->_maxMapV + 1) * xBottom );
                if ( posIter > ( (this->_maxMapU+1) * (this->_maxMapV+1) * (this->_maxMapW+1) ) ) { this->_shellMappedData[sh][static_cast<unsigned int> ( phIt * this->_thetaAngle + thIt)] = 0.0; continue; }
                c010[0]                       = (this->_densityMapCorCoords[posIter])[0] * this->_xSamplingRate;
                c010[1]                       = (this->_densityMapCorCoords[posIter])[1] * this->_ySamplingRate;
                c010[2]                       = (this->_densityMapCorCoords[posIter])[2] * this->_zSamplingRate;
                c010[3]                       = this->_densityMapCor[posIter];
                
                posIter                       = zTop  + (this->_maxMapW + 1) * ( yTop  + (this->_maxMapV + 1) * xBottom );
                if ( posIter > ( (this->_maxMapU+1) * (this->_maxMapV+1) * (this->_maxMapW+1) ) ) { this->_shellMappedData[sh][static_cast<unsigned int> ( phIt * this->_thetaAngle + thIt)] = 0.0; continue; }
                c011[0]                       = (this->_densityMapCorCoords[posIter])[0] * this->_xSamplingRate;
                c011[1]                       = (this->_densityMapCorCoords[posIter])[1] * this->_ySamplingRate;
                c011[2]                       = (this->_densityMapCorCoords[posIter])[2] * this->_zSamplingRate;
                c011[3]                       = this->_densityMapCor[posIter];
                
                posIter                       = zBottom  + (this->_maxMapW + 1) * ( yBottom  + (this->_maxMapV + 1) * xTop );
                if ( posIter > ( (this->_maxMapU+1) * (this->_maxMapV+1) * (this->_maxMapW+1) ) ) { this->_shellMappedData[sh][static_cast<unsigned int> ( phIt * this->_thetaAngle + thIt)] = 0.0; continue; }
                c100[0]                       = (this->_densityMapCorCoords[posIter])[0] * this->_xSamplingRate;
                c100[1]                       = (this->_densityMapCorCoords[posIter])[1] * this->_ySamplingRate;
                c100[2]                       = (this->_densityMapCorCoords[posIter])[2] * this->_zSamplingRate;
                c100[3]                       = this->_densityMapCor[posIter];
                
                posIter                       = zTop  + (this->_maxMapW + 1) * ( yBottom  + (this->_maxMapV + 1) * xTop );
                if ( posIter > ( (this->_maxMapU+1) * (this->_maxMapV+1) * (this->_maxMapW+1) ) ) { this->_shellMappedData[sh][static_cast<unsigned int> ( phIt * this->_thetaAngle + thIt)] = 0.0; continue; }
                c101[0]                       = (this->_densityMapCorCoords[posIter])[0] * this->_xSamplingRate;
                c101[1]                       = (this->_densityMapCorCoords[posIter])[1] * this->_ySamplingRate;
                c101[2]                       = (this->_densityMapCorCoords[posIter])[2] * this->_zSamplingRate;
                c101[3]                       = this->_densityMapCor[posIter];
                
                posIter                       = zBottom  + (this->_maxMapW + 1) * ( yTop  + (this->_maxMapV + 1) * xTop );
                if ( posIter > ( (this->_maxMapU+1) * (this->_maxMapV+1) * (this->_maxMapW+1) ) ) { this->_shellMappedData[sh][static_cast<unsigned int> ( phIt * this->_thetaAngle + thIt)] = 0.0; continue; }
                c110[0]                       = (this->_densityMapCorCoords[posIter])[0] * this->_xSamplingRate;
                c110[1]                       = (this->_densityMapCorCoords[posIter])[1] * this->_ySamplingRate;
                c110[2]                       = (this->_densityMapCorCoords[posIter])[2] * this->_zSamplingRate;
                c110[3]                       = this->_densityMapCor[posIter];
                
                posIter                       = zTop + (this->_maxMapW + 1) * ( yTop  + (this->_maxMapV + 1) * xTop );
                if ( posIter > ( (this->_maxMapU+1) * (this->_maxMapV+1) * (this->_maxMapW+1) ) ) { this->_shellMappedData[sh][static_cast<unsigned int> ( phIt * this->_thetaAngle + thIt)] = 0.0; continue; }
                c111[0]                       = (this->_densityMapCorCoords[posIter])[0] * this->_xSamplingRate;
                c111[1]                       = (this->_densityMapCorCoords[posIter])[1] * this->_ySamplingRate;
                c111[2]                       = (this->_densityMapCorCoords[posIter])[2] * this->_zSamplingRate;
                c111[3]                       = this->_densityMapCor[posIter];
                
                //============================ Interpolate along x
                xd                            = (x - ( ( xBottom - static_cast<int> ( ((this->_maxMapU+1)/2) ) ) * this->_xSamplingRate ) ) / this->_xSamplingRate;
                c00[0]                        = (this->_xSamplingRate * xd) + c000[0]; c00[1] = c000[1]; c00[2] = c000[2]; c00[3] = ( c000[3] * ( 1.0 - xd ) ) + ( c100[3] * xd );
                c01[0]                        = (this->_xSamplingRate * xd) + c001[0]; c01[1] = c001[1]; c01[2] = c001[2]; c01[3] = ( c001[3] * ( 1.0 - xd ) ) + ( c101[3] * xd );
                c10[0]                        = (this->_xSamplingRate * xd) + c010[0]; c10[1] = c010[1]; c10[2] = c010[2]; c10[3] = ( c010[3] * ( 1.0 - xd ) ) + ( c110[3] * xd );
                c11[0]                        = (this->_xSamplingRate * xd) + c011[0]; c11[1] = c011[1]; c11[2] = c011[2]; c11[3] = ( c011[3] * ( 1.0 - xd ) ) + ( c111[3] * xd );
                
                //============================ Interpolate along y
                yd                            = (y - ( ( yBottom - static_cast<int> ( ((this->_maxMapV+1)/2) ) ) * this->_ySamplingRate ) ) / this->_ySamplingRate;
                c0[0]                         = c00[0]; c0[1] = (this->_ySamplingRate * yd) + c00[1]; c0[2] = c00[2]; c0[3] = ( c00[3] * ( 1.0 - yd ) ) + ( c10[3] * yd );
                c1[0]                         = c01[0]; c1[1] = (this->_ySamplingRate * yd) + c01[1]; c1[2] = c01[2]; c1[3] = ( c01[3] * ( 1.0 - yd ) ) + ( c11[3] * yd );
                
                //============================ Interpolate along z
                zd                            = (z - ( ( zBottom - static_cast<int> ( ((this->_maxMapW+1)/2) ) ) * this->_zSamplingRate ) ) / this->_zSamplingRate;

                //============================ Save results
                this->_shellMappedData[sh][static_cast<unsigned int> ( phIt * this->_thetaAngle + thIt)] = ( c0[3] * ( 1.0 - zd ) ) + ( c1[3] * zd );
            }
        }
    }
    
    //======================================== Free unnecessary memory
    delete[] this->_densityMapCorCoords;
    this->_densityMapCorCoords                = nullptr;
    
    if ( !keepInMemory )
    {
        delete[] this->_densityMapCor;
        this->_densityMapCor                  = nullptr;
    }
    
    //======================================== Done
    this->_sphereMapped                       = true;
    
}

void ProSHADE_internal::ProSHADE_data::getSphericalHarmonicsCoeffs ( unsigned int bandwidth, ProSHADE::ProSHADE_settings* settings )
{
    //======================================== Sanity check
    if ( !this->_sphereMapped )
    {
        std::cerr << "!!! ProSHADE ERROR !!! Error in file " << this->_inputFileName << " !!! Attempted to obtain spherical harmonics before mapping the map onto spheres. Call the mapPhaselessToSphere function BEFORE the getSphericalHarmonicsCoeffs function." << std::endl;
        
        if ( settings->htmlReport )
        {
            std::stringstream hlpSS;
            hlpSS << "<font color=\"red\">" << "Attempted to obtain spherical harmonics before mapping the map onto spheres. This looks like an internal bug, please report this case." << "</font>";
            rvapi_set_text                    ( hlpSS.str().c_str(),
                                               "ProgressSection",
                                               settings->htmlReportLineProgress,
                                               1,
                                               1,
                                               1 );
            settings->htmlReportLineProgress += 1;
            rvapi_flush                       ( );
        }
        
        exit ( -1 );
    }
    
    //======================================== Initialise internal values
    this->_bandwidthLimit                     = bandwidth;
    this->_oneDimmension                      = this->_bandwidthLimit * 2;
    
    //======================================== Allocate Memory
    double *inputReal                         = new double [static_cast<unsigned int> ( this->_oneDimmension  * this->_oneDimmension )];
    double *inputZeroes                       = new double [static_cast<unsigned int> ( this->_oneDimmension  * this->_oneDimmension )];
    this->_realSHCoeffs                       = new double*[static_cast<unsigned int> ( this->_noShellsWithData )];
    this->_imagSHCoeffs                       = new double*[static_cast<unsigned int> ( this->_noShellsWithData )];
    this->_sphericalHarmonicsWeights          = new double [static_cast<unsigned int> ( this->_bandwidthLimit * 4 )];
    this->_semiNaiveTableSpace                = new double [static_cast<unsigned int> ( Reduced_Naive_TableSize     ( this->_bandwidthLimit, this->_bandwidthLimit ) +
                                                                              Reduced_SpharmonicTableSize ( this->_bandwidthLimit, this->_bandwidthLimit ) )];
    
    this->_shWorkspace                        = (fftw_complex *) fftw_malloc ( sizeof(fftw_complex) * ( (  8 * this->_bandwidthLimit * this->_bandwidthLimit ) +
                                                                                              ( 10 * this->_bandwidthLimit ) ) );
    
    for ( unsigned int shIt = 0; shIt < this->_noShellsWithData; shIt++ )
    {
        this->_realSHCoeffs[shIt]             = new double [static_cast<unsigned int> ( this->_oneDimmension * this->_oneDimmension )];
        this->_imagSHCoeffs[shIt]             = new double [static_cast<unsigned int> ( this->_oneDimmension * this->_oneDimmension )];
    }
    
    //======================================== Check Memory Allocation
    if ( this->_shWorkspace == nullptr )
    {
        std::cerr << "!!! ProSHADE ERROR !!! Error in file " << this->_inputFileName << " !!! malloc (Memory Allocation) failed for _shWorkspace object in function getSphericalHarmonicsCoeffs." << std::endl;
        
        if ( settings->htmlReport )
        {
            std::stringstream hlpSS;
            hlpSS << "<font color=\"red\">" << "Cannot allocate memory for spherical harmonics decomposition. Could you have run out of memory?" << "</font>";
            rvapi_set_text                    ( hlpSS.str().c_str(),
                                               "ProgressSection",
                                               settings->htmlReportLineProgress,
                                               1,
                                               1,
                                               1 );
            settings->htmlReportLineProgress += 1;
            rvapi_flush                       ( );
        }
        
        exit( -1 ) ;
    }
    
    //======================================== Generate Seminaive and naive table
    this->_semiNaiveTable                     = SemiNaive_Naive_Pml_Table ( this->_bandwidthLimit,
                                                                            this->_bandwidthLimit,
                                                                            this->_semiNaiveTableSpace,
                                                                            reinterpret_cast<double*> ( this->_shWorkspace ) );
    
    //======================================== Make weights for spherical transform
    makeweights                               ( this->_bandwidthLimit,
                                                this->_sphericalHarmonicsWeights );
    
    //======================================== Declare pointers to within workspace
    double *rres;
    double *ires;
    double *fltres;
    double *scratchpad;
    double *rdataptr;
    double *idataptr;
    
    //======================================== Initialize pointers to within workspace
    rres                                      = reinterpret_cast<double*> ( this->_shWorkspace );
    ires                                      = rres + ( this->_oneDimmension * this->_oneDimmension );
    fltres                                    = ires + ( this->_oneDimmension * this->_oneDimmension );
    scratchpad                                = fltres + ( this->_bandwidthLimit );
    
    //======================================== Initialize fft plan along phi angles
    fftw_plan fftPlan;
    fftw_plan dctPlan;
    
    fftw_iodim dims[1];
    fftw_iodim howmany_dims[1];
    
    int rank                                  = 1;
    int howmany_rank                          = 1;
    double normCoeff                          = ( 1.0 / ( static_cast<double> ( this->_oneDimmension ) ) ) * sqrt( 2.0 * M_PI );
    double powerOne                           = 1.0;
    unsigned int hlp1                         = 0;
    unsigned int hlp2                         = 0;
    
    dims[0].n                                 = this->_oneDimmension;
    dims[0].is                                = 1;
    dims[0].os                                = this->_oneDimmension;
    
    howmany_dims[0].n                         = this->_oneDimmension;
    howmany_dims[0].is                        = this->_oneDimmension;
    howmany_dims[0].os                        = 1;
    
    //======================================== Plan fft transform
    fftPlan                                   = fftw_plan_guru_split_dft ( rank,
                                                                           dims,
                                                                           howmany_rank,
                                                                           howmany_dims,
                                                                           inputReal,
                                                                           inputZeroes,
                                                                           rres,
                                                                           ires,
                                                                           FFTW_ESTIMATE  );
    
    //======================================== Initialize dct plan for SHT
    dctPlan                                   = fftw_plan_r2r_1d ( this->_oneDimmension,
                                                                   scratchpad,
                                                                   scratchpad + this->_oneDimmension,
                                                                   FFTW_REDFT10,
                                                                   FFTW_ESTIMATE ) ;

    //======================================== For each shell
    for ( unsigned int shIter = 0; shIter < static_cast<unsigned int> ( this->_noShellsWithData ); shIter++ )
    {
        //==================================== Load pixel hits to signal array
        for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( this->_oneDimmension * this->_oneDimmension ); iter++ )
        {
            inputReal[iter]                   = this->_shellMappedData[shIter][iter];
            inputZeroes[iter]                 = 0.0;
        }

        //==================================== Execute fft plan along phi
        fftw_execute_split_dft                ( fftPlan, inputReal, inputZeroes, rres, ires ) ;
        
        //==================================== Normalize
        for ( unsigned int iter = 0; iter < ( this->_oneDimmension * this->_oneDimmension ); iter++ )
        {
            rres[iter]                       *= normCoeff;
            ires[iter]                       *= normCoeff;
        }
        
        //==================================== Calculate the coefficients for each band
        rdataptr                              = this->_realSHCoeffs[shIter];
        idataptr                              = this->_imagSHCoeffs[shIter];
        for ( unsigned int bandIter = 0; bandIter < this->_bandwidthLimit; bandIter++ )
        {
            //================================ *real* part calculation
            SemiNaiveReduced                  ( rres + ( bandIter * this->_oneDimmension ),
                                                this->_bandwidthLimit,
                                                bandIter,
                                                fltres,
                                                scratchpad,
                                                this->_semiNaiveTable[bandIter],
                                                this->_sphericalHarmonicsWeights,
                                               &dctPlan);
            
            //================================ Save the *real* results
            memcpy                            ( rdataptr, fltres, sizeof(double) * ( this->_bandwidthLimit - bandIter ) );
            rdataptr                         += this->_bandwidthLimit - bandIter;
            
            //================================ *imaginary* part calculation
            SemiNaiveReduced                  ( ires + ( bandIter * this->_oneDimmension ),
                                                this->_bandwidthLimit,
                                                bandIter,
                                                fltres,
                                                scratchpad,
                                                this->_semiNaiveTable[bandIter],
                                                this->_sphericalHarmonicsWeights,
                                               &dctPlan );
            
            //================================ Save the *imaginary* results
            memcpy                            ( idataptr, fltres, sizeof(double) * ( this->_bandwidthLimit - bandIter ) );
            idataptr                         += this->_bandwidthLimit - bandIter;
        }
        
        //==================================== Since the data are real, there is no need to calculate negative order data, as there is the equivalence
        powerOne                              = 1.0;
        hlp1                                  = 0;
        hlp2                                  = 0;
        for( unsigned int order = 1; order < this->_bandwidthLimit; order++)
        {
            powerOne                         *= -1.0;
            for( unsigned int bandIter = order; bandIter < this->_bandwidthLimit; bandIter++){
                
                hlp1                          = seanindex (  order, bandIter, this->_bandwidthLimit );
                hlp2                          = seanindex ( -order, bandIter, this->_bandwidthLimit );
                
                this->_realSHCoeffs[shIter][hlp2] =  powerOne * this->_realSHCoeffs[shIter][hlp1];
                this->_imagSHCoeffs[shIter][hlp2] = -powerOne * this->_imagSHCoeffs[shIter][hlp1];
            }
        }
    }
    
    delete[] inputReal;
    delete[] inputZeroes;
    
    delete[] this->_semiNaiveTable;
    delete[] this->_semiNaiveTableSpace;
    delete[] this->_sphericalHarmonicsWeights;
    
    fftw_free                                 ( this->_shWorkspace );
    
    this->_semiNaiveTable                     = nullptr;
    this->_semiNaiveTableSpace                = nullptr;
    this->_sphericalHarmonicsWeights          = nullptr;
    this->_shWorkspace                        = nullptr;
    
    fftw_destroy_plan                         ( dctPlan );
    fftw_destroy_plan                         ( fftPlan );
    
    //======================================== Done
    this->_sphericalCoefficientsComputed      = true;
    
}