From 200496c4216ab08824851c15632b1a1db6e76b31 Mon Sep 17 00:00:00 2001 From: Ryan Baumann Date: Fri, 29 Jul 2016 15:34:38 -0400 Subject: Add files/code from https://github.com/manuelruder/artistic-videos --- consistencyChecker/CFilter.h | 2210 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 2210 insertions(+) create mode 100644 consistencyChecker/CFilter.h (limited to 'consistencyChecker/CFilter.h') diff --git a/consistencyChecker/CFilter.h b/consistencyChecker/CFilter.h new file mode 100644 index 0000000..29bef9e --- /dev/null +++ b/consistencyChecker/CFilter.h @@ -0,0 +1,2210 @@ +// - 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