include/zpolbasist.hh

Go to the documentation of this file.

#ifndef ZSIG_ZPOLBASIST_HH
#define ZSIG_ZPOLBASIST_HH

//== INCLUDES =================================================================

#include <cmath>
#include <vector>
#include <complex>
#include <iostream>

//== NAMESPACES ===============================================================

namespace zsig {

//=== IMPLEMENTATION ==========================================================

template< class T >
T fac( const int& _n ) {
        T f = (T)1;
        for (int i = 2; i <= _n; ++i) f *= i;
        return f;
}

template< class T >
T compute_R( const unsigned& _p,
             const int& _q,
             const T& _r ) {

        int a, b, c, d;
        T sum = (T)0;

        for (int k = (_q<0?-_q:_q); k <= _p; k += 2) {

                a = (_p - k) / 2; b = (_p + k) / 2;
                c = (k - _q) / 2; d = (k + _q) / 2;

                sum += ( ( pow( (T)-1.0, (T)a ) * fac<T>(b) ) /
                         ( fac<T>(a) * fac<T>(c) * fac<T>(d) ) ) * pow( _r, (T)k );

        }

        return sum;

}

template< class T >
std::complex< T > compute_V( const unsigned& _p,
                             const int& _q,
                             const T& _r,
                             const T& _t ) {

        return std::polar( compute_R<T>( _p, _q, _r), _q * _t );

}

//== CLASS DEFINITION =========================================================

template< unsigned Order = 8, class T = long double >
class ZernikePolynomialsBasisT {

public:

        typedef ZernikePolynomialsBasisT< Order, T > zpolbasis_type; 
        typedef std::complex< T > value_type; 
        typedef std::vector< value_type > radial_polynomial; 

        ZernikePolynomialsBasisT() : pol(Order+1) {
                for (unsigned p = 0; p <= Order; ++p) pol[p].resize( 1+p/2 );
        }

        void project( const T **_fxy,
                      const zpolbasis_type **_zpb,
                      const unsigned int& _szx,
                      const unsigned int& _szy ) {
                T x, y, r; // (x, y) cartesian coordinates and radius r
                for (unsigned p = 0; p <= Order; ++p) {
                        T p1pi = (p + 1) / M_PI; // frac outside summation
                        for (int qi = 0; qi <= p/2; ++qi) {
                                pol[p][qi] = value_type();
                                for (int gx = 0; gx < _szx; ++gx) {
                                        x = ( (T)2 * gx + (T)1 ) / (T)_szx - (T)1;
                                        for (int gy = 0; gy < _szy; ++gy) {
                                                y = ( (T)2 * gy + (T)1 ) / (T)_szy - (T)1;
                                                r = (T)sqrt( x*x + y*y );
                                                if( r > (T)1 ) continue; // outside unit circle is zero
                                                pol[p][qi] += std::conj( _zpb[gx][gy][p][qi] ) * _fxy[gx][gy];
                                        } // gy
                                } // gx
                                pol[p][qi] *= p1pi;
                        } // qi
                } // p
        }

        void project( const T **_fxy,
                      zpolbasis_type **_zpb,
                      const unsigned int& _szx,
                      const unsigned int& _szy ) {
                project( _fxy, (const zpolbasis_type **)_zpb, _szx, _szy );
        }

        void project( T **_fxy,
                      const zpolbasis_type **_zpb,
                      const unsigned int& _szx,
                      const unsigned int& _szy ) {
                project( (const T **)_fxy, (const zpolbasis_type **)_zpb, _szx, _szy );
        }

        void project( T **_fxy,
                      zpolbasis_type **_zpb,
                      const unsigned int& _szx,
                      const unsigned int& _szy ) {
                project( (const T **)_fxy, (const zpolbasis_type **)_zpb, _szx, _szy );
        }

        void project( const T **_fxy,
                      const unsigned int& _szx,
                      const unsigned int& _szy ) {
                T x, y, r, t; // (x, y) cartesian and (r, t) polar coordinates
                for (unsigned p = 0; p <= Order; ++p) {
                        T p1pi = (p + 1) / M_PI; // frac outside summation
                        int qi = 0; // q index: q[0, 2, 4] -> qi[0, 1, 2]
                        for (int q = p % 2; q <= p; q += 2, qi++) {
                                pol[p][qi] = value_type();
                                for (int gx = 0; gx < _szx; ++gx) {
                                        x = ( (T)2 * gx + (T)1 ) / (T)_szx - (T)1;
                                        for (int gy = 0; gy < _szy; ++gy) {
                                                y = ( (T)2 * gy + (T)1 ) / (T)_szy - (T)1;
                                                r = (T)sqrt( x*x + y*y );
                                                if( r > (T)1 ) continue; // outside unit circle is zero
                                                t = (T)atan2( y, x );
                                                pol[p][qi] += std::conj( compute_V(p, q, r, t) ) * _fxy[gx][gy];
                                        } // gy
                                } // gx
                                pol[p][qi] *= p1pi;
                        } // q
                } // p
        }

        void project( T **_fxy,
                      const unsigned int& _szx,
                      const unsigned int& _szy ) {
                project( (const T **)_fxy, _szx, _szy );
        }

        void reconstruct( T **_fxy,
                          const zpolbasis_type **_zpb,
                          const unsigned int& _szx,
                          const unsigned int& _szy ) const {
                value_type rv; // reconstructed value
                for (int gx = 0; gx < _szx; ++gx) {
                        for (int gy = 0; gy < _szy; ++gy) {
                                rv = value_type();
                                for (unsigned p = 0; p <= Order; ++p) {
                                        for (int qi = 0; qi <= p/2; ++qi) {
                                                rv += pol[p][qi] * _zpb[gx][gy][p][qi]; // using +q
                                                if( p % 2 != 0 or qi != 0 ) // for p even and q zero z does not repeat itself
                                                        rv += std::conj( pol[p][qi] ) * std::conj( _zpb[gx][gy][p][qi] ); // using -q
                                        } // qi
                                } // p
                                _fxy[gx][gy] = rv.real();
                        } // gy
                } // gx
        }

        void reconstruct( T **_fxy,
                          zpolbasis_type **_zpb,
                          const unsigned int& _szx,
                          const unsigned int& _szy ) const {
                reconstruct( _fxy, (const zpolbasis_type **)_zpb, _szx, _szy );
        }

        void reconstruct( T **_fxy,
                          const unsigned int& _szx,
                          const unsigned int& _szy ) const {
                T x, y, r, t; // (x, y) cartesian and (r, t) polar coordinates
                value_type rv, zpV; // reconstructed value and Zernike Polynomial V
                for (int gx = 0; gx < _szx; ++gx) {
                        x = ( (T)2 * gx + (T)1 ) / (T)_szx - (T)1;
                        for (int gy = 0; gy < _szy; ++gy) {
                                y = ( (T)2 * gy + (T)1 ) / (T)_szy - (T)1;
                                r = (T)sqrt( x*x + y*y );
                                rv = value_type();
                                if( r <= (T)1 ) { // reconstruct inside unit circle only
                                        for (unsigned p = 0; p <= Order; ++p) {
                                                int qi = 0; // q index: q[0, 2, 4] -> qi[0, 1, 2]
                                                for (int q = p % 2; q <= p; q += 2, qi++) {
                                                        t = (T)atan2( y, x );
                                                        zpV = compute_V(p, q, r, t);
                                                        rv += pol[p][qi] * zpV; // using +q
                                                        if( p % 2 != 0 or qi != 0 ) // for p even and q zero z does not repeat itself
                                                                rv += std::conj( pol[p][qi] ) * std::conj( zpV ); // using -q
                                                } // q
                                        } // p
                                } // if
                                _fxy[gx][gy] = rv.real();
                        } // gy
                } // gx
        }

        T compare( const zpolbasis_type& _zp ) const {
                T dist = (T)0; // distance
                T modz[2]; // modulus of z
                for (unsigned p = 0; p <= Order; ++p) {
                        for (unsigned qi = 0; qi <= p/2; ++qi) {
                                modz[0] = std::abs( pol[p][qi] );
                                modz[1] = std::abs( _zp[p][qi] );
                                dist += (modz[0] - modz[1] ) * (modz[0] - modz[1] );
                        } // qi
                } // p
                return sqrt( dist );
        }

        radial_polynomial& operator[] ( const unsigned& _p ) { return this->pol[_p]; }

        const radial_polynomial& operator[] ( const unsigned& _p ) const { return this->pol[_p]; }

        friend std::ostream& operator << ( std::ostream& _out,
                                           const zpolbasis_type& _z ) {
                for (unsigned p = 0; p <= Order; ++p) {
                        _out << _z[p][0];
                        for (unsigned qi = 1; qi <= p/2; ++qi)
                                _out << " " << _z[p][qi];
                        if( p < Order ) _out << "\n";
                }
                return _out;
        }

        friend std::istream& operator >> ( std::istream& _in,
                                           zpolbasis_type& _z ) {
                for (unsigned p = 0; p <= Order; ++p)
                        for (unsigned qi = 0; qi <= p/2; ++qi)
                                _in >> _z[p][qi];
                return _in;
        }

private:

        std::vector< radial_polynomial > pol; 

};
//=== IMPLEMENTATION ==========================================================

template< unsigned Order, class T >
void compute_basis( ZernikePolynomialsBasisT< Order, T > **_zpb,
                    const unsigned int& _szx,
                    const unsigned int& _szy ) {

        T x, y, r, t; // (x, y) cartesian and (r, t) polar coordinates

        for (unsigned gx = 0; gx < _szx; ++gx) {

                x = ( (T)2 * gx + (T)1 ) / (T)_szx - (T)1;

                for (unsigned gy = 0; gy < _szy; ++gy) {

                        y = ( (T)2 * gy + (T)1 ) / (T)_szy - (T)1;

                        r = (T)sqrt( x*x + y*y );
                        if( r > (T)1 ) continue; // outside unit circle is zero

                        t = (T)atan2( y, x );

                        for (int p = 0; p <= Order; ++p) {

                                int qi = 0; // q index: q[0, 2, 4] -> qi[0, 1, 2]

                                for (int q = p % 2; q <= p; q += 2, qi++) {

                                        _zpb[gx][gy][p][qi] = compute_V(p, q, r, t);

                                } // q

                        } // p

                } // gy

        } // gx

}
//=============================================================================
} // namespace zsig
//=============================================================================
#endif // ZSIG_ZPOLBASIST_HH
//=============================================================================