📄 nddirectionalfilterbank.cpp
字号:
void NdDirectionalFilterBank::FilterDownsampleF(int nDims, int pDims[], double *pData, int K, int M, int Level)
{
// Some validity checks
assert((nDims >= 2) && (Level >= 1) && (K >= 0) && (K < nDims) && (M >= 0) && (M < nDims) && (K != M));
assert((FilterSlice0.nx == pDims[M]) && (FilterSlice0.ny / 4 + 1 == pDims[K]));
int nChannels = (int)pow(2.0, (double)Level);
int nRepeat = nChannels / 2;
// We will divide this dimension into nChannel pieces
assert(pDims[M] % nChannels == 0);
// Record the downsampling dimension
int downsample_dim = M;
// We always make M < K
if (M > K)
{
// Note: we are transposing a complex-valued matrix
FilterSlice0.Transpose(false);
FilterSlice1.Transpose(false);
// Swap the two dimensions
int tmp;
tmp = K;
K = M;
M = tmp;
}
int dimA, dimB, dimC, dimD, dimE;
register int s, m, n, k, j, i;
dimA = 1;
for (m = 0; m < M; m++)
dimA *= pDims[m];
dimB = pDims[M];
dimC = 1;
for (m = M + 1; m < K; m++)
dimC *= pDims[m];
dimD = pDims[K];
dimE = 1;
for (m = K + 1; m < nDims; m++)
dimE *= pDims[m];
double *pFilter0_base, *pFilter1_base, *pFilter0_up, *pFilter0_down, *pFilter1_up, *pFilter1_down;
pFilter0_base = FilterSlice0.GetPointer();
pFilter1_base = FilterSlice1.GetPointer();
double f0_up_real, f0_up_imag, f0_down_real, f0_down_imag;
double f1_up_real, f1_up_imag, f1_down_real, f1_down_imag;
double *pArray_up, *pArray_down;
double a_up_real, a_up_imag, a_down_real, a_down_imag;
int filter_block_skip, filter_line_skip, filter_line_retreat;
int array_block_skip;
// There are two cases:
// (1) downsample_dim == M
// (2) downsample_dim == K
if (downsample_dim == M)
{
// preparation
dimB /= nChannels;
filter_block_skip = dimB * FilterSlice0.ny;
filter_line_retreat = 2 * dimD;
filter_line_skip = FilterSlice0.ny;
array_block_skip = dimB * 2;
for (i = M + 1; i < nDims; i++)
array_block_skip *= pDims[i];
pArray_up = pData;
pArray_down = pArray_up + array_block_skip;
for (i = 0; i < dimA; i++)
{
pFilter0_up = pFilter0_base;
pFilter0_down = pFilter0_up + filter_block_skip;
pFilter1_up = pFilter1_base;
pFilter1_down = pFilter1_up + filter_block_skip;
for (j = 0; j < nRepeat; j++)
{
for (k = 0; k < dimB; k++)
{
for (m = 0; m < dimC; m++)
{
for (n = 0; n < dimD; n++)
{
f0_up_real = *(pFilter0_up++); f0_up_imag = *(pFilter0_up++);
f0_down_real = *(pFilter0_down++); f0_down_imag = *(pFilter0_down++);
f1_up_real = *(pFilter1_up++); f1_up_imag = *(pFilter1_up++);
f1_down_real = *(pFilter1_down++); f1_down_imag = *(pFilter1_down++);
for (s = 0; s < dimE; s++)
{
a_up_real = *(pArray_up++);
a_up_imag = *(pArray_up++);
a_down_real = *(pArray_down++);
a_down_imag = *(pArray_down++);
*(pArray_up - 2) = a_up_real * f0_up_real - a_up_imag * f0_up_imag
+ a_down_real * f0_down_real - a_down_imag * f0_down_imag;
*(pArray_up - 1) = a_up_real * f0_up_imag + a_up_imag * f0_up_real
+ a_down_real * f0_down_imag + a_down_imag * f0_down_real;
*(pArray_down - 2) = a_up_real * f1_up_real - a_up_imag * f1_up_imag
+ a_down_real * f1_down_real - a_down_imag * f1_down_imag;
*(pArray_down - 1) = a_up_real * f1_up_imag + a_up_imag * f1_up_real
+ a_down_real * f1_down_imag + a_down_imag * f1_down_real;
} // s = 0 : dime
} // n = 0 : dimD
pFilter0_up -= filter_line_retreat;
pFilter0_down -= filter_line_retreat;
pFilter1_up -= filter_line_retreat;
pFilter1_down -= filter_line_retreat;
} // m = 0 : dimC
pFilter0_up += filter_line_skip;
pFilter0_down += filter_line_skip;
pFilter1_up += filter_line_skip;
pFilter1_down += filter_line_skip;
} // k = 0 : dimB
pFilter0_up = pFilter0_down;
pFilter0_down = pFilter0_up + filter_block_skip;
pFilter1_up = pFilter1_down;
pFilter1_down = pFilter1_up + filter_block_skip;
pArray_up = pArray_down;
pArray_down = pArray_up + array_block_skip;
} // j = 0 : nRepeat
} // i = 0 : dimA
} // if downsample_dim == M
else
{
// preparation
filter_block_skip = FilterSlice0.ny / nChannels;
filter_line_skip = FilterSlice0.ny;
dimD /= nChannels;
array_block_skip = dimD * 2;
for (i = K + 1; i < nDims; i++)
array_block_skip *= pDims[i];
pArray_up = pData;
pArray_down = pArray_up + array_block_skip;
for (i = 0; i < dimA; i++)
{
pFilter0_up = pFilter0_base;
pFilter0_down = pFilter0_up + filter_block_skip;
pFilter1_up = pFilter1_base;
pFilter1_down = pFilter1_up + filter_block_skip;
for (j = 0; j < dimB; j++)
{
for (k = 0; k < dimC; k++)
{
for (m = 0; m < nRepeat; m++)
{
for (n = 0; n < dimD; n++)
{
f0_up_real = *(pFilter0_up++); f0_up_imag = *(pFilter0_up++);
f0_down_real = *(pFilter0_down++); f0_down_imag = *(pFilter0_down++);
f1_up_real = *(pFilter1_up++); f1_up_imag = *(pFilter1_up++);
f1_down_real = *(pFilter1_down++); f1_down_imag = *(pFilter1_down++);
for (s = 0; s < dimE; s++)
{
a_up_real = *(pArray_up++);
a_up_imag = *(pArray_up++);
a_down_real = *(pArray_down++);
a_down_imag = *(pArray_down++);
*(pArray_up - 2) = a_up_real * f0_up_real - a_up_imag * f0_up_imag
+ a_down_real * f0_down_real - a_down_imag * f0_down_imag;
*(pArray_up - 1) = a_up_real * f0_up_imag + a_up_imag * f0_up_real
+ a_down_real * f0_down_imag + a_down_imag * f0_down_real;
*(pArray_down - 2) = a_up_real * f1_up_real - a_up_imag * f1_up_imag
+ a_down_real * f1_down_real - a_down_imag * f1_down_imag;
*(pArray_down - 1) = a_up_real * f1_up_imag + a_up_imag * f1_up_real
+ a_down_real * f1_down_imag + a_down_imag * f1_down_real;
} // s = 0 : dimE
} // n = 0 : dimD
pFilter0_up = pFilter0_down;
pFilter0_down = pFilter0_up + filter_block_skip;
pFilter1_up = pFilter1_down;
pFilter1_down = pFilter1_up + filter_block_skip;
pArray_up = pArray_down;
pArray_down = pArray_up + array_block_skip;
} // m = 0 : nRepeat
pFilter0_up -= filter_line_skip;
pFilter0_down -= filter_line_skip;
pFilter1_up -= filter_line_skip;
pFilter1_down -= filter_line_skip;
} // k = 0 : dimC
pFilter0_up += filter_line_skip;
pFilter0_down = pFilter0_up + filter_block_skip;
pFilter1_up += filter_line_skip;
pFilter1_down = pFilter1_up + filter_block_skip;
} // j = 0 : dimB
} // i = 0 : dimA
}
}
//////////////////////////////////////////////////////////////////////////
// Filtering followed by upsampling
//////////////////////////////////////////////////////////////////////////
//
// PARAMETERS:
//
// nDims
// Number of dimensions of the input data
//
// pDims
// Dimensions of the input data
//
// pData
// pointer to the input data
//
// K
// The major dimension
//
// M
// The minor dimension
//
// Level
// The level of multiscale pyramid
void NdDirectionalFilterBank::FilterUpsampleF(int nDims, int pDims[], double *pData, int K, int M, int Level)
{
// Some validity checks
assert((nDims >= 2) && (Level >= 1) && (K >= 0) && (K < nDims) && (M >= 0) && (M < nDims) && (K != M));
assert((FilterSlice0.nx == pDims[M]) && (FilterSlice0.ny / 4 + 1 == pDims[K]));
int nChannels = (int)pow(2.0, (double)Level);
int nRepeat = nChannels / 2;
// We will divide this dimension into nChannel pieces
assert(pDims[M] % nChannels == 0);
// Record the upsampling dimension
int upsample_dim = M;
// We always make M < K
if (M > K)
{
// Note: we are transposing a complex-valued matrix
FilterSlice0.Transpose(false);
FilterSlice1.Transpose(false);
// Swap the two dimensions
int tmp;
tmp = K;
K = M;
M = tmp;
}
int dimA, dimB, dimC, dimD, dimE;
register int s, m, n, k, j, i;
dimA = 1;
for (m = 0; m < M; m++)
dimA *= pDims[m];
dimB = pDims[M];
dimC = 1;
for (m = M + 1; m < K; m++)
dimC *= pDims[m];
dimD = pDims[K];
dimE = 1;
for (m = K + 1; m < nDims; m++)
dimE *= pDims[m];
double *pFilter0_base, *pFilter1_base, *pFilter0_up, *pFilter0_down, *pFilter1_up, *pFilter1_down;
pFilter0_base = FilterSlice0.GetPointer();
pFilter1_base = FilterSlice1.GetPointer();
double f0_up_real, f0_up_imag, f0_down_real, f0_down_imag;
double f1_up_real, f1_up_imag, f1_down_real, f1_down_imag;
double *pArray_up, *pArray_down;
double a_up_real, a_up_imag, a_down_real, a_down_imag;
int filter_block_skip, filter_line_skip, filter_line_retreat;
int array_block_skip;
// There are two cases:
// (1) downsample_dim == M
// (2) downsample_dim == K
if (upsample_dim == M)
{
// preparation
dimB /= nChannels;
filter_block_skip = dimB * FilterSlice0.ny;
filter_line_retreat = 2 * dimD;
filter_line_skip = FilterSlice0.ny;
array_block_skip = dimB * 2;
for (i = M + 1; i < nDims; i++)
array_block_skip *= pDims[i];
pArray_up = pData;
pArray_down = pArray_up + array_block_skip;
for (i = 0; i < dimA; i++)
{
pFilter0_up = pFilter0_base;
pFilter0_down = pFilter0_up + filter_block_skip;
pFilter1_up = pFilter1_base;
pFilter1_down = pFilter1_up + filter_block_skip;
for (j = 0; j < nRepeat; j++)
{
for (k = 0; k < dimB; k++)
{
for (m = 0; m < dimC; m++)
{
for (n = 0; n < dimD; n++)
{
f0_up_real = *(pFilter0_up++); f0_up_imag = *(pFilter0_up++);
f0_down_real = *(pFilter0_down++); f0_down_imag = *(pFilter0_down++);
f1_up_real = *(pFilter1_up++); f1_up_imag = *(pFilter1_up++);
f1_down_real = *(pFilter1_down++); f1_down_imag = *(pFilter1_down++);
for (s = 0; s < dimE; s++)
{
a_up_real = *(pArray_up++);
a_up_imag = *(pArray_up++);
a_down_real = *(pArray_down++);
a_down_imag = *(pArray_down++);
*(pArray_up - 2) = a_up_real * f0_up_real - a_up_imag * f0_up_imag
+ a_down_real * f1_up_real - a_down_imag * f1_up_imag;
*(pArray_up - 1) = a_up_real * f0_up_imag + a_up_imag * f0_up_real
+ a_down_real * f1_up_imag + a_down_imag * f1_up_real;
*(pArray_down - 2) = a_up_real * f0_down_real - a_up_imag * f0_down_imag
+ a_down_real * f1_down_real - a_down_imag * f1_down_imag;
*(pArray_down - 1) = a_up_real * f0_down_imag + a_up_imag * f0_down_real
+ a_down_real * f1_down_imag + a_down_imag * f1_down_real;
} // s = 0 : dime
} // n = 0 : dimD
pFilter0_up -= filter_line_retreat;
pFilter0_down -= filter_line_retreat;
pFilter1_up -= filter_line_retreat;
pFilter1_down -= filter_line_retreat;
} // m = 0 : dimC
pFilter0_up += filter_line_skip;
pFilter0_down += filter_line_skip;
pFilter1_up += filter_line_skip;
pFilter1_down += filter_line_skip;
} // k = 0 : dimB
pFilter0_up = pFilter0_down;
pFilter0_down = pFilter0_up + filter_block_skip;
pFilter1_up = pFilter1_down;
pFilter1_down = pFilter1_up + filter_block_skip;
pArray_up = pArray_down;
pArray_down = pArray_up + array_block_skip;
} // j = 0 : nRepeat
} // i = 0 : dimA
} // if downsample_dim == M
else
{
// preparation
filter_block_skip = FilterSlice0.ny / nChannels;
filter_line_skip = FilterSlice0.ny;
dimD /= nChannels;
array_block_skip = dimD * 2;
for (i = K + 1; i < nDims; i++)
array_block_skip *= pDims[i];
pArray_up = pData;
pArray_down = pArray_up + array_block_skip;
for (i = 0; i < dimA; i++)
{
pFilter0_up = pFilter0_base;
pFilter0_down = pFilter0_up + filter_block_skip;
pFilter1_up = pFilter1_base;
pFilter1_down = pFilter1_up + filter_block_skip;
for (j = 0; j < dimB; j++)
{
for (k = 0; k < dimC; k++)
{
for (m = 0; m < nRepeat; m++)
{
for (n = 0; n < dimD; n++)
{
f0_up_real = *(pFilter0_up++); f0_up_imag = *(pFilter0_up++);
f0_down_real = *(pFilter0_down++); f0_down_imag = *(pFilter0_down++);
f1_up_real = *(pFilter1_up++); f1_up_imag = *(pFilter1_up++);
f1_down_real = *(pFilter1_down++); f1_down_imag = *(pFilter1_down++);
for (s = 0; s < dimE; s++)
{
a_up_real = *(pArray_up++);
a_up_imag = *(pArray_up++);
a_down_real = *(pArray_down++);
a_down_imag = *(pArray_down++);
*(pArray_up - 2) = a_up_real * f0_up_real - a_up_imag * f0_up_imag
+ a_down_real * f1_up_real - a_down_imag * f1_up_imag;
*(pArray_up - 1) = a_up_real * f0_up_imag + a_up_imag * f0_up_real
+ a_down_real * f1_up_imag + a_down_imag * f1_up_real;
*(pArray_down - 2) = a_up_real * f0_down_real - a_up_imag * f0_down_imag
+ a_down_real * f1_down_real - a_down_imag * f1_down_imag;
*(pArray_down - 1) = a_up_real * f0_down_imag + a_up_imag * f0_down_real
+ a_down_real * f1_down_imag + a_down_imag * f1_down_real;
} // s = 0 : dimE
} // n = 0 : dimD
pFilter0_up = pFilter0_down;
pFilter0_down = pFilter0_up + filter_block_skip;
pFilter1_up = pFilter1_down;
pFilter1_down = pFilter1_up + filter_block_skip;
pArray_up = pArray_down;
pArray_down = pArray_up + array_block_skip;
} // m = 0 : nRepeat
pFilter0_up -= filter_line_skip;
pFilter0_down -= filter_line_skip;
pFilter1_up -= filter_line_skip;
pFilter1_down -= filter_line_skip;
} // k = 0 : dimC
pFilter0_up += filter_line_skip;
pFilter0_down = pFilter0_up + filter_block_skip;
pFilter1_up += filter_line_skip;
pFilter1_down = pFilter1_up + filter_block_skip;
} // j = 0 : dimB
} // i = 0 : dimA
}
}
// This software is provided "as-is", without any express or implied
// warranty. In no event will the authors be held liable for any
// damages arising from the use of this software.⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -