RegularCoordinateSlice.hpp

#ifndef __FIBEROPERATIONS_SLICE_SELECTION
#define __FIBEROPERATIONS_SLICE_SELECTION
 
#include <vector/MultiArray.hpp>
#include <vector/MultiCopy.hpp>
 
#include <field/MemBase.hpp>
#include <field/Field.hpp>
 
#include <eagle/FixedArray.hpp>
 
#ifdef  USE_OPENMP
#include <omp.h>
#endif
 
 
#include "gridopDllApi.h"
 
namespace Fiber
{
 
struct  gridop_API      RegularCoordinateSlice
{
        int     X, Y, Z;
 
        void    setup(int orient);
 
        RegularCoordinateSlice(int o)
        {
                setup(o);
        }
 
        MultiIndex<3>   getIndex(index_t x,
                                 index_t y,
                                 index_t z) const
        {
        MultiIndex<3> I;
 
                I[ X ] = x;
                I[ Y ] = y;
                I[ Z ] = z;
 
                return I;
        }
 
        double  getComponent(const Eagle::FixedArray<double,3>&P) const
        {
                return P[ Z ];
        }
 
 
        bool    getRectilinearIndex(double CoordinateValue,
                                    const MemCore::RefPtr<MemBase>&CoordValues,
                                    double&CoordinateIndex);
 
        bool    getRectilinearIndex(double SliceCoordinateValue,
                                    const RefPtr<Field>&Coords,
                                    const RefPtr<FragmentID>&fid,
                                    double&CoordinateIndex);
 
 
        bool    getSliceCornerIndex(MultiIndex<3>&I,
                                    const MultiIndex<3>&DataDims,
                                    int CornerX, int CornerY,
                                    index_t sliceIndex) const
        {
//              if (sliceIndex < 0) //(mr) unsigned cant be <0
//                      return false;
 
                if (sliceIndex >= DataDims[ Z ] )
                        return false;
 
                I[ X ] = CornerX*(DataDims[X]-1);
                I[ Y ] = CornerY*(DataDims[Y]-1);
                I[ Z ] = sliceIndex;
 
                return true;
        }
 
        MemCore::RefPtr<MemBase>
                getCoordinates(index_t sliceIndex,
                               const MultiIndex<3>&DataDims,
                               const RefPtr<Field>&Coords, const RefPtr<FragmentID>&fid);
 
 
        template <class Dest, class Source, class Converter>
        void    extractData(MultiArray<2, Dest>&DestArray,
                            const MultiArray<3, Source>&SrcArray,
                            const index_t&SliceNumber,
                            const Converter&C,
                            const MultiIndex<3>&SourceOffset = MIndex(0,0,0) )
        {
        MultiIndex<3> Dims = SrcArray.Size();
 
        index_t DimsX = SrcArray.Size()[X],
                DimsY = SrcArray.Size()[Y];
 
                assert( DestArray.Size()[0] >= DimsX );
                assert( DestArray.Size()[1] >= DimsY );
 
#ifdef  USE_OPENMP
#pragma omp parallel for
#endif
                for(index_t y = 0; y<DimsY; y++)
                for(index_t x = 0; x<DimsX; x++)
                {
                MultiIndex<3> P = SourceOffset;
                        P[X] += x; P[Y] += y; P[Z] += SliceNumber;
 
                MultiIndex<2> Pdst = MIndex(x,y);
                        C.copy( DestArray[ Pdst ], SrcArray[ P ]);
                }
        }
 
        bool    validDestinationSize(const MultiIndex<3>&Dst,
                                     index_t DestSliceNumber,
                                     const MultiIndex<3>&Src)
        {
                // ensure the destination array is large enough to accomodate the data
                if (Src[X] > Dst[X] ) return false;
                if (Src[Y] > Dst[Y] ) return false;
                if (DestSliceNumber >= Dst[Z] ) return false;
 
                assert( Src[X] <= Dst[X] );
                assert( Src[Y] <= Dst[Y] );
                assert( DestSliceNumber < Dst[Z] );
 
                return true;
        }
 
        template <class Dest, class Source, class Converter>
        bool    extractDataToOriented3D(MultiArray<3, Dest>&DestArray,
                                        index_t DestSliceNumber,
                                        const MultiArray<3, Source>&SrcArray,
                                        index_t SliceNumber,
                                        const Converter&C)
        {
#ifdef  VERBOSE
                printf("extractDataTo3D: src (%dx%dx%d) ==> dst (%dx%dx%d) \n",
                       int(SrcArray .Size()[0]), int(SrcArray .Size()[1]), int(SrcArray .Size()[2]),
                       int(DestArray.Size()[0]), int(DestArray.Size()[1]), int(DestArray.Size()[2]) );
                fflush(stdout);
#endif
                assert( DestArray.count()>0 );
                assert( SrcArray.count() >0 );
 
                if (!validDestinationSize(DestArray.Size(),
                                          DestSliceNumber,
                                          SrcArray.Size()) )
                {
                        printf("RegularCoordinateSlice::extractDataToOriented3D(): "
                               "Invalid (too small) Destination in slicing src (%dx%dx%d) ==> dst (%dx%dx%d) \n",
                               int(SrcArray .Size()[0]), int(SrcArray .Size()[1]), int(SrcArray .Size()[2]),
                               int(DestArray.Size()[0]), int(DestArray.Size()[1]), int(DestArray.Size()[2]) );
                        fflush(stdout);
                        return false;
                }
 
        index_t DimsX = SrcArray.Size()[X],
                DimsY = SrcArray.Size()[Y];
 
                // ensure the destination array is large enough to accomodate the data
                assert( DimsX <= DestArray.Size()[X] );
                assert( DimsY <= DestArray.Size()[Y] );
                assert( DestSliceNumber < DestArray.Size()[Z] );
 
#ifdef  USE_OPENMP
#pragma omp parallel for
#endif
                for(index_t y = 0; y<DimsY; y++)
                for(index_t x = 0; x<DimsX; x++)
                {
                MultiIndex<3>   P;
                        P[X] = x; P[Y] = y; P[Z] = SliceNumber;
 
                MultiIndex<3> Pdst = P;
                        Pdst[Z] = DestSliceNumber;
 
#if     0
                        printf("extractDataTo3D: %dx%dx%d ==> %dx%dx%d\n",
                               int(P[0]), int(P[1]), int(P[2]),
                               int(Pdst[0]), int(Pdst[1]), int(Pdst[2]) ); fflush(stdout);
#endif
 
                        C.copy( DestArray[ Pdst ], SrcArray[ P ]);
                }
                return true;
        }
};
 
 
} // namespace Fiber
 
#endif // __FIBEROPERATIONS_SLICE_SELECTION