diff options
Diffstat (limited to 'applications/rview/rviewMDD.cpp')
| -rw-r--r-- | applications/rview/rviewMDD.cpp | 1006 |
1 files changed, 1006 insertions, 0 deletions
diff --git a/applications/rview/rviewMDD.cpp b/applications/rview/rviewMDD.cpp new file mode 100644 index 0000000..4f175fc --- /dev/null +++ b/applications/rview/rviewMDD.cpp @@ -0,0 +1,1006 @@ +/* +* This file is part of rasdaman community. +* +* Rasdaman community is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* Rasdaman community is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with rasdaman community. If not, see <http://www.gnu.org/licenses/>. +* +* Copyright 2003, 2004, 2005, 2006, 2007, 2008, 2009 Peter Baumann / +rasdaman GmbH. +* +* For more information please see <http://www.rasdaman.org> +* or contact Peter Baumann via <baumann@rasdaman.com>. +/ + +/** + * PURPOSE: + * + * The rviewMDD module encapsulates various generic operations on + * MDD objects; it relies heavily on templates for this. All functions + * provided here are a) completely independent from the rest of rView + * and b) work on any kind of MDD, i.e. no restrictions whatsoever on + * dimensionality or base type. Most functions utilize + * mdd_objectFunctionType() for this which iterates over the MDD object + * basetype first and calls functions provided in the mdd_function_pointers + * structure which perform the required operation on each of the primitive + * basetypes. + * + * Functions provided are: + * + * - mdd_objectRange(): return the minimum and maximum cell values. + * - mdd_createSubcube(): creates a new MDD object by copying a subcube + * (or all of) a source MDD object; optionally persistent. + * - mdd_objectScaleInter(): creates new object by resampling (part of) + * the source object using n-linear interpolation. Use for upsampling. + * - mdd_objectScaleAverage(): creates new object by resampling (part of) + * the source object using averaging. Use for downsampling. + * - mdd_objectScaleSimple(): creates new object by scaling (part of) + * the source object using simple nearest neighbour. Fast but blocky. + * - mdd_objectChangeEndianness(): creates a new object where the + * endianness is inverted compared to the source object, or changes + * the object itself. + * + * COMMENTS: + * No comments + */ + + + +// For early templates +#ifndef _RVIEW_MDD_CPP_ +#define _RVIEW_MDD_CPP_ + + + +#ifdef __GNUG__ +#pragma implementation +#endif + + + +#include <limits.h> +#include <float.h> +#include <string.h> + + + +// Was this file included from rviewMDD.hh? +#ifdef EARLY_TEMPLATE +#ifndef __EXECUTABLE__ +#define __EXECUTABLE__ +#define _RVIEW_MDD_EXECUTABLE_ +#endif +#endif + +#include "rasodmg/transaction.hh" +#include "raslib/primitivetype.hh" +#include "raslib/structuretype.hh" +#include "raslib/endian.hh" +#include "rasodmg/alignedtiling.hh" +#include "rasodmg/storagelayout.hh" +#include "rviewMDD.hh" +#include "rviewTypes.hh" + +#ifdef _RVIEW_MDD_EXECUTABLE_ +#undef __EXECUTABLE__ +#endif + + + + +#define MDD_FIXPREC 16 + + + + + +// auxData interpreted as double[2] containing min, max, initialised externally +template<class T> +int objectRange_temp(T *dest, const T *src, const mdd_function_desc *mfd, int dim, int datastep, void *auxData) +{ + long *pos; + char *srcdata; + T **srcptr; + int i, j; + double min, max, val; + + pos = new long[dim]; srcptr = new T *[dim]; + + srcdata = objectCalcStart((const char *)src, mfd, dim); + + for (i=0; i<dim; i++) + { + pos[i] = mfd[i].useLow; srcptr[i] = (T*)srcdata; + } + + min = ((double*)auxData)[0]; max = ((double*)auxData)[1]; + + do + { + val = (double)(*srcptr[dim-1]); + if (val < min) min = val; + if (val > max) max = val; + i = dim-1; + while (i >= 0) + { + pos[i]++; srcptr[i] = (T*)(((char*)(srcptr[i])) + mfd[i].srcoff); + if (pos[i] <= mfd[i].useHigh) break; + pos[i] = mfd[i].useLow; i--; + } + if (i >= 0) + { + for (j=i+1; j<dim; j++) srcptr[j] = srcptr[i]; + } + } + while (i >= 0); + + ((double*)auxData)[0] = min; ((double*)auxData)[1] = max; + + delete [] srcptr; delete [] pos; + + return 1; +} + + +// Additional template code for objectScaleInter template +inline void objectScaleInterICore(double *srcf, double *destf, double pos, int idx1, int idx2) +{ + destf[idx1] = srcf[idx1] + pos * (srcf[idx2] - srcf[idx1]); +} + +template<class T> +int objectScaleInter_temp(T *dest, const T *src, const mdd_function_desc *mfd, int dim, int datastep, void *auxData) +{ + double *pos, *f, *fwork; + long *lastpos; + long *ipos; + long *dp; + char **srcptr; + char *srcbase; + long offset; + int i, j, k; + int slowOffset; + T *srcdata, *destdata; + long *cubeoff = (long*)auxData; + + f = new double[(1<<dim)]; fwork = new double[(1<<dim)]; + pos = new double[dim]; ipos = new long[dim]; lastpos = new long[dim]; srcptr = new char*[dim]; + dp = new long[dim]; + + // Calculate the starting point + srcbase = objectCalcStart((const char*)src, mfd, dim); + + for (i=0; i<dim; i++) + { + pos[i] = (double)(mfd[i].useLow); lastpos[i] = mfd[i].useLow; + srcptr[i] = srcbase; dp[i] = mfd[i].newLow; + } + + destdata = dest; + + // Main loop: iterate over all dimensions except for the last which is handled separately + // in its own loop for efficiency reasons. + do + { + pos[dim-1] = (double)(mfd[dim-1].useLow); dp[dim-1] = mfd[dim-1].newLow; + + // slowOffset determines whether the fast table-lookup offset or the slower dimension looping + // offset calculation can be used. Table-lookup is always possible as long as the current + // position in the cube is not at the cube's boundary. In the innermost loop slowOffset only + // has to be updated, not recalculated, because there is only a (new) overflow possible in the + // innermost loop-dimension. + slowOffset = 0; + for (i=0; i<dim; i++) + { + ipos[i] = (long)(pos[i]); + // Check for high rather than useHigh because we only get problems if its > high. + if (ipos[i] >= mfd[i].high) {ipos[i] = mfd[i].useHigh; slowOffset = 1;} + } + + srcbase = srcptr[dim-1]; + + if (slowOffset != 0) + { + for (i=0; i<(1<<dim); i++) + { + offset = 0; + for (j=0; j<dim; j++) + { + // Check against high, not useHigh (see above) + if (((i & (1<<j)) != 0) && (ipos[j] < mfd[j].high)) offset += mfd[j].srcoff; + } + srcdata = (T*)(srcbase + offset); + f[i] = (double)(*srcdata); + //if (srcdata+offset >= srchigh) cout << "Overflow1" << endl; + } + } + else + { + for (i=0; i<(1<<dim); i++) + { + srcdata = (T*)(srcbase + cubeoff[i]); + f[i] = (double)(*srcdata); + } + } + + lastpos[dim-1] = (long)(pos[dim-1]); + // Innermost loop + for (k = mfd[dim-1].newHigh - mfd[dim-1].newLow + 1; k>0; k--) + { + // n-linear interpolation + for (j=0; j<(1<<dim); j+=2) + { + objectScaleInterICore(f, fwork, (double)(pos[0] - ipos[0]), j, j+1); + } + for (i=1; i<dim; i++) + { + for (j=0; j<(1<<dim); j+=(1<<(i+1))) + { + objectScaleInterICore(fwork, fwork, (double)(pos[i] - ipos[i]), j, j + (1<<i)); + } + } + // We go in sequential order (highest coordinate first) + *destdata = (T)(fwork[0]); + destdata = (T*)(((char*)destdata) + datastep); + + // end of innermost loop + pos[dim-1] += mfd[dim-1].step; offset = (long)(pos[dim-1]); + if (offset != lastpos[dim-1]) + { + offset -= lastpos[dim-1]; lastpos[dim-1] = (long)(pos[dim-1]); + // Check against high, not useHigh + if (lastpos[dim-1] <= mfd[dim-1].high) + { + srcbase += offset*mfd[dim-1].srcoff; ipos[dim-1] = lastpos[dim-1]; + } + // high, not useHigh + if (lastpos[dim-1] >= mfd[dim-1].high) slowOffset = 1; + offset = (1<<(dim-1)); + memcpy(f, f + offset, offset*sizeof(double)); + + if (slowOffset != 0) + { + for (i=(1<<(dim-1)); i<(1<<dim); i++) + { + offset = 0; + for (j=0; j<dim-1; j++) + { + // high, not useHigh + if (((i & (1<<j)) != 0) && (ipos[j] < mfd[j].high)) offset += mfd[j].srcoff; + } + srcdata = (T*)(srcbase + offset); + f[i] = (double)(*srcdata); + //if (srcdata+offset >= srchigh) cout << "Overflow2" << endl; + } + } + else + { + for (i=(1<<(dim-1)); i<(1<<dim); i++) + { + srcdata = (T*)(srcbase + cubeoff[i - (1<<(dim-1))]); + f[i] = (double)(*srcdata); + } + } + } + } + i = dim-2; + while (i >= 0) + { + dp[i]++; + pos[i] += mfd[i].step; offset = (long)(pos[i]); + if (offset != lastpos[i]) + { + offset -= lastpos[i]; lastpos[i] = (long)(pos[i]); + // high, not useHigh + if (lastpos[i] <= mfd[i].high) + { + srcptr[i] += offset*mfd[i].srcoff; ipos[i] = lastpos[i]; + } + } + if (dp[i] <= mfd[i].newHigh) break; + dp[i] = mfd[i].newLow; pos[i] = (double)(mfd[i].useLow); lastpos[i] = mfd[i].useLow; i--; + } + if (i >= 0) + { + for (j=i+1; j<dim; j++) srcptr[j] = srcptr[i]; + } + } + while (i >= 0); + + delete [] pos; delete [] f; delete [] fwork; + delete [] ipos; delete [] lastpos; delete [] dp; delete [] srcptr; + + return 1; +} + + + + +template<class T> +int objectScaleAverage_temp(T *dest, const T *src, const mdd_function_desc *mfd, int dim, int datastep, void *aux) +{ + long *pos, *subpos, *dp; + T **srcptr, **subptr; + long lastpos, number; + double avg; + int i, j; + T *destdata; + char *srcdata; + + pos = new long[dim]; subpos = new long[dim]; dp = new long[dim]; + srcptr = new T*[dim]; subptr = new T*[dim]; + + srcdata = objectCalcStart((const char*)src, mfd, dim); + + for (i=0; i<dim; i++) + { + pos[i] = (mfd[i].useLow << MDD_FIXPREC); + dp[i] = mfd[i].newLow; srcptr[i] = (T*)srcdata; + } + + destdata = dest; + // Iterate over all cells in the destination object + do + { + //for (i=0; i<dim; i++) cout << (pos[i]>>MDD_FIXPREC) << ':' << (void*)(srcptr[i]) << ' '; cout << endl; + // start averaging over a subcube + avg = 0; number = 0; + // init subcube iteration + for (i=0; i<dim; i++) + { + subpos[i] = pos[i]; subptr[i] = srcptr[dim-1]; // no mistake! + } + do + { + avg += (double)(*(subptr[dim-1])); number++; + i = dim-1; + while (i >= 0) + { + subpos[i] += (1<<MDD_FIXPREC); + subptr[i] = (T*)(((char*)(subptr[i])) + mfd[i].srcoff); + if ((subpos[i] >> MDD_FIXPREC) < ((pos[i] + mfd[i].lstep) >> MDD_FIXPREC)) break; + subpos[i] = pos[i]; + i--; + } + if (i >= 0) + { + for (j=i+1; j<dim; j++) subptr[j] = subptr[i]; + } + } + while (i >= 0); + + *destdata = (T)(avg/number); + destdata = (T*)(((char*)destdata) + datastep); + + i = dim-1; + do + { + dp[i]++; lastpos = (pos[i] >> MDD_FIXPREC); pos[i] += mfd[i].lstep; + srcptr[i] = (T*)(((char*)(srcptr[i])) + ((pos[i] >> MDD_FIXPREC) - lastpos) * mfd[i].srcoff); + if (dp[i] <= mfd[i].newHigh) break; + dp[i] = mfd[i].newLow; pos[i] = (mfd[i].useLow << MDD_FIXPREC); + i--; + } + while (i >= 0); + if (i >= 0) + { + for (j=i+1; j<dim; j++) srcptr[j] = srcptr[i]; + } + } + while (i >= 0); + + delete [] subptr; delete [] srcptr; + delete [] dp; delete [] subpos; delete [] pos; + + return 1; +} + + + +template<class T> +int objectScaleSimple_temp(T *dest, const T *src, const mdd_function_desc *mfd, int dim, int datastep, void *aux) +{ + long *pos, *dp; + long lastpos; + T **srcptr; + T *destdata; + char *srcdata; + int i, j; + + pos = new long[dim]; dp = new long[dim]; srcptr = new T*[dim]; + + srcdata = objectCalcStart((const char*)src, mfd, dim); + + for (i=0; i<dim; i++) + { + pos[i] = (mfd[i].useLow << MDD_FIXPREC); dp[i] = mfd[i].newLow; + srcptr[i] = (T*)srcdata; + } + destdata = dest; + + do + { + *destdata = *(srcptr[dim-1]); + destdata = (T*)(((char*)destdata) + datastep); + i=dim-1; + while (i >= 0) + { + dp[i]++; lastpos = (pos[i] >> MDD_FIXPREC); pos[i] += mfd[i].lstep; + srcptr[i] = (T*)(((char*)srcptr[i]) + ((pos[i] >> MDD_FIXPREC) - lastpos)*mfd[i].srcoff); + if (dp[i] <= mfd[i].newHigh) break; + dp[i] = mfd[i].newLow; pos[i] = (mfd[i].useLow << MDD_FIXPREC); + i--; + } + if (i >= 0) + { + for (j=i+1; j<dim; j++) srcptr[j] = srcptr[i]; + } + } + while (i >= 0); + + delete [] srcptr; delete [] dp; delete [] pos; + + return 1; +} + + + +/* Must not compile non-template code when this file was just included */ +#ifndef __EXECUTABLE__ + +// in raslib/odmgtypes.hh there are the following typedefs: +// typedef unsigned char r_Boolean; +// typedef unsigned char r_Char; +// so the template instantiation is done twice for the same type. gcc 2.95.2 doesn't +// like this. + +// Sadly, this is necessary with VisualC... +#define DECLARE_MDD_FUNC(f) \ +// template int f<r_Boolean> MDD_FUNCTION_SIGNATURE(r_Boolean); \ + template int f<r_Char> MDD_FUNCTION_SIGNATURE(r_Char); \ + template int f<r_Octet> MDD_FUNCTION_SIGNATURE(r_Octet); \ + template int f<r_Short> MDD_FUNCTION_SIGNATURE(r_Short); \ + template int f<r_UShort> MDD_FUNCTION_SIGNATURE(r_UShort); \ + template int f<r_Long> MDD_FUNCTION_SIGNATURE(r_Long); \ + template int f<r_ULong> MDD_FUNCTION_SIGNATURE(r_ULong); \ + template int f<r_Float> MDD_FUNCTION_SIGNATURE(r_Float); \ + template int f<r_Double> MDD_FUNCTION_SIGNATURE(r_Double); + +DECLARE_MDD_FUNC(objectScaleInter_temp); +DECLARE_MDD_FUNC(objectScaleAverage_temp); +DECLARE_MDD_FUNC(objectScaleSimple_temp); +DECLARE_MDD_FUNC(objectRange_temp); + + + +// Calculate the starting point in the source cube +char *objectCalcStart(const char *src, const mdd_function_desc *mfd, int dim) +{ + int i; + char *d; + + d = (char*)src; + for (i=0; i<dim; i++) + { + d += (mfd[i].useLow - mfd[i].low) * mfd[i].srcoff; + } + return d; +} + + +/* + * Wrapper functions for those templates + */ + +#define CAST_MDDPTR(t, src) \ + r_Ref<r_Marray<t> > mddPtr = (r_Ref<r_Marray<t> >)(src) + + +int mdd_createSubcube(r_Ref<r_GMarray> srcMdd, r_Ref<r_GMarray> &newMdd, r_Minterval *domain, r_Database *db) +{ + r_Dimension dim; + char *destdata, *srcdata; + r_Ref<r_GMarray> newMddPtr; + r_Minterval interv; + r_Minterval newInterv; + int tpsize; + int i, j; + long length; + mdd_function_desc *mfd; + long *pos; + char **srcptr; + + interv = srcMdd->spatial_domain(); dim = interv.dimension(); + + if (dim < 1) return 0; + + tpsize = srcMdd->get_type_length(); + + newInterv = r_Minterval(dim); + for (i=0; i<(int)dim; i++) + { + if (domain == NULL) + { + newInterv << interv[i]; + } + else + { + newInterv << (*domain)[i]; + } + } + + if ((destdata = mdd_objectFunctionInitMdd(srcMdd, newMddPtr, newInterv, tpsize, dim, db)) == NULL) + { + cerr << "Couldn't create new MDD object" << endl; + return 0; + } + + mfd = mdd_objectFunctionInitData(interv, newInterv, newInterv, tpsize, MDD_OBJECT_INIT_NEWIV); + + srcdata = objectCalcStart(srcMdd->get_array(), mfd, (int)dim); + + pos = new long[dim]; srcptr = new char *[dim]; + + for (i=0; i<(int)dim; i++) + { + pos[i] = mfd[i].useLow; srcptr[i] = srcdata; + } + + // Lenght of one strip of data to copy + length = (mfd[dim-1].useHigh - mfd[dim-1].useLow + 1) * tpsize; + do + { + memcpy(destdata, srcptr[dim-1], length); + destdata += length; + i = (int)dim-2; + while (i >= 0) + { + pos[i]++; srcptr[i] += mfd[i].srcoff; + if (pos[i] <= mfd[i].useHigh) break; + pos[i] = mfd[i].useLow; i--; + } + if (i >= 0) + { + for (j=i+1; j<(int)dim; j++) srcptr[j] = srcptr[i]; + } + } + while (i >= 0); + + delete [] srcptr; delete [] pos; + delete [] mfd; + + newMdd = newMddPtr; + + return 1; +} + + + +// Support functions for (recursive) member-by-member resampling of a data cube via interpolation +int mdd_objectFunctionPrim(const mdd_function_pointers *mfp, r_Primitive_Type *primType, const char *src, char *dest, const mdd_function_desc *mfd, int dim, int tpsize, void *auxData) +{ + switch (primType->type_id()) + { + case r_Primitive_Type::BOOL: + if (mfp->mddf_bool != NULL) + mfp->mddf_bool((r_Boolean*)dest, (const r_Boolean *)src, mfd, dim, tpsize, auxData); + break; + case r_Primitive_Type::CHAR: + if (mfp->mddf_char != NULL) + mfp->mddf_char((r_Char*)dest, (const r_Char *)src, mfd, dim, tpsize, auxData); + break; + case r_Primitive_Type::OCTET: + if (mfp->mddf_octet != NULL) + mfp->mddf_octet((r_Octet *)dest, (const r_Octet*)src, mfd, dim, tpsize, auxData); + break; + case r_Primitive_Type::SHORT: + if (mfp->mddf_short != NULL) + mfp->mddf_short((r_Short *)dest, (const r_Short*)src, mfd, dim, tpsize, auxData); + break; + case r_Primitive_Type::USHORT: + if (mfp->mddf_ushort != NULL) + mfp->mddf_ushort((r_UShort *)dest, (const r_UShort*)src, mfd, dim, tpsize, auxData); + break; + case r_Primitive_Type::LONG: + if (mfp->mddf_long != NULL) + mfp->mddf_long((r_Long *)dest, (const r_Long*)src, mfd, dim, tpsize, auxData); + break; + case r_Primitive_Type::ULONG: + if (mfp->mddf_ulong != NULL) + mfp->mddf_ulong((r_ULong *)dest, (const r_ULong*)src, mfd, dim, tpsize, auxData); + break; + case r_Primitive_Type::FLOAT: + if (mfp->mddf_float != NULL) + mfp->mddf_float((r_Float *)dest, (const r_Float*)src, mfd, dim, tpsize, auxData); + break; + case r_Primitive_Type::DOUBLE: + if (mfp->mddf_double != NULL) + mfp->mddf_double((r_Double *)dest, (const r_Double*)src, mfd, dim, tpsize, auxData); + break; + default: return 0; + } + return 1; +} + + +// Support functions for (recursive) member-by-member resampling of a data cube via interpolation +int mdd_objectFunctionStruct(const mdd_function_pointers *mfp, r_Structure_Type *structType, const char *src, char *dest, const mdd_function_desc *mfd, int dim, int tpsize, void *auxData) +{ + r_Type *newType; + unsigned long offset; + char *inc, *outc; + + r_Structure_Type::attribute_iterator iter(structType->defines_attribute_begin()); + while (iter != structType->defines_attribute_end()) + { + newType = (*iter).type_of().clone(); + offset = (*iter).offset(); + inc = (char*)src + offset; outc = (char*)dest + offset; + + if (newType->isStructType()) + { + r_Structure_Type *newStructType = (r_Structure_Type*)newType; + mdd_objectFunctionStruct(mfp, newStructType, inc, outc, mfd, dim, tpsize, auxData); + } + else + { + r_Primitive_Type *newPrimType = (r_Primitive_Type*)newType; + mdd_objectFunctionPrim(mfp, newPrimType, inc, outc, mfd, dim, tpsize, auxData); + } + delete newType; + iter++; + } + return 1; +} + + +int mdd_objectFunctionType(const mdd_function_pointers *mfp, const r_Type *baseType, const char *src, char *dest, const mdd_function_desc *mfd, int dim, int tpsize, void *auxData) +{ + if (((r_Type*)baseType)->isStructType()) + { + r_Structure_Type *structType = (r_Structure_Type*)baseType; + return mdd_objectFunctionStruct(mfp, structType, src, dest, mfd, dim, tpsize, auxData); + } + else + { + r_Primitive_Type *primType = (r_Primitive_Type*)baseType; + return mdd_objectFunctionPrim(mfp, primType, src, dest, mfd, dim, tpsize, auxData); + } +} + + +char *mdd_objectFunctionInitMdd(r_Ref<r_GMarray> mddPtr, r_Ref<r_GMarray> &newMddPtr, r_Minterval &newInterv, int tpsize, r_Dimension dim, r_Database *db) +{ + char *destdata; + const char *b; + r_Storage_Layout *store; + + if (mddPtr->get_storage_layout() != NULL) + { + store = new r_Storage_Layout(*(mddPtr->get_storage_layout())); + } + else + { + r_Minterval tileInterv(dim); + r_Aligned_Tiling *dfltTiling; + int dim, i; + + dim = newInterv.dimension(); + for (i=0; i<dim; i++) + { + tileInterv << r_Sinterval((r_Range)0, (r_Range)256); + } + dfltTiling = new r_Aligned_Tiling(tileInterv, tpsize * tileInterv.cell_count()); + store = new r_Storage_Layout(dfltTiling); + } + + if (db == NULL) + { + newMddPtr = new r_GMarray(newInterv, tpsize, store); + } + else + { + newMddPtr = new (db) r_GMarray(newInterv, tpsize, store); + } + + destdata = newMddPtr->get_array(); + newMddPtr->set_current_format(mddPtr->get_current_format()); + if ((b = mddPtr->get_type_name()) != NULL) newMddPtr->set_type_by_name(b); + if ((b = mddPtr->get_type_structure()) != NULL) newMddPtr->set_type_structure(b); + + return destdata; +} + + +mdd_function_desc *mdd_objectFunctionInitData(r_Minterval &interv, r_Minterval &useInterv, r_Minterval &newInterv, int tpsize, unsigned int flags) +{ + int i, offset; + r_Dimension dim; + mdd_function_desc *mfd; + + dim = interv.dimension(); + mfd = new mdd_function_desc[dim]; + + // Calculate scaling factors + for (i=0; i<(int)dim; i++) + { + mfd[i].low = interv[i].low(); mfd[i].high = interv[i].high(); + mfd[i].useLow = useInterv[i].low(); mfd[i].useHigh = useInterv[i].high(); + if ((flags & MDD_OBJECT_INIT_NEWIV) != 0) + { + mfd[i].newLow = newInterv[i].low(); mfd[i].newHigh = newInterv[i].high(); + if ((flags & MDD_OBJECT_INIT_FPSTEP) == MDD_OBJECT_INIT_FPSTEP) + { + // In order to avoid doubling right-boundary pixels we have to use (h-l), not (h-l+1) + // as nominator + mfd[i].step = ((double)(mfd[i].useHigh - mfd[i].useLow)) / (mfd[i].newHigh - mfd[i].newLow + 1); + } + else + { + double h = ((double)(mfd[i].useHigh - mfd[i].useLow + 1)) / (mfd[i].newHigh - mfd[i].newLow + 1); + mfd[i].lstep = (long)(h * (1<<MDD_FIXPREC)); + } + } + } + + offset = tpsize; + for (i=(int)dim-1; i>=0; i--) + { + mfd[i].srcoff = offset; + offset *= (mfd[i].high - mfd[i].low + 1); + } + + return mfd; +} + + + +int mdd_objectRange(r_Ref<r_GMarray> mddPtr, r_Minterval &useInterv, double &min, double &max) +{ + r_Dimension dim; + double minmax[2]; + r_Ref<r_GMarray> newMddPtr; // Dummy... + r_Minterval newInterv; // Dummy... + r_Minterval interv; + const r_Type *baseType; + int tpsize; + mdd_function_desc *mfd; + mdd_function_pointers mfp; + + interv = mddPtr->spatial_domain(); dim = interv.dimension(); + + if (dim < 1) return 0; + + if ((baseType = mddPtr->get_base_type_schema()) == NULL) + { + cerr << "No schema information available!" << endl; + return 0; + } + + tpsize = mddPtr->get_type_length(); + + MDD_INIT_FUNCTIONS(mfp, objectRange_temp); + + mfd = mdd_objectFunctionInitData(interv, useInterv, newInterv, tpsize, 0); + + minmax[0] = DBL_MAX; minmax[1] = -DBL_MAX; // these are unsigned!!! + + mdd_objectFunctionType(&mfp, baseType, mddPtr->get_array(), NULL, mfd, (int)dim, tpsize, (void*)minmax); + + min = minmax[0]; max = minmax[1]; + + delete [] mfd; + + return 1; +} + + + +// Scale an object using n-linear interpolation (recommended for magnification) +int mdd_objectScaleInter(r_Ref<r_GMarray> mddPtr, r_Minterval &useInterv, r_Ref<r_GMarray> &newMdd, r_Minterval &newInterv) +{ + r_Dimension dim; + r_Minterval interv; + r_Ref<r_GMarray> newMddPtr; + mdd_function_desc *mfd; + int i, j; + long *cubeoff; + char *srcdata, *destdata; + int tpsize; + long offset; + const r_Type *baseType; + mdd_function_pointers mfp; + + interv = mddPtr->spatial_domain(); + dim = interv.dimension(); + if ((dim < 1) || (dim != newInterv.dimension())) return 0; + + // Type information available? + if ((baseType = mddPtr->get_base_type_schema()) == NULL) + { + cerr << "No schema information available!" << endl; + return 0; + } + + tpsize = mddPtr->get_type_length(); + + if ((destdata = mdd_objectFunctionInitMdd(mddPtr, newMddPtr, newInterv, tpsize, dim)) == NULL) + return 0; + + MDD_INIT_FUNCTIONS(mfp, objectScaleInter_temp); + + // create temporary arrays + mfd = mdd_objectFunctionInitData(interv, useInterv, newInterv, tpsize, MDD_OBJECT_INIT_FPSTEP); + cubeoff = new long[(1<<dim)]; + + // Debugging + /*T *srchigh; + offset=1; + for (i=0; i<dim; i++) offset *= (interv[i].high() - interv[i].low() + 1); + srchigh = (T*)(mddPtr->get_array()) + offset;*/ + + srcdata = (char*)(mddPtr->get_array()); + + // Precalculate cube-offsets for more speed in inner loop + for (i=0; i<(1<<dim); i++) + { + offset = 0; + for (j=0; j<(int)dim; j++) + { + if ((i & (1<<j)) != 0) offset += mfd[j].srcoff; + } + cubeoff[i] = offset; + } + + mdd_objectFunctionType(&mfp, baseType, srcdata, destdata, mfd, (int)dim, tpsize, (void*)cubeoff); + + delete [] cubeoff; delete [] mfd; + + newMdd = newMddPtr; + + return 1; +} + + +// Scale an object by averaging cell values -- recommended for scaling down +int mdd_objectScaleAverage(r_Ref<r_GMarray> mddPtr, r_Minterval &useInterv, r_Ref<r_GMarray> &newMdd, r_Minterval &newInterv) +{ + r_Dimension dim; + r_Minterval interv; + r_Ref<r_GMarray> newMddPtr; + mdd_function_desc *mfd; + char *srcdata, *destdata; + int tpsize; + const r_Type *baseType; + mdd_function_pointers mfp; + + interv = mddPtr->spatial_domain(); + dim = interv.dimension(); + + if ((dim < 1) || (dim != newInterv.dimension())) return 0; + + if ((baseType = mddPtr->get_base_type_schema()) == NULL) + { + cerr << "No schema information available!" << endl; + return 0; + } + + tpsize = mddPtr->get_type_length(); + + if ((destdata = mdd_objectFunctionInitMdd(mddPtr, newMddPtr, newInterv, tpsize, dim)) == NULL) + return 0; + + MDD_INIT_FUNCTIONS(mfp, objectScaleAverage_temp); + + mfd = mdd_objectFunctionInitData(interv, useInterv, newInterv, tpsize, MDD_OBJECT_INIT_NEWIV); + + srcdata = mddPtr->get_array(); + + mdd_objectFunctionType(&mfp, baseType, srcdata, destdata, mfd, (int)dim, tpsize, NULL); + + delete [] mfd; + + newMdd = (r_Ref<r_GMarray>)newMddPtr; + + return 1; +} + + +int mdd_objectScaleSimple(r_Ref<r_GMarray> mddPtr, r_Minterval &useInterv, r_Ref<r_GMarray> &newMdd, r_Minterval &newInterv) +{ + r_Dimension dim; + r_Minterval interv; + r_Ref<r_GMarray> newMddPtr; + mdd_function_desc *mfd; + char *srcdata, *destdata; + int tpsize; + const r_Type *baseType; + mdd_function_pointers mfp; + + interv = mddPtr->spatial_domain(); + dim = interv.dimension(); + + if ((dim < 1) || (dim != newInterv.dimension())) return 0; + + if ((baseType = mddPtr->get_base_type_schema()) == NULL) + { + cerr << "No schema information available!" << endl; + return 0; + } + + tpsize = mddPtr->get_type_length(); + + if ((destdata = mdd_objectFunctionInitMdd(mddPtr, newMddPtr, newInterv, tpsize, dim)) == NULL) + return 0; + + MDD_INIT_FUNCTIONS(mfp, objectScaleSimple_temp); + + mfd = mdd_objectFunctionInitData(interv, useInterv, newInterv, tpsize, MDD_OBJECT_INIT_NEWIV); + + srcdata = mddPtr->get_array(); + + mdd_objectFunctionType(&mfp, baseType, srcdata, destdata, mfd, (int)dim, tpsize, NULL); + + newMdd = (r_Ref<r_GMarray>)newMddPtr; + + delete [] mfd; + + return 1; +} + + +int mdd_objectChangeEndianness(r_Ref<r_GMarray> mddPtr, r_Minterval &useInterv, r_Ref<r_GMarray> *newMdd, r_Minterval *newInterv) +{ + char *srcdata; + const r_Base_Type *baseType; + r_Ref<r_GMarray> newMddPtr; // only used in "new" mode + r_Minterval dom, *iterDom; + int tpsize; + + if ((baseType = (r_Base_Type*)(mddPtr->get_base_type_schema())) == NULL) + { + cerr << "No schema information available!" << endl; + return 0; + } + + tpsize = mddPtr->get_type_length(); + + srcdata = mddPtr->get_array(); + + dom = mddPtr->spatial_domain(); + + if (newMdd == NULL) + { + r_Endian::swap_array( baseType, dom, dom, srcdata, srcdata ); + } + else + { + void *destdata; + r_Dimension dim = dom.dimension(); + + if ((destdata = mdd_objectFunctionInitMdd(mddPtr, newMddPtr, *newInterv, tpsize, dim)) == NULL) + return 0; + + r_Endian::swap_array( baseType, dom, *newInterv, *newInterv, *newInterv, srcdata, destdata ); + } + + if (newMdd != NULL) + { + *newMdd = newMddPtr; + } + + return 1; +} + +#endif + +#endif |