|| || || || || || || ||

---

FuncGeometricOp.h

Go to the documentation of this file.

#ifndef Impala_Core_Array_Pattern_FuncGeometricOp_h
#define Impala_Core_Array_Pattern_FuncGeometricOp_h

#include "Core/Array/Pattern/ArrayFunc.h"
#include "Core/Array/Pattern/PtrFunc.h"
#include "Core/Matrix/MatMulVec.h"
#include "Core/Array/Element/E1Cast.h"

namespace Impala
{
namespace Core
{
namespace Array
{
namespace Pattern
{


template <class DstArrayT, class SrcArrayT>
void
FuncGeometricOp(DstArrayT* dst, SrcArrayT* src, Matrix::Mat* func,
                Geometry::GeoIntType gi, Element::Vec3Real64 translate,
                typename DstArrayT::ArithType background)
{
    typedef typename DstArrayT::StorType DstStorT;
    typedef typename DstArrayT::ArithType DstArithT;
    typedef typename SrcArrayT::StorType SrcStorT;
    typedef typename SrcArrayT::ArithType SrcArithT;

    Real64 srcWidth = ArrayCW(src);
    Real64 srcHeight = ArrayCH(src);
    Real64 maxSrcX = (gi == Geometry::LINEAR) ? srcWidth - 1 : srcWidth - 0.5;
    Real64 maxSrcY = (gi == Geometry::LINEAR) ? srcHeight - 1 : srcHeight - 0.5;
    int dstWidth = ArrayCW(dst);
    int dstHeight = ArrayCH(dst);

    for (int y=0 ; y<dstHeight ; y++)
    {
        DstStorT* dstPtr = ArrayCPB(dst, 0, y);
        for (int x=0 ; x<dstWidth ; x++)
        {
            Element::Vec3Real64 vDst(x, y, 1);
            Element::Vec3Real64 vSrc = Matrix::MatMulVec(func, vDst + translate);
            Real64 srcXf = vSrc.X() / vSrc.Z(); // homogeneous coordinates
            Real64 srcYf = vSrc.Y() / vSrc.Z();

            DstArithT pixVal;
            if ((srcXf<0) || (srcXf>=maxSrcX) || (srcYf<0) || (srcYf>=maxSrcY))
                pixVal = background; // outside image
            else
            {
                if (gi == Geometry::LINEAR)
                {
                    Int32 srcX = Int32(srcXf); // truncate
                    Int32 srcY = Int32(srcYf);
                    DstArithT alpha = Element::E1Cast(Real64(srcXf - srcX), DstArithT());
                    DstArithT beta  = Element::E1Cast(Real64(srcYf - srcY), DstArithT());
                    SrcStorT* sPtr;
                    sPtr = ArrayCPB(src, srcX, srcY);
                    DstArithT v1 = PtrRead(sPtr, SrcArithT());
                    sPtr += SrcArrayT::ElemSize();
                    DstArithT v2 = PtrRead(sPtr, SrcArithT());
                    sPtr = ArrayCPB(src, srcX, srcY+1);
                    DstArithT v3 = PtrRead(sPtr, SrcArithT());
                    sPtr += SrcArrayT::ElemSize();
                    DstArithT v4 = PtrRead(sPtr, SrcArithT());
                    pixVal = v1 + alpha*(v2-v1) + beta*(v3-v1)
                                + alpha*beta*(v1-v2-v3+v4);
                }
                else
                { // gi == NEAREST
                    SrcStorT* sPtr = ArrayCPB(src, Int32(srcXf + 0.5),
                                                Int32(srcYf + 0.5));
                    pixVal = PtrRead(sPtr, SrcArithT());
                }
            }
            PtrWrite(dstPtr, pixVal);
            dstPtr += DstArrayT::ElemSize();
        }
    }
}

} // namespace Pattern
} // namespace Array
} // namespace Core
} // namespace Impala

#endif

---