Rescale.hpp

//
// $Id: Rescale.hpp,v 1.1 2008/04/16 21:24:28 werner Exp $
//
#ifndef __BASEOP_RESCALE_HPP
#define __BASEOP_RESCALE_HPP

#include "gridopDllApi.h"
#include <grid/Representation.hpp>

#include <vector/LinearIpol.hpp>
#include <vector/Interpolate.hpp>

namespace Fiber
{

template <int Dims, typename Type>
        struct  ResampleArray
{

static  RefPtr<MemBase> compute(const MultiIndex<Dims>&DestDims, const RefPtr<MemBase>&DataSrc)
        {
        RefPtr<ReferencingMemArray<Dims, Type> > Src = DataSrc;
                if (!Src)
                        return NullPtr();

        MultiArray<Dims,Type> & SourceData = *Src;
        MultiIndex<Dims> DimsSrc = SourceData.Size();

        RefPtr<MemArray<Dims, Type> > Dst = new MemArray<Dims, Type>( DestDims );
        MultiArray<Dims,Type> & DestData = *Dst;


        FixedArray<double, Dims> IndexScale;
                for(int i=0; i<Dims; i++)
                {
                        IndexScale[i] = double(DimsSrc[i])/DestDims[i];
                }

        MultiIndex<Dims> MI; 
        FixedArray<double,Dims> EvalPoint;
        Interpolate<Dims, Type, LinearIpol<Type> > Ipol(SourceData, EvalPoint);
                do
                {
                        for(int i=0; i<Dims; i++)
                                EvalPoint[i] = MI[i]*IndexScale[i];

                        DestData[ MI ] = Ipol.eval();
                }
                while( MI.inc(DestDims) );

                return Dst;
        }

static  RefPtr<MemBase> rescale(double scale, const RefPtr<MemBase>&DataSrc, const MemBase::Creator_t&Crec)
        {
        RefPtr<ReferencingMemArray<Dims, Type> > Src = DataSrc;
                if (!Src)
                        return NullPtr();

        MultiArray<Dims,Type> & SourceData = *Src;
        MultiIndex<Dims> DimsSrc = SourceData.Size();

        MultiIndex<Dims> DestDims;
                for(int i=0; i<Dims; i++)
                {
                        DestDims[i] = index_t( scale*DimsSrc[i] );
                } 

        RefPtr<MemArray<Dims, Type> > Dst = new MemArray<Dims, Type>( DestDims, Crec );
        MultiArray<Dims,Type> & DestData = *Dst;

        MultiIndex<Dims> MI;
        FixedArray<double,Dims> EvalPoint;
        Interpolate<Dims, Type, LinearIpol<Type> > Ipol(SourceData, EvalPoint);
                do
                {
                        for(int i=0; i<Dims; i++)
                                EvalPoint[i] = MI[i]/scale;

                        DestData[ MI ] = Ipol.eval();
                }
                while( MI.inc(DestDims) ); 

                return Dst;
        }
};


extern gridop_API RefPtr<MemBase> rescale(double scale, const RefPtr<MemBase>&Data, const MemBase::Creator_t&Crec = NullPtr() );

} // namespace Fiber

#endif