gramSchmidt.hpp

Go to the documentation of this file.

/** \file gramSchmidt.hpp
 * \brief Procedures to orthogonalize vector basis sets
 *
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 * \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 gramSchmidt_hpp
#define gramSchmidt_hpp
 
namespace mx
{
namespace sigproc
{
 
/// Perform Gram-Schmidt ortogonalization of a basis set, and normalize the result.
/** Performs the stabilized Gram-Schmidt procedure on the input basis set, which
 * is the columns of the array.  Optionally the output is normalized.
 *
 * \tparam progress if true, then the loop index is printed for progress reporting
 * \tparam eigenTout is the Eigen array type of the desired output
 * \tparam eigenTin is the Eigen array type of the input
 *
 * \ingroup signal_processing
 */
template <int progress = 0, typename eigenTout, typename eigenTin>
void gramSchmidt( eigenTout &out,       ///< [out] the orthonormal basis set constructed from the input
                  const eigenTin &in,   ///< [in] a basis set, where each column represents one vector.
                  bool normalize = true ///< [in] [optional] whether or not to normalize the output
)
{
    out.resize( in.rows(), in.cols() );
 
    out.col( 0 ) = in.col( 0 );
 
    for( int i = 1; i < in.cols(); ++i )
    {
        if( progress )
        {
            std::cout << i + 1 << "/" << in.cols() << "\n";
        }
 
        // out.col(i) = in.col(i);
 
        out.col( i ) = in.col( i ) - ( ( in.col( i ).matrix().dot( out.col( 0 ).matrix() ) ) /
                                       ( out.col( 0 ).matrix().dot( out.col( 0 ).matrix() ) ) ) *
                                         out.col( 0 );
 
        for( int j = 1; j < i; ++j )
        {
            out.col( i ) = out.col( i ) - ( ( out.col( i ).matrix().dot( out.col( j ).matrix() ) ) /
                                            ( out.col( j ).matrix().dot( out.col( j ).matrix() ) ) ) *
                                              out.col( j );
        }
    }
 
    if( normalize )
    {
        for( int i = 0; i < out.cols(); ++i )
        {
            out.col( i ) = out.col( i ) / out.col( i ).matrix().norm();
        }
    }
}
 
/// Perform Gram-Schmidt ortogonalization of a basis set on a window, and normalize the result.
/** Performs the stabilized Gram-Schmidt procedure on the input basis set over a window (or
 * weight function), followed by normalization of the result.
 *
 * \param out [out] is the orthonormal basis set constructed from the input
 * \param in [in] is a basis set, where each column represents one vector.
 * \param window [in] is the window, or weighting function
 *
 * \tparam progress if true, then the loop index is printed for progress reporting
 * \tparam eigenTout is the Eigen array type of the desired output
 * \tparam eigenTin is the Eigen array type of the input
 * \tparam eigenTWin is the Eigen array type of the window
 *
 * \ingroup signal_processing
 */
template <int progress = 0, typename eigenTout, typename eigenTin, typename eigenTWin>
void gramSchmidt( eigenTout &out, const eigenTin &in, const eigenTWin &window )
{
    // out.resize(in.rows(), in.cols());
 
    out.col( 0 ) = in.col( 0 );
 
    for( int i = 1; i < in.cols(); ++i )
    {
        if( progress )
        {
            std::cout << i + 1 << "/" << in.cols() << "\n";
        }
 
        // out.col(i) = in.col(i);
 
        out.col( i ) = in.col( i ) - ( ( ( in.col( i ) * window ).matrix().dot( out.col( 0 ).matrix() ) ) /
                                       ( ( out.col( 0 ) * window ).matrix().dot( out.col( 0 ).matrix() ) ) ) *
                                         out.col( 0 );
 
        for( int j = 1; j < i; ++j )
        {
            out.col( i ) = out.col( i ) - ( ( ( out.col( i ) * window ).matrix().dot( out.col( j ).matrix() ) ) /
                                            ( ( out.col( j ) * window ).matrix().dot( out.col( j ).matrix() ) ) ) *
                                              out.col( j );
        }
    }
 
    for( int i = 0; i < out.cols(); ++i )
    {
        out.col( i ) = out.col( i ) / ( out.col( i ) * window.sqrt() ).matrix().norm();
    }
}
 
// Unwraps the gram schmidt coefficients to give a spectrum in terms of the original basis set
// Helper function for gramSchmidtSpectrum (below)
template <int progress = 0, typename eigenT>
void baseSpectrum( eigenT &bspect, eigenT &gsspect )
{
    bspect.resize( gsspect.rows(), gsspect.cols() );
    bspect.setZero();
 
    // #pragma omp parallel for
    for( int i = 0; i < gsspect.rows(); ++i )
    {
        bspect( i, i ) = gsspect( i, i );
 
        for( int j = i - 1; j >= 0; --j )
        {
            bspect.row( i ) -= gsspect( i, j ) * bspect.row( j );
        }
    }
}
 
/// Perform Gram-Schmidt ortogonalization of a basis set, and normalize the result, while recording the spectrum.
/** Performs the stabilized Gram-Schmidt procedure on the input basis set, followed
 * by normalization of the result.  Also records the spectrum, that is the coefficients of the linear expansion
 * in the orginal basis set for the resultant basis set.
 *
 *
 * \tparam progress if true, then the loop index is printed for progress reporting
 * \tparam eigenTout is the Eigen array type of the output orthogonalized array
 * \tparam eigenTout2 is the Eigen array type of the spectrum
 * \tparam eigenTin is the Eigen array type of the input
 *
 * \ingroup signal_processing
 */
template <int progress = 0, typename eigenTout, typename eigenTout2, typename eigenTin>
void gramSchmidtSpectrum(
    eigenTout &out,                         ///< [out] the orthonormal basis set constructed from the input
    eigenTout2 &spect,                      ///< [out] the spectrum
    const eigenTin &in,                     ///< [in] a basis set, where each column represents one vector
    typename eigenTin::Scalar normPix = 0.0 /**< [in] [optional] area of (usually number of pixels in) the orthogonal
                                                       region for normalization.  If 0 the basis is not renormalized */
)
{
    typedef typename eigenTout::Scalar Scalar;
 
    out.resize( in.rows(), in.cols() );
 
    eigenTout2 gsspect;
 
    gsspect.resize( in.cols(), in.cols() );
    gsspect.setZero();
 
    out.col( 0 ) = in.col( 0 );
    gsspect( 0, 0 ) = 1;
 
    for( int i = 1; i < in.cols(); ++i )
    {
        if( progress )
        {
            std::cout << i + 1 << "/" << in.cols() << "\n";
        }
 
        gsspect( i, i ) = 1;
 
        gsspect( i, 0 ) = ( ( in.col( i ).matrix().dot( out.col( 0 ).matrix() ) ) /
                            ( out.col( 0 ).matrix().dot( out.col( 0 ).matrix() ) ) );
        out.col( i ) = in.col( i ) - gsspect( i, 0 ) * out.col( 0 );
 
        for( int j = 1; j < i; ++j )
        {
            gsspect( i, j ) = ( ( out.col( i ).matrix().dot( out.col( j ).matrix() ) ) /
                                ( out.col( j ).matrix().dot( out.col( j ).matrix() ) ) );
            out.col( i ) = out.col( i ) - gsspect( i, j ) * out.col( j );
        }
    }
 
    // Here we unwrap the gram schmidt coefficients, giving us the coefficients in the original basis
    baseSpectrum( spect, gsspect );
 
    if( normPix > 0 )
    {
        Scalar norm;
 
        for( int i = 0; i < out.cols(); ++i )
        {
            norm = sqrt( out.col( i ).square().sum() / normPix );
 
            out.col( i ) /= norm;
            // spect.row(i) /= norm;
        }
    }
}
 
} // namespace sigproc
} // namespace mx
 
#endif // gramSchmidt_hpp