Merge pull request #20103 from thezane:make-div-spectrums-public

* Make divSpectrums public

* Add unit test

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

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();

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<vector<Size> >& sizes, vector<vector<int> >& 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<vector<Size> >& sizes, vector<vector<int> >& 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<typename depth_t>
+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<depth_t> *dst_data = dst.ptr<Complex<depth_t> >();
+    const depth_t* src0_data = src0.ptr<depth_t>();
+    const depth_t* src1_data = src1.ptr<depth_t>();
+    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<float>( src0, src1, dst );
+        }
+        else
+        {
+            convert_from_ccs_helper<double>( src0, src1, dst );
+        }
+    }
+}
+
+/// Helper function to compute complex number (nu_re + nu_im * i) / (de_re + de_im * i).
+static std::pair<double, double> 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<double, double>(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 <typename depth_t>
+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<depth_t>(i);
+        const depth_t *src2_data = src2.ptr<depth_t>(i);
+        depth_t *dst_data = dst.ptr<depth_t>(i);
+        for( int j = 0; j < cols; j += 2 )
+        {
+            std::pair<double, double> 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<float>( src1, src2, dst, flags );
+    }
+    else
+    {
+        return div_complex_helper<double>( 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