Merge pull request #3384 from ilya-lavrenov:neon_new

modules/core/include/opencv2/core/base.hpp

@@ -605,6 +605,48 @@ inline uint32x4_t cv_vrndq_u32_f32(float32x4_t v)
     return vcvtq_u32_f32(vaddq_f32(v, v_05));
 }
 
+inline float32x4_t cv_vrecpq_f32(float32x4_t val)
+{
+    float32x4_t reciprocal = vrecpeq_f32(val);
+    reciprocal = vmulq_f32(vrecpsq_f32(val, reciprocal), reciprocal);
+    reciprocal = vmulq_f32(vrecpsq_f32(val, reciprocal), reciprocal);
+    return reciprocal;
+}
+
+inline float32x2_t cv_vrecp_f32(float32x2_t val)
+{
+    float32x2_t reciprocal = vrecpe_f32(val);
+    reciprocal = vmul_f32(vrecps_f32(val, reciprocal), reciprocal);
+    reciprocal = vmul_f32(vrecps_f32(val, reciprocal), reciprocal);
+    return reciprocal;
+}
+
+inline float32x4_t cv_vrsqrtq_f32(float32x4_t val)
+{
+    float32x4_t e = vrsqrteq_f32(val);
+    e = vmulq_f32(vrsqrtsq_f32(vmulq_f32(e, e), val), e);
+    e = vmulq_f32(vrsqrtsq_f32(vmulq_f32(e, e), val), e);
+    return e;
+}
+
+inline float32x2_t cv_vrsqrt_f32(float32x2_t val)
+{
+    float32x2_t e = vrsqrte_f32(val);
+    e = vmul_f32(vrsqrts_f32(vmul_f32(e, e), val), e);
+    e = vmul_f32(vrsqrts_f32(vmul_f32(e, e), val), e);
+    return e;
+}
+
+inline float32x4_t cv_vsqrtq_f32(float32x4_t val)
+{
+    return cv_vrecpq_f32(cv_vrsqrtq_f32(val));
+}
+
+inline float32x2_t cv_vsqrt_f32(float32x2_t val)
+{
+    return cv_vrecp_f32(cv_vrsqrt_f32(val));
+}
+
 #endif
 
 } // cv

modules/core/src/mathfuncs.cpp

@@ -168,6 +168,31 @@ static void FastAtan2_32f(const float *Y, const float *X, float *angle, int len,
             _mm_storeu_ps(angle + i, a);
         }
     }
+#elif CV_NEON
+    float32x4_t eps = vdupq_n_f32((float)DBL_EPSILON);
+    float32x4_t _90 = vdupq_n_f32(90.f), _180 = vdupq_n_f32(180.f), _360 = vdupq_n_f32(360.f);
+    float32x4_t z = vdupq_n_f32(0.0f), scale4 = vdupq_n_f32(scale);
+    float32x4_t p1 = vdupq_n_f32(atan2_p1), p3 = vdupq_n_f32(atan2_p3);
+    float32x4_t p5 = vdupq_n_f32(atan2_p5), p7 = vdupq_n_f32(atan2_p7);
+
+    for( ; i <= len - 4; i += 4 )
+    {
+        float32x4_t x = vld1q_f32(X + i), y = vld1q_f32(Y + i);
+        float32x4_t ax = vabsq_f32(x), ay = vabsq_f32(y);
+        float32x4_t tmin = vminq_f32(ax, ay), tmax = vmaxq_f32(ax, ay);
+        float32x4_t c = vmulq_f32(tmin, cv_vrecpq_f32(vaddq_f32(tmax, eps)));
+        float32x4_t c2 = vmulq_f32(c, c);
+        float32x4_t a = vmulq_f32(c2, p7);
+        a = vmulq_f32(vaddq_f32(a, p5), c2);
+        a = vmulq_f32(vaddq_f32(a, p3), c2);
+        a = vmulq_f32(vaddq_f32(a, p1), c);
+
+        a = vbslq_f32(vcgeq_f32(ax, ay), a, vsubq_f32(_90, a));
+        a = vbslq_f32(vcltq_f32(x, z), vsubq_f32(_180, a), a);
+        a = vbslq_f32(vcltq_f32(y, z), vsubq_f32(_360, a), a);
+
+        vst1q_f32(angle + i, vmulq_f32(a, scale4));
+    }
 #endif
 
     for( ; i < len; i++ )
@@ -268,17 +293,15 @@ static void Magnitude_32f(const float* x, const float* y, float* mag, int len)
         }
     }
 #elif CV_NEON
-    float CV_DECL_ALIGNED(16) m[4];
-
     for( ; i <= len - 4; i += 4 )
     {
         float32x4_t v_x = vld1q_f32(x + i), v_y = vld1q_f32(y + i);
-        vst1q_f32(m, vaddq_f32(vmulq_f32(v_x, v_x), vmulq_f32(v_y, v_y)));
-
-        mag[i] = std::sqrt(m[0]);
-        mag[i+1] = std::sqrt(m[1]);
-        mag[i+2] = std::sqrt(m[2]);
-        mag[i+3] = std::sqrt(m[3]);
+        vst1q_f32(mag + i, cv_vsqrtq_f32(vmlaq_f32(vmulq_f32(v_x, v_x), v_y, v_y)));
+    }
+    for( ; i <= len - 2; i += 2 )
+    {
+        float32x2_t v_x = vld1_f32(x + i), v_y = vld1_f32(y + i);
+        vst1_f32(mag + i, cv_vsqrt_f32(vmla_f32(vmul_f32(v_x, v_x), v_y, v_y)));
     }
 #endif
 
@@ -370,6 +393,12 @@ static void InvSqrt_32f(const float* src, float* dst, int len)
                 _mm_storeu_ps(dst + i, t0); _mm_storeu_ps(dst + i + 4, t1);
             }
     }
+#elif CV_NEON
+    for ( ; i <= len - 8; i += 8)
+    {
+        vst1q_f32(dst + i, cv_vrsqrtq_f32(vld1q_f32(src + i)));
+        vst1q_f32(dst + i + 4, cv_vrsqrtq_f32(vld1q_f32(src + i + 4)));
+    }
 #endif
 
     for( ; i < len; i++ )
@@ -428,6 +457,12 @@ static void Sqrt_32f(const float* src, float* dst, int len)
                 _mm_storeu_ps(dst + i, t0); _mm_storeu_ps(dst + i + 4, t1);
             }
     }
+#elif CV_NEON
+    for ( ; i <= len - 8; i += 8)
+    {
+        vst1q_f32(dst + i, cv_vsqrtq_f32(vld1q_f32(src + i)));
+        vst1q_f32(dst + i + 4, cv_vsqrtq_f32(vld1q_f32(src + i + 4)));
+    }
 #endif
 
     for( ; i < len; i++ )
@@ -869,11 +904,24 @@ void polarToCart( InputArray src1, InputArray src2,
 
                 SinCos_32f( angle, y, x, len, angleInDegrees );
                 if( mag )
-                    for( k = 0; k < len; k++ )
+                {
+                    k = 0;
+
+                    #if CV_NEON
+                    for( ; k <= len - 4; k += 4 )
+                    {
+                        float32x4_t v_m = vld1q_f32(mag + k);
+                        vst1q_f32(x + k, vmulq_f32(vld1q_f32(x + k), v_m));
+                        vst1q_f32(y + k, vmulq_f32(vld1q_f32(y + k), v_m));
+                    }
+                    #endif
+
+                    for( ; k < len; k++ )
                     {
                         float m = mag[k];
                         x[k] *= m; y[k] *= m;
                     }
+                }
             }
             else
             {
@@ -2121,12 +2169,230 @@ void log( InputArray _src, OutputArray _dst )
 *                                    P O W E R                                           *
 \****************************************************************************************/
 
+template <typename T, typename WT>
+struct iPow_SIMD
+{
+    int operator() ( const T *, T *, int, int)
+    {
+        return 0;
+    }
+};
+
+#if CV_NEON
+
+template <>
+struct iPow_SIMD<uchar, int>
+{
+    int operator() ( const uchar * src, uchar * dst, int len, int power)
+    {
+        int i = 0;
+        uint32x4_t v_1 = vdupq_n_u32(1u);
+
+        for ( ; i <= len - 8; i += 8)
+        {
+            uint32x4_t v_a1 = v_1, v_a2 = v_1;
+            uint16x8_t v_src = vmovl_u8(vld1_u8(src + i));
+            uint32x4_t v_b1 = vmovl_u16(vget_low_u16(v_src)), v_b2 = vmovl_u16(vget_high_u16(v_src));
+            int p = power;
+
+            while( p > 1 )
+            {
+                if (p & 1)
+                {
+                    v_a1 = vmulq_u32(v_a1, v_b1);
+                    v_a2 = vmulq_u32(v_a2, v_b2);
+                }
+                v_b1 = vmulq_u32(v_b1, v_b1);
+                v_b2 = vmulq_u32(v_b2, v_b2);
+                p >>= 1;
+            }
+
+            v_a1 = vmulq_u32(v_a1, v_b1);
+            v_a2 = vmulq_u32(v_a2, v_b2);
+            vst1_u8(dst + i, vqmovn_u16(vcombine_u16(vqmovn_u32(v_a1), vqmovn_u32(v_a2))));
+        }
+
+        return i;
+    }
+};
+
+template <>
+struct iPow_SIMD<schar, int>
+{
+    int operator() ( const schar * src, schar * dst, int len, int power)
+    {
+        int i = 0;
+        int32x4_t v_1 = vdupq_n_s32(1);
+
+        for ( ; i <= len - 8; i += 8)
+        {
+            int32x4_t v_a1 = v_1, v_a2 = v_1;
+            int16x8_t v_src = vmovl_s8(vld1_s8(src + i));
+            int32x4_t v_b1 = vmovl_s16(vget_low_s16(v_src)), v_b2 = vmovl_s16(vget_high_s16(v_src));
+            int p = power;
+
+            while( p > 1 )
+            {
+                if (p & 1)
+                {
+                    v_a1 = vmulq_s32(v_a1, v_b1);
+                    v_a2 = vmulq_s32(v_a2, v_b2);
+                }
+                v_b1 = vmulq_s32(v_b1, v_b1);
+                v_b2 = vmulq_s32(v_b2, v_b2);
+                p >>= 1;
+            }
+
+            v_a1 = vmulq_s32(v_a1, v_b1);
+            v_a2 = vmulq_s32(v_a2, v_b2);
+            vst1_s8(dst + i, vqmovn_s16(vcombine_s16(vqmovn_s32(v_a1), vqmovn_s32(v_a2))));
+        }
+
+        return i;
+    }
+};
+
+template <>
+struct iPow_SIMD<ushort, int>
+{
+    int operator() ( const ushort * src, ushort * dst, int len, int power)
+    {
+        int i = 0;
+        uint32x4_t v_1 = vdupq_n_u32(1u);
+
+        for ( ; i <= len - 8; i += 8)
+        {
+            uint32x4_t v_a1 = v_1, v_a2 = v_1;
+            uint16x8_t v_src = vld1q_u16(src + i);
+            uint32x4_t v_b1 = vmovl_u16(vget_low_u16(v_src)), v_b2 = vmovl_u16(vget_high_u16(v_src));
+            int p = power;
+
+            while( p > 1 )
+            {
+                if (p & 1)
+                {
+                    v_a1 = vmulq_u32(v_a1, v_b1);
+                    v_a2 = vmulq_u32(v_a2, v_b2);
+                }
+                v_b1 = vmulq_u32(v_b1, v_b1);
+                v_b2 = vmulq_u32(v_b2, v_b2);
+                p >>= 1;
+            }
+
+            v_a1 = vmulq_u32(v_a1, v_b1);
+            v_a2 = vmulq_u32(v_a2, v_b2);
+            vst1q_u16(dst + i, vcombine_u16(vqmovn_u32(v_a1), vqmovn_u32(v_a2)));
+        }
+
+        return i;
+    }
+};
+
+template <>
+struct iPow_SIMD<short, int>
+{
+    int operator() ( const short * src, short * dst, int len, int power)
+    {
+        int i = 0;
+        int32x4_t v_1 = vdupq_n_s32(1);
+
+        for ( ; i <= len - 8; i += 8)
+        {
+            int32x4_t v_a1 = v_1, v_a2 = v_1;
+            int16x8_t v_src = vld1q_s16(src + i);
+            int32x4_t v_b1 = vmovl_s16(vget_low_s16(v_src)), v_b2 = vmovl_s16(vget_high_s16(v_src));
+            int p = power;
+
+            while( p > 1 )
+            {
+                if (p & 1)
+                {
+                    v_a1 = vmulq_s32(v_a1, v_b1);
+                    v_a2 = vmulq_s32(v_a2, v_b2);
+                }
+                v_b1 = vmulq_s32(v_b1, v_b1);
+                v_b2 = vmulq_s32(v_b2, v_b2);
+                p >>= 1;
+            }
+
+            v_a1 = vmulq_s32(v_a1, v_b1);
+            v_a2 = vmulq_s32(v_a2, v_b2);
+            vst1q_s16(dst + i, vcombine_s16(vqmovn_s32(v_a1), vqmovn_s32(v_a2)));
+        }
+
+        return i;
+    }
+};
+
+
+template <>
+struct iPow_SIMD<int, int>
+{
+    int operator() ( const int * src, int * dst, int len, int power)
+    {
+        int i = 0;
+        int32x4_t v_1 = vdupq_n_s32(1);
+
+        for ( ; i <= len - 4; i += 4)
+        {
+            int32x4_t v_b = vld1q_s32(src + i), v_a = v_1;
+            int p = power;
+
+            while( p > 1 )
+            {
+                if (p & 1)
+                    v_a = vmulq_s32(v_a, v_b);
+                v_b = vmulq_s32(v_b, v_b);
+                p >>= 1;
+            }
+
+            v_a = vmulq_s32(v_a, v_b);
+            vst1q_s32(dst + i, v_a);
+        }
+
+        return i;
+    }
+};
+
+template <>
+struct iPow_SIMD<float, float>
+{
+    int operator() ( const float * src, float * dst, int len, int power)
+    {
+        int i = 0;
+        float32x4_t v_1 = vdupq_n_f32(1.0f);
+
+        for ( ; i <= len - 4; i += 4)
+        {
+            float32x4_t v_b = vld1q_f32(src + i), v_a = v_1;
+            int p = power;
+
+            while( p > 1 )
+            {
+                if (p & 1)
+                    v_a = vmulq_f32(v_a, v_b);
+                v_b = vmulq_f32(v_b, v_b);
+                p >>= 1;
+            }
+
+            v_a = vmulq_f32(v_a, v_b);
+            vst1q_f32(dst + i, v_a);
+        }
+
+        return i;
+    }
+};
+
+#endif
+
 template<typename T, typename WT>
 static void
 iPow_( const T* src, T* dst, int len, int power )
 {
-    int i;
-    for( i = 0; i < len; i++ )
+    iPow_SIMD<T, WT> vop;
+    int i = vop(src, dst, len, power);
+
+    for( ; i < len; i++ )
     {
         WT a = 1, b = src[i];
         int p = power;

modules/imgproc/src/corner.cpp

@@ -69,7 +69,7 @@ static void calcMinEigenVal( const Mat& _cov, Mat& _dst )
         if( simd )
         {
             __m128 half = _mm_set1_ps(0.5f);
-            for( ; j <= size.width - 5; j += 4 )
+            for( ; j <= size.width - 4; j += 4 )
             {
                 __m128 t0 = _mm_loadu_ps(cov + j*3); // a0 b0 c0 x
                 __m128 t1 = _mm_loadu_ps(cov + j*3 + 3); // a1 b1 c1 x
@@ -90,6 +90,19 @@ static void calcMinEigenVal( const Mat& _cov, Mat& _dst )
                 _mm_storeu_ps(dst + j, a);
             }
         }
+    #elif CV_NEON
+        float32x4_t v_half = vdupq_n_f32(0.5f);
+        for( ; j <= size.width - 4; j += 4 )
+        {
+            float32x4x3_t v_src = vld3q_f32(cov + j * 3);
+            float32x4_t v_a = vmulq_f32(v_src.val[0], v_half);
+            float32x4_t v_b = v_src.val[1];
+            float32x4_t v_c = vmulq_f32(v_src.val[2], v_half);
+
+            float32x4_t v_t = vsubq_f32(v_a, v_c);
+            v_t = vmlaq_f32(vmulq_f32(v_t, v_t), v_b, v_b);
+            vst1q_f32(dst + j, vsubq_f32(vaddq_f32(v_a, v_c), cv_vsqrtq_f32(v_t)));
+        }
     #endif
         for( ; j < size.width; j++ )
         {
@@ -290,8 +303,24 @@ cornerEigenValsVecs( const Mat& src, Mat& eigenv, int block_size,
         float* cov_data = cov.ptr<float>(i);
         const float* dxdata = Dx.ptr<float>(i);
         const float* dydata = Dy.ptr<float>(i);
+        j = 0;
 
-        for( j = 0; j < size.width; j++ )
+        #if CV_NEON
+        for( ; j <= size.width - 4; j += 4 )
+        {
+            float32x4_t v_dx = vld1q_f32(dxdata + j);
+            float32x4_t v_dy = vld1q_f32(dydata + j);
+
+            float32x4x3_t v_dst;
+            v_dst.val[0] = vmulq_f32(v_dx, v_dx);
+            v_dst.val[1] = vmulq_f32(v_dx, v_dy);
+            v_dst.val[2] = vmulq_f32(v_dy, v_dy);
+
+            vst3q_f32(cov_data + j * 3, v_dst);
+        }
+        #endif
+
+        for( ; j < size.width; j++ )
         {
             float dx = dxdata[j];
             float dy = dydata[j];

modules/imgproc/src/histogram.cpp

@@ -2316,7 +2316,16 @@ double cv::compareHist( InputArray _H1, InputArray _H2, int method )
         }
         else if( method == CV_COMP_INTERSECT )
         {
-            for( j = 0; j < len; j++ )
+            j = 0;
+            #if CV_NEON
+            float32x4_t v_result = vdupq_n_f32(0.0f);
+            for( ; j <= len - 4; j += 4 )
+                v_result = vaddq_f32(v_result, vminq_f32(vld1q_f32(h1 + j), vld1q_f32(h2 + j)));
+            float CV_DECL_ALIGNED(16) ar[4];
+            vst1q_f32(ar, v_result);
+            result += ar[0] + ar[1] + ar[2] + ar[3];
+            #endif
+            for( ; j < len; j++ )
                 result += std::min(h1[j], h2[j]);
         }
         else if( method == CV_COMP_BHATTACHARYYA )

modules/imgproc/test/test_histograms.cpp

@@ -1100,7 +1100,7 @@ int CV_CompareHistTest::validate_test_results( int /*test_case_idx*/ )
             code = cvtest::TS::FAIL_INVALID_OUTPUT;
             break;
         }
-        else if( fabs(v0 - v) > FLT_EPSILON*10*MAX(fabs(v0),0.1) )
+        else if( fabs(v0 - v) > FLT_EPSILON*14*MAX(fabs(v0),0.1) )
         {
             ts->printf( cvtest::TS::LOG, "The comparison result using the method #%d (%s)\n\tis inaccurate (=%g, should be =%g)\n",
                 i, method_name, v, v0 );

modules/imgproc/test/test_imgwarp.cpp

@@ -1548,9 +1548,28 @@ TEST(Imgproc_GetQuadSubPix, accuracy) { CV_GetQuadSubPixTest test; test.safe_run
 //////////////////////////////////////////////////////////////////////////
 
 template <typename T, typename WT>
+struct IntCast
+{
+    T operator() (WT val) const
+    {
+        return cv::saturate_cast<T>(val >> 2);
+    }
+};
+
+template <typename T, typename WT>
+struct FltCast
+{
+    T operator() (WT val) const
+    {
+        return cv::saturate_cast<T>(val * 0.25);
+    }
+};
+
+template <typename T, typename WT, int one, typename CastOp>
 void resizeArea(const cv::Mat & src, cv::Mat & dst)
 {
     int cn = src.channels();
+    CastOp castOp;
 
     for (int y = 0; y < dst.rows; ++y)
     {
@@ -1565,9 +1584,9 @@ void resizeArea(const cv::Mat & src, cv::Mat & dst)
             for (int c = 0; c < cn; ++c)
             {
                 WT sum = WT(sptr0[x1 + c]) + WT(sptr0[x1 + c + cn]);
-                sum += WT(sptr1[x1 + c]) + WT(sptr1[x1 + c + cn]) + (WT)(2);
+                sum += WT(sptr1[x1 + c]) + WT(sptr1[x1 + c + cn]) + (WT)(one);
 
-                dptr[x + c] = cv::saturate_cast<T>(sum >> 2);
+                dptr[x + c] = castOp(sum);
             }
         }
     }
@@ -1575,32 +1594,38 @@ void resizeArea(const cv::Mat & src, cv::Mat & dst)
 
 TEST(Resize, Area_half)
 {
-    const int size = 10;
-    int types[] = { CV_8UC1, CV_8UC4, CV_16UC1, CV_16UC4 };
+    const int size = 1000;
+    int types[] = { CV_8UC1, CV_8UC4, CV_16UC1, CV_16UC4, CV_16SC1, CV_16SC4, CV_32FC1, CV_32FC4 };
 
     cv::RNG rng(17);
 
     for (int i = 0, _size = sizeof(types) / sizeof(types[0]); i < _size; ++i)
     {
-        int type = types[i], depth = CV_MAT_DEPTH(type);
+        int type = types[i], depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
+        const float eps = depth <= CV_32S ? 0 : 7e-5f;
 
         SCOPED_TRACE(depth);
+        SCOPED_TRACE(cn);
 
         cv::Mat src(size, size, type), dst_actual(size >> 1, size >> 1, type),
             dst_reference(size >> 1, size >> 1, type);
 
-        rng.fill(src, cv::RNG::UNIFORM, 0, 1000, true);
+        rng.fill(src, cv::RNG::UNIFORM, -1000, 1000, true);
 
         if (depth == CV_8U)
-            resizeArea<uchar, ushort>(src, dst_reference);
+            resizeArea<uchar, ushort, 2, IntCast<uchar, ushort> >(src, dst_reference);
         else if (depth == CV_16U)
-            resizeArea<ushort, int>(src, dst_reference);
+            resizeArea<ushort, uint, 2, IntCast<ushort, uint> >(src, dst_reference);
+        else if (depth == CV_16S)
+            resizeArea<short, int, 2, IntCast<short, int> >(src, dst_reference);
+        else if (depth == CV_32F)
+            resizeArea<float, float, 0, FltCast<float, float> >(src, dst_reference);
         else
             CV_Assert(0);
 
         cv::resize(src, dst_actual, dst_actual.size(), 0, 0, cv::INTER_AREA);
 
-        ASSERT_EQ(0, cvtest::norm(dst_reference, dst_actual, cv::NORM_INF));
+        ASSERT_GE(eps, cvtest::norm(dst_reference, dst_actual, cv::NORM_INF));
     }
 }