From af50ae73cfe6f1533e5390e5131e25c6d322dc8a Mon Sep 17 00:00:00 2001 From: cam Date: Wed, 4 Jan 2017 23:09:30 -0700 Subject: removed .cpp and .h files --- video_input/consistencyChecker/CFilter.h | 2210 ------------------------------ 1 file changed, 2210 deletions(-) delete mode 100644 video_input/consistencyChecker/CFilter.h (limited to 'video_input/consistencyChecker/CFilter.h') diff --git a/video_input/consistencyChecker/CFilter.h b/video_input/consistencyChecker/CFilter.h deleted file mode 100644 index 29bef9e..0000000 --- a/video_input/consistencyChecker/CFilter.h +++ /dev/null @@ -1,2210 +0,0 @@ -// - Classes for 1D and 2D convolution stencils -// - Pre-defined convolution stencils for binomial filters -// - Pre-defined convolution stencils for 1st, 2nd, 3rd and 4th derivatives up to order 10 -// - Functions for convolution -// -// Author: Thomas Brox - -#ifndef CFILTER -#define CFILTER - -#include -#include -#include -#include -#include -#include - -// CFilter is an extention of CVector. It has an additional property Delta -// which shifts the data to the left (a vector always begins with index 0). -// This enables a filter's range to go from A to B where A can also -// be less than zero. -// -// Example: -// CFilter filter(3,1); -// filter = 1.0; -// cout << filter(-1) << ", " << filter(0) << ", " << filter(1) << endl; -// -// CFilter2D behaves the same way as CFilter but is an extension of CMatrix - -template -class CFilter : public CVector { -public: - // constructor - inline CFilter(const int aSize, const int aDelta = 0); - // copy constructor - CFilter(const CFilter& aCopyFrom); - // constructor initialized by a vector - CFilter(const CVector& aCopyFrom, const int aDelta = 0); - - // Access to the filter's values - inline T& operator()(const int aIndex) const; - inline T& operator[](const int aIndex) const; - // Copies a filter into this filter - CFilter& operator=(const CFilter& aCopyFrom); - - // Access to the filter's delta - inline int delta() const; - // Access to the filter's range A<=i -class CFilter2D : public CMatrix { -public: - // constructor - inline CFilter2D(); - inline CFilter2D(const int aXSize, const int aYSize, const int aXDelta = 0, const int aYDelta = 0); - // copy contructor - CFilter2D(const CFilter2D& aCopyFrom); - // constructor initialized by a matrix - CFilter2D(const CMatrix& aCopyFrom, const int aXDelta = 0, const int aYDelta = 0); - // Normalize sum of values to 1.0 - void normalizeSum(); - // Moves the filter's center - void shift(int aXDelta, int aYDelta); - - // Access to filter's values - inline T& operator()(const int ax, const int ay) const; - // Copies a filter into this filter - CFilter2D& operator=(const CFilter2D& aCopyFrom); - - // Access to the filter's delta - inline int deltaX() const; - inline int deltaY() const; - // Access to the filter's range A<=i inline void filter(CVector& aVector, const CFilter& aFilter); - // Convolution of the vector aVector with aFilter, the initial values of aVector will persist. - template void filter(const CVector& aVector, CVector& aResult, const CFilter& aFilter); - - // Convolution with a rectangle -> approximation of Gaussian - template inline void boxFilter(CVector& aVector, int aWidth); - template void boxFilter(const CVector& aVector, CVector& aResult, int aWidth); - - // Linear 2D filtering - - // Convolution of the matrix aMatrix with aFilter, aFilter should be a separable filter - // The result will be written into aMatrix, so its initial values will get lost - template inline void filter(CMatrix& aMatrix, const CFilter& aFilterX, const CFilter& aFilterY); - template inline void filterMin(CMatrix& aMatrix, const CFilter& aFilterX, const CFilter& aFilterY); - // Convolution of the matrix aMatrix with aFilter, aFilter must be separable - // The initial values of aMatrix will persist. - template inline void filter(const CMatrix& aMatrix, CMatrix& aResult, const CFilter& aFilterX, const CFilter& aFilterY); - template inline void filterMin(const CMatrix& aMatrix, CMatrix& aResult, const CFilter& aFilterX, const CFilter& aFilterY); - - // Convolution of the matrix aMatrix with aFilter only in x-direction, aDummy can be set to 1 - // The result will be written into aMatrix, so its initial values will get lost - template inline void filter(CMatrix& aMatrix, const CFilter& aFilter, const int aDummy); - template inline void filterMin(CMatrix& aMatrix, const CFilter& aFilter, const int aDummy); - // Convolution of the matrix aMatrix with aFilter only in x-direction, aDummy can be set to 1 - // The initial values of aMatrix will persist. - template void filter(const CMatrix& aMatrix, CMatrix& aResult, const CFilter& aFilter, const int aDummy); - template void filterMin(const CMatrix& aMatrix, const CMatrix& aOrig, CMatrix& aResult, const CFilter& aFilter, const int aDummy); - // Convolution of the matrix aMatrix with aFilter only in y-direction, aDummy can be set to 1 - // The result will be written into aMatrix, so its initial values will get lost - template inline void filter(CMatrix& aMatrix, const int aDummy, const CFilter& aFilter); - template inline void filterMin(CMatrix& aMatrix, const int aDummy, const CFilter& aFilter); - // Convolution of the matrix aMatrix with aFilter only in y-direction, aDummy can be set to 1 - // The initial values of aMatrix will persist. - template void filter(const CMatrix& aMatrix, CMatrix& aResult, const int aDummy, const CFilter& aFilter); - template void filterMin(const CMatrix& aMatrix, const CMatrix& aOrig, CMatrix& aResult, const int aDummy, const CFilter& aFilter); - - // Convolution of the matrix aMatrix with aFilter - // The result will be written to aMatrix, so its initial values will get lost - template inline void filter(CMatrix& aMatrix, const CFilter2D& aFilter); - // Convolution of the matrix aMatrix with aFilter, the initial values of aMatrix will persist - template void filter(const CMatrix& aMatrix, CMatrix& aResult, const CFilter2D& aFilter); - - // Convolution with a rectangle -> approximation of Gaussian - template inline void boxFilterX(CMatrix& aMatrix, int aWidth); - template void boxFilterX(const CMatrix& aMatrix, CMatrix& aResult, int aWidth); - template inline void boxFilterY(CMatrix& aMatrix, int aWidth); - template void boxFilterY(const CMatrix& aMatrix, CMatrix& aResult, int aWidth); - - // Recursive filter -> approximation of Gaussian - template void recursiveSmoothX(CMatrix& aMatrix, float aSigma); - template void recursiveSmoothY(CMatrix& aMatrix, float aSigma); - template inline void recursiveSmooth(CMatrix& aMatrix, float aSigma); - - // Linear 3D filtering - - // Convolution of the 3D Tensor aTensor with aFilter, aFilter must be separable - // The result will be written back to aTensor so its initial values will get lost - template inline void filter(CTensor& aTensor, const CFilter& aFilterX, const CFilter& aFilterY, const CFilter& aFilterZ); - // Convolution of the 3D Tensor aTensor with aFilter, aFilter must be separable - // The initial values of aTensor will persist - template inline void filter(const CTensor& aTensor, CTensor& aResult, const CFilter& aFilterX, const CFilter& aFilterY, const CFilter& aFilterZ); - - // Convolution of the 3D Tensor aTensor with aFilter only in x-Direction - template inline void filter(CTensor& aTensor, const CFilter& aFilter, const int aDummy1, const int aDummy2); - template void filter(const CTensor& aTensor, CTensor& aResult, const CFilter& aFilter, const int aDummy1, const int aDummy2); - // Convolution of the 3D Tensor aTensor with aFilter only in y-Direction - template inline void filter(CTensor& aTensor, const int aDummy1, const CFilter& aFilter, const int aDummy2); - template void filter(const CTensor& aTensor, CTensor& aResult, const int aDummy1, const CFilter& aFilter, const int aDummy2); - // Convolution of the 3D Tensor aTensor with aFilter only in z-Direction - template inline void filter(CTensor& aTensor, const int aDummy1, const int aDummy2, const CFilter& aFilter); - template void filter(const CTensor& aTensor, CTensor& aResult, const int aDummy1, const int aDummy2, const CFilter& aFilter); - - // Convolution with a rectangle -> approximation of Gaussian - template inline void boxFilterX(CTensor& aTensor, int aWidth); - template void boxFilterX(const CTensor& aTensor, CTensor& aResult, int aWidth); - template inline void boxFilterY(CTensor& aTensor, int aWidth); - template void boxFilterY(const CTensor& aTensor, CTensor& aResult, int aWidth); - template inline void boxFilterZ(CTensor& aTensor, int aWidth); - template void boxFilterZ(const CTensor& aTensor, CTensor& aResult, int aWidth); - - // Recursive filter -> approximation of Gaussian - template void recursiveSmoothX(CTensor& aTensor, float aSigma); - template void recursiveSmoothY(CTensor& aTensor, float aSigma); - template void recursiveSmoothZ(CTensor& aTensor, float aSigma); - - // Linear 4D filtering - - // Convolution of the 4D Tensor aTensor with aFilter, aFilter must be separable - // The result will be written back to aTensor so its initial values will get lost - template inline void filter(CTensor4D& aTensor, const CFilter& aFilterX, const CFilter& aFilterY, const CFilter& aFilterZ, const CFilter& aFilterA); - - // Convolution of the 4D Tensor aTensor with aFilter only in x-Direction - template inline void filter(CTensor4D& aTensor, const CFilter& aFilter, const int aDummy1, const int aDummy2, const int aDummy3); - template void filter(const CTensor4D& aTensor, CTensor4D& aResult, const CFilter& aFilter, const int aDummy1, const int aDummy2, const int aDummy3); - // Convolution of the 4D Tensor aTensor with aFilter only in y-Direction - template inline void filter(CTensor4D& aTensor, const int aDummy1, const CFilter& aFilter, const int aDummy2, const int aDummy3); - template void filter(const CTensor4D& aTensor, CTensor4D& aResult, const int aDummy1, const CFilter& aFilter, const int aDummy2, const int aDummy3); - // Convolution of the 4D Tensor aTensor with aFilter only in z-Direction - template inline void filter(CTensor4D& aTensor, const int aDummy1, const int aDummy2, const CFilter& aFilter, const int aDummy3); - template void filter(const CTensor4D& aTensor, CTensor4D& aResult, const int aDummy1, const int aDummy2, const CFilter& aFilter, const int aDummy3); - // Convolution of the 4D Tensor aTensor with aFilter only in a-Direction - template inline void filter(CTensor4D& aTensor, const int aDummy1, const int aDummy2, const int aDummy3, const CFilter& aFilter); - template void filter(const CTensor4D& aTensor, CTensor4D& aResult, const int aDummy1, const int aDummy2, const int aDummy3, const CFilter& aFilter); - - // Recursive filter -> approximation of Gaussian - template void recursiveSmoothX(CTensor4D& aTensor, float aSigma); - template void recursiveSmoothY(CTensor4D& aTensor, float aSigma); - template void recursiveSmoothZ(CTensor4D& aTensor, float aSigma); - template void recursiveSmoothA(CTensor4D& aTensor, float aSigma); - - // Nonlinear filtering: Osher shock filter - template void osher(CMatrix& aData, int aIterations = 20); - template inline void osher(const CMatrix& aData, CMatrix& aResult, int aIterations = 20); -} - -// Common filters - -template -class CGauss : public CFilter { -public: - CGauss(const int aSize, const int aDegreeOfDerivative); -}; - -template -class CSmooth : public CFilter { -public: - CSmooth(float aSigma, float aPrecision); -}; - -template -class CGaussianFirstDerivative : public CFilter { -public: - CGaussianFirstDerivative(float aSigma, float aPrecision); -}; - -template -class CGaussianSecondDerivative : public CFilter { -public: - CGaussianSecondDerivative(float aSigma, float aPrecision); -}; - -template -class CDerivative : public CFilter { -public: - CDerivative(const int aSize); -}; - -template -class CHighOrderDerivative : public CFilter { -public: - CHighOrderDerivative(int aOrder, int aSize); -}; - -template -class CGaborReal : public CFilter2D { -public: - CGaborReal(float aFrequency, float aAngle, float aSigma1 = 3.0, float aSigma2 = 3.0); -}; - -template -class CGaborImaginary : public CFilter2D { -public: - CGaborImaginary(float aFrequency, float aAngle, float aSigma1 = 3.0, float aSigma2 = 3.0); -}; - -// Exceptions ----------------------------------------------------------------- - -// Thrown if one tries to access an element of a filter which is out of the filter's bounds -struct EFilterRangeOverflow { - EFilterRangeOverflow(const int aIndex, const int aA, const int aB) { - using namespace std; - cerr << "Exception EFilterRangeOverflow: i = " << aIndex; - cerr << " Allowed Range: " << aA << " <= i < " << aB << endl; - } - EFilterRangeOverflow(const int ax, const int ay, const int aAX, const int aBX, const int aAY, const int aBY) { - using namespace std; - cerr << "Exception EFilterRangeOverflow: (x,y) = (" << ax << "," << ay << ") "; - cerr << "Allowed Range: " << aAX << " <= x < " << aBX << " " << aAY << " <= y < " << aBY << endl; - } -}; - -// Thrown if the resulting container has not the same size as the initial container -struct EFilterIncompatibleSize { - EFilterIncompatibleSize(const int aSize1, const int aSize2) { - using namespace std; - cerr << "Exception EFilterIncompatibleSize: Initial container size: " << aSize1; - cerr << " Resulting container size: " << aSize2 << endl; - } -}; - -// Thrown if the demanded filter is not available -struct EFilterNotAvailable { - EFilterNotAvailable(int aSize, int aOrder) { - using namespace std; - cerr << "Exception EFilterNotAvailable: Mask size: " << aSize; - if (aOrder >= 0) cerr << " Derivative order: " << aOrder; - cerr << endl; - } -}; - -// I M P L E M E N T A T I O N ------------------------------------------------ -// -// You might wonder why there is implementation code in a header file. -// The reason is that not all C++ compilers yet manage separate compilation -// of templates. Inline functions cannot be compiled separately anyway. -// So in this case the whole implementation code is added to the header -// file. -// Users should ignore everything that's beyond this line :) -// ---------------------------------------------------------------------------- - -// C F I L T E R -------------------------------------------------------------- -// P U B L I C ---------------------------------------------------------------- -// constructor -template -inline CFilter::CFilter(const int aSize, const int aDelta) - : CVector(aSize),mDelta(aDelta) { -} - -// copy constructor -template -CFilter::CFilter(const CFilter& aCopyFrom) - : CVector(aCopyFrom.mSize),mDelta(aCopyFrom.mDelta) { - for (register int i = 0; i < this->mSize; i++) - this->mData[i] = aCopyFrom.mData[i]; -} - -// constructor initialized by a vector -template -CFilter::CFilter(const CVector& aCopyFrom, const int aDelta) - : CVector(aCopyFrom.size()),mDelta(aDelta) { - for (register int i = 0; i < this->mSize; i++) - this->mData[i] = aCopyFrom(i); -} - -// operator() -template -inline T& CFilter::operator()(const int aIndex) const { - #ifdef DEBUG - if (aIndex < A() || aIndex >= B()) - throw EFilterRangeOverflow(aIndex,A(),B()); - #endif - return this->mData[aIndex+mDelta]; -} - -// operator[] -template -inline T& CFilter::operator[](const int aIndex) const { - return operator()(aIndex); -} - -// operator= -template -CFilter& CFilter::operator=(const CFilter& aCopyFrom) { - if (this != &aCopyFrom) { - delete[] this->mData; - this->mSize = aCopyFrom.mSize; - mDelta = aCopyFrom.mDelta; - this->mData = new T[this->mSize]; - for (register int i = 0; i < this->mSize; i++) - this->mData[i] = aCopyFrom.mData[i]; - } - return *this; -} - -// delta -template -inline int CFilter::delta() const { - return mDelta; -} - -// A -template -inline int CFilter::A() const { - return -mDelta; -} - -// B -template -inline int CFilter::B() const { - return this->mSize-mDelta; -} - -// sum -template -T CFilter::sum() const { - T aResult = 0; - for (int i = 0; i < this->mSize; i++) - aResult += fabs(this->mData[i]); - return aResult; -} - -// shift -template -inline void CFilter::shift(int aDelta) { - mDelta += aDelta; -} - -// C F I L T E R 2 D ----------------------------------------------------------- -// P U B L I C ---------------------------------------------------------------- -// constructor -template -inline CFilter2D::CFilter2D() - : CMatrix(),mDeltaX(0),mDeltaY(0) { -} - -template -inline CFilter2D::CFilter2D(const int aXSize, const int aYSize, const int aDeltaX, const int aDeltaY) - : CMatrix(aXSize,aYSize),mDeltaX(aDeltaX),mDeltaY(aDeltaY) { -} - -// copy constructor -template -CFilter2D::CFilter2D(const CFilter2D& aCopyFrom) - : CMatrix(aCopyFrom.mXSize,aCopyFrom.mYSize),mDeltaX(aCopyFrom.mDeltaX,aCopyFrom.mDeltaY) { - for (int i = 0; i < this->mXSize*this->mYSize; i++) - this->mData[i] = aCopyFrom.mData[i]; -} - -// constructor initialized by a matrix -template -CFilter2D::CFilter2D(const CMatrix& aCopyFrom, const int aDeltaX, const int aDeltaY) - : CMatrix(aCopyFrom.xSize(),aCopyFrom.ySize()),mDeltaX(aDeltaX),mDeltaY(aDeltaY) { - for (register int i = 0; i < this->mXSize*this->mYSize; i++) - this->mData[i] = aCopyFrom.data()[i]; -} - -// normalizeSum -template -void CFilter2D::normalizeSum() { - int aSize = this->size(); - T aSum = 0; - for (int i = 0; i < aSize; i++) - aSum += this->mData[i]; - T invSum = 1.0/aSum; - for (int i = 0; i < aSize; i++) - this->mData[i] *= invSum; -} - -// shift -template -void CFilter2D::shift(int aXDelta, int aYDelta) { - mDeltaX = aXDelta; - mDeltaY = aYDelta; -} - -// operator() -template -inline T& CFilter2D::operator()(const int ax, const int ay) const { - #ifdef DEBUG - if (ax < AX() || ax >= BX() || ay < AY() || ay >= BY) - throw EFilterRangeOverflow(ax,ay,AX(),BX(),AY(),BY()); - #endif - return this->mData[(ay+mDeltaY)*this->mXSize+ax+mDeltaX]; -} - -// operator= -template -CFilter2D& CFilter2D::operator=(const CFilter2D& aCopyFrom) { - if (this != &aCopyFrom) { - delete[] this->mData; - this->mXSize = aCopyFrom.mXSize; - this->mYSize = aCopyFrom.mYSize; - mDeltaX = aCopyFrom.mDeltaX; - mDeltaY = aCopyFrom.mDeltaY; - this->mData = new T[this->mXSize*this->mYSize]; - for (register int i = 0; i < this->mXSize*this->mYSize; i++) - this->mData[i] = aCopyFrom.mData[i]; - } - return *this; -} - -// deltaX -template -inline int CFilter2D::deltaX() const { - return mDeltaX; -} - -// deltaY -template -inline int CFilter2D::deltaY() const { - return mDeltaY; -} - -// AX -template -inline int CFilter2D::AX() const { - return -mDeltaX; -} - -// AY -template -inline int CFilter2D::AY() const { - return -mDeltaY; -} - -// BX -template -inline int CFilter2D::BX() const { - return this->mXSize-mDeltaX; -} - -// BY -template -inline int CFilter2D::BY() const { - return this->mYSize-mDeltaY; -} - -// sum -template -T CFilter2D::sum() const { - T aResult = 0; - for (int i = 0; i < this->mXSize*this->mYSize; i++) - aResult += abs(this->mData[i]); - return aResult; -} - -// C G A U S S ----------------------------------------------------------------- -template -CGauss::CGauss(const int aSize, const int aDegreeOfDerivative) - : CFilter(aSize,aSize >> 1) { - CVector *oldData; - CVector *newData; - CVector *temp; - oldData = new CVector(aSize); - newData = new CVector(aSize); - - (*oldData)(0) = 1; - (*oldData)(1) = 1; - - for (int i = 2; i < aSize-aDegreeOfDerivative; i++) { - (*newData)(0) = 1; - for (int j = 1; j < i; j++) - (*newData)(j) = (*oldData)(j)+(*oldData)(j-1); - (*newData)(i) = 1; - temp = oldData; - oldData = newData; - newData = temp; - } - for (int i = aSize-aDegreeOfDerivative; i < aSize; i++) { - (*newData)(0) = 1; - for (int j = 1; j < i; j++) - (*newData)(j) = (*oldData)(j)-(*oldData)(j-1); - (*newData)(i) = -(*oldData)(i-1); - temp = oldData; - oldData = newData; - newData = temp; - } - - int aSum = 0; - for (int i = 0; i < aSize; i++) - aSum += abs((*oldData)(i)); - double aInvSum = 1.0/aSum; - for (int i = 0; i < aSize; i++) - this->mData[aSize-1-i] = (*oldData)(i)*aInvSum; - - delete newData; - delete oldData; -} - -// C S M O O T H --------------------------------------------------------------- -template -CSmooth::CSmooth(float aSigma, float aPrecision) - : CFilter(2*(int)ceil(aPrecision*aSigma)+1,(int)ceil(aPrecision*aSigma)) { - float aSqrSigma = aSigma*aSigma; - for (int i = 0; i <= (this->mSize >> 1); i++) { - T aTemp = exp(i*i/(-2.0*aSqrSigma))/(aSigma*sqrt(2.0*NMath::Pi)); - this->operator()(i) = aTemp; - this->operator()(-i) = aTemp; - } - T invSum = 1.0/this->sum(); - for (int i = 0; i < this->mSize; i++) - this->mData[i] *= invSum; -} - -template -CGaussianFirstDerivative::CGaussianFirstDerivative(float aSigma, float aPrecision) - : CFilter(2*(int)ceil(aPrecision*aSigma)+1,(int)ceil(aPrecision*aSigma)) { - float aSqrSigma = aSigma*aSigma; - float aPreFactor = 1.0/(aSqrSigma*aSigma*sqrt(2.0*NMath::Pi)); - for (int i = 0; i <= (this->mSize >> 1); i++) { - T aTemp = exp(i*i/(-2.0*aSqrSigma))*i*aPreFactor; - this->operator()(i) = aTemp; - this->operator()(-i) = -aTemp; - } -} - -template -CGaussianSecondDerivative::CGaussianSecondDerivative(float aSigma, float aPrecision) - : CFilter(2*(int)ceil(aPrecision*aSigma)+1,(int)ceil(aPrecision*aSigma)) { - float aSqrSigma = aSigma*aSigma; - float aPreFactor = 1.0/(aSqrSigma*aSigma*sqrt(2.0*NMath::Pi)); - for (int i = 0; i <= (this->mSize >> 1); i++) { - T aTemp = exp(i*i/(-2.0*aSqrSigma))*(i*i/aSqrSigma-1.0)*aPreFactor; - this->operator()(i) = aTemp; - this->operator()(-i) = aTemp; - } -} - -// C D E R I V A T I V E ------------------------------------------------------- -template -CDerivative::CDerivative(const int aSize) - : CFilter(aSize,(aSize-1) >> 1) { - switch (aSize) { - case 2: - this->mData[0] = -1; - this->mData[1] = 1; - break; - case 3: - this->mData[0] = -0.5; - this->mData[1] = 0; - this->mData[2] = 0.5; - break; - case 4: - this->mData[0] = 0.041666666666666666666666666666667; - this->mData[1] = -1.125; - this->mData[2] = 1.125; - this->mData[3] = -0.041666666666666666666666666666667; - break; - case 5: - this->mData[0] = 0.083333333333; - this->mData[1] = -0.66666666666; - this->mData[2] = 0; - this->mData[3] = 0.66666666666; - this->mData[4] = -0.083333333333; - break; - case 6: - this->mData[0] = -0.0046875; - this->mData[1] = 0.0651041666666666666666666666666667; - this->mData[2] = -1.171875; - this->mData[3] = 1.171875; - this->mData[4] = -0.0651041666666666666666666666666667; - this->mData[5] = 0.0046875; - break; - case 7: - this->mData[0] = -0.016666666666666666666666666666667; - this->mData[1] = 0.15; - this->mData[2] = -0.75; - this->mData[3] = 0; - this->mData[4] = 0.75; - this->mData[5] = -0.15; - this->mData[6] = 0.016666666666666666666666666666667; - break; - case 8: - this->mData[0] = 6.9754464285714285714285714285714e-4; - this->mData[1] = -0.0095703125; - this->mData[2] = 0.079752604166666666666666666666667; - this->mData[3] = -1.1962890625; - this->mData[4] = 1.1962890625; - this->mData[5] = -0.079752604166666666666666666666667; - this->mData[6] = 0.0095703125; - this->mData[7] = -6.9754464285714285714285714285714e-4; - break; - case 9: - this->mData[0] = 0.0035714285714285714285714285714286; - this->mData[1] = -0.038095238095238095238095238095238; - this->mData[2] = 0.2; - this->mData[3] = -0.8; - this->mData[4] = 0; - this->mData[5] = 0.8; - this->mData[6] = -0.2; - this->mData[7] = 0.038095238095238095238095238095238; - - this->mData[8] = -0.0035714285714285714285714285714286; - break; - case 10: - this->mData[0] = -1.1867947048611111111111111111111e-4; - this->mData[1] = 0.0017656598772321428571428571428571; - this->mData[2] = -0.0138427734375; - this->mData[3] = 0.0897216796875; - this->mData[4] = -1.21124267578125; - this->mData[5] = 1.21124267578125; - this->mData[6] = -0.0897216796875; - this->mData[7] = 0.0138427734375; - this->mData[8] = -0.0017656598772321428571428571428571; - this->mData[9] = 1.1867947048611111111111111111111e-4; - break; - default: - throw EFilterNotAvailable(aSize,-1); - } -} - -// C H I G H O R D E R D E R I V A T I V E ------------------------------------- -template -CHighOrderDerivative::CHighOrderDerivative(int aOrder, int aSize) - : CFilter(aSize,(aSize-1) >> 1) { - switch (aSize) { - case 3: - switch (aOrder) { - case 2: - this->mData[0] = 1; - this->mData[1] = -2; - this->mData[2] = 1; - break; - default: - throw EFilterNotAvailable(aSize,aOrder); - } - break; - case 4: - switch (aOrder) { - case 2: - this->mData[0] = 0.25; - this->mData[1] = -0.25; - this->mData[2] = -0.25; - this->mData[3] = 0.25; - break; - case 3: - this->mData[0] = -0.25; - this->mData[1] = 0.75; - this->mData[2] = -0.75; - this->mData[3] = 0.25; - break; - default: - throw EFilterNotAvailable(aSize,aOrder); - } - break; - case 5: - switch (aOrder) { - case 2: - this->mData[0] = -0.083333333333333333333333333333333; - this->mData[1] = 1.3333333333333333333333333333333; - this->mData[2] = -2.5; - this->mData[3] = 1.3333333333333333333333333333333; - this->mData[4] = -0.083333333333333333333333333333333; - break; - case 3: - this->mData[0] = -0.5; - this->mData[1] = 1; - this->mData[2] = 0; - this->mData[3] = -1; - this->mData[4] = 0.5; - break; - case 4: - this->mData[0] = 1; - this->mData[1] = -4; - this->mData[2] = 6; - this->mData[3] = -4; - this->mData[4] = 1; - break; - default: - throw EFilterNotAvailable(aSize,aOrder); - } - break; - case 6: - switch (aOrder) { - case 2: - this->mData[0] = -0.052083333333333333333333333333333; - this->mData[1] = 0.40625; - this->mData[2] = -0.35416666666666666666666666666667; - this->mData[3] = -0.35416666666666666666666666666667; - this->mData[4] = 0.40625; - this->mData[5] = -0.052083333333333333333333333333333; - break; - case 3: - this->mData[0] = 0.03125; - this->mData[1] = -0.40625; - this->mData[2] = 1.0625; - this->mData[3] = -1.0625; - this->mData[4] = 0.40625; - this->mData[5] = -0.03125; - break; - case 4: - this->mData[0] = 0.0625; - this->mData[1] = -0.1875; - this->mData[2] = 0.125; - this->mData[3] = 0.125; - this->mData[4] = -0.1875; - this->mData[5] = 0.0625; - break; - default: - throw EFilterNotAvailable(aSize,aOrder); - } - break; - case 7: - switch (aOrder) { - case 2: - this->mData[0] = 0.011111111111111111111111111111111; - this->mData[1] = -0.15; - this->mData[2] = 1.5; - this->mData[3] = -2.6666666666666666666666666666667; - this->mData[4] = 1.5; - this->mData[5] = -0.15; - this->mData[6] = 0.011111111111111111111111111111111; - break; - case 3: - this->mData[0] = 0.125; - this->mData[1] = -1; - this->mData[2] = 1.625; - this->mData[3] = 0; - this->mData[4] = -1.625; - this->mData[5] = 1; - this->mData[6] = -0.125; - break; - case 4: - this->mData[0] = -0.16666666666666666666666666666667; - this->mData[1] = 2; - this->mData[2] = -6.5; - this->mData[3] = 9.3333333333333333333333333333333; - this->mData[4] = -6.5; - this->mData[5] = 2; - this->mData[6] = -0.16666666666666666666666666666667; - break; - default: - throw EFilterNotAvailable(aSize,aOrder); - } - break; - case 8: - switch (aOrder) { - case 2: - this->mData[0] = 0.011241319444444444444444444444444; - this->mData[1] = -0.10828993055555555555555555555556; - this->mData[2] = 0.507421875; - this->mData[3] = -0.41037326388888888888888888888889; - this->mData[4] = -0.41037326388888888888888888888889; - this->mData[5] = 0.507421875; - this->mData[6] = -0.10828993055555555555555555555556; - this->mData[7] = 0.011241319444444444444444444444444; - break; - case 3: - this->mData[0] = -0.0048177083333333333333333333333333; - this->mData[1] = 0.064973958333333333333333333333333; - this->mData[2] = -0.507421875; - this->mData[3] = 1.2311197916666666666666666666667; - this->mData[4] = -1.2311197916666666666666666666667; - this->mData[5] = 0.507421875; - this->mData[6] = -0.064973958333333333333333333333333; - this->mData[7] = 0.0048177083333333333333333333333333; - break; - case 4: - this->mData[0] = -0.018229166666666666666666666666667; - this->mData[1] = 0.15364583333333333333333333333333; - this->mData[2] = -0.3515625; - this->mData[3] = 0.21614583333333333333333333333333; - this->mData[4] = 0.21614583333333333333333333333333; - this->mData[5] = -0.3515625; - this->mData[6] = 0.15364583333333333333333333333333; - this->mData[7] = -0.018229166666666666666666666666667; - break; - default: - throw EFilterNotAvailable(aSize,aOrder); - } - break; - case 9: - switch (aOrder) { - case 2: - this->mData[0] = -0.0017857142857142857142857142857143; - this->mData[1] = 0.025396825396825396825396825396825; - this->mData[2] = -0.2; - this->mData[3] = 1.6; - this->mData[4] = -2.8472222222222222222222222222222; - this->mData[5] = 1.6; - this->mData[6] = -0.2; - this->mData[7] = 0.025396825396825396825396825396825; - this->mData[8] = -0.0017857142857142857142857142857143; - break; - case 3: - this->mData[0] = -0.029166666666666666666666666666667; - this->mData[1] = 0.3; - this->mData[2] = -1.4083333333333333333333333333333; - this->mData[3] = 2.0333333333333333333333333333333; - this->mData[4] = 0; - this->mData[5] = -2.0333333333333333333333333333333; - this->mData[6] = 1.4083333333333333333333333333333; - this->mData[7] = -0.3; - this->mData[8] = 0.029166666666666666666666666666667; - break; - case 4: - this->mData[0] = 0.029166666666666666666666666666667; - this->mData[1] = -0.4; - this->mData[2] = 2.8166666666666666666666666666667; - this->mData[3] = -8.1333333333333333333333333333333; - this->mData[4] = 11.375; - this->mData[5] = -8.1333333333333333333333333333333; - this->mData[6] = 2.8166666666666666666666666666667; - this->mData[7] = -0.4; - this->mData[8] = 0.029166666666666666666666666666667; - break; - default: - throw EFilterNotAvailable(aSize,aOrder); - } - break; - case 10: - switch (aOrder) { - case 2: - this->mData[0] = -0.0025026351686507936507936507936508; - this->mData[1] = 0.028759765625; - this->mData[2] = -0.15834263392857142857142857142857; - this->mData[3] = 0.57749565972222222222222222222222; - this->mData[4] = -0.44541015625; - this->mData[5] = -0.44541015625; - this->mData[6] = 0.57749565972222222222222222222222; - this->mData[7] = -0.15834263392857142857142857142857; - this->mData[8] = 0.028759765625; - this->mData[9] = -0.0025026351686507936507936507936508; - break; - case 3: - this->mData[0] = 0.0008342117228835978835978835978836; - this->mData[1] = -0.012325613839285714285714285714286; - this->mData[2] = 0.095005580357142857142857142857143; - this->mData[3] = -0.57749565972222222222222222222222; - this->mData[4] = 1.33623046875; - this->mData[5] = -1.33623046875; - this->mData[6] = 0.57749565972222222222222222222222; - this->mData[7] = -0.095005580357142857142857142857143; - this->mData[8] = 0.012325613839285714285714285714286; - this->mData[9] = -0.0008342117228835978835978835978836; - break; - case 4: - this->mData[0] = 0.00458984375; - this->mData[1] = -0.050358072916666666666666666666667; - this->mData[2] = 0.24544270833333333333333333333333; - this->mData[3] = -0.480078125; - this->mData[4] = 0.28040364583333333333333333333333; - this->mData[5] = 0.28040364583333333333333333333333; - this->mData[6] = -0.480078125; - this->mData[7] = 0.24544270833333333333333333333333; - this->mData[8] = -0.050358072916666666666666666666667; - this->mData[9] = 0.00458984375; - break; - default: - throw EFilterNotAvailable(aSize,aOrder); - } - break; - default: - throw EFilterNotAvailable(aSize,aOrder); - } -} - -// C G A B O R ----------------------------------------------------------------- -template -CGaborReal::CGaborReal(float aFrequency, float aAngle, float aSigma1, float aSigma2) - : CFilter2D() { - // sqrt(2.0*log(2.0))/(2.0*NMath::Pi) = 0.18739 - float sigma1Sqr2 = aSigma1*0.18739/aFrequency; - sigma1Sqr2 = 0.5/(sigma1Sqr2*sigma1Sqr2); - float sigma2Sqr2 = aSigma2*0.18739/aFrequency; - sigma2Sqr2 = 0.5/(sigma2Sqr2*sigma2Sqr2); - float aCos = cos(aAngle); - float aSin = sin(aAngle); - float a = 0.6*aSigma1/aFrequency; - float b = 0.6*aSigma2/aFrequency; - float aXSize = fabs(a*aCos)+fabs(b*aSin); - float aYSize = fabs(b*aCos)+fabs(a*aSin); - this->setSize(1+2.0*floor(aXSize),1+2.0*floor(aYSize)); - this->shift(floor(aXSize),floor(aYSize)); - for (int y = this->AY(); y < this->BY(); y++) - for (int x = this->AX(); x < this->BX(); x++) { - float a = x*aCos+y*aSin; - float b = y*aCos-x*aSin; - float aGauss = exp(-sigma1Sqr2*a*a-sigma2Sqr2*b*b); - float aHelp = 2.0*NMath::Pi*aFrequency*(x*aCos+y*aSin); - this->operator()(x,y) = aGauss*cos(aHelp); - } -} - -template -CGaborImaginary::CGaborImaginary(float aFrequency, float aAngle, float aSigma1, float aSigma2) - : CFilter2D() { - // sqrt(2.0*log(2.0))/(2.0*NMath::Pi) = 0.18739 - float sigma1Sqr2 = aSigma1*0.18739/aFrequency; - sigma1Sqr2 = 0.5/(sigma1Sqr2*sigma1Sqr2); - float sigma2Sqr2 = aSigma2*0.18739/aFrequency; - sigma2Sqr2 = 0.5/(sigma2Sqr2*sigma2Sqr2); - float aCos = cos(aAngle); - float aSin = sin(aAngle); - float a = 0.6*aSigma1/aFrequency; - float b = 0.6*aSigma2/aFrequency; - float aXSize = fabs(a*aCos)+fabs(b*aSin); - float aYSize = fabs(b*aCos)+fabs(a*aSin); - this->setSize(1+2.0*floor(aXSize),1+2.0*floor(aYSize)); - this->shift(floor(aXSize),floor(aYSize)); - for (int y = this->AY(); y < this->BY(); y++) - for (int x = this->AX(); x < this->BX(); x++) { - float a = x*aCos+y*aSin; - float b = y*aCos-x*aSin; - float aGauss = exp(-sigma1Sqr2*a*a-sigma2Sqr2*b*b); - float aHelp = 2.0*NMath::Pi*aFrequency*(x*aCos+y*aSin); - this->operator()(x,y) = aGauss*sin(aHelp); - } -} - -// F I L T E R ----------------------------------------------------------------- - -namespace NFilter { - -// 1D linear filtering --------------------------------------------------------- - -template -inline void filter(CVector& aVector, const CFilter& aFilter) { - CVector oldVector(aVector); - filter(oldVector,aVector,aFilter); -} - -template -void filter(const CVector& aVector, CVector& aResult, const CFilter& aFilter) { - if (aResult.size() != aVector.size()) throw EFilterIncompatibleSize(aVector.size(),aResult.size()); - int x1 = -aFilter.A(); - int x2 = aVector.size()-aFilter.B(); - int a2Size = 2*aVector.size()-1; - // Left rim - for (int i = 0; i < x1; i++) { - aResult[i] = 0; - for (int j = aFilter.A(); j < aFilter.B(); j++) - if (j+i < 0) aResult(i) += aFilter(j)*aVector(-1-j-i); - else aResult(i) += aFilter(j)*aVector(j+i); - } - // Middle - for (int i = x1; i < x2; i++) { - aResult[i] = 0; - for (int j = aFilter.A(); j < aFilter.B(); j++) - aResult(i) += aFilter(j)*aVector(j+i); - } - // Right rim - for (int i = x2; i < aResult.size(); i++) { - aResult[i] = 0; - for (int j = aFilter.A(); j < aFilter.B(); j++) - if (j+i >= aVector.size()) aResult(i) += aFilter(j)*aVector(a2Size-j-i); - else aResult(i) += aFilter(j)*aVector(j+i); - } -} - -// boxfilter -template -inline void boxFilter(CVector& aVector, int aWidth) { - CVector aTemp(aVector); - boxFilter(aTemp,aVector,aWidth); -} - -template -void boxFilter(const CVector& aVector, CVector& aResult, int aWidth) { - if (aWidth % 2 == 0) aWidth += 1; - T* invWidth = new T[aWidth+1]; - invWidth[0] = 1.0f; - for (int i = 1; i <= aWidth; i++) - invWidth[i] = 1.0/i; - int halfWidth = (aWidth >> 1); - int aRight = halfWidth; - if (aRight >= aVector.size()) aRight = aVector.size()-1; - // Initialize - T aSum = 0.0f; - for (int x = 0; x <= aRight; x++) - aSum += aVector(x); - int aNum = aRight+1; - // Shift - for (int x = 0; x < aVector.size(); x++) { - aResult(x) = aSum*invWidth[aNum]; - if (x-halfWidth >= 0) { - aSum -= aVector(x-halfWidth); aNum--; - } - if (x+halfWidth+1 < aVector.size()) { - aSum += aVector(x+halfWidth+1); aNum++; - } - } - delete[] invWidth; -} - -// 2D linear filtering --------------------------------------------------------- - -template -inline void filter(CMatrix& aMatrix, const CFilter& aFilterX, const CFilter& aFilterY) { - CMatrix tempMatrix(aMatrix.xSize(),aMatrix.ySize()); - filter(aMatrix,tempMatrix,aFilterX,1); - filter(tempMatrix,aMatrix,1,aFilterY); -} - -template -inline void filter(const CMatrix& aMatrix, CMatrix& aResult, const CFilter& aFilterX, const CFilter& aFilterY) { - CMatrix tempMatrix(aMatrix.xSize(),aMatrix.ySize()); - filter(aMatrix,tempMatrix,aFilterX,1); - filter(tempMatrix,aResult,1,aFilterY); -} - -template -inline void filter(CMatrix& aMatrix, const CFilter& aFilter, const int aDummy) { - CMatrix tempMatrix(aMatrix.xSize(),aMatrix.ySize()); - filter(aMatrix,tempMatrix,aFilter,1); - aMatrix = tempMatrix; -} - -template -void filter(const CMatrix& aMatrix, CMatrix& aResult, const CFilter& aFilter, const int aDummy) { - if (aResult.xSize() != aMatrix.xSize() || aResult.ySize() != aMatrix.ySize()) - throw EFilterIncompatibleSize(aMatrix.xSize()*aMatrix.ySize(),aResult.xSize()*aResult.ySize()); - int x1 = -aFilter.A(); - int x2 = aMatrix.xSize()-aFilter.B(); - int a2Size = 2*aMatrix.xSize()-1; - aResult = 0; - for (int y = 0; y < aMatrix.ySize(); y++) { - int aOffset = y*aMatrix.xSize(); - // Left rim - for (int x = 0; x < x1; x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) { - if (x+i < 0) aResult.data()[aOffset+x] += aFilter[i]*aMatrix.data()[aOffset-1-x-i]; - else if (x+i >= aMatrix.xSize()) aResult.data()[aOffset+x] += aFilter[i]*aMatrix.data()[aOffset+a2Size-x-i]; - else aResult.data()[aOffset+x] += aFilter[i]*aMatrix.data()[aOffset+x+i]; - } - // Center - for (int x = x1; x < x2; x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) - aResult.data()[aOffset+x] += aFilter[i]*aMatrix.data()[aOffset+x+i]; - // Right rim - for (int x = x2; x < aMatrix.xSize(); x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) { - if (x+i < 0) aResult.data()[aOffset+x] += aFilter[i]*aMatrix.data()[aOffset-1-x-i]; - else if (x+i >= aMatrix.xSize()) aResult.data()[aOffset+x] += aFilter[i]*aMatrix.data()[aOffset+a2Size-x-i]; - else aResult.data()[aOffset+x] += aFilter[i]*aMatrix.data()[aOffset+x+i]; - } - } -} - -template -inline void filter(CMatrix& aMatrix, const int aDummy, const CFilter& aFilter) { - CMatrix tempMatrix(aMatrix.xSize(),aMatrix.ySize()); - filter(aMatrix,tempMatrix,1,aFilter); - aMatrix = tempMatrix; -} - -template -void filter(const CMatrix& aMatrix, CMatrix& aResult, const int aDummy, const CFilter& aFilter) { - if (aResult.xSize() != aMatrix.xSize() || aResult.ySize() != aMatrix.ySize()) - throw EFilterIncompatibleSize(aMatrix.xSize()*aMatrix.ySize(),aResult.xSize()*aResult.ySize()); - int y1 = -aFilter.A(); - int y2 = aMatrix.ySize()-aFilter.B(); - int a2Size = 2*aMatrix.ySize()-1; - // Upper rim - for (int y = 0; y < y1; y++) - for (int x = 0; x < aMatrix.xSize(); x++) { - aResult(x,y) = 0; - for (int j = aFilter.A(); j < aFilter.B(); j++) { - if (y+j < 0) aResult(x,y) += aFilter[j]*aMatrix(x,-1-y-j); - else if (y+j >= aMatrix.ySize()) aResult(x,y) += aFilter[j]*aMatrix(x,a2Size-y-j); - else aResult(x,y) += aFilter[j]*aMatrix(x,y+j); - } - } - // Lower rim - for (int y = y2; y < aMatrix.ySize(); y++) - for (int x = 0; x < aMatrix.xSize(); x++) { - aResult(x,y) = 0; - for (int j = aFilter.A(); j < aFilter.B(); j++) { - if (y+j < 0) aResult(x,y) += aFilter[j]*aMatrix(x,-1-y-j); - else if (y+j >= aMatrix.ySize()) aResult(x,y) += aFilter[j]*aMatrix(x,a2Size-y-j); - else aResult(x,y) += aFilter[j]*aMatrix(x,y+j); - } - } - // Center - for (int y = y1; y < y2; y++) - for (int x = 0; x < aMatrix.xSize(); x++) { - aResult(x,y) = 0; - for (int j = aFilter.A(); j < aFilter.B(); j++) - aResult(x,y) += aFilter[j]*aMatrix(x,y+j); - } -} - -template -inline void filter(CMatrix& aMatrix, const CFilter2D& aFilter) { - CMatrix tempMatrix(aMatrix.xSize(),aMatrix.ySize()); - filter(aMatrix,tempMatrix,aFilter); - aMatrix = tempMatrix; -} - -template -void filter(const CMatrix& aMatrix, CMatrix& aResult, const CFilter2D& aFilter) { - if (aResult.xSize() != aMatrix.xSize() || aResult.ySize() != aMatrix.ySize()) - throw EFilterIncompatibleSize(aMatrix.xSize()*aMatrix.ySize(),aResult.xSize()*aResult.ySize()); - int x1 = -aFilter.AX(); - int y1 = -aFilter.AY(); - int x2 = aMatrix.xSize()-aFilter.BX(); - int y2 = aMatrix.ySize()-aFilter.BY(); - int a2XSize = 2*aMatrix.xSize()-1; - int a2YSize = 2*aMatrix.ySize()-1; - // Upper rim - for (int y = 0; y < y1; y++) - for (int x = 0; x < aMatrix.xSize(); x++) { - aResult(x,y) = 0; - for (int j = aFilter.AY(); j < aFilter.BY(); j++) { - int tempY; - if (y+j < 0) tempY = -1-y-j; - else if (y+j >= aMatrix.ySize()) tempY = a2YSize-y-j; - else tempY = y+j; - for (int i = aFilter.AX(); i < aFilter.BX(); i++) { - if (x+i < 0) aResult(x,y) += aFilter(i,j)*aMatrix(-1-x-i,tempY); - else if (x+i >= aMatrix.xSize()) aResult(x,y) += aFilter(i,j)*aMatrix(a2XSize-x-i,tempY); - else aResult(x,y) += aFilter(i,j)*aMatrix(x+i,tempY); - } - } - } - // Lower rim - for (int y = y2; y < aMatrix.ySize(); y++) - for (int x = 0; x < aMatrix.xSize(); x++) { - aResult(x,y) = 0; - for (int j = aFilter.AY(); j < aFilter.BY(); j++) { - int tempY; - if (y+j < 0) tempY = -1-y-j; - else if (y+j >= aMatrix.ySize()) tempY = a2YSize-y-j; - else tempY = y+j; - for (int i = aFilter.AX(); i < aFilter.BX(); i++) { - if (x+i < 0) aResult(x,y) += aFilter(i,j)*aMatrix(-1-x-i,tempY); - else if (x+i >= aMatrix.xSize()) aResult(x,y) += aFilter(i,j)*aMatrix(a2XSize-x-i,tempY); - else aResult(x,y) += aFilter(i,j)*aMatrix(x+i,tempY); - } - } - } - for (int y = y1; y < y2; y++) { - // Left rim - for (int x = 0; x < x1; x++) { - aResult(x,y) = 0; - for (int j = aFilter.AY(); j < aFilter.BY(); j++) { - for (int i = aFilter.AX(); i < aFilter.BX(); i++) { - if (x+i < 0) aResult(x,y) += aFilter(i,j)*aMatrix(-1-x-i,y+j); - else if (x+i >= aMatrix.xSize()) aResult(x,y) += aFilter(i,j)*aMatrix(a2XSize-x-i,y+j); - else aResult(x,y) += aFilter(i,j)*aMatrix(x+i,y+j); - } - } - } - // Right rim - for (int x = x2; x < aMatrix.xSize(); x++) { - aResult(x,y) = 0; - for (int j = aFilter.AY(); j < aFilter.BY(); j++) { - for (int i = aFilter.AX(); i < aFilter.BX(); i++) { - if (x+i < 0) aResult(x,y) += aFilter(i,j)*aMatrix(-1-x-i,y+j); - else if (x+i >= aMatrix.xSize()) aResult(x,y) += aFilter(i,j)*aMatrix(a2XSize-x-i,y+j); - else aResult(x,y) += aFilter(i,j)*aMatrix(x+i,y+j); - } - } - } - } - // Center - for (int y = y1; y < y2; y++) - for (int x = x1; x < x2; x++) { - aResult(x,y) = 0; - for (int j = aFilter.AY(); j < aFilter.BY(); j++) - for (int i = aFilter.AX(); i < aFilter.BX(); i++) - aResult(x,y) += aFilter(i,j)*aMatrix(x+i,y+j); - } -} - - - -template -inline void filterMin(CMatrix& aMatrix, const CFilter& aFilterX, const CFilter& aFilterY) { - CMatrix tempMatrix(aMatrix.xSize(),aMatrix.ySize()); - filterMin(aMatrix,aMatrix,tempMatrix,aFilterX,1); - CMatrix tmp(aMatrix); - filterMin(tempMatrix,tmp,aMatrix,1,aFilterY); -} - -template -inline void filterMin(const CMatrix& aMatrix, CMatrix& aResult, const CFilter& aFilterX, const CFilter& aFilterY) { - CMatrix tempMatrix(aMatrix.xSize(),aMatrix.ySize()); - filterMin(aMatrix,aMatrix,tempMatrix,aFilterX,1); - filterMin(tempMatrix,aMatrix,aResult,1,aFilterY); -} - -template -inline void filterMin(CMatrix& aMatrix, const CFilter& aFilter, const int aDummy) { - CMatrix tempMatrix(aMatrix.xSize(),aMatrix.ySize()); - filterMin(aMatrix,aMatrix,tempMatrix,aFilter,1); - aMatrix = tempMatrix; -} - -template -void filterMin(const CMatrix& aMatrix, const CMatrix& aOrig, CMatrix& aResult, const CFilter& aFilter, const int aDummy) { - if (aResult.xSize() != aMatrix.xSize() || aResult.ySize() != aMatrix.ySize()) - throw EFilterIncompatibleSize(aMatrix.xSize()*aMatrix.ySize(),aResult.xSize()*aResult.ySize()); - int x1 = -aFilter.A(); - int x2 = aMatrix.xSize()-aFilter.B(); - int a2Size = 2*aMatrix.xSize()-1; - aResult = 0; - for (int y = 0; y < aMatrix.ySize(); y++) { - int aOffset = y*aMatrix.xSize(); - // Left rim - for (int x = 0; x < x1; x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) { - int matrixIdx; - if (x+i < 0) matrixIdx = aOffset-1-x-i; - else if (x+i >= aMatrix.xSize()) matrixIdx = aOffset+a2Size-x-i; - else matrixIdx = aOffset+x+i; - if (matrixIdx == aOffset+x || aOrig.data()[matrixIdx] - 1e-5 <= aOrig.data()[aOffset+x]) - aResult.data()[aOffset+x] += aFilter[i]*aMatrix.data()[matrixIdx]; - } - // Center - for (int x = x1; x < x2; x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) - if (i == 0 || aOrig.data()[aOffset+x+i] - 1e-5 <= aOrig.data()[aOffset+x]) - aResult.data()[aOffset+x] += aFilter[i]*aMatrix.data()[aOffset+x+i]; - // Right rim - for (int x = x2; x < aMatrix.xSize(); x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) { - int matrixIdx; - if (x+i < 0) matrixIdx = aOffset-1-x-i; - else if (x+i >= aMatrix.xSize()) matrixIdx = aOffset+a2Size-x-i; - else matrixIdx = aOffset+x+i; - if (matrixIdx == aOffset+x || aOrig.data()[matrixIdx] - 1e-5 <= aOrig.data()[aOffset+x]) - aResult.data()[aOffset+x] += aFilter[i]*aMatrix.data()[matrixIdx]; - } - } -} - -template -inline void filterMin(CMatrix& aMatrix, const int aDummy, const CFilter& aFilter) { - CMatrix tempMatrix(aMatrix.xSize(),aMatrix.ySize()); - filterMin(aMatrix, aMatrix,tempMatrix,1,aFilter); - aMatrix = tempMatrix; -} - -template -void filterMin(const CMatrix& aMatrix, const CMatrix& aOrig, CMatrix& aResult, const int aDummy, const CFilter& aFilter) { - if (aResult.xSize() != aMatrix.xSize() || aResult.ySize() != aMatrix.ySize()) - throw EFilterIncompatibleSize(aMatrix.xSize()*aMatrix.ySize(),aResult.xSize()*aResult.ySize()); - int y1 = -aFilter.A(); - int y2 = aMatrix.ySize()-aFilter.B(); - int a2Size = 2*aMatrix.ySize()-1; - // Upper rim - for (int y = 0; y < y1; y++) - for (int x = 0; x < aMatrix.xSize(); x++) { - aResult(x,y) = 0; - for (int j = aFilter.A(); j < aFilter.B(); j++) { - int matrixIdx; - if (y+j < 0) matrixIdx = -1-y-j; - else if (y+j >= aMatrix.ySize()) matrixIdx = a2Size-y-j; - else matrixIdx = y+j; - if (matrixIdx == y || aOrig(x, matrixIdx) - 1e-5 <= aOrig(x, y)) - aResult(x,y) += aFilter[j]*aMatrix(x,matrixIdx); - } - } - // Lower rim - for (int y = y2; y < aMatrix.ySize(); y++) - for (int x = 0; x < aMatrix.xSize(); x++) { - aResult(x,y) = 0; - for (int j = aFilter.A(); j < aFilter.B(); j++) { - int matrixIdx; - if (y+j < 0) matrixIdx = -1-y-j; - else if (y+j >= aMatrix.ySize()) matrixIdx = a2Size-y-j; - else matrixIdx = y+j; - if (matrixIdx == y || aOrig(x, matrixIdx) - 1e-5 <= aOrig(x, y)) - aResult(x,y) += aFilter[j]*aMatrix(x,matrixIdx); - } - } - // Center - for (int y = y1; y < y2; y++) - for (int x = 0; x < aMatrix.xSize(); x++) { - aResult(x,y) = 0; - for (int j = aFilter.A(); j < aFilter.B(); j++) - if (j == 0 || aOrig(x,y+j) - 1e-5 <= aOrig(x,y)) - aResult(x,y) += aFilter[j]*aMatrix(x,y+j); - } -} - - - -// boxfilterX -template -inline void boxFilterX(CMatrix& aMatrix, int aWidth) { - CMatrix aTemp(aMatrix); - boxFilterX(aTemp,aMatrix,aWidth); -} - -template -void boxFilterX(const CMatrix& aMatrix, CMatrix& aResult, int aWidth) { - if (aWidth & 1 == 0) aWidth += 1; - T invWidth = 1.0/aWidth; - int halfWidth = (aWidth >> 1); - int aRight = halfWidth; - if (aRight >= aMatrix.xSize()) aRight = aMatrix.xSize()-1; - for (int y = 0; y < aMatrix.ySize(); y++) { - int aOffset = y*aMatrix.xSize(); - // Initialize - T aSum = 0.0f; - for (int x = aRight-aWidth+1; x <= aRight; x++) - if (x < 0) aSum += aMatrix.data()[aOffset-x-1]; - else aSum += aMatrix.data()[aOffset+x]; - // Shift - int xm = -halfWidth; - int xp = halfWidth+1; - for (int x = 0; x < aMatrix.xSize(); x++,xm++,xp++) { - aResult.data()[aOffset+x] = aSum*invWidth; - if (xm < 0) aSum -= aMatrix.data()[aOffset-xm-1]; - else aSum -= aMatrix.data()[aOffset+xm]; - if (xp >= aMatrix.xSize()) aSum += aMatrix.data()[aOffset+2*aMatrix.xSize()-1-xp]; - else aSum += aMatrix.data()[aOffset+xp]; - } - } -} - -// boxfilterY -template -inline void boxFilterY(CMatrix& aMatrix, int aWidth) { - CMatrix aTemp(aMatrix); - boxFilterY(aTemp,aMatrix,aWidth); -} - -template -void boxFilterY(const CMatrix& aMatrix, CMatrix& aResult, int aWidth) { - if (aWidth & 1 == 0) aWidth += 1; - T invWidth = 1.0/aWidth; - int halfWidth = (aWidth >> 1); - int aBottom = halfWidth; - if (aBottom >= aMatrix.ySize()) aBottom = aMatrix.xSize()-1; - for (int x = 0; x < aMatrix.xSize(); x++) { - // Initialize - T aSum = 0.0f; - for (int y = aBottom-aWidth+1; y <= aBottom; y++) - if (y < 0) aSum += aMatrix(x,-1-y); - else aSum += aMatrix(x,y); - // Shift - int ym = -halfWidth; - int yp = halfWidth+1; - for (int y = 0; y < aMatrix.ySize(); y++,ym++,yp++) { - aResult(x,y) = aSum*invWidth; - if (ym < 0) aSum -= aMatrix(x,-1-ym); - else aSum -= aMatrix(x,ym); - if (yp >= aMatrix.ySize()) aSum += aMatrix(x,2*aMatrix.ySize()-1-yp); - else aSum += aMatrix(x,yp); - } - } -} - -template -void recursiveSmoothX(CMatrix& aMatrix, float aSigma) { - CVector aVals1(aMatrix.xSize()); - CVector aVals2(aMatrix.xSize()); - float aAlpha = 2.5/(sqrt(NMath::Pi)*aSigma); - float aExp = exp(-aAlpha); - float aExpSqr = aExp*aExp; - float a2Exp = 2.0*aExp; - float k = (1.0-aExp)*(1.0-aExp)/(1.0+2.0*aAlpha*aExp-aExpSqr); - float aPreMinus = aExp*(aAlpha-1.0); - float aPrePlus = aExp*(aAlpha+1.0); - for (int y = 0; y < aMatrix.ySize(); y++) { - aVals1(0) = (0.5f-k*aPreMinus)*aMatrix(0,y); - aVals1(1) = k*(aMatrix(1,y)+aPreMinus*aMatrix(0,y))+(a2Exp-aExpSqr)*aVals1(0); - for (int x = 2; x < aMatrix.xSize(); x++) - aVals1(x) = k*(aMatrix(x,y)+aPreMinus*aMatrix(x-1,y))+a2Exp*aVals1(x-1)-aExpSqr*aVals1(x-2); - aVals2(aMatrix.xSize()-1) = (0.5f+k*aPreMinus)*aMatrix(aMatrix.xSize()-1,y); - aVals2(aMatrix.xSize()-2) = k*((aPrePlus-aExpSqr)*aMatrix(aMatrix.xSize()-1,y))+(a2Exp-aExpSqr)*aVals2(aMatrix.xSize()-1); - for (int x = aMatrix.xSize()-3; x >= 0; x--) - aVals2(x) = k*(aPrePlus*aMatrix(x+1,y)-aExpSqr*aMatrix(x+2,y))+a2Exp*aVals2(x+1)-aExpSqr*aVals2(x+2); - for (int x = 0; x < aMatrix.xSize(); x++) - aMatrix(x,y) = aVals1(x)+aVals2(x); - } -} - -template -void recursiveSmoothY(CMatrix& aMatrix, float aSigma) { - CVector aVals1(aMatrix.ySize()); - CVector aVals2(aMatrix.ySize()); - float aAlpha = 2.5/(sqrt(NMath::Pi)*aSigma); - float aExp = exp(-aAlpha); - float aExpSqr = aExp*aExp; - float a2Exp = 2.0*aExp; - float k = (1.0-aExp)*(1.0-aExp)/(1.0+2.0*aAlpha*aExp-aExpSqr); - float aPreMinus = aExp*(aAlpha-1.0); - float aPrePlus = aExp*(aAlpha+1.0); - for (int x = 0; x < aMatrix.xSize(); x++) { - aVals1(0) = (0.5f-k*aPreMinus)*aMatrix(x,0); - aVals1(1) = k*(aMatrix(x,1)+aPreMinus*aMatrix(x,0))+(a2Exp-aExpSqr)*aVals1(0); - for (int y = 2; y < aMatrix.ySize(); y++) - aVals1(y) = k*(aMatrix(x,y)+aPreMinus*aMatrix(x,y-1))+a2Exp*aVals1(y-1)-aExpSqr*aVals1(y-2); - aVals2(aMatrix.ySize()-1) = (0.5f+k*aPreMinus)*aMatrix(x,aMatrix.ySize()-1); - aVals2(aMatrix.ySize()-2) = k*((aPrePlus-aExpSqr)*aMatrix(x,aMatrix.ySize()-1))+(a2Exp-aExpSqr)*aVals2(aMatrix.ySize()-1); - for (int y = aMatrix.ySize()-3; y >= 0; y--) - aVals2(y) = k*(aPrePlus*aMatrix(x,y+1)-aExpSqr*aMatrix(x,y+2))+a2Exp*aVals2(y+1)-aExpSqr*aVals2(y+2); - for (int y = 0; y < aMatrix.ySize(); y++) - aMatrix(x,y) = aVals1(y)+aVals2(y); - } -} - -template -inline void recursiveSmooth(CMatrix& aMatrix, float aSigma) { - recursiveSmoothX(aMatrix,aSigma); - recursiveSmoothY(aMatrix,aSigma); -} - -// Linear 3D filtering --------------------------------------------------------- - -template -inline void filter(CTensor& aTensor, const CFilter& aFilterX, const CFilter& aFilterY, const CFilter& aFilterZ) { - CTensor tempTensor(aTensor.xSize(),aTensor.ySize(),aTensor.zSize()); - filter(aTensor,tempTensor,aFilterX,1,1); - filter(tempTensor,aTensor,1,aFilterY,1); - filter(aTensor,tempTensor,1,1,aFilterZ); - aTensor = tempTensor; -} - -template -inline void filter(const CTensor& aTensor, CTensor& aResult, const CFilter& aFilterX, const CFilter& aFilterY, const CFilter& aFilterZ) { - CTensor tempTensor(aTensor.xSize(),aTensor.ySize(),aTensor.zSize()); - filter(aTensor,aResult,aFilterX,1,1); - filter(aResult,tempTensor,1,aFilterY,1); - filter(tempTensor,aResult,1,1,aFilterZ); -} - -template -inline void filter(CTensor& aTensor, const CFilter& aFilter, const int aDummy1, const int aDummy2) { - CTensor tempTensor(aTensor.xSize(),aTensor.ySize(),aTensor.zSize()); - filter(aTensor,tempTensor,aFilter,1,1); - aTensor = tempTensor; -} - -template -void filter(const CTensor& aTensor, CTensor& aResult, const CFilter& aFilter, const int aDummy1, const int aDummy2) { - if (aResult.xSize() != aTensor.xSize() || aResult.ySize() != aTensor.ySize() || aResult.zSize() != aTensor.zSize()) - throw EFilterIncompatibleSize(aTensor.xSize()*aTensor.ySize()*aTensor.zSize(),aResult.xSize()*aResult.ySize()*aResult.zSize()); - int x1 = -aFilter.A(); - int x2 = aTensor.xSize()-aFilter.B(); - int a2Size = 2*aTensor.xSize()-1; - for (int z = 0; z < aTensor.zSize(); z++) - for (int y = 0; y < aTensor.ySize(); y++) { - // Left rim - for (int x = 0; x < x1; x++) { - aResult(x,y,z) = 0; - for (int i = aFilter.A(); i < aFilter.B(); i++) { - if (x+i < 0) aResult(x,y,z) += aFilter[i]*aTensor(-1-x-i,y,z); - else if (x+i >= aTensor.xSize()) aResult(x,y,z) += aFilter[i]*aTensor(a2Size-x-i,y,z); - else aResult(x,y,z) += aFilter[i]*aTensor(x+i,y,z); - } - } - // Center - for (int x = x1; x < x2; x++) { - aResult(x,y,z) = 0; - for (int i = aFilter.A(); i < aFilter.B(); i++) - aResult(x,y,z) += aFilter[i]*aTensor(x+i,y,z); - } - // Right rim - for (int x = x2; x < aTensor.xSize(); x++) { - aResult(x,y,z) = 0; - for (int i = aFilter.A(); i < aFilter.B(); i++) { - if (x+i < 0) aResult(x,y,z) += aFilter[i]*aTensor(-1-x-i,y,z); - else if (x+i >= aTensor.xSize()) aResult(x,y,z) += aFilter[i]*aTensor(a2Size-x-i,y,z); - else aResult(x,y,z) += aFilter[i]*aTensor(x+i,y,z); - } - } - } -} - -template -inline void filter(CTensor& aTensor, const int aDummy1, const CFilter& aFilter, const int aDummy2) { - CTensor tempTensor(aTensor.xSize(),aTensor.ySize(),aTensor.zSize()); - filter(aTensor,tempTensor,1,aFilter,1); - aTensor = tempTensor; -} - -template -void filter(const CTensor& aTensor, CTensor& aResult, const int aDummy1, const CFilter& aFilter, const int aDummy2) { - if (aResult.xSize() != aTensor.xSize() || aResult.ySize() != aTensor.ySize() || aResult.zSize() != aTensor.zSize()) - throw EFilterIncompatibleSize(aTensor.xSize()*aTensor.ySize()*aTensor.zSize(),aResult.xSize()*aResult.ySize()*aResult.zSize()); - int y1 = -aFilter.A(); - int y2 = aTensor.ySize()-aFilter.B(); - int a2Size = 2*aTensor.ySize()-1; - for (int z = 0; z < aTensor.zSize(); z++) { - // Upper rim - for (int y = 0; y < y1; y++) - for (int x = 0; x < aTensor.xSize(); x++) { - aResult(x,y,z) = 0; - for (int i = aFilter.A(); i < aFilter.B(); i++) { - if (y+i < 0) aResult(x,y,z) += aFilter[i]*aTensor(x,-1-y-i,z); - else if (y+i >= aTensor.ySize()) aResult(x,y,z) += aFilter[i]*aTensor(x,a2Size-y-i,z); - else aResult(x,y,z) += aFilter[i]*aTensor(x,y+i,z); - } - } - // Lower rim - for (int y = y2; y < aTensor.ySize(); y++) - for (int x = 0; x < aTensor.xSize(); x++) { - aResult(x,y,z) = 0; - for (int i = aFilter.A(); i < aFilter.B(); i++) { - if (y+i < 0) aResult(x,y,z) += aFilter[i]*aTensor(x,-1-y-i,z); - else if (y+i >= aTensor.ySize()) aResult(x,y,z) += aFilter[i]*aTensor(x,a2Size-y-i,z); - else aResult(x,y,z) += aFilter[i]*aTensor(x,y+i,z); - } - } - } - // Center - for (int z = 0; z < aTensor.zSize(); z++) - for (int y = y1; y < y2; y++) - for (int x = 0; x < aTensor.xSize(); x++) { - aResult(x,y,z) = 0; - for (int i = aFilter.A(); i < aFilter.B(); i++) - aResult(x,y,z) += aFilter[i]*aTensor(x,y+i,z); - } -} - -template -inline void filter(CTensor& aTensor, const int aDummy1, const int aDummy2, const CFilter& aFilter) { - CTensor tempTensor(aTensor.xSize(),aTensor.ySize(),aTensor.zSize()); - filter(aTensor,tempTensor,1,1,aFilter); - aTensor = tempTensor; -} - -template -void filter(const CTensor& aTensor, CTensor& aResult, const int aDummy1, const int aDummy2, const CFilter& aFilter) { - if (aResult.xSize() != aTensor.xSize() || aResult.ySize() != aTensor.ySize() || aResult.zSize() != aTensor.zSize()) - throw EFilterIncompatibleSize(aTensor.xSize()*aTensor.ySize()*aTensor.zSize(),aResult.xSize()*aResult.ySize()*aResult.zSize()); - int z1 = -aFilter.A(); - int z2 = aTensor.zSize()-aFilter.B(); - if (z2 < 0) z2 = 0; - int a2Size = 2*aTensor.zSize()-1; - // Front rim - for (int z = 0; z < z1; z++) - for (int y = 0; y < aTensor.ySize(); y++) - for (int x = 0; x < aTensor.xSize(); x++) { - aResult(x,y,z) = 0; - for (int i = aFilter.A(); i < aFilter.B(); i++) { - if (z+i < 0) aResult(x,y,z) += aFilter[i]*aTensor(x,y,-1-z-i); - else if (z+i >= aTensor.zSize()) aResult(x,y,z) += aFilter[i]*aTensor(x,y,a2Size-z-i); - else aResult(x,y,z) += aFilter[i]*aTensor(x,y,z+i); - } - } - // Back rim - for (int z = z2; z < aTensor.zSize(); z++) - for (int y = 0; y < aTensor.ySize(); y++) - for (int x = 0; x < aTensor.xSize(); x++) { - aResult(x,y,z) = 0; - for (int i = aFilter.A(); i < aFilter.B(); i++) { - if (z+i < 0) aResult(x,y,z) += aFilter[i]*aTensor(x,y,-1-z-i); - else if (z+i >= aTensor.zSize()) aResult(x,y,z) += aFilter[i]*aTensor(x,y,a2Size-z-i); - else aResult(x,y,z) += aFilter[i]*aTensor(x,y,z+i); - } - } - // Center - for (int z = z1; z < z2; z++) - for (int y = 0; y < aTensor.ySize(); y++) - for (int x = 0; x < aTensor.xSize(); x++) { - aResult(x,y,z) = 0; - for (int i = aFilter.A(); i < aFilter.B(); i++) - aResult(x,y,z) += aFilter[i]*aTensor(x,y,z+i); - } -} - -// boxfilterX -template -inline void boxFilterX(CTensor& aTensor, int aWidth) { - CTensor aTemp(aTensor); - boxFilterX(aTemp,aTensor,aWidth); -} - -template -void boxFilterX(const CTensor& aTensor, CTensor& aResult, int aWidth) { - if (aWidth % 2 == 0) aWidth += 1; - T* invWidth = new T[aWidth+1]; - invWidth[0] = 1.0f; - for (int i = 1; i <= aWidth; i++) - invWidth[i] = 1.0/i; - int halfWidth = (aWidth >> 1); - int aRight = halfWidth; - if (aRight >= aTensor.xSize()) aRight = aTensor.xSize()-1; - for (int z = 0; z < aTensor.zSize(); z++) - for (int y = 0; y < aTensor.ySize(); y++) { - int aOffset = (z*aTensor.ySize()+y)*aTensor.xSize(); - // Initialize - int aNum = 0; - T aSum = 0.0f; - for (int x = 0; x <= aRight; x++) { - aSum += aTensor.data()[aOffset+x]; aNum++; - } - // Shift - for (int x = 0; x < aTensor.xSize(); x++) { - aResult.data()[aOffset+x] = aSum*invWidth[aNum]; - if (x-halfWidth >= 0) { - aSum -= aTensor.data()[aOffset+x-halfWidth]; aNum--; - } - if (x+halfWidth+1 < aTensor.xSize()) { - aSum += aTensor.data()[aOffset+x+halfWidth+1]; aNum++; - } - } - } - delete[] invWidth; -} - -// boxfilterY -template -inline void boxFilterY(CTensor& aTensor, int aWidth) { - CTensor aTemp(aTensor); - boxFilterY(aTemp,aTensor,aWidth); -} - -template -void boxFilterY(const CTensor& aTensor, CTensor& aResult, int aWidth) { - if (aWidth % 2 == 0) aWidth += 1; - T* invWidth = new T[aWidth+1]; - invWidth[0] = 1.0f; - for (int i = 1; i <= aWidth; i++) - invWidth[i] = 1.0/i; - int halfWidth = (aWidth >> 1); - int aBottom = halfWidth; - if (aBottom >= aTensor.ySize()) aBottom = aTensor.ySize()-1; - for (int z = 0; z < aTensor.zSize(); z++) - for (int x = 0; x < aTensor.xSize(); x++) { - // Initialize - int aNum = 0; - T aSum = 0.0f; - for (int y = 0; y <= aBottom; y++) { - aSum += aTensor(x,y,z); aNum++; - } - // Shift - for (int y = 0; y < aTensor.ySize(); y++) { - aResult(x,y,z) = aSum*invWidth[aNum]; - if (y-halfWidth >= 0) { - aSum -= aTensor(x,y-halfWidth,z); aNum--; - } - if (y+halfWidth+1 < aTensor.ySize()) { - aSum += aTensor(x,y+halfWidth+1,z); aNum++; - } - } - } - delete[] invWidth; -} - -// boxfilterZ -template -inline void boxFilterZ(CTensor& aTensor, int aWidth) { - CTensor aTemp(aTensor); - boxFilterZ(aTemp,aTensor,aWidth); -} - -template -void boxFilterZ(const CTensor& aTensor, CTensor& aResult, int aWidth) { - if (aWidth % 2 == 0) aWidth += 1; - T* invWidth = new T[aWidth+1]; - invWidth[0] = 1.0f; - for (int i = 1; i <= aWidth; i++) - invWidth[i] = 1.0/i; - int halfWidth = (aWidth >> 1); - int aBottom = halfWidth; - if (aBottom >= aTensor.zSize()) aBottom = aTensor.zSize()-1; - for (int y = 0; y < aTensor.ySize(); y++) - for (int x = 0; x < aTensor.xSize(); x++) { - // Initialize - int aNum = 0; - T aSum = 0.0f; - for (int z = 0; z <= aBottom; z++) { - aSum += aTensor(x,y,z); aNum++; - } - // Shift - for (int z = 0; z < aTensor.zSize(); z++) { - aResult(x,y,z) = aSum*invWidth[aNum]; - if (z-halfWidth >= 0) { - aSum -= aTensor(x,y,z-halfWidth); aNum--; - } - if (z+halfWidth+1 < aTensor.zSize()) { - aSum += aTensor(x,y,z+halfWidth+1); aNum++; - } - } - } - delete[] invWidth; -} - -template -void recursiveSmoothX(CTensor& aTensor, float aSigma) { - CVector aVals1(aTensor.xSize()); - CVector aVals2(aTensor.xSize()); - float aAlpha = 2.5/(sqrt(NMath::Pi)*aSigma); - float aExp = exp(-aAlpha); - float aExpSqr = aExp*aExp; - float a2Exp = 2.0*aExp; - float k = (1.0-aExp)*(1.0-aExp)/(1.0+2.0*aAlpha*aExp-aExpSqr); - float aPreMinus = aExp*(aAlpha-1.0); - float aPrePlus = aExp*(aAlpha+1.0); - for (int z = 0; z < aTensor.zSize(); z++) - for (int y = 0; y < aTensor.ySize(); y++) { - int aOffset = (z*aTensor.ySize()+y)*aTensor.xSize(); - aVals1(0) = (0.5-k*aPreMinus)*aTensor.data()[aOffset]; - aVals1(1) = k*(aTensor.data()[aOffset+1]+aPreMinus*aTensor.data()[aOffset])+(2.0*aExp-aExpSqr)*aVals1(0); - for (int x = 2; x < aTensor.xSize(); x++) - aVals1(x) = k*(aTensor.data()[aOffset+x]+aPreMinus*aTensor.data()[aOffset+x-1])+a2Exp*aVals1(x-1)-aExpSqr*aVals1(x-2); - aVals2(aTensor.xSize()-1) = (0.5+k*aPreMinus)*aTensor.data()[aOffset+aTensor.xSize()-1]; - aVals2(aTensor.xSize()-2) = k*((aPrePlus-aExpSqr)*aTensor.data()[aOffset+aTensor.xSize()-1])+(a2Exp-aExpSqr)*aVals2(aTensor.xSize()-1); - for (int x = aTensor.xSize()-3; x >= 0; x--) - aVals2(x) = k*(aPrePlus*aTensor.data()[aOffset+x+1]-aExpSqr*aTensor.data()[aOffset+x+2])+a2Exp*aVals2(x+1)-aExpSqr*aVals2(x+2); - for (int x = 0; x < aTensor.xSize(); x++) - aTensor.data()[aOffset+x] = aVals1(x)+aVals2(x); - } -} - -template -void recursiveSmoothY(CTensor& aTensor, float aSigma) { - CVector aVals1(aTensor.ySize()); - CVector aVals2(aTensor.ySize()); - float aAlpha = 2.5/(sqrt(NMath::Pi)*aSigma); - float aExp = exp(-aAlpha); - float aExpSqr = aExp*aExp; - float a2Exp = 2.0*aExp; - float k = (1.0-aExp)*(1.0-aExp)/(1.0+2.0*aAlpha*aExp-aExpSqr); - float aPreMinus = aExp*(aAlpha-1.0); - float aPrePlus = aExp*(aAlpha+1.0); - for (int z = 0; z < aTensor.zSize(); z++) - for (int x = 0; x < aTensor.xSize(); x++) { - aVals1(0) = (0.5-k*aPreMinus)*aTensor(x,0,z); - aVals1(1) = k*(aTensor(x,1,z)+aPreMinus*aTensor(x,0,z))+(2.0*aExp-aExpSqr)*aVals1(0); - for (int y = 2; y < aTensor.ySize(); y++) - aVals1(y) = k*(aTensor(x,y,z)+aPreMinus*aTensor(x,y-1,z))+a2Exp*aVals1(y-1)-aExpSqr*aVals1(y-2); - aVals2(aTensor.ySize()-1) = (0.5+k*aPreMinus)*aTensor(x,aTensor.ySize()-1,z); - aVals2(aTensor.ySize()-2) = k*((aPrePlus-aExpSqr)*aTensor(x,aTensor.ySize()-1,z))+(a2Exp-aExpSqr)*aVals2(aTensor.ySize()-1); - for (int y = aTensor.ySize()-3; y >= 0; y--) - aVals2(y) = k*(aPrePlus*aTensor(x,y+1,z)-aExpSqr*aTensor(x,y+2,z))+a2Exp*aVals2(y+1)-aExpSqr*aVals2(y+2); - for (int y = 0; y < aTensor.ySize(); y++) - aTensor(x,y,z) = aVals1(y)+aVals2(y); - } -} - -template -void recursiveSmoothZ(CTensor& aTensor, float aSigma) { - CVector aVals1(aTensor.zSize()); - CVector aVals2(aTensor.zSize()); - float aAlpha = 2.5/(sqrt(NMath::Pi)*aSigma); - float aExp = exp(-aAlpha); - float aExpSqr = aExp*aExp; - float a2Exp = 2.0*aExp; - float k = (1.0-aExp)*(1.0-aExp)/(1.0+2.0*aAlpha*aExp-aExpSqr); - float aPreMinus = aExp*(aAlpha-1.0); - float aPrePlus = aExp*(aAlpha+1.0); - for (int y = 0; y < aTensor.ySize(); y++) - for (int x = 0; x < aTensor.xSize(); x++) { - aVals1(0) = (0.5-k*aPreMinus)*aTensor(x,y,0); - aVals1(1) = k*(aTensor(x,y,1)+aPreMinus*aTensor(x,y,0))+(2.0*aExp-aExpSqr)*aVals1(0); - for (int z = 2; z < aTensor.zSize(); z++) - aVals1(z) = k*(aTensor(x,y,z)+aPreMinus*aTensor(x,y,z-1))+a2Exp*aVals1(z-1)-aExpSqr*aVals1(z-2); - aVals2(aTensor.zSize()-1) = (0.5+k*aPreMinus)*aTensor(x,y,aTensor.zSize()-1); - aVals2(aTensor.zSize()-2) = k*((aPrePlus-aExpSqr)*aTensor(x,y,aTensor.zSize()-1))+(a2Exp-aExpSqr)*aVals2(aTensor.zSize()-1); - for (int z = aTensor.zSize()-3; z >= 0; z--) - aVals2(z) = k*(aPrePlus*aTensor(x,y,z+1)-aExpSqr*aTensor(x,y,z+2))+a2Exp*aVals2(z+1)-aExpSqr*aVals2(z+2); - for (int z = 0; z < aTensor.zSize(); z++) - aTensor(x,y,z) = aVals1(z)+aVals2(z); - } -} - -// Linear 4D filtering --------------------------------------------------------- - -template -inline void filter(CTensor4D& aTensor, const CFilter& aFilterX, const CFilter& aFilterY, const CFilter& aFilterZ, const CFilter& aFilterA) { - CTensor4D tempTensor(aTensor.xSize(),aTensor.ySize(),aTensor.zSize()); - filter(aTensor,tempTensor,aFilterX,1,1,1); - filter(tempTensor,aTensor,1,aFilterY,1,1); - filter(aTensor,tempTensor,1,1,aFilterZ,1); - filter(tempTensor,aTensor,1,1,1,aFilterA); -} - -template -inline void filter(CTensor4D& aTensor, const CFilter& aFilter, const int aDummy1, const int aDummy2, const int aDummy3) { - CTensor4D tempTensor(aTensor.xSize(),aTensor.ySize(),aTensor.zSize(),aTensor.aSize()); - filter(aTensor,tempTensor,aFilter,1,1,1); - aTensor = tempTensor; -} - -template -void filter(const CTensor4D& aTensor, CTensor4D& aResult, const CFilter& aFilter, const int aDummy1, const int aDummy2, const int aDummy3) { - if (aResult.xSize() != aTensor.xSize() || aResult.ySize() != aTensor.ySize() || aResult.zSize() != aTensor.zSize() || aResult.aSize() != aTensor.aSize()) - throw EFilterIncompatibleSize(aTensor.xSize()*aTensor.ySize()*aTensor.zSize()*aTensor.aSize(),aResult.xSize()*aResult.ySize()*aResult.zSize()*aResult.aSize()); - int x1 = -aFilter.A(); - int x2 = aTensor.xSize()-aFilter.B(); - int a2Size = 2*aTensor.xSize()-1; - aResult = 0; - for (int a = 0; a < aTensor.aSize(); a++) - for (int z = 0; z < aTensor.zSize(); z++) - for (int y = 0; y < aTensor.ySize(); y++) { - int aOffset = aTensor.xSize()*(y+aTensor.ySize()*(z+aTensor.zSize()*a)); - // Left rim - for (int x = 0; x < x1; x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) { - if (x+i < 0) aResult.data()[aOffset+x] += aFilter[i]*aTensor.data()[aOffset-1-x-i]; - else if (x+i >= aTensor.xSize()) aResult.data()[aOffset+x] += aFilter[i]*aTensor.data()[aOffset+a2Size-x-i]; - else aResult.data()[aOffset+x] += aFilter[i]*aTensor.data()[x+i+aOffset]; - } - // Center - for (int x = x1; x < x2; x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) - aResult.data()[aOffset+x] += aFilter[i]*aTensor.data()[aOffset+x+i]; - // Right rim - for (int x = x2; x < aTensor.xSize(); x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) { - if (x+i < 0) aResult.data()[aOffset+x] += aFilter[i]*aTensor.data()[aOffset-1-x-i]; - else if (x+i >= aTensor.xSize()) aResult.data()[aOffset+x] += aFilter[i]*aTensor.data()[aOffset+a2Size-x-i]; - else aResult.data()[aOffset+x] += aFilter[i]*aTensor.data()[x+i+aOffset]; - } - } -} - -template -inline void filter(CTensor4D& aTensor, const int aDummy1, const CFilter& aFilter, const int aDummy2, const int aDummy3) { - CTensor4D tempTensor(aTensor.xSize(),aTensor.ySize(),aTensor.zSize(),aTensor.aSize()); - filter(aTensor,tempTensor,1,aFilter,1,1); - aTensor = tempTensor; -} - -template -void filter(const CTensor4D& aTensor, CTensor4D& aResult, const int aDummy1, const CFilter& aFilter, const int aDummy2, const int aDummy3) { - if (aResult.xSize() != aTensor.xSize() || aResult.ySize() != aTensor.ySize() || aResult.zSize() != aTensor.zSize() || aResult.aSize() != aTensor.aSize()) - throw EFilterIncompatibleSize(aTensor.xSize()*aTensor.ySize()*aTensor.zSize()*aTensor.aSize(),aResult.xSize()*aResult.ySize()*aResult.zSize()*aResult.aSize()); - int y1 = -aFilter.A(); - int y2 = aTensor.ySize()-aFilter.B(); - int a2Size = 2*aTensor.ySize()-1; - aResult = 0; - for (int a = 0; a < aTensor.aSize(); a++) { - for (int z = 0; z < aTensor.zSize(); z++) { - // Upper rim - for (int y = 0; y < y1; y++) - for (int x = 0; x < aTensor.xSize(); x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) { - if (y+i < 0) aResult(x,y,z,a) += aFilter[i]*aTensor(x,-1-y-i,z,a); - else if (y+i >= aTensor.ySize()) aResult(x,y,z,a) += aFilter[i]*aTensor(x,a2Size-y-i,z,a); - else aResult(x,y,z,a) += aFilter[i]*aTensor(x,y+i,z,a); - } - // Lower rim - for (int y = y2; y < aTensor.ySize(); y++) - for (int x = 0; x < aTensor.xSize(); x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) { - if (y+i < 0) aResult(x,y,z,a) += aFilter[i]*aTensor(x,-1-y-i,z,a); - else if (y+i >= aTensor.ySize()) aResult(x,y,z,a) += aFilter[i]*aTensor(x,a2Size-y-i,z,a); - else aResult(x,y,z,a) += aFilter[i]*aTensor(x,y+i,z,a); - } - } - // Center - for (int z = 0; z < aTensor.zSize(); z++) - for (int y = y1; y < y2; y++) - for (int x = 0; x < aTensor.xSize(); x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) - aResult(x,y,z,a) += aFilter[i]*aTensor(x,y+i,z,a); - } -} - -template -inline void filter(CTensor4D& aTensor, const int aDummy1, const int aDummy2, const CFilter& aFilter, const int aDummy3) { - CTensor4D tempTensor(aTensor.xSize(),aTensor.ySize(),aTensor.zSize(),aTensor.aSize()); - filter(aTensor,tempTensor,1,1,aFilter,1); - aTensor = tempTensor; -} - -template -void filter(const CTensor4D& aTensor, CTensor4D& aResult, const int aDummy1, const int aDummy2, const CFilter& aFilter, const int aDummy3) { - if (aResult.xSize() != aTensor.xSize() || aResult.ySize() != aTensor.ySize() || aResult.zSize() != aTensor.zSize() || aResult.aSize() != aTensor.aSize()) - throw EFilterIncompatibleSize(aTensor.xSize()*aTensor.ySize()*aTensor.zSize()*aTensor.aSize(),aResult.xSize()*aResult.ySize()*aResult.zSize()*aResult.aSize()); - int z1 = -aFilter.A(); - int z2 = aTensor.zSize()-aFilter.B(); - int a2Size = 2*aTensor.zSize()-1; - aResult = 0; - for (int a = 0; a < aTensor.aSize(); a++) { - // Front rim - for (int z = 0; z < z1; z++) - for (int y = 0; y < aTensor.ySize(); y++) - for (int x = 0; x < aTensor.xSize(); x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) { - if (z+i < 0) aResult(x,y,z,a) += aFilter[i]*aTensor(x,y,-1-z-i,a); - else if (z+i >= aTensor.zSize()) aResult(x,y,z,a) += aFilter[i]*aTensor(x,y,a2Size-z-i,a); - else aResult(x,y,z,a) += aFilter[i]*aTensor(x,y,z+i,a); - } - // Back rim - for (int z = z2; z < aTensor.zSize(); z++) - for (int y = 0; y < aTensor.ySize(); y++) - for (int x = 0; x < aTensor.xSize(); x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) { - if (z+i < 0) aResult(x,y,z,a) += aFilter[i]*aTensor(x,y,-1-z-i,a); - else if (z+i >= aTensor.zSize()) aResult(x,y,z,a) += aFilter[i]*aTensor(x,y,a2Size-z-i,a); - else aResult(x,y,z,a) += aFilter[i]*aTensor(x,y,z+i,a); - } - // Center - for (int z = z1; z < z2; z++) - for (int y = 0; y < aTensor.ySize(); y++) - for (int x = 0; x < aTensor.xSize(); x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) - aResult(x,y,z,a) += aFilter[i]*aTensor(x,y,z+i,a); - } -} - -template -inline void filter(CTensor4D& aTensor, const int aDummy1, const int aDummy2, const int aDummy3, const CFilter& aFilter) { - CTensor4D tempTensor(aTensor.xSize(),aTensor.ySize(),aTensor.zSize(),aTensor.aSize()); - filter(aTensor,tempTensor,1,1,1,aFilter); - aTensor = tempTensor; -} - -template -void filter(const CTensor4D& aTensor, CTensor4D& aResult, const int aDummy1, const int aDummy2, const int aDummy3, const CFilter& aFilter) { - if (aResult.xSize() != aTensor.xSize() || aResult.ySize() != aTensor.ySize() || aResult.zSize() != aTensor.zSize() || aResult.aSize() != aTensor.aSize()) - throw EFilterIncompatibleSize(aTensor.xSize()*aTensor.ySize()*aTensor.zSize()*aTensor.aSize(),aResult.xSize()*aResult.ySize()*aResult.zSize()*aResult.aSize()); - int a1 = -aFilter.A(); - int a2 = aTensor.aSize()-aFilter.B(); - int a2Size = 2*aTensor.aSize()-1; - aResult = 0; - // Front rim - for (int a = 0; a < a1; a++) - for (int z = 0; z < aTensor.zSize(); z++) - for (int y = 0; y < aTensor.ySize(); y++) - for (int x = 0; x < aTensor.xSize(); x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) { - if (a+i < 0) aResult(x,y,z,a) += aFilter[i]*aTensor(x,y,z,-1-a-i); - else if (a+i >= aTensor.aSize()) aResult(x,y,z,a) += aFilter[i]*aTensor(x,y,z,a2Size-a-i); - else aResult(x,y,z,a) += aFilter[i]*aTensor(x,y,z,a+i); - } - // Back rim - for (int a = a2; a < aTensor.aSize(); a++) - for (int z = 0; z < aTensor.zSize(); z++) - for (int y = 0; y < aTensor.ySize(); y++) - for (int x = 0; x < aTensor.xSize(); x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) { - if (a+i < 0) aResult(x,y,z,a) += aFilter[i]*aTensor(x,y,z,-1-a-i); - else if (a+i >= aTensor.aSize()) aResult(x,y,z,a) += aFilter[i]*aTensor(x,y,z,a2Size-a-i); - else aResult(x,y,z,a) += aFilter[i]*aTensor(x,y,z,a+i); - } - // Center - for (int a = a1; a < a2; a++) - for (int z = 0; z < aTensor.zSize(); z++) - for (int y = 0; y < aTensor.ySize(); y++) - for (int x = 0; x < aTensor.xSize(); x++) - for (int i = aFilter.A(); i < aFilter.B(); i++) - aResult(x,y,z,a) += aFilter[i]*aTensor(x,y,z,a+i); -} - -template -void recursiveSmoothX(CTensor4D& aTensor, float aSigma) { - CVector aVals1(aTensor.xSize()); - CVector aVals2(aTensor.xSize()); - float aAlpha = 2.5/(sqrt(NMath::Pi)*aSigma); - float aExp = exp(-aAlpha); - float aExpSqr = aExp*aExp; - float a2Exp = 2.0*aExp; - float k = (1.0-aExp)*(1.0-aExp)/(1.0+2.0*aAlpha*aExp-aExpSqr); - float aPreMinus = aExp*(aAlpha-1.0); - float aPrePlus = aExp*(aAlpha+1.0); - for (int a = 0; a < aTensor.aSize(); a++) - for (int z = 0; z < aTensor.zSize(); z++) - for (int y = 0; y < aTensor.ySize(); y++) { - int aOffset = ((a*aTensor.zSize()+z)*aTensor.ySize()+y)*aTensor.xSize(); - aVals1(0) = (0.5-k*aPreMinus)*aTensor.data()[aOffset]; - aVals1(1) = k*(aTensor.data()[aOffset+1]+aPreMinus*aTensor.data()[aOffset])+(2.0*aExp-aExpSqr)*aVals1(0); - for (int x = 2; x < aTensor.xSize(); x++) - aVals1(x) = k*(aTensor.data()[aOffset+x]+aPreMinus*aTensor.data()[aOffset+x-1])+a2Exp*aVals1(x-1)-aExpSqr*aVals1(x-2); - aVals2(aTensor.xSize()-1) = (0.5+k*aPreMinus)*aTensor.data()[aOffset+aTensor.xSize()-1]; - aVals2(aTensor.xSize()-2) = k*((aPrePlus-aExpSqr)*aTensor.data()[aOffset+aTensor.xSize()-1])+(a2Exp-aExpSqr)*aVals2(aTensor.xSize()-1); - for (int x = aTensor.xSize()-3; x >= 0; x--) - aVals2(x) = k*(aPrePlus*aTensor.data()[aOffset+x+1]-aExpSqr*aTensor.data()[aOffset+x+2])+a2Exp*aVals2(x+1)-aExpSqr*aVals2(x+2); - for (int x = 0; x < aTensor.xSize(); x++) - aTensor.data()[aOffset+x] = aVals1(x)+aVals2(x); - } -} - -template -void recursiveSmoothY(CTensor4D& aTensor, float aSigma) { - CVector aVals1(aTensor.ySize()); - CVector aVals2(aTensor.ySize()); - CVector aVals3(aTensor.ySize()); - float aAlpha = 2.5/(sqrt(NMath::Pi)*aSigma); - float aExp = exp(-aAlpha); - float aExpSqr = aExp*aExp; - float a2Exp = 2.0*aExp; - float k = (1.0-aExp)*(1.0-aExp)/(1.0+2.0*aAlpha*aExp-aExpSqr); - float aPreMinus = aExp*(aAlpha-1.0); - float aPrePlus = aExp*(aAlpha+1.0); - for (int a = 0; a < aTensor.aSize(); a++) - for (int z = 0; z < aTensor.zSize(); z++) - for (int x = 0; x < aTensor.xSize(); x++) { - for (int y = 0; y < aTensor.ySize(); y++) - aVals3(y) = aTensor(x,y,z,a); - aVals1(0) = (0.5-k*aPreMinus)*aVals3(0); - aVals1(1) = k*(aVals3(1)+aPreMinus*aVals3(0))+(2.0*aExp-aExpSqr)*aVals1(0); - for (int y = 2; y < aTensor.ySize(); y++) - aVals1(y) = k*(aVals3(y)+aPreMinus*aVals3(y-1))+a2Exp*aVals1(y-1)-aExpSqr*aVals1(y-2); - aVals2(aTensor.ySize()-1) = (0.5+k*aPreMinus)*aVals3(aTensor.ySize()-1); - aVals2(aTensor.ySize()-2) = k*((aPrePlus-aExpSqr)*aVals3(aTensor.ySize()-1))+(a2Exp-aExpSqr)*aVals2(aTensor.ySize()-1); - for (int y = aTensor.ySize()-3; y >= 0; y--) - aVals2(y) = k*(aPrePlus*aVals3(y+1)-aExpSqr*aVals3(y+2))+a2Exp*aVals2(y+1)-aExpSqr*aVals2(y+2); - for (int y = 0; y < aTensor.ySize(); y++) - aTensor(x,y,z,a) = aVals1(y)+aVals2(y); - } -} - -template -void recursiveSmoothZ(CTensor4D& aTensor, float aSigma) { - CVector aVals1(aTensor.zSize()); - CVector aVals2(aTensor.zSize()); - CVector aVals3(aTensor.zSize()); - float aAlpha = 2.5/(sqrt(NMath::Pi)*aSigma); - float aExp = exp(-aAlpha); - float aExpSqr = aExp*aExp; - float a2Exp = 2.0*aExp; - float k = (1.0-aExp)*(1.0-aExp)/(1.0+2.0*aAlpha*aExp-aExpSqr); - float aPreMinus = aExp*(aAlpha-1.0); - float aPrePlus = aExp*(aAlpha+1.0); - for (int a = 0; a < aTensor.aSize(); a++) - for (int y = 0; y < aTensor.ySize(); y++) - for (int x = 0; x < aTensor.xSize(); x++) { - for (int z = 0; z < aTensor.zSize(); z++) - aVals3(z) = aTensor(x,y,z,a); - aVals1(0) = (0.5-k*aPreMinus)*aVals3(0); - aVals1(1) = k*(aVals3(1)+aPreMinus*aVals3(0))+(2.0*aExp-aExpSqr)*aVals1(0); - for (int z = 2; z < aTensor.zSize(); z++) - aVals1(z) = k*(aVals3(z)+aPreMinus*aVals3(z-1))+a2Exp*aVals1(z-1)-aExpSqr*aVals1(z-2); - aVals2(aTensor.zSize()-1) = (0.5+k*aPreMinus)*aVals3(aTensor.zSize()-1); - aVals2(aTensor.zSize()-2) = k*((aPrePlus-aExpSqr)*aVals3(aTensor.zSize()-1))+(a2Exp-aExpSqr)*aVals2(aTensor.zSize()-1); - for (int z = aTensor.zSize()-3; z >= 0; z--) - aVals2(z) = k*(aPrePlus*aVals3(z+1)-aExpSqr*aVals3(z+2))+a2Exp*aVals2(z+1)-aExpSqr*aVals2(z+2); - for (int z = 0; z < aTensor.zSize(); z++) - aTensor(x,y,z,a) = aVals1(z)+aVals2(z); - } -} - -template -void recursiveSmoothA(CTensor4D& aTensor, float aSigma) { - CVector aVals1(aTensor.aSize()); - CVector aVals2(aTensor.aSize()); - CVector aVals3(aTensor.aSize()); - float aAlpha = 2.5/(sqrt(NMath::Pi)*aSigma); - float aExp = exp(-aAlpha); - float aExpSqr = aExp*aExp; - float a2Exp = 2.0*aExp; - float k = (1.0-aExp)*(1.0-aExp)/(1.0+2.0*aAlpha*aExp-aExpSqr); - float aPreMinus = aExp*(aAlpha-1.0); - float aPrePlus = aExp*(aAlpha+1.0); - for (int z = 0; z < aTensor.zSize(); z++) - for (int y = 0; y < aTensor.ySize(); y++) - for (int x = 0; x < aTensor.xSize(); x++) { - for (int a = 0; a < aTensor.aSize(); a++) - aVals3(a) = aTensor(x,y,z,a); - aVals1(0) = (0.5-k*aPreMinus)*aVals3(0); - aVals1(1) = k*(aVals3(1)+aPreMinus*aVals3(0))+(2.0*aExp-aExpSqr)*aVals1(0); - for (int a = 2; a < aTensor.aSize(); a++) - aVals1(a) = k*(aVals3(a)+aPreMinus*aVals3(a-1))+a2Exp*aVals1(a-1)-aExpSqr*aVals1(a-2); - aVals2(aTensor.aSize()-1) = (0.5+k*aPreMinus)*aVals3(aTensor.aSize()-1); - aVals2(aTensor.aSize()-2) = k*((aPrePlus-aExpSqr)*aVals3(aTensor.aSize()-1))+(a2Exp-aExpSqr)*aVals2(aTensor.aSize()-1); - for (int a = aTensor.aSize()-3; a >= 0; a--) - aVals2(a) = k*(aPrePlus*aVals3(a+1)-aExpSqr*aVals3(a+2))+a2Exp*aVals2(a+1)-aExpSqr*aVals2(a+2); - for (int a = 0; a < aTensor.aSize(); a++) - aTensor(x,y,z,a) = aVals1(a)+aVals2(a); - } -} - -// Nonlinear filters ----------------------------------------------------------- - -// osher (2D) -template -void osher(CMatrix& aData, int aIterations) { - CMatrix aDiff(aData.xSize(),aData.ySize()); - for (int t = 0; t < aIterations; t++) { - for (int y = 0; y < aData.ySize(); y++) - for (int x = 0; x < aData.xSize(); x++) { - T u00,u01,u02,u10,u11,u12,u20,u21,u22; - if (x > 0) { - if (y > 0) u00 = aData(x-1,y-1); - else u00 = aData(x-1,0); - u01 = aData(x-1,y); - if (y < aData.ySize()-1) u02 = aData(x-1,y+1); - else u02 = aData(x-1,y); - } - else { - if (y > 0) u00 = aData(0,y-1); - else u00 = aData(0,0); - u01 = aData(0,y); - if (y < aData.ySize()-1) u02 = aData(0,y+1); - else u02 = aData(0,y); - } - if (y > 0) u10 = aData(x,y-1); - else u10 = aData(x,y); - u11 = aData(x,y); - if (y < aData.ySize()-1) u12 = aData(x,y+1); - else u12 = aData(x,y); - if (x < aData.xSize()-1) { - if (y > 0) u20 = aData(x+1,y-1); - else u20 = aData(x+1,y); - u21 = aData(x+1,y); - if (y < aData.ySize()-1) u22 = aData(x+1,y+1); - else u22 = aData(x+1,y); - } - else { - if (y > 0) u20 = aData(x,y-1); - else u20 = aData(x,y); - u21 = aData(x,y); - if (y < aData.ySize()-1) u22 = aData(x,y+1); - else u22 = aData(x,y); - } - T ux = 0.5*(u21-u01); - T uy = 0.5*(u12-u10); - T uxuy = ux*uy; - T uxx = u01-2.0*u11+u21; - T uyy = u10-2.0*u11+u12; - T uxy; - if (uxuy < 0) uxy = 2.0*u11+u00+u22-u10-u12-u01-u21; - else uxy = u10+u12+u01+u21-2.0*u11-u02-u20; - T laPlace = uyy*uy*uy+uxy*uxuy+uxx*ux*ux; - T uxLeft = u11-u01; - T uxRight = u21-u11; - T uyUp = u11-u10; - T uyDown = u12-u11; - if (laPlace < 0) { - T aSum = 0; - if (uxRight > 0) aSum += uxRight*uxRight; - if (uxLeft < 0) aSum += uxLeft*uxLeft; - if (uyDown > 0) aSum += uyDown*uyDown; - if (uyUp < 0) aSum += uyUp*uyUp; - aDiff(x,y) = sqrt(aSum); - } - else if (laPlace > 0) { - T aSum = 0; - if (uxRight < 0) aSum += uxRight*uxRight; - if (uxLeft > 0) aSum += uxLeft*uxLeft; - if (uyDown < 0) aSum += uyDown*uyDown; - if (uyUp > 0) aSum += uyUp*uyUp; - aDiff(x,y) = -sqrt(aSum); - } - } - for (int i = 0; i < aData.size(); i++) - aData.data()[i] += 0.25*aDiff.data()[i]; - } -} - -template -inline void osher(const CMatrix& aData, CMatrix& aResult, int aIterations) { - aResult = aData; - osher(aResult,aIterations); -} - -} - -#endif - -- cgit v1.2.3-70-g09d2