From 6946f510fe134b55989cf0268b6cbfad6ea6d2e4 Mon Sep 17 00:00:00 2001 From: Alexander Alekhin Date: Thu, 8 Dec 2016 13:20:38 +0300 Subject: [PATCH] mulSpectrums: refactor code --- modules/core/src/dxt.cpp | 271 ++++++++++++++++++++------------------- 1 file changed, 138 insertions(+), 133 deletions(-) diff --git a/modules/core/src/dxt.cpp b/modules/core/src/dxt.cpp index b94f65a1e9..ea587b1870 100644 --- a/modules/core/src/dxt.cpp +++ b/modules/core/src/dxt.cpp @@ -1891,13 +1891,131 @@ void cv::idft( InputArray src, OutputArray dst, int flags, int nonzero_rows ) dft( src, dst, flags | DFT_INVERSE, nonzero_rows ); } +namespace { + +#define VAL(buf, elem) (((T*)((char*)data ## buf + (step ## buf * (elem))))[0]) +#define MUL_SPECTRUMS_COL(A, B, C) \ + VAL(C, 0) = VAL(A, 0) * VAL(B, 0); \ + for (size_t j = 1; j <= rows - 2; j += 2) \ + { \ + double a_re = VAL(A, j), a_im = VAL(A, j + 1); \ + double b_re = VAL(B, j), b_im = VAL(B, j + 1); \ + if (conjB) b_im = -b_im; \ + double c_re = a_re * b_re - a_im * b_im; \ + double c_im = a_re * b_im + a_im * b_re; \ + VAL(C, j) = (T)c_re; VAL(C, j + 1) = (T)c_im; \ + } \ + if ((rows & 1) == 0) \ + VAL(C, rows-1) = VAL(A, rows-1) * VAL(B, rows-1) + +template static inline +void mulSpectrums_processCol_noinplace(const T* dataA, const T* dataB, T* dataC, size_t stepA, size_t stepB, size_t stepC, size_t rows) +{ + MUL_SPECTRUMS_COL(A, B, C); +} + +template static inline +void mulSpectrums_processCol_inplaceA(const T* dataB, T* dataAC, size_t stepB, size_t stepAC, size_t rows) +{ + MUL_SPECTRUMS_COL(AC, B, AC); +} +template static inline +void mulSpectrums_processCol(const T* dataA, const T* dataB, T* dataC, size_t stepA, size_t stepB, size_t stepC, size_t rows) +{ + if (inplaceA) + mulSpectrums_processCol_inplaceA(dataB, dataC, stepB, stepC, rows); + else + mulSpectrums_processCol_noinplace(dataA, dataB, dataC, stepA, stepB, stepC, rows); +} +#undef MUL_SPECTRUMS_COL +#undef VAL + +template static inline +void mulSpectrums_processCols(const T* dataA, const T* dataB, T* dataC, size_t stepA, size_t stepB, size_t stepC, size_t rows, size_t cols) +{ + mulSpectrums_processCol(dataA, dataB, dataC, stepA, stepB, stepC, rows); + if ((cols & 1) == 0) + { + mulSpectrums_processCol(dataA + cols - 1, dataB + cols - 1, dataC + cols - 1, stepA, stepB, stepC, rows); + } +} + +#define VAL(buf, elem) (data ## buf[(elem)]) +#define MUL_SPECTRUMS_ROW(A, B, C) \ + for (size_t j = j0; j < j1; j += 2) \ + { \ + double a_re = VAL(A, j), a_im = VAL(A, j + 1); \ + double b_re = VAL(B, j), b_im = VAL(B, j + 1); \ + if (conjB) b_im = -b_im; \ + double c_re = a_re * b_re - a_im * b_im; \ + double c_im = a_re * b_im + a_im * b_re; \ + VAL(C, j) = (T)c_re; VAL(C, j + 1) = (T)c_im; \ + } +template static inline +void mulSpectrums_processRow_noinplace(const T* dataA, const T* dataB, T* dataC, size_t j0, size_t j1) +{ + MUL_SPECTRUMS_ROW(A, B, C); +} +template static inline +void mulSpectrums_processRow_inplaceA(const T* dataB, T* dataAC, size_t j0, size_t j1) +{ + MUL_SPECTRUMS_ROW(AC, B, AC); +} +template static inline +void mulSpectrums_processRow(const T* dataA, const T* dataB, T* dataC, size_t j0, size_t j1) +{ + if (inplaceA) + mulSpectrums_processRow_inplaceA(dataB, dataC, j0, j1); + else + mulSpectrums_processRow_noinplace(dataA, dataB, dataC, j0, j1); +} +#undef MUL_SPECTRUMS_ROW +#undef VAL + +template static inline +void mulSpectrums_processRows(const T* dataA, const T* dataB, T* dataC, size_t stepA, size_t stepB, size_t stepC, size_t rows, size_t cols, size_t j0, size_t j1, bool is_1d_CN1) +{ + while (rows-- > 0) + { + if (is_1d_CN1) + dataC[0] = dataA[0]*dataB[0]; + mulSpectrums_processRow(dataA, dataB, dataC, j0, j1); + if (is_1d_CN1 && (cols & 1) == 0) + dataC[j1] = dataA[j1]*dataB[j1]; + + dataA = (const T*)(((char*)dataA) + stepA); + dataB = (const T*)(((char*)dataB) + stepB); + dataC = (T*)(((char*)dataC) + stepC); + } +} + + +template static inline +void mulSpectrums_Impl_(const T* dataA, const T* dataB, T* dataC, size_t stepA, size_t stepB, size_t stepC, size_t rows, size_t cols, size_t j0, size_t j1, bool is_1d, bool isCN1) +{ + if (!is_1d && isCN1) + { + mulSpectrums_processCols(dataA, dataB, dataC, stepA, stepB, stepC, rows, cols); + } + mulSpectrums_processRows(dataA, dataB, dataC, stepA, stepB, stepC, rows, cols, j0, j1, is_1d && isCN1); +} +template static inline +void mulSpectrums_Impl(const T* dataA, const T* dataB, T* dataC, size_t stepA, size_t stepB, size_t stepC, size_t rows, size_t cols, size_t j0, size_t j1, bool is_1d, bool isCN1) +{ + if (dataA == dataC) + mulSpectrums_Impl_(dataA, dataB, dataC, stepA, stepB, stepC, rows, cols, j0, j1, is_1d, isCN1); + else + mulSpectrums_Impl_(dataA, dataB, dataC, stepA, stepB, stepC, rows, cols, j0, j1, is_1d, isCN1); +} + +} // namespace + void cv::mulSpectrums( InputArray _srcA, InputArray _srcB, OutputArray _dst, int flags, bool conjB ) { Mat srcA = _srcA.getMat(), srcB = _srcB.getMat(); int depth = srcA.depth(), cn = srcA.channels(), type = srcA.type(); - int rows = srcA.rows, cols = srcA.cols; - int j, k; + size_t rows = srcA.rows, cols = srcA.cols; CV_Assert( type == srcB.type() && srcA.size() == srcB.size() ); CV_Assert( type == CV_32FC1 || type == CV_32FC2 || type == CV_64FC1 || type == CV_64FC2 ); @@ -1906,154 +2024,41 @@ void cv::mulSpectrums( InputArray _srcA, InputArray _srcB, Mat dst = _dst.getMat(); // correct inplace support - if (dst.data == srcA.data) - srcA = srcA.clone(); + // Case 'dst.data == srcA.data' is handled by implementation, + // because it is used frequently (filter2D, matchTemplate) if (dst.data == srcB.data) - srcB = srcB.clone(); + srcB = srcB.clone(); // workaround for B only - bool is_1d = (flags & DFT_ROWS) || (rows == 1 || (cols == 1 && - srcA.isContinuous() && srcB.isContinuous() && dst.isContinuous())); + bool is_1d = (flags & DFT_ROWS) + || (rows == 1) + || (cols == 1 && srcA.isContinuous() && srcB.isContinuous() && dst.isContinuous()); if( is_1d && !(flags & DFT_ROWS) ) cols = cols + rows - 1, rows = 1; - int ncols = cols*cn; - int j0 = cn == 1; - int j1 = ncols - (cols % 2 == 0 && cn == 1); + bool isCN1 = cn == 1; + size_t j0 = isCN1 ? 1 : 0; + size_t j1 = cols*cn - (((cols & 1) == 0 && cn == 1) ? 1 : 0); - if( depth == CV_32F ) + if (depth == CV_32F) { const float* dataA = (const float*)srcA.data; const float* dataB = (const float*)srcB.data; float* dataC = (float*)dst.data; - - size_t stepA = srcA.step/sizeof(dataA[0]); - size_t stepB = srcB.step/sizeof(dataB[0]); - size_t stepC = dst.step/sizeof(dataC[0]); - - if( !is_1d && cn == 1 ) - { - for( k = 0; k < (cols % 2 ? 1 : 2); k++ ) - { - if( k == 1 ) - dataA += cols - 1, dataB += cols - 1, dataC += cols - 1; - dataC[0] = dataA[0]*dataB[0]; - if( rows % 2 == 0 ) - dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA]*dataB[(rows-1)*stepB]; - if( !conjB ) - for( j = 1; j <= rows - 2; j += 2 ) - { - double re = (double)dataA[j*stepA]*dataB[j*stepB] - - (double)dataA[(j+1)*stepA]*dataB[(j+1)*stepB]; - double im = (double)dataA[j*stepA]*dataB[(j+1)*stepB] + - (double)dataA[(j+1)*stepA]*dataB[j*stepB]; - dataC[j*stepC] = (float)re; dataC[(j+1)*stepC] = (float)im; - } - else - for( j = 1; j <= rows - 2; j += 2 ) - { - double re = (double)dataA[j*stepA]*dataB[j*stepB] + - (double)dataA[(j+1)*stepA]*dataB[(j+1)*stepB]; - double im = (double)dataA[(j+1)*stepA]*dataB[j*stepB] - - (double)dataA[j*stepA]*dataB[(j+1)*stepB]; - dataC[j*stepC] = (float)re; dataC[(j+1)*stepC] = (float)im; - } - if( k == 1 ) - dataA -= cols - 1, dataB -= cols - 1, dataC -= cols - 1; - } - } - - for( ; rows--; dataA += stepA, dataB += stepB, dataC += stepC ) - { - if( is_1d && cn == 1 ) - { - dataC[0] = dataA[0]*dataB[0]; - if( cols % 2 == 0 ) - dataC[j1] = dataA[j1]*dataB[j1]; - } - - if( !conjB ) - for( j = j0; j < j1; j += 2 ) - { - double re = (double)dataA[j]*dataB[j] - (double)dataA[j+1]*dataB[j+1]; - double im = (double)dataA[j+1]*dataB[j] + (double)dataA[j]*dataB[j+1]; - dataC[j] = (float)re; dataC[j+1] = (float)im; - } - else - for( j = j0; j < j1; j += 2 ) - { - double re = (double)dataA[j]*dataB[j] + (double)dataA[j+1]*dataB[j+1]; - double im = (double)dataA[j+1]*dataB[j] - (double)dataA[j]*dataB[j+1]; - dataC[j] = (float)re; dataC[j+1] = (float)im; - } - } + if (!conjB) + mulSpectrums_Impl(dataA, dataB, dataC, srcA.step, srcB.step, dst.step, rows, cols, j0, j1, is_1d, isCN1); + else + mulSpectrums_Impl(dataA, dataB, dataC, srcA.step, srcB.step, dst.step, rows, cols, j0, j1, is_1d, isCN1); } else { const double* dataA = (const double*)srcA.data; const double* dataB = (const double*)srcB.data; double* dataC = (double*)dst.data; - - size_t stepA = srcA.step/sizeof(dataA[0]); - size_t stepB = srcB.step/sizeof(dataB[0]); - size_t stepC = dst.step/sizeof(dataC[0]); - - if( !is_1d && cn == 1 ) - { - for( k = 0; k < (cols % 2 ? 1 : 2); k++ ) - { - if( k == 1 ) - dataA += cols - 1, dataB += cols - 1, dataC += cols - 1; - dataC[0] = dataA[0]*dataB[0]; - if( rows % 2 == 0 ) - dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA]*dataB[(rows-1)*stepB]; - if( !conjB ) - for( j = 1; j <= rows - 2; j += 2 ) - { - double re = dataA[j*stepA]*dataB[j*stepB] - - dataA[(j+1)*stepA]*dataB[(j+1)*stepB]; - double im = dataA[j*stepA]*dataB[(j+1)*stepB] + - dataA[(j+1)*stepA]*dataB[j*stepB]; - dataC[j*stepC] = re; dataC[(j+1)*stepC] = im; - } - else - for( j = 1; j <= rows - 2; j += 2 ) - { - double re = dataA[j*stepA]*dataB[j*stepB] + - dataA[(j+1)*stepA]*dataB[(j+1)*stepB]; - double im = dataA[(j+1)*stepA]*dataB[j*stepB] - - dataA[j*stepA]*dataB[(j+1)*stepB]; - dataC[j*stepC] = re; dataC[(j+1)*stepC] = im; - } - if( k == 1 ) - dataA -= cols - 1, dataB -= cols - 1, dataC -= cols - 1; - } - } - - for( ; rows--; dataA += stepA, dataB += stepB, dataC += stepC ) - { - if( is_1d && cn == 1 ) - { - dataC[0] = dataA[0]*dataB[0]; - if( cols % 2 == 0 ) - dataC[j1] = dataA[j1]*dataB[j1]; - } - - if( !conjB ) - for( j = j0; j < j1; j += 2 ) - { - double re = dataA[j]*dataB[j] - dataA[j+1]*dataB[j+1]; - double im = dataA[j+1]*dataB[j] + dataA[j]*dataB[j+1]; - dataC[j] = re; dataC[j+1] = im; - } - else - for( j = j0; j < j1; j += 2 ) - { - double re = dataA[j]*dataB[j] + dataA[j+1]*dataB[j+1]; - double im = dataA[j+1]*dataB[j] - dataA[j]*dataB[j+1]; - dataC[j] = re; dataC[j+1] = im; - } - } + if (!conjB) + mulSpectrums_Impl(dataA, dataB, dataC, srcA.step, srcB.step, dst.step, rows, cols, j0, j1, is_1d, isCN1); + else + mulSpectrums_Impl(dataA, dataB, dataC, srcA.step, srcB.step, dst.step, rows, cols, j0, j1, is_1d, isCN1); } }