Merge pull request #4139 from swook:spatial_gradient

modules/imgproc/include/opencv2/imgproc.hpp

@@ -1369,6 +1369,28 @@ CV_EXPORTS_W void Sobel( InputArray src, OutputArray dst, int ddepth,
                          double scale = 1, double delta = 0,
                          int borderType = BORDER_DEFAULT );
 
+/** @brief Calculates the first order image derivative in both x and y using a Sobel operator
+
+Equivalent to calling:
+
+@code
+Sobel( src, dx, CV_16SC1, 1, 0, 3 );
+Sobel( src, dy, CV_16SC1, 0, 1, 3 );
+@endcode
+
+@param src input image.
+@param dx output image with first-order derivative in x.
+@param dy output image with first-order derivative in y.
+@param ksize size of Sobel kernel. It must be 3.
+@param borderType pixel extrapolation method, see cv::BorderTypes
+
+@sa Sobel
+ */
+
+CV_EXPORTS_W void spatialGradient( InputArray src, OutputArray dx,
+                                   OutputArray dy, int ksize = 3,
+                                   int borderType = BORDER_DEFAULT );
+
 /** @brief Calculates the first x- or y- image derivative using Scharr operator.
 
 The function computes the first x- or y- spatial image derivative using the Scharr operator. The

modules/imgproc/perf/perf_spatialgradient.cpp

+#include "perf_precomp.hpp"
+
+using namespace std;
+using namespace cv;
+using namespace perf;
+using namespace testing;
+using std::tr1::make_tuple;
+using std::tr1::get;
+
+typedef std::tr1::tuple<Size, int, int> Size_Ksize_BorderType_t;
+typedef perf::TestBaseWithParam<Size_Ksize_BorderType_t> Size_Ksize_BorderType;
+
+PERF_TEST_P( Size_Ksize_BorderType, spatialGradient,
+    Combine(
+        SZ_ALL_HD,
+        Values( 3 ),
+        Values( BORDER_DEFAULT, BORDER_REPLICATE )
+    )
+)
+{
+    Size size = std::tr1::get<0>(GetParam());
+    int ksize = std::tr1::get<1>(GetParam());
+    int borderType = std::tr1::get<2>(GetParam());
+
+    Mat src(size, CV_8UC1);
+    Mat dx(size, CV_16SC1);
+    Mat dy(size, CV_16SC1);
+
+    declare.in(src, WARMUP_RNG).out(dx, dy);
+
+    TEST_CYCLE() spatialGradient(src, dx, dy, ksize, borderType);
+
+    SANITY_CHECK(dx);
+    SANITY_CHECK(dy);
+}

modules/imgproc/src/precomp.hpp

@@ -49,6 +49,7 @@
 #include "opencv2/imgproc/imgproc_c.h"
 #include "opencv2/core/private.hpp"
 #include "opencv2/core/ocl.hpp"
+#include "opencv2/hal.hpp"
 
 #include <math.h>
 #include <assert.h>

modules/imgproc/src/spatialgradient.cpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "precomp.hpp"
+#include "opencv2/hal/intrin.hpp"
+
+#include <iostream>
+namespace cv
+{
+
+/* NOTE:
+ *
+ * Sobel-x: -1  0  1
+ *          -2  0  2
+ *          -1  0  1
+ *
+ * Sobel-y: -1 -2 -1
+ *           0  0  0
+ *           1  2  1
+ */
+template <typename T>
+static inline void spatialGradientKernel( T& vx, T& vy,
+                                          const T& v00, const T& v01, const T& v02,
+                                          const T& v10,               const T& v12,
+                                          const T& v20, const T& v21, const T& v22 )
+{
+    // vx = (v22 - v00) + (v02 - v20) + 2 * (v12 - v10)
+    // vy = (v22 - v00) + (v20 - v02) + 2 * (v21 - v01)
+
+    T tmp_add = v22 - v00,
+      tmp_sub = v02 - v20,
+      tmp_x   = v12 - v10,
+      tmp_y   = v21 - v01;
+
+    vx = tmp_add + tmp_sub + tmp_x + tmp_x;
+    vy = tmp_add - tmp_sub + tmp_y + tmp_y;
+}
+
+void spatialGradient( InputArray _src, OutputArray _dx, OutputArray _dy,
+                      int ksize, int borderType )
+{
+
+    // Prepare InputArray src
+    Mat src = _src.getMat();
+    CV_Assert( !src.empty() );
+    CV_Assert( src.type() == CV_8UC1 );
+    CV_Assert( borderType == BORDER_DEFAULT || borderType == BORDER_REPLICATE );
+
+    // Prepare OutputArrays dx, dy
+    _dx.create( src.size(), CV_16SC1 );
+    _dy.create( src.size(), CV_16SC1 );
+    Mat dx = _dx.getMat(),
+        dy = _dy.getMat();
+
+    // TODO: Allow for other kernel sizes
+    CV_Assert(ksize == 3);
+
+    // Get dimensions
+    const int H = src.rows,
+              W = src.cols;
+
+    // Row, column indices
+    int i = 0,
+        j = 0;
+
+    // Handle border types
+    int i_top    = 0,     // Case for H == 1 && W == 1 && BORDER_REPLICATE
+        i_bottom = H - 1,
+        j_offl   = 0,     // j offset from 0th   pixel to reach -1st pixel
+        j_offr   = 0;     // j offset from W-1th pixel to reach Wth  pixel
+
+    if ( borderType == BORDER_DEFAULT ) // Equiv. to BORDER_REFLECT_101
+    {
+        if ( H > 1 )
+        {
+            i_top    = 1;
+            i_bottom = H - 2;
+        }
+        if ( W > 1 )
+        {
+            j_offl = 1;
+            j_offr = -1;
+        }
+    }
+
+    // Pointer to row vectors
+    uchar *p_src, *c_src, *n_src; // previous, current, next row
+    short *c_dx,  *c_dy;
+
+    int i_start = 0;
+    int j_start = 0;
+#if CV_SIMD128 && CV_SSE2
+    uchar *m_src;
+    short *n_dx, *n_dy;
+
+    // Characters in variable names have the following meanings:
+    // u: unsigned char
+    // s: signed int
+    //
+    // [row][column]
+    // m: offset -1
+    // n: offset  0
+    // p: offset  1
+    // Example: umn is offset -1 in row and offset 0 in column
+    for ( i = 0; i < H - 1; i += 2 )
+    {
+        if   ( i == 0 ) p_src = src.ptr<uchar>(i_top);
+        else            p_src = src.ptr<uchar>(i-1);
+
+        c_src = src.ptr<uchar>(i);
+        n_src = src.ptr<uchar>(i+1);
+
+        if ( i == H - 2 ) m_src = src.ptr<uchar>(i_bottom);
+        else              m_src = src.ptr<uchar>(i+2);
+
+        c_dx = dx.ptr<short>(i);
+        c_dy = dy.ptr<short>(i);
+        n_dx = dx.ptr<short>(i+1);
+        n_dy = dy.ptr<short>(i+1);
+
+        v_uint8x16 v_select_m = v_uint8x16(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                           0, 0, 0, 0xFF);
+
+        // Process rest of columns 16-column chunks at a time
+        for ( j = 1; j < W - 16; j += 16 )
+        {
+            // Load top row for 3x3 Sobel filter
+            v_uint8x16 v_um = v_load(&p_src[j-1]);
+            v_uint8x16 v_up = v_load(&p_src[j+1]);
+            // TODO: Replace _mm_slli_si128 with hal method
+            v_uint8x16 v_un = v_select(v_select_m, v_uint8x16(_mm_slli_si128(v_up.val, 1)),
+                                                   v_uint8x16(_mm_srli_si128(v_um.val, 1)));
+            v_uint16x8 v_um1, v_um2, v_un1, v_un2, v_up1, v_up2;
+            v_expand(v_um, v_um1, v_um2);
+            v_expand(v_un, v_un1, v_un2);
+            v_expand(v_up, v_up1, v_up2);
+            v_int16x8 v_s1m1 = v_reinterpret_as_s16(v_um1);
+            v_int16x8 v_s1m2 = v_reinterpret_as_s16(v_um2);
+            v_int16x8 v_s1n1 = v_reinterpret_as_s16(v_un1);
+            v_int16x8 v_s1n2 = v_reinterpret_as_s16(v_un2);
+            v_int16x8 v_s1p1 = v_reinterpret_as_s16(v_up1);
+            v_int16x8 v_s1p2 = v_reinterpret_as_s16(v_up2);
+
+            // Load second row for 3x3 Sobel filter
+            v_um = v_load(&c_src[j-1]);
+            v_up = v_load(&c_src[j+1]);
+            // TODO: Replace _mm_slli_si128 with hal method
+            v_un = v_select(v_select_m, v_uint8x16(_mm_slli_si128(v_up.val, 1)),
+                                        v_uint8x16(_mm_srli_si128(v_um.val, 1)));
+            v_expand(v_um, v_um1, v_um2);
+            v_expand(v_un, v_un1, v_un2);
+            v_expand(v_up, v_up1, v_up2);
+            v_int16x8 v_s2m1 = v_reinterpret_as_s16(v_um1);
+            v_int16x8 v_s2m2 = v_reinterpret_as_s16(v_um2);
+            v_int16x8 v_s2n1 = v_reinterpret_as_s16(v_un1);
+            v_int16x8 v_s2n2 = v_reinterpret_as_s16(v_un2);
+            v_int16x8 v_s2p1 = v_reinterpret_as_s16(v_up1);
+            v_int16x8 v_s2p2 = v_reinterpret_as_s16(v_up2);
+
+            // Load third row for 3x3 Sobel filter
+            v_um = v_load(&n_src[j-1]);
+            v_up = v_load(&n_src[j+1]);
+            // TODO: Replace _mm_slli_si128 with hal method
+            v_un = v_select(v_select_m, v_uint8x16(_mm_slli_si128(v_up.val, 1)),
+                                        v_uint8x16(_mm_srli_si128(v_um.val, 1)));
+            v_expand(v_um, v_um1, v_um2);
+            v_expand(v_un, v_un1, v_un2);
+            v_expand(v_up, v_up1, v_up2);
+            v_int16x8 v_s3m1 = v_reinterpret_as_s16(v_um1);
+            v_int16x8 v_s3m2 = v_reinterpret_as_s16(v_um2);
+            v_int16x8 v_s3n1 = v_reinterpret_as_s16(v_un1);
+            v_int16x8 v_s3n2 = v_reinterpret_as_s16(v_un2);
+            v_int16x8 v_s3p1 = v_reinterpret_as_s16(v_up1);
+            v_int16x8 v_s3p2 = v_reinterpret_as_s16(v_up2);
+
+            // dx & dy for rows 1, 2, 3
+            v_int16x8 v_sdx1, v_sdy1;
+            spatialGradientKernel<v_int16x8>( v_sdx1, v_sdy1,
+                                              v_s1m1, v_s1n1, v_s1p1,
+                                              v_s2m1,         v_s2p1,
+                                              v_s3m1, v_s3n1, v_s3p1 );
+
+            v_int16x8 v_sdx2, v_sdy2;
+            spatialGradientKernel<v_int16x8>( v_sdx2, v_sdy2,
+                                              v_s1m2, v_s1n2, v_s1p2,
+                                              v_s2m2,         v_s2p2,
+                                              v_s3m2, v_s3n2, v_s3p2 );
+
+            // Store
+            v_store(&c_dx[j],   v_sdx1);
+            v_store(&c_dx[j+8], v_sdx2);
+            v_store(&c_dy[j],   v_sdy1);
+            v_store(&c_dy[j+8], v_sdy2);
+
+            // Load fourth row for 3x3 Sobel filter
+            v_um = v_load(&m_src[j-1]);
+            v_up = v_load(&m_src[j+1]);
+            // TODO: Replace _mm_slli_si128 with hal method
+            v_un = v_select(v_select_m, v_uint8x16(_mm_slli_si128(v_up.val, 1)),
+                                        v_uint8x16(_mm_srli_si128(v_um.val, 1)));
+            v_expand(v_um, v_um1, v_um2);
+            v_expand(v_un, v_un1, v_un2);
+            v_expand(v_up, v_up1, v_up2);
+            v_int16x8 v_s4m1 = v_reinterpret_as_s16(v_um1);
+            v_int16x8 v_s4m2 = v_reinterpret_as_s16(v_um2);
+            v_int16x8 v_s4n1 = v_reinterpret_as_s16(v_un1);
+            v_int16x8 v_s4n2 = v_reinterpret_as_s16(v_un2);
+            v_int16x8 v_s4p1 = v_reinterpret_as_s16(v_up1);
+            v_int16x8 v_s4p2 = v_reinterpret_as_s16(v_up2);
+
+            // dx & dy for rows 2, 3, 4
+            spatialGradientKernel<v_int16x8>( v_sdx1, v_sdy1,
+                                              v_s2m1, v_s2n1, v_s2p1,
+                                              v_s3m1,         v_s3p1,
+                                              v_s4m1, v_s4n1, v_s4p1 );
+
+            spatialGradientKernel<v_int16x8>( v_sdx2, v_sdy2,
+                                              v_s2m2, v_s2n2, v_s2p2,
+                                              v_s3m2,         v_s3p2,
+                                              v_s4m2, v_s4n2, v_s4p2 );
+
+            // Store
+            v_store(&n_dx[j],   v_sdx1);
+            v_store(&n_dx[j+8], v_sdx2);
+            v_store(&n_dy[j],   v_sdy1);
+            v_store(&n_dy[j+8], v_sdy2);
+        }
+    }
+    i_start = i;
+    j_start = j;
+#endif
+    int j_p, j_n;
+    uchar v00, v01, v02, v10, v11, v12, v20, v21, v22;
+    for ( i = 0; i < H; i++ )
+    {
+        if   ( i == 0 ) p_src = src.ptr<uchar>(i_top);
+        else            p_src = src.ptr<uchar>(i-1);
+
+        c_src = src.ptr<uchar>(i);
+
+        if ( i == H - 1 ) n_src = src.ptr<uchar>(i_bottom);
+        else              n_src = src.ptr<uchar>(i+1);
+
+        c_dx = dx.ptr<short>(i);
+        c_dy = dy.ptr<short>(i);
+
+        // Process left-most column
+        j = 0;
+        j_p = j + j_offl;
+        j_n = 1;
+        if ( j_n >= W ) j_n = j + j_offr;
+        v00 = p_src[j_p]; v01 = p_src[j]; v02 = p_src[j_n];
+        v10 = c_src[j_p]; v11 = c_src[j]; v12 = c_src[j_n];
+        v20 = n_src[j_p]; v21 = n_src[j]; v22 = n_src[j_n];
+        spatialGradientKernel<short>( c_dx[0], c_dy[0], v00, v01, v02, v10,
+                                      v12, v20, v21, v22 );
+        v00 = v01; v10 = v11; v20 = v21;
+        v01 = v02; v11 = v12; v21 = v22;
+
+        // Process middle columns
+        j = i >= i_start ? 1 : j_start;
+        j_p = j - 1;
+        v00 = p_src[j_p]; v01 = p_src[j];
+        v10 = c_src[j_p]; v11 = c_src[j];
+        v20 = n_src[j_p]; v21 = n_src[j];
+
+        for ( ; j < W - 1; j++ )
+        {
+            // Get values for next column
+            j_n = j + 1; v02 = p_src[j_n]; v12 = c_src[j_n]; v22 = n_src[j_n];
+            spatialGradientKernel<short>( c_dx[j], c_dy[j], v00, v01, v02, v10,
+                                          v12, v20, v21, v22 );
+
+            // Move values back one column for next iteration
+            v00 = v01; v10 = v11; v20 = v21;
+            v01 = v02; v11 = v12; v21 = v22;
+        }
+
+        // Process right-most column
+        if ( j < W )
+        {
+            j_n = j + j_offr; v02 = p_src[j_n]; v12 = c_src[j_n]; v22 = n_src[j_n];
+            spatialGradientKernel<short>( c_dx[j], c_dy[j], v00, v01, v02, v10,
+                                          v12, v20, v21, v22 );
+        }
+    }
+
+}
+
+}

modules/imgproc/test/test_filter.cpp

@@ -552,6 +552,68 @@ void CV_SobelTest::prepare_to_validation( int /*test_case_idx*/ )
 }
 
 
+/////////////// spatialGradient ///////////////
+
+class CV_SpatialGradientTest : public CV_DerivBaseTest
+{
+public:
+    CV_SpatialGradientTest();
+
+protected:
+    void prepare_to_validation( int test_case_idx );
+    void run_func();
+    void get_test_array_types_and_sizes( int test_case_idx,
+        vector<vector<Size> >& sizes, vector<vector<int> >& types );
+    int ksize;
+};
+
+CV_SpatialGradientTest::CV_SpatialGradientTest() {
+    test_array[OUTPUT].push_back(NULL);
+    test_array[REF_OUTPUT].push_back(NULL);
+    inplace = false;
+}
+
+
+void CV_SpatialGradientTest::get_test_array_types_and_sizes( int test_case_idx,
+                                                             vector<vector<Size> >& sizes,
+                                                             vector<vector<int> >& types )
+{
+    CV_DerivBaseTest::get_test_array_types_and_sizes( test_case_idx, sizes, types );
+
+    sizes[OUTPUT][1] = sizes[REF_OUTPUT][1] = sizes[OUTPUT][0];
+
+    // Inputs are only CV_8UC1 for now
+    types[INPUT][0] = CV_8UC1;
+
+    // Outputs are only CV_16SC1 for now
+    types[OUTPUT][0] = types[OUTPUT][1] = types[REF_OUTPUT][0]
+                     = types[REF_OUTPUT][1] = CV_16SC1;
+
+    ksize = 3;
+    border = BORDER_DEFAULT; // TODO: Add BORDER_REPLICATE
+}
+
+
+void CV_SpatialGradientTest::run_func()
+{
+    spatialGradient( test_mat[INPUT][0], test_mat[OUTPUT][0],
+                     test_mat[OUTPUT][1], ksize, border );
+}
+
+void CV_SpatialGradientTest::prepare_to_validation( int /*test_case_idx*/ )
+{
+    int dx, dy;
+
+    dx = 1; dy = 0;
+    Sobel( test_mat[INPUT][0], test_mat[REF_OUTPUT][0], CV_16SC1, dx, dy, ksize,
+           1, 0, border );
+
+    dx = 0; dy = 1;
+    Sobel( test_mat[INPUT][0], test_mat[REF_OUTPUT][1], CV_16SC1, dx, dy, ksize,
+           1, 0, border );
+}
+
+
 /////////////// laplace ///////////////
 
 class CV_LaplaceTest : public CV_DerivBaseTest
@@ -1773,6 +1835,7 @@ TEST(Imgproc_Dilate, accuracy) { CV_DilateTest test; test.safe_run(); }
 TEST(Imgproc_MorphologyEx, accuracy) { CV_MorphExTest test; test.safe_run(); }
 TEST(Imgproc_Filter2D, accuracy) { CV_FilterTest test; test.safe_run(); }
 TEST(Imgproc_Sobel, accuracy) { CV_SobelTest test; test.safe_run(); }
+TEST(Imgproc_SpatialGradient, accuracy) { CV_SpatialGradientTest test; test.safe_run(); }
 TEST(Imgproc_Laplace, accuracy) { CV_LaplaceTest test; test.safe_run(); }
 TEST(Imgproc_Blur, accuracy) { CV_BlurTest test; test.safe_run(); }
 TEST(Imgproc_GaussianBlur, accuracy) { CV_GaussianBlurTest test; test.safe_run(); }