00001 #ifndef Impala_Core_Array_Pattern_PatGeometricOp_h
00002 #define Impala_Core_Array_Pattern_PatGeometricOp_h
00003
00004 #include "Core/Array/Element/ArithTypes.h"
00005 #include "Core/Matrix/MatCopy.h"
00006 #include "Core/Matrix/MatInverse.h"
00007 #include "Core/Matrix/MatMulVec.h"
00008 #include "Core/Array/Element/E2Inf.h"
00009 #include "Core/Array/Element/E2Sup.h"
00010 #include "Core/Geometry/GeoIntType.h"
00011 #include "Core/Array/Pattern/FuncGeometricOp.h"
00012
00013 #ifdef PX_HORUS_USED
00014 #include "Core/Array/Pattern/PxArrayFunc.h"
00015 #include "Core/Array/Pattern/PxStateTrans.h"
00016 #endif
00017
00018 namespace Impala
00019 {
00020 namespace Core
00021 {
00022 namespace Array
00023 {
00024 namespace Pattern
00025 {
00026
00027
00028 template<class DstArrayT, class SrcArrayT>
00029 inline void
00030 PatGeometricOp(DstArrayT*& dst, SrcArrayT* src, Matrix::Mat* func,
00031 bool isForward, Geometry::GeoIntType gi, bool adjustSize,
00032 typename DstArrayT::ArithType background)
00033 {
00034 if ((Matrix::MatNrRow(func) != 3) || (Matrix::MatNrCol(func) != 3))
00035 {
00036 std::cout << "PatGeometricOp: matrix needs to be 3x3" << std::endl;
00037 return;
00038 }
00039
00040 Matrix::Mat* fForw = (isForward) ? Matrix::MatCopy(func) : Matrix::MatInverse(func);
00041 Matrix::Mat* fBack = (isForward) ? Matrix::MatInverse(func) : Matrix::MatCopy(func);
00042
00043 int sW = ArrayCW(src);
00044 int sH = ArrayCH(src);
00045 int reqW = (dst) ? ArrayCW(dst) : sW;
00046 int reqH = (dst) ? ArrayCH(dst) : sH;
00047
00048 Vec3Real64 translate(0, 0, 0);
00049
00050 if (adjustSize)
00051 {
00052 Vec3Real64 ulP = Matrix::MatMulVec(fForw, Vec3Real64(0, 0, 1));
00053 Vec3Real64 llP = Matrix::MatMulVec(fForw, Vec3Real64(0, sH, 1));
00054 Vec3Real64 urP = Matrix::MatMulVec(fForw, Vec3Real64(sW, 0, 1));
00055 Vec3Real64 lrP = Matrix::MatMulVec(fForw, Vec3Real64(sW, sH, 1));
00056
00057
00058 ulP /= Vec3Real64(ulP.Z(), ulP.Z(), ulP.Z());
00059 llP /= Vec3Real64(llP.Z(), ulP.Z(), ulP.Z());
00060 urP /= Vec3Real64(urP.Z(), ulP.Z(), ulP.Z());
00061 lrP /= Vec3Real64(lrP.Z(), ulP.Z(), ulP.Z());
00062
00063 Vec3Real64 maerB = E2Inf(E2Inf(E2Inf(ulP,llP),urP),lrP);
00064 Vec3Real64 maerE = E2Sup(E2Sup(E2Sup(ulP,llP),urP),lrP);
00065 reqW = int(maerE.X() - maerB.X() + 0.5);
00066 reqH = int(maerE.Y() - maerB.Y() + 0.5);
00067 translate = Vec3Real64(maerB.X(), maerB.Y(), 0);
00068 }
00069
00070 if (dst == 0)
00071 {
00072 dst = ArrayCreate<DstArrayT>(reqW, reqH,
00073 ArrayBW(src), ArrayBH(src));
00074 }
00075 else
00076 if ((ArrayCW(dst) != reqW) || (ArrayCH(dst) != reqH))
00077 {
00078 delete dst;
00079 dst = ArrayCreate<DstArrayT>(reqW, reqH,
00080 ArrayBW(src), ArrayBH(src));
00081 }
00082
00083 #ifdef PX_HORUS_USED
00084 if (!PxRunParallel()) {
00085 #endif
00086 FuncGeometricOp(dst, src, fBack, gi, translate, background);
00087
00088 #ifdef PX_HORUS_USED
00089 } else {
00090 PxArrayPreStateTrans(src, PAR_FULL, STRONG);
00091 PxArrayPreStateTrans(dst, PAR_PART, WEAK);
00092 translate += Vec3Real64(PxLclStartX(PxMyCPU()),
00093 PxLclStartY(PxMyCPU()), PxLclStartZ(PxMyCPU()));
00094 FuncGeometricOp(PxArrayPD(dst), PxArrayPD(src),
00095 fBack, gi, translate, background);
00096 PxArrayPostStateTrans(dst);
00097
00098 if (!PxRunLazyParallel()) {
00099 PxArrayForceNonDistributed(src);
00100 PxArrayForceNonDistributed(dst);
00101 }
00102 }
00103 #endif
00104
00105 delete fForw;
00106 delete fBack;
00107 }
00108
00109 }
00110 }
00111 }
00112 }
00113
00114 #endif