Resample.hpp

#ifndef __BASEOP_RESAMPLE_HPP
#define __BASEOP_RESAMPLE_HPP
 
#include <memcore/RefPtr.hpp>
#include <memcore/Verbose.hpp>
 
#include "gridopDllApi.h"
#include <grid/Representation.hpp>
 
#include <vector/LinearIpol.hpp>
#include <vector/Interpolate.hpp>
#include <vector/MultiCopy.hpp>
 
 
namespace Fiber
{
 
template <Dims_t Dims, typename DestType, typename SrcType, class Converter = PlainCopyConvert>
MultiArray<Dims,DestType>&Resample(MultiArray<Dims,DestType>&DestData
                                   , const MultiArray<Dims,SrcType>&SourceData
                                   , const MultiIndex<Dims>&DestOffset = MultiIndex<Dims>()
                                   , const MultiIndex<Dims>&SourceOffset = MultiIndex<Dims>()
                                   , const Converter&C = Converter()
                                   )
{
const MultiIndex<Dims>&SrcDims  = SourceData.Size();
const MultiIndex<Dims>&DstDims  = DestData.Size();
 
/*
Given the destination index we get the float-point source index via:
 
SrcIndex = DstIndex*SrcDims/DstDims;
 
   Rescaling is always specified via integer dimensions to ensure
   consistency of the dimensions, though the interpolation could
   handly arbitrary floating point mappings. Do we have a need
   for exact index scaling without rounding to integer factors?
 */
FixedArray<double, Dims> IndexScale;
        for(int i=0; i<Dims; i++)
        {
                IndexScale[i] = SrcDims[i] / double(DstDims[i]);
        }
 
#pragma message "TODO: Check whether the index scale is an integer in all dimensions, and if so, use MultiCopy.hpp's Copy operation"
 
FixedArray<double,Dims> SrcIndex;
#pragma message "TODO: Pass Interpolation type - linear/cubic - as template argument"
Interpolate<Dims, SrcType, LinearIpol<SrcType> > Ipol(SourceData, SrcIndex);
 
        Verbose(30) << " Resampling over " << DstDims << " destination indices with scale " << IndexScale;
 
        for(auto DstIndex : DstDims)
        {
                for(int i=0; i<Dims; i++)
                        SrcIndex[i] = DstIndex[i]*IndexScale[i];
 
                C.copy( DestData[ DstIndex], Ipol.eval() );
        }
 
        return DestData;
}
 
} // namespace Fiber
 
 
#include <field/RegularlyFragmentedField.hpp>
 
namespace Fiber
{
 
template <Dims_t Dims>
struct  CreativeFieldResampler : FieldFragmentIterator
{
        FixedArray<double, Dims> IndexScale;
        Field::CreationFailAlternative CreationMode = Field::CreationFailAlternative::Discardable;
 
        RefPtr<Field>   OutputField;
        RefPtr<Field>   InputField;
 
        CreativeFieldResampler()
        {
                IndexScale.set(1.0);
        }
 
        RefPtr<Field>&createEmptyField(const RefPtr<Field>&theInputField,
                                       const RefPtr<FragmentIDCollection>&theOutputFragmentIDCollection)
        {
                InputField = theInputField;
 
                if (RefPtr< RegularlyFragmentedField<Dims> > RegularInputField = InputField)
                {
#pragma warning "reuse fragment ID collection!?"
//              const RefPtr<CreatorConstructionContext> CT;
 
                        OutputField = new RegularlyFragmentedField<Dims>( RegularInputField->getNumberOfFragments() );
 
                        Verbose(20) << "CreativeFieldResampler created empty field with regular fragmentation, using "
                                    << RegularInputField->getNumberOfFragments() << " fragments";
 
                }
                else
                {
                        OutputField = new Field( theOutputFragmentIDCollection );
 
                        Verbose(20) << "CreativeFieldResampler created empty field with no predefined fragmentation";
                }
                return OutputField;
        }
 
        RefPtr<Field>&compute()
        {
                if (OutputField)
                {
                        InputField->iterate( *this );
                }
 
                return  OutputField;
        }
};
 
 
template <int Dims, typename DestType, typename SrcType = DestType, class Converter = PlainCopyConvert>
struct  CreativeResampleFragmentIterator : CreativeFieldResampler<Dims>
{
        bool    apply(const RefPtr<FragmentID>&f, const RefPtr<CreativeArrayBase>&SourceFragment) override
        {
                Assert(SourceFragment);
 
                Verbose(20) << "CreativeResampleFragmentIterator: Iterating over a source fragment of type "
                            << typeid(SrcType) << " to be put into destination type " << typeid(DestType);
 
                if (RefPtr<ReferencingMemArray<Dims, SrcType>> SrcData = SourceFragment->create())
                {
                MultiIndex<Dims> DestDims;
                        for(int i=0; i<Dims; i++)
                        {
//                              DestDims[i] = index_t( SrcData->Size()[i] * this->IndexScale[i] );
                                DestDims[i] = index_t( SrcData->Size()[i] * this->IndexScale[i] ) -1;
                        }
 
                        Verbose(20) << "CreativeResampleFragmentIterator: Creating a destination fragment of size " << DestDims;
 
                RefPtr<MemArray<Dims, DestType> > DestData = new MemArray<Dims, DestType>( DestDims );
 
                        Resample( *DestData, *SrcData);
 
                        if (!this->OutputField->createCreator( DestData, f, this->CreationMode))
                        {
                                Verbose(0) << "CreativeResampleFragmentIterator: Could not create a Creator, which is not so good.";
                        }
                }
                else
                {
                        if (RefPtr<MemBase> SourceData = SourceFragment->create() )
                                Verbose(15) << "SourceData is " << SourceData;
                        else
                                Verbose(15) << "No acceptable input data type was found.";
                }
 
                return true;
        }
};
 
 
struct  gridop_API CreativeResampleUniformCartesian3DIterator : CreativeFieldResampler<3>
{
        bool    apply(const RefPtr<FragmentID>&f, const RefPtr<CreativeArrayBase>&SourceFragment) override;
};
 
} // namespace Fiber
 
 
#include <grid/Grid.hpp>
 
 
namespace Fiber
{
 
extern gridop_API
RefPtr<Field> ResampleRegularGrid3D(Grid&OutputGrid, const Grid&InputGrid,
                                    const string&Fieldname,
                                    const FixedArray<double, 3>&IndexScale,
                                    const std::string&ChartName = std::string() );
 
 
} // namespace Fiber
 
#endif