convolve2d.h

Go to the documentation of this file.

#pragma once

/*@+
________________________________________________________________________

 (C) dGB Beheer B.V.; (LICENSE) http://opendtect.org/OpendTect_license.txt
 Author:        Kristofer Tingdahl
 Date:          07-10-1999
________________________________________________________________________


*/

#include "odcomplex.h"
#include "arraynd.h"
#include "paralleltask.h"
#include "rowcol.h"

#if defined(__msvc__) && (defined(ALGO_EXPORTS) || defined(Algo_EXPORTS))
# define mAlgoClass     class dll_export
#else
# define mAlgoClass     class
#endif

namespace Fourier { class CC; }

template <class T>
mAlgoClass Convolver2D : public ParallelTask
{
public:
                                Convolver2D();
                                ~Convolver2D();

    inline void                 setX(const Array2D<T>&,bool hasudfs);
    inline const Array2D<T>*    getX() const            { return x_; }
    inline void                 setY(const Array2D<T>&,bool hasudfs);
    inline const Array2D<T>*    getY() const            { return x_; }
    inline void                 setZ(Array2D<T>& z )    { z_ = &z; }
    void                        setNormalize( bool n )  { normalize_ = n; }
    void                        setCorrelate( bool yn ) { correlate_ = yn; }
    od_int64    nrIterations() const;

protected:
    bool                doWork(od_int64,od_int64,int);
    bool                doNonFFTWork(od_int64,od_int64,int);
    bool                doPrepare(int);
    bool                shouldFFT() const;

    const Array2D<T>*   x_;
    const Array2D<T>*   y_;
    Array2D<T>*         z_;
    bool                normalize_;
    bool                correlate_;
    bool                xhasudfs_;
    bool                yhasudfs_;

    float_complex*      xf_;
    float_complex*      yf_;
    float_complex*      zf_;

    bool                updatexf_;
    bool                updateyf_;

    Fourier::CC*        fft_;
};




template <> inline
bool Convolver2D<float>::shouldFFT() const
{
    if ( xhasudfs_ || yhasudfs_ || x_->info()!=y_->info() ||
         x_->info()!=z_->info() )
        return false;

    return true;
}


template <> inline
od_int64 Convolver2D<float>::nrIterations() const
{
    return shouldFFT()
        ? 1
        : z_->info().getTotalSz();
}


template <class T> inline
bool Convolver2D<T>::shouldFFT() const
{ return false; }


template <class T> inline
Convolver2D<T>::Convolver2D()
    : x_( 0 )
    , y_( 0 )
    , z_( 0 )
    , normalize_( false )
    , correlate_( false )
    , xhasudfs_( false )
    , yhasudfs_( false )
    , updatexf_( true )
    , updateyf_( true )
    , xf_( 0 )
    , yf_( 0 )
    , zf_( 0 )
    , fft_( 0 )
{}


template <>
Convolver2D<float>::~Convolver2D();


template <class T>
Convolver2D<T>::~Convolver2D()
{}

template <class T> inline
void Convolver2D<T>::setX( const Array2D<T>& x, bool hasudfs )
{
    x_ = &x;
    updatexf_ = true;
    xhasudfs_ = hasudfs;
}


template <class T> inline
void Convolver2D<T>::setY( const Array2D<T>& y, bool hasudfs )
{
    y_ = &y;
    updateyf_ = true;
    yhasudfs_ = hasudfs;
}




#define mConvolver2DSetIndex( dim ) \
const char y##dim##inc = correlate_ ? 1 : -1; \
const int idy##dim = firsty##dim+idx##dim*y##dim##inc; \
if ( idy##dim<0 ) \
{ \
    if ( correlate_ ) \
    { \
        idx##dim += (-idy##dim)-1; \
        continue; \
    } \
 \
    break; \
} \
 \
if ( idy##dim>=ysz##dim ) \
{ \
    if ( correlate_ ) \
        break; \
 \
    const int diff = idy##dim-(ysz##dim-1); \
    idx##dim += diff-1; \
    continue; \
}


template <class T> inline
od_int64 Convolver2D<T>::nrIterations() const
{
    return shouldFFT()
        ? 1
        : z_->info().getTotalSz();
}


template <>
bool Convolver2D<float>::doWork( od_int64 start, od_int64 stop, int );


template <class T> inline
bool Convolver2D<T>::doWork( od_int64 start, od_int64 stop, int thread )
{
    return doNonFFTWork( start, stop, thread );
}


template <class T> inline
bool Convolver2D<T>::doNonFFTWork( od_int64 start, od_int64 stop, int )
{
    const int xsz0 = x_->info().getSize( 0 );
    const int xsz1 = x_->info().getSize( 1 );
    const int ysz0 = y_->info().getSize( 0 );
    const int ysz1 = y_->info().getSize( 1 );

    int startpos[2];

    if ( !z_->info().getArrayPos( start, startpos ) )
        return false;

    ArrayNDIter iterator( z_->info() );
    iterator.setPos( startpos );

    const ValueSeries<T>* xstor_ = x_->getStorage();
    const T* xptr_ = x_->getData();

    const ValueSeries<T>* ystor_ = y_->getStorage();
    const T* yptr_ = y_->getData();

    for ( od_int64 idx=start; idx<=stop; idx++ )
    {
        const int* zvar = iterator.getPos();
        const int firsty0 = correlate_ ? -zvar[0] : zvar[0];
        const int firsty1 = correlate_ ? -zvar[1] : zvar[1];
        T sum = 0;
        T ysum = 0;
        int nrsamples = 0;

        for ( int idx0=0; idx0<xsz0 && shouldContinue(); idx0++ )
        {
            mConvolver2DSetIndex( 0 );

            const od_int64 yoffset = ystor_ || yptr_ ?
                y_->info().getOffset( idy0, 0 ) : 0;

            const od_int64 xoffset = xstor_ || xptr_ ?
                x_->info().getOffset( idx0, 0 ) : 0;

            for ( int idx1=0; idx1<xsz1; idx1++ )
            {
                mConvolver2DSetIndex( 1 );

                const T yval = yptr_
                    ? yptr_[yoffset+idy1]
                    : ystor_
                        ? ystor_->value( yoffset+idy1 )
                        : y_->get( idy0, idy1 );

                if ( mIsUdf(yval) )
                    continue;

                const T xval = xptr_
                    ? xptr_[xoffset+idx1]
                    : xstor_
                        ? xstor_->value( xoffset+idx1 )
                        : x_->get( idx0, idx1 );

                if ( mIsUdf(xval) )
                    continue;

                sum += xval * yval;
                ysum += yval;
                nrsamples++;
            }
        }

        if ( !nrsamples ) z_->setND( zvar, 0 );
        else if ( normalize_ && !mIsZero(ysum,1e-8) )
            z_->setND( zvar, sum/ysum );
        else z_->setND( zvar, sum );

        addToNrDone( 1 );

        if ( !iterator.next() && idx!=stop )
            return false;
    }

    return true;
}