diff --git a/modules/imgproc/include/opencv2/imgproc.hpp b/modules/imgproc/include/opencv2/imgproc.hpp index aba7821e3e..a1cfff991d 100644 --- a/modules/imgproc/include/opencv2/imgproc.hpp +++ b/modules/imgproc/include/opencv2/imgproc.hpp @@ -2835,6 +2835,22 @@ An example is shown below: */ CV_EXPORTS_W void createHanningWindow(OutputArray dst, Size winSize, int type); +/** @brief Performs the per-element division of the first Fourier spectrum by the second Fourier spectrum. + +The function cv::divSpectrums performs the per-element division of the first array by the second array. +The arrays are CCS-packed or complex matrices that are results of a real or complex Fourier transform. + +@param a first input array. +@param b second input array of the same size and type as src1 . +@param c output array of the same size and type as src1 . +@param flags operation flags; currently, the only supported flag is cv::DFT_ROWS, which indicates that +each row of src1 and src2 is an independent 1D Fourier spectrum. If you do not want to use this flag, then simply add a `0` as value. +@param conjB optional flag that conjugates the second input array before the multiplication (true) +or not (false). +*/ +CV_EXPORTS_W void divSpectrums(InputArray a, InputArray b, OutputArray c, + int flags, bool conjB = false); + //! @} imgproc_motion //! @addtogroup imgproc_misc diff --git a/modules/imgproc/src/phasecorr.cpp b/modules/imgproc/src/phasecorr.cpp index a2ba79ee30..da2aa5c8d1 100644 --- a/modules/imgproc/src/phasecorr.cpp +++ b/modules/imgproc/src/phasecorr.cpp @@ -154,7 +154,7 @@ static void magSpectrums( InputArray _src, OutputArray _dst) } } -static void divSpectrums( InputArray _srcA, InputArray _srcB, OutputArray _dst, int flags, bool conjB) +void divSpectrums( 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(); diff --git a/modules/imgproc/test/test_pc.cpp b/modules/imgproc/test/test_pc.cpp index 22c4bb5d76..edfe0701e5 100644 --- a/modules/imgproc/test/test_pc.cpp +++ b/modules/imgproc/test/test_pc.cpp @@ -120,4 +120,278 @@ TEST(Imgproc_PhaseCorrelatorTest, accuracy_1d_odd_fft) { ASSERT_NEAR(phaseShift.x, (double)xShift, 1.); } +////////////////////// DivSpectrums //////////////////////// +class CV_DivSpectrumsTest : public cvtest::ArrayTest +{ +public: + CV_DivSpectrumsTest(); +protected: + void run_func(); + void get_test_array_types_and_sizes( int, vector >& sizes, vector >& types ); + void prepare_to_validation( int test_case_idx ); + int flags; +}; + + +CV_DivSpectrumsTest::CV_DivSpectrumsTest() : flags(0) +{ + // Allocate test matrices. + test_array[INPUT].push_back(NULL); // first input DFT as a CCS-packed array or complex matrix. + test_array[INPUT].push_back(NULL); // second input DFT as a CCS-packed array or complex matrix. + test_array[OUTPUT].push_back(NULL); // output DFT as a complex matrix. + test_array[REF_OUTPUT].push_back(NULL); // reference output DFT as a complex matrix. + test_array[TEMP].push_back(NULL); // first input DFT converted to a complex matrix. + test_array[TEMP].push_back(NULL); // second input DFT converted to a complex matrix. + test_array[TEMP].push_back(NULL); // output DFT as a CCV-packed array. +} + +void CV_DivSpectrumsTest::get_test_array_types_and_sizes( int test_case_idx, vector >& sizes, vector >& types ) +{ + cvtest::ArrayTest::get_test_array_types_and_sizes(test_case_idx, sizes, types); + RNG& rng = ts->get_rng(); + + // Get the flag of the input. + const int rand_int_flags = cvtest::randInt(rng); + flags = rand_int_flags & (CV_DXT_MUL_CONJ | CV_DXT_ROWS); + + // Get input type. + const int rand_int_type = cvtest::randInt(rng); + int type; + + if (rand_int_type % 4) + { + type = CV_32FC1; + } + else if (rand_int_type % 4 == 1) + { + type = CV_32FC2; + } + else if (rand_int_type % 4 == 2) + { + type = CV_64FC1; + } + else + { + type = CV_64FC2; + } + + for( size_t i = 0; i < types.size(); i++ ) + { + for( size_t j = 0; j < types[i].size(); j++ ) + { + types[i][j] = type; + } + } + + // Inputs are CCS-packed arrays. Prepare outputs and temporary inputs as complex matrices. + if( type == CV_32FC1 || type == CV_64FC1 ) + { + types[OUTPUT][0] += 8; + types[REF_OUTPUT][0] += 8; + types[TEMP][0] += 8; + types[TEMP][1] += 8; + } +} + +/// Helper function to convert a ccs array of depth_t into a complex matrix. +template +static void convert_from_ccs_helper( const Mat& src0, const Mat& src1, Mat& dst ) +{ + const int cn = src0.channels(); + int srcstep = cn; + int dststep = 1; + + if( !dst.isContinuous() ) + dststep = (int)(dst.step/dst.elemSize()); + + if( !src0.isContinuous() ) + srcstep = (int)(src0.step/src0.elemSize1()); + + Complex *dst_data = dst.ptr >(); + const depth_t* src0_data = src0.ptr(); + const depth_t* src1_data = src1.ptr(); + dst_data->re = src0_data[0]; + dst_data->im = 0; + const int n = dst.cols + dst.rows - 1; + const int n2 = (n+1) >> 1; + + if( (n & 1) == 0 ) + { + dst_data[n2*dststep].re = src0_data[(cn == 1 ? n-1 : n2)*srcstep]; + dst_data[n2*dststep].im = 0; + } + + int delta0 = srcstep; + int delta1 = delta0 + (cn == 1 ? srcstep : 1); + + if( cn == 1 ) + srcstep *= 2; + + for( int i = 1; i < n2; i++, delta0 += srcstep, delta1 += srcstep ) + { + depth_t t0 = src0_data[delta0]; + depth_t t1 = src0_data[delta1]; + + dst_data[i*dststep].re = t0; + dst_data[i*dststep].im = t1; + + t0 = src1_data[delta0]; + t1 = -src1_data[delta1]; + + dst_data[(n-i)*dststep].re = t0; + dst_data[(n-i)*dststep].im = t1; + } +} + +/// Helper function to convert a ccs array into a complex matrix. +static void convert_from_ccs( const Mat& src0, const Mat& src1, Mat& dst, const int flags ) +{ + if( dst.rows > 1 && (dst.cols > 1 || (flags & DFT_ROWS)) ) + { + const int count = dst.rows; + const int len = dst.cols; + const bool is2d = (flags & DFT_ROWS) == 0; + for( int i = 0; i < count; i++ ) + { + const int j = !is2d || i == 0 ? i : count - i; + const Mat& src0row = src0.row(i); + const Mat& src1row = src1.row(j); + Mat dstrow = dst.row(i); + convert_from_ccs( src0row, src1row, dstrow, 0 ); + } + + if( is2d ) + { + const Mat& src0row = src0.col(0); + Mat dstrow = dst.col(0); + convert_from_ccs( src0row, src0row, dstrow, 0 ); + + if( (len & 1) == 0 ) + { + const Mat& src0row_even = src0.col(src0.cols - 1); + Mat dstrow_even = dst.col(len/2); + convert_from_ccs( src0row_even, src0row_even, dstrow_even, 0 ); + } + } + } + else + { + if( dst.depth() == CV_32F ) + { + convert_from_ccs_helper( src0, src1, dst ); + } + else + { + convert_from_ccs_helper( src0, src1, dst ); + } + } +} + +/// Helper function to compute complex number (nu_re + nu_im * i) / (de_re + de_im * i). +static std::pair divide_complex_numbers( const double nu_re, const double nu_im, + const double de_re, const double de_im, + const bool conj_de ) +{ + if ( conj_de ) + { + return divide_complex_numbers( nu_re, nu_im, de_re, -de_im, false /* conj_de */ ); + } + + const double result_de = de_re * de_re + de_im * de_im + DBL_EPSILON; + const double result_re = nu_re * de_re + nu_im * de_im; + const double result_im = nu_re * (-de_im) + nu_im * de_re; + return std::pair(result_re / result_de, result_im / result_de); +}; + +/// Helper function to divide a DFT in src1 by a DFT in src2 with depths depth_t. The DFTs are +/// complex matrices. +template +static void div_complex_helper( const Mat& src1, const Mat& src2, Mat& dst, int flags ) +{ + CV_Assert( src1.size == src2.size && src1.type() == src2.type() ); + dst.create( src1.rows, src1.cols, src1.type() ); + const int cn = src1.channels(); + int cols = src1.cols * cn; + + for( int i = 0; i < dst.rows; i++ ) + { + const depth_t *src1_data = src1.ptr(i); + const depth_t *src2_data = src2.ptr(i); + depth_t *dst_data = dst.ptr(i); + for( int j = 0; j < cols; j += 2 ) + { + std::pair result = + divide_complex_numbers( src1_data[j], src1_data[j + 1], + src2_data[j], src2_data[j + 1], + (flags & CV_DXT_MUL_CONJ) != 0 ); + dst_data[j] = (depth_t)result.first; + dst_data[j + 1] = (depth_t)result.second; + } + } +} + +/// Helper function to divide a DFT in src1 by a DFT in src2. The DFTs are complex matrices. +static void div_complex( const Mat& src1, const Mat& src2, Mat& dst, const int flags ) +{ + const int type = src1.type(); + CV_Assert( type == CV_32FC2 || type == CV_64FC2 ); + + if ( src1.depth() == CV_32F ) + { + return div_complex_helper( src1, src2, dst, flags ); + } + else + { + return div_complex_helper( src1, src2, dst, flags ); + } +} + +void CV_DivSpectrumsTest::prepare_to_validation( int /* test_case_idx */ ) +{ + Mat &src1 = test_mat[INPUT][0]; + Mat &src2 = test_mat[INPUT][1]; + Mat &ref_dst = test_mat[REF_OUTPUT][0]; + const int cn = src1.channels(); + // Inputs are CCS-packed arrays. Convert them to complex matrices and get the expected output + // as a complex matrix. + if( cn == 1 ) + { + Mat &converted_src1 = test_mat[TEMP][0]; + Mat &converted_src2 = test_mat[TEMP][1]; + convert_from_ccs( src1, src1, converted_src1, flags ); + convert_from_ccs( src2, src2, converted_src2, flags ); + div_complex( converted_src1, converted_src2, ref_dst, flags ); + } + // Inputs are complex matrices. Get the expected output as a complex matrix. + else + { + div_complex( src1, src2, ref_dst, flags ); + } +} + +void CV_DivSpectrumsTest::run_func() +{ + const Mat &src1 = test_mat[INPUT][0]; + const Mat &src2 = test_mat[INPUT][1]; + const int cn = src1.channels(); + + // Inputs are CCS-packed arrays. Get the output as a CCS-packed array and convert it to a + // complex matrix. + if ( cn == 1 ) + { + Mat &dst = test_mat[TEMP][2]; + cv::divSpectrums( src1, src2, dst, flags, (flags & CV_DXT_MUL_CONJ) != 0 ); + Mat &converted_dst = test_mat[OUTPUT][0]; + convert_from_ccs( dst, dst, converted_dst, flags ); + } + // Inputs are complex matrices. Get the output as a complex matrix. + else + { + Mat &dst = test_mat[OUTPUT][0]; + cv::divSpectrums( src1, src2, dst, flags, (flags & CV_DXT_MUL_CONJ) != 0 ); + } +} + +TEST(Imgproc_DivSpectrums, accuracy) { CV_DivSpectrumsTest test; test.safe_run(); } + }} // namespace