fourierModes.hpp

Go to the documentation of this file.

/** \file fourierModes.hpp
 * \author Jared R. Males
 * \brief Functions for generating 2D Fourier modes
 * \ingroup signal_processing_files
 *
 */
 
//***********************************************************************//
// Copyright 2015, 2016, 2017 Jared R. Males (jaredmales@gmail.com)
//
// This file is part of mxlib.
//
// mxlib is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// mxlib is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with mxlib.  If not, see <http://www.gnu.org/licenses/>.
//***********************************************************************//
 
#ifndef fourierModes_hpp
#define fourierModes_hpp
 
#include <vector>
#include <algorithm>
 
#include "../mxError.hpp"
 
#include "../math/constants.hpp"
#include "../math/geo.hpp"
 
namespace mx
{
 
namespace sigproc
{
 
/**
 * \ingroup fourier_basis
 *
 * \todo need to handle references and r-values so this works directly with Eigen arrays.
 *
 * @{
 */
 
/** \def MX_FOURIER_BASIC
 * \brief Signifies the basic sine/cosine Fourier basis.
 */
#define MX_FOURIER_BASIC 0
 
/** \def MX_FOURIER_MODIFIED
 * \brief Signifies the modified Fourier basis.
 */
#define MX_FOURIER_MODIFIED 1
 
/// Populate a single cosine mode
/** Makes a 2D image of the cosine mode:
    \f[
    M_c(\vec{q}) = \cos\left( 2\pi\frac{m}{D}u + 2\pi\frac{n}{D}v \right)
    \f]
  * where
    \f[
    \vec{q} = u \hat{u} + v \hat{v}
    \f]
  * and  \f$ D \f$ is taken to be the maximum of the number of columns and rows in the image.
  *
  * \retval 0 on success
  *
  * \tparam typeN is an Eigen-like reference type.
  */
template <class typeN>
int makeFourierModeC( typeN im,                 ///<  [out] is an Eigen-like image
                      typename typeN::Scalar m, ///<  [in] specifies the spatial frequency in the u direction
                      typename typeN::Scalar n  ///<  [in] n specifies the spatial frequency in the v direction
)
{
    typedef typename typeN::Scalar realT;
 
    int dim1 = im.cols();
    int dim2 = im.rows();
 
    realT uc, vc, u, v, D;
 
    uc = 0.5 * ( dim1 - 1.0 );
    vc = 0.5 * ( dim2 - 1.0 );
 
    D = std::max( dim1, dim2 );
 
    for( int i = 0; i < dim1; ++i )
    {
        u = i - uc;
        for( int j = 0; j < dim2; ++j )
        {
            v = j - vc;
 
            im( i, j ) = cos( math::two_pi<realT>() / D * ( m * u + n * v ) );
        }
    }
 
    return 0;
}
 
/// Populate a single sine mode
/** Makes a 2D image of the cosine mode:
    \f[
    M_s(\vec{q}) = \sin\left( 2\pi\frac{m}{D}u + 2\pi\frac{n}{D}v \right)
    \f]
  * where
    \f[
    \vec{q} = u \hat{u} + v \hat{v}
    \f]
  * and  \f$ D \f$ is taken to be the maximum of the number of columns and rows in the image.
  * \param [out] im is an Eigen-like image
  * \param [in] m specifies the spatial frequency in the u direction
  * \param [in] n specifies the spatial frequency in the v direction
  *
  * \retval 0 on success
  *
  * \tparam typeN is an Eigen-like reference type.
  */
template <class typeN>
int makeFourierModeS( typeN im, typename typeN::Scalar m, typename typeN::Scalar n )
{
    typedef typename typeN::Scalar realT;
 
    int dim1 = im.cols();
    int dim2 = im.rows();
 
    realT uc, vc, u, v, D;
 
    uc = 0.5 * ( dim1 - 1.0 );
    vc = 0.5 * ( dim2 - 1.0 );
 
    D = std::max( dim1, dim2 );
 
    for( int i = 0; i < dim1; ++i )
    {
        u = i - uc;
        for( int j = 0; j < dim2; ++j )
        {
            v = j - vc;
 
            im( i, j ) = sin( math::two_pi<realT>() / D * ( m * u + n * v ) );
        }
    }
 
    return 0;
}
 
/// Populate a single basic Fourier mode
/** Makes a 2D image of the basic sine or cosine mode.  Calls \ref makeFourierModeC and \ref makeFourierModeS.
 *
 * \param [out] im is an Eigen-like image
 * \param [in] m specifies the spatial frequency in the u direction
 * \param [in] n specifies the spatial frequency in the v direction
 * \param [in] p species sine (-1) or cosine (+1)
 *
 * \retval 0 on success
 *
 * \tparam typeN is an Eigen-like reference type.
 */
template <class typeN>
int makeFourierMode( typeN im, typename typeN::Scalar m, typename typeN::Scalar n, int p )
{
    if( p == -1 )
        return makeFourierModeS( im, m, n );
 
    if( p == 1 )
        return makeFourierModeC( im, m, n );
 
    mxError( "makeFourierMode", MXE_INVALIDARG, "p must be +1 (cosine) or -1 (sine)." );
 
    return -1;
}
 
/// Populate a single modified Fourier mode
/** Makes a 2D image of the modified fourier mode:
    \f[
    M_p(\vec{q}) = \cos\left( 2\pi \frac{m}{D}u + 2\pi\frac{n}{D}v \right) + p  \sin\left( 2\pi\frac{m}{D}u +
  2\pi\frac{n}{D}v \right) \f]
  * where \f$ p = \pm 1 \f$, \f$ \vec{q} = u \hat{u} + v \hat{v} \f$,
  * and  \f$ D \f$ is taken to be the maximum of the number of columns in the image.
  *
  * \param [out] im is an Eigen-like image
  * \param [in] m specifies the spatial frequency in the u direction
  * \param [in] n specifies the spatial frequency in the v direction
  * \param [in] p is +/- 1 specifying which modified Fourier mode
  *
  * \retval 0 on success
  *
  * \tparam typeN is an Eigen-like reference type.
  */
template <class typeN>
int makeModifiedFourierMode(
    typeN im, typename typeN::Scalar m, typename typeN::Scalar n, int p, typename typeN::Scalar ang = 0 )
{
    typedef typename typeN::Scalar realT;
 
    if( p != 1 && p != -1 )
    {
        mxError( "makeModifiedFourierMode", MXE_INVALIDARG, "p must be +1 or -1." );
    }
 
    int dim1 = im.cols();
    int dim2 = im.rows();
 
    realT xc, yc, x, y, arg, D;
 
    xc = 0.5 * ( dim1 - 1.0 );
    yc = 0.5 * ( dim2 - 1.0 );
 
    D = std::max( dim1, dim2 );
 
    realT c_ang, s_ang;
    if( ang != 0 )
    {
        c_ang = cos( math::dtor( ang ) );
        s_ang = sin( math::dtor( ang ) );
    }
 
    for( int i = 0; i < dim1; ++i )
    {
 
        for( int j = 0; j < dim2; ++j )
        {
            x = i - xc;
            y = j - yc;
 
            if( ang != 0 )
            {
                realT x0 = x * c_ang - y * s_ang;
                realT y0 = x * s_ang + y * c_ang;
 
                x = x0;
                y = y0;
            }
 
            arg = math::two_pi<realT>() / D * ( m * x + n * y );
 
            im( i, j ) = cos( arg ) + p * sin( arg );
        }
    }
 
    return 0;
}
 
/// Populate a single modified Fourier +1 mode
/** See \ref makeModifiedFourierMode.
 *
 * \param [out] im is an Eigen-like image
 * \param [in] m specifies the spatial frequency in the u direction
 * \param [in] n specifies the spatial frequency in the v direction
 *
 * \tparam typeN is an Eigen-like reference type.
 */
template <class typeN>
int makeModifiedFourierModeP( typeN im,
                              typename typeN::Scalar m,
                              typename typeN::Scalar n,
                              typename typeN::Scalar ang = 0 )
{
    return makeModifiedFourierMode( im, m, n, 1, ang );
}
 
/// Populate a single modified Fourier -1 mode
/** See \ref makeModifiedFourierMode.
 *
 * \param [out] im is an Eigen-like image
 * \param [in] m specifies the spatial frequency in the u direction
 * \param [in] n specifies the spatial frequency in the v direction
 *
 * \tparam typeN is an Eigen-like reference type.
 */
template <class typeN>
int makeModifiedFourierModeM( typeN im,
                              typename typeN::Scalar m,
                              typename typeN::Scalar n,
                              typename typeN::Scalar ang = 0 )
{
    return makeModifiedFourierMode( im, m, n, -1, ang );
}
 
/// Structure to contain the parameters of a Fourier mode.
struct fourierModeDef
{
    int m; ///< Spatial frequency m index
    int n; ///< Spatial frequency n index
    int p; ///< +/- 1 for modified Fourier modes, -1==>sine, +1==>cosine for basic Fourier modes.
 
    /// Constructor
    fourierModeDef()
    {
        m = 0;
        n = 0;
        p = 0;
    }
};
 
/// Sorting functor for modes according to the Poyner and Veran (2005) standard.
/** Follows the conventions of  Poyneer & Veran (2005) \cite poyneer_and_veran_2005
 *
 * \param[in] spf1 is an object of type \ref fourierModeDef to compare with spf2
 * \param[in] spf2 is an object of type \ref fourierModeDef to compare with spf1
 *
 * \returns the result of (spf1 < spf2) according to the above rules.
 */
bool comp_fourierModeDef_PandV( const fourierModeDef &spf1, const fourierModeDef &spf2 )
{
    double fr1 = spf1.m * spf1.m + spf1.n * spf1.n;
    double fr2 = spf2.m * spf2.m + spf2.n * spf2.n;
 
    if( fr1 == fr2 )
    {
        if( spf1.n == spf2.n )
        {
            if( spf1.m == spf2.m )
                return ( spf1.p < spf2.p );
            else
                return ( spf1.m < spf2.m );
        }
        return ( spf1.n < spf2.n );
    }
 
    // Otherwise sort by lowest absolute value
    return ( fr1 < fr2 );
}
 
/// Generate a Poyneer and Veran spatial frequency grid.
/** Follows the convention of Poyneer and Veran (2005) \cite poyneer_and_veran_2005
 * \todo use push_back instead of resize
 * \param spf
 * \param N is the linear number of degrees of freedom.
 */
int makeFourierModeFreqs_PandV( std::vector<fourierModeDef> &spf, int N )
{
    int Nmodes = N * N;
 
    spf.resize( Nmodes );
 
    int _k, _l;
    int modeCount = 0;
    for( int ll = 0; ll <= N; ++ll )
    {
        for( int l = ll, k = 0; l >= 0; --l, ++k )
        {
            if( modeCount >= Nmodes )
            {
                mxError( "makeFourierModeFreqs_PandV", MXE_SIZEERR, "mode count exceeded expected number of modes" );
                return -1;
            }
 
            if( k == 0 && l > .5 * N )
                continue;
 
            if( k > .5 * N && ( l == 0 || l >= N - k ) )
                continue;
 
            if( k > 0.5 * N )
                _k = k - N;
            else
                _k = k;
 
            if( l > 0.5 * N )
                _l = l - N;
            else
                _l = l;
 
            // There is always a cosine mode
            spf[modeCount].m = _k;
            spf[modeCount].n = _l;
            spf[modeCount].p = 1;
 
            ++modeCount;
 
            // Make sure it isn't one of the corners, then populate the sine mode
            if( !( ( k == 0 && l == 0 ) || ( k == 0.5 * N && l == 0 ) || ( k == 0 && l == 0.5 * N ) ||
                   ( k == 0.5 * N && l == 0.5 * N ) ) )
            {
                spf[modeCount].m = _k;
                spf[modeCount].n = _l;
                spf[modeCount].p = -1;
                ++modeCount;
            }
        }
    }
 
    std::sort( spf.begin(), spf.end(), comp_fourierModeDef_PandV );
 
    return 0;
}
 
/// Sorting functor for modes according to the mx::AO standard.
/** Given a spatial frequency vector \f$ \vec{k} = k_u \hat{u} + k_v \hat{v} \f$, sorts first by
 * \f$ | \vec{k} | \f$, then by the angle from the u axis.  Results in mode numbers increasing radially
 * and counter-clockwise from the \f$ u \f$ axis.
 *
 * \param[in] spf1 is an object of type \ref fourierModeDef to compare with spf2
 * \param[in] spf2 is an object of type \ref fourierModeDef to compare with spf1
 *
 * \returns the result of (spf1 < spf2) according to the above rules.
 */
bool comp_fourierModeDef( const fourierModeDef &spf1, const fourierModeDef &spf2 )
{
    double k1 = spf1.m * spf1.m + spf1.n * spf1.n;
    double k2 = spf2.m * spf2.m + spf2.n * spf2.n;
 
    if( k1 == k2 )
    {
        // Now compare by angle
        double a1 = atan2( spf1.n, spf1.m );
        if( a1 < 0 )
            a1 += math::two_pi<double>();
 
        double a2 = atan2( spf2.n, spf2.m );
        if( a2 < 0 )
            a2 += math::two_pi<double>();
 
        if( a1 == a2 )
            return ( spf1.p < spf2.p );
 
        return ( a1 < a2 );
    }
 
    // Otherwise sort by lowest absolute value
    return ( k1 < k2 );
}
 
/// Generate a circular spatial frequency grid.
/** \todo use push_back instead of resize
 *
 * \param spf
 * \param N is the linear number of degrees of freedom.  The number of modes is ((N+1)(N+1) - 1).
 */
int makeFourierModeFreqs_Circ( std::vector<fourierModeDef> &spf, int N )
{
    int krad;
    int Nmodes;
 
    krad = floor( 0.5 * N );
 
    Nmodes = 0.25 * math::pi<double>() * N * N * 1.1; // The 1.1 is a buffer
 
    spf.resize( Nmodes );
 
    int modeCount = 0;
 
    for( int m = -krad; m <= krad; ++m )
    {
        int nmax = sqrt( krad * krad - m * m );
 
        for( int n = 0; n <= nmax; ++n )
        {
            if( n == 0 && m <= 0 )
                continue;
 
            for( int p = -1; p <= 1; p += 2 )
            {
                if( modeCount >= Nmodes )
                {
                    mxError( "makeFourierModeFreqs_Circ", MXE_SIZEERR, "mode count exceeded expected number of modes" );
                    return -1;
                }
 
                spf[modeCount].m = m;
                spf[modeCount].n = n;
                spf[modeCount].p = p;
                ++modeCount;
            }
        }
    }
 
    // Erase any extra modes (there will bedue to the buffer).
    spf.erase( spf.begin() + modeCount, spf.end() );
 
    // And now sort it
    std::sort( spf.begin(), spf.end(), comp_fourierModeDef );
 
    return 0;
}
 
/// Calculate the index of a Fourier mode given its (m,n,p) coordinates.
/** The index increases counter-clockwise from the m axis in squares, with p==-1 being even and p==1 being odd.
 *
 * \param [in] m is the m spatial frequency coordinate
 * \param [in] n is the n spatial frequency coordinage
 * \param [in] p is the parity, +/-1, for sine or cosine or specifiying the modified Fourier mode.
 *
 * \returns the index of the Fourier mode on success
 * \returns -1 on an error.
 */
int fourierModeNumber( int m, int n, int p )
{
 
    if( p != -1 && p != 1 )
    {
        mxError( "fourierModeCoordinates", MXE_INVALIDARG, "p must +/-1" );
        return -1;
    }
 
    if( m == 0 && n == 0 )
    {
        mxError( "fourierModeCoordinates", MXE_INVALIDARG, "piston is ignored" );
        return -1;
    }
 
    int m0 = std::max( abs( m ), abs( n ) );
 
    int i0 = 2 * ( m0 * m0 - m0 );
 
    // int di0= 4*m0;
 
    int di;
 
    if( m == m0 )
    {
        di = n;
    }
    else if( m == -m0 )
    {
        di = 4 * m0 - n;
    }
    else
    {
        di = 2 * m0 - m;
    }
 
    return 2 * ( i0 + di ) + ( 1 + p ) / 2;
}
 
/// Calculate the (m,n,p) coordinates of a Fourier mode given its index.
/** The index increases counter-clockwise from the m axis in squares, with p==-1 being even and p==1 being odd.
 *
 * \returns 0 on success, a negative number otherwise
 */
int fourierModeCoordinates(
    int &m, ///< [out] is the m spatial frequency coordinate
    int &n, ///< [out] is the n spatial frequency coordinage
    int &p, ///< [out] is the parity, +/-1, for sine or cosine or specifiying the modified Fourier mode.
    int i ) ///< [in] is the mode index.
{
    if( i < 0 )
    {
        mxError( "fourierModeCoordinates", MXE_INVALIDARG, "i must be >= 0" );
        m = 0;
        n = 0;
        p = 0;
        return -1; /// \retval -1 if i < 0
    }
 
    int ii = ( i ) / 2;
 
    int m0 = 1;
 
    while( 2 * ( m0 * m0 - m0 ) <= ii )
        ++m0;
    --m0;
 
    int di = ii - 2 * ( m0 * m0 - m0 );
 
    if( di <= m0 )
    {
        m = m0;
        n = di;
    }
    else if( di <= 3 * m0 )
    {
        m = 2 * m0 - di;
        n = m0;
    }
    else
    {
        m = -m0;
        n = 4 * m0 - di;
    }
 
    p = -1 + 2 * ( i % 2 );
 
    return 0;
}
 
/// Generate a rectangular spatial frequency grid.
/**
 * \param spf is a vector of fourierModeDef structures, which will be resized and populated.
 * \param N is the linear number of degrees of freedom.  The number of modes is ((N+1)(N+1) - 1).
 */
int makeFourierModeFreqs_Rect( std::vector<fourierModeDef> &spf, int N )
{
    int Ndx = floor( 0.5 * ( N ) );
 
    int Nmodes = 2 * ( ( 2 * Ndx + 1 ) * ( Ndx + 1 ) - ( Ndx + 1 ) );
 
    spf.resize( Nmodes );
 
    int m, n, p;
 
    for( int i = 0; i < Nmodes; ++i )
    {
        fourierModeCoordinates( m, n, p, i );
 
        spf[i].m = m;
        spf[i].n = n;
        spf[i].p = p;
    }
 
    return 0;
}
 
/// Fill in a cube with a Fourier basis.
/** Fills the cube with either the basic or modified Fourier basis for the modes specified.
 *
 * \param[out] cube will be allocated to hold and will be filled with the modes
 * \param[in] dim is the linear size of the maps, each is
 * \param[in] spf is a vector of mode definitions to use for each mode
 * \param[in] basisType is either MX_FOURIER_MODIFIED or MX_FOURIER_BASIC
 *
 * \tparam cubeT is an eigen-like cube, e.g. \ref mx::eigenCube.
 */
template <typename cubeT>
int fillFourierBasis(
    cubeT &cube, int dim, std::vector<fourierModeDef> spf, int basisType, typename cubeT::Scalar ang = 0 )
{
    int Nmodes = spf.size();
 
    cube.resize( dim, dim, Nmodes );
 
    for( int j = 0; j < Nmodes; ++j )
    {
        if( basisType == MX_FOURIER_MODIFIED )
        {
            if( makeModifiedFourierMode( cube.image( j ), spf[j].m, spf[j].n, spf[j].p, ang ) < 0 )
                return -1;
        }
        else
        {
            if( makeFourierMode( cube.image( j ), spf[j].m, spf[j].n, spf[j].p ) < 0 )
                return -1;
        }
    }
 
    return 0;
}
 
/// Generate a rectangular modified Fourier basis with the Poyneer and Veran ordering.
/** Fills the cube with modes generated by \ref makeFourierMode or \ref makeModifiedFourierMode using the mode grid
 * generated by \ref makeFourierModeFreqs_PandV. This follows the convention of Poyneer and Veran (2005) \cite
 * poyneer_and_veran_2005 for mode numbering.
 *
 * \param[out] cube will be allocated to hold and will be filled with the modes
 * \param[in] dim is the linear size of the maps, each is
 * \param[in] N is the linear number of degrees of freedom.  The number of modes is ((N+1)(N+1) - 1).
 * \param[in] basisType is either MX_FOURIER_MODIFIED or MX_FOURIER_BASIC
 *
 * \tparam cubeT is an eigen-like cube, e.g. \ref mx::eigenCube.
 */
template <typename cubeT>
int makeFourierBasis_PandV( cubeT &cube, int dim, int N, int basisType )
{
    std::vector<fourierModeDef> spf;
 
    if( makeFourierModeFreqs_PandV( spf, N ) < 0 )
    {
        return -1;
    }
 
    return fillFourierBasis( cube, dim, spf, basisType );
}
 
/// Generate a circular Fourier basis.
/** Fills the cube with modes generated by \ref makeFourierMode or \ref makeModifiedFourierMode using the mode grid
 * generated by \ref makeFourierModeFreqs_Circ.
 *
 * \param[out] cube will be allocated to hold and will be filled with the modes
 * \param[in] dim is the linear size of the maps, each is
 * \param[in] N is the linear number of degrees of freedom.  The number of modes is ((N+1)(N+1) - 1).
 * \param[in] basisType is either MX_FOURIER_MODIFIED or MX_FOURIER_BASIC
 *
 * \tparam cubeT is an eigen-like cube, e.g. \ref mx::eigenCube.
 */
template <typename cubeT>
int makeFourierBasis_Circ( cubeT &cube, int dim, int N, int basisType )
{
    std::vector<fourierModeDef> spf;
 
    if( makeFourierModeFreqs_Circ( spf, N ) < 0 )
    {
        return -1;
    }
 
    return fillFourierBasis( cube, dim, spf, basisType );
}
 
/// Generate a rectangular Fourier basis.
/** Fills the cube with modes generated by \ref makeFourierMode or \ref makeModifiedFourierMode using the mode grid
 * generated by \ref makeFourierModeFreqs_Rect.
 *
 * \param[out] cube will be allocated to hold and will be filled with the modes
 * \param[in] dim is the linear size of the maps, each is
 * \param[in] N is the linear number of degrees of freedom.  The number of modes is ((N+1)(N+1) - 1).
 * \param[in] basisType is either MX_FOURIER_MODIFIED or MX_FOURIER_BASIC
 *
 * \tparam cubeT is an eigen-like cube, e.g. \ref mx::eigenCube.
 */
template <typename cubeT>
int makeFourierBasis_Rect( cubeT &cube, int dim, int N, int basisType, typename cubeT::Scalar ang = 0 )
{
    std::vector<fourierModeDef> spf;
 
    if( makeFourierModeFreqs_Rect( spf, N ) < 0 )
    {
        return -1;
    }
 
    return fillFourierBasis( cube, dim, spf, basisType, ang );
}
 
///@}
 
} // namespace sigproc
} // namespace mx
 
#endif // fourierModes_hpp