modules/imgproc/src/corner.cpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 //
  10 //                           License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
  14 // Copyright (C) 2009, Willow Garage Inc., all rights reserved.
  15 // Third party copyrights are property of their respective owners.
  16 //
  17 // Redistribution and use in source and binary forms, with or without modification,
  18 // are permitted provided that the following conditions are met:
  19 //
  20 //   * Redistribution's of source code must retain the above copyright notice,
  21 //     this list of conditions and the following disclaimer.
  22 //
  23 //   * Redistribution's in binary form must reproduce the above copyright notice,
  24 //     this list of conditions and the following disclaimer in the documentation
  25 //     and/or other materials provided with the distribution.
  26 //
  27 //   * The name of the copyright holders may not be used to endorse or promote products
  28 //     derived from this software without specific prior written permission.
  29 //
  30 // This software is provided by the copyright holders and contributors "as is" and
  31 // any express or implied warranties, including, but not limited to, the implied
  32 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  33 // In no event shall the Intel Corporation or contributors be liable for any direct,
  34 // indirect, incidental, special, exemplary, or consequential damages
  35 // (including, but not limited to, procurement of substitute goods or services;
  36 // loss of use, data, or profits; or business interruption) however caused
  37 // and on any theory of liability, whether in contract, strict liability,
  38 // or tort (including negligence or otherwise) arising in any way out of
  39 // the use of this software, even if advised of the possibility of such damage.
  40 //
  41 //M*/
  42
  43 #include "precomp.hpp"
  44 #include "opencl_kernels_imgproc.hpp"
  45
  46 namespace cv
  47 {
  48
  49 static void calcMinEigenVal( const Mat& _cov, Mat& _dst )
  50 {
  51     int i, j;
  52     Size size = _cov.size();
  53 #if CV_SSE
  54     volatile bool simd = checkHardwareSupport(CV_CPU_SSE);
  55 #endif
  56
  57     if( _cov.isContinuous() && _dst.isContinuous() )
  58     {
  59         size.width *= size.height;
  60         size.height = 1;
  61     }
  62
  63     for( i = 0; i < size.height; i++ )
  64     {
  65         const float* cov = _cov.ptr<float>(i);
  66         float* dst = _dst.ptr<float>(i);
  67         j = 0;
  68     #if CV_SSE
  69         if( simd )
  70         {
  71             __m128 half = _mm_set1_ps(0.5f);
  72             for( ; j <= size.width - 5; j += 4 )
  73             {
  74                 __m128 t0 = _mm_loadu_ps(cov + j*3); // a0 b0 c0 x
  75                 __m128 t1 = _mm_loadu_ps(cov + j*3 + 3); // a1 b1 c1 x
  76                 __m128 t2 = _mm_loadu_ps(cov + j*3 + 6); // a2 b2 c2 x
  77                 __m128 t3 = _mm_loadu_ps(cov + j*3 + 9); // a3 b3 c3 x
  78                 __m128 a, b, c, t;
  79                 t = _mm_unpacklo_ps(t0, t1); // a0 a1 b0 b1
  80                 c = _mm_unpackhi_ps(t0, t1); // c0 c1 x x
  81                 b = _mm_unpacklo_ps(t2, t3); // a2 a3 b2 b3
  82                 c = _mm_movelh_ps(c, _mm_unpackhi_ps(t2, t3)); // c0 c1 c2 c3
  83                 a = _mm_movelh_ps(t, b);
  84                 b = _mm_movehl_ps(b, t);
  85                 a = _mm_mul_ps(a, half);
  86                 c = _mm_mul_ps(c, half);
  87                 t = _mm_sub_ps(a, c);
  88                 t = _mm_add_ps(_mm_mul_ps(t, t), _mm_mul_ps(b,b));
  89                 a = _mm_sub_ps(_mm_add_ps(a, c), _mm_sqrt_ps(t));
  90                 _mm_storeu_ps(dst + j, a);
  91             }
  92         }
  93     #endif
  94         for( ; j < size.width; j++ )
  95         {
  96             float a = cov[j*3]*0.5f;
  97             float b = cov[j*3+1];
  98             float c = cov[j*3+2]*0.5f;
  99             dst[j] = (float)((a + c) - std::sqrt((a - c)*(a - c) + b*b));
 100         }
 101     }
 102 }
 103
 104
 105 static void calcHarris( const Mat& _cov, Mat& _dst, double k )
 106 {
 107     int i, j;
 108     Size size = _cov.size();
 109 #if CV_SSE
 110     volatile bool simd = checkHardwareSupport(CV_CPU_SSE);
 111 #endif
 112
 113     if( _cov.isContinuous() && _dst.isContinuous() )
 114     {
 115         size.width *= size.height;
 116         size.height = 1;
 117     }
 118
 119     for( i = 0; i < size.height; i++ )
 120     {
 121         const float* cov = _cov.ptr<float>(i);
 122         float* dst = _dst.ptr<float>(i);
 123         j = 0;
 124
 125     #if CV_SSE
 126         if( simd )
 127         {
 128             __m128 k4 = _mm_set1_ps((float)k);
 129             for( ; j <= size.width - 4; j += 4 )
 130             {
 131                 __m128 t0 = _mm_loadu_ps(cov + j*3); // a0 b0 c0 x
 132                 __m128 t1 = _mm_loadu_ps(cov + j*3 + 3); // a1 b1 c1 x
 133                 __m128 t2 = _mm_loadu_ps(cov + j*3 + 6); // a2 b2 c2 x
 134                 __m128 t3 = _mm_loadu_ps(cov + j*3 + 9); // a3 b3 c3 x
 135                 __m128 a, b, c, t;
 136                 t = _mm_unpacklo_ps(t0, t1); // a0 a1 b0 b1
 137                 c = _mm_unpackhi_ps(t0, t1); // c0 c1 x x
 138                 b = _mm_unpacklo_ps(t2, t3); // a2 a3 b2 b3
 139                 c = _mm_movelh_ps(c, _mm_unpackhi_ps(t2, t3)); // c0 c1 c2 c3
 140                 a = _mm_movelh_ps(t, b);
 141                 b = _mm_movehl_ps(b, t);
 142                 t = _mm_add_ps(a, c);
 143                 a = _mm_sub_ps(_mm_mul_ps(a, c), _mm_mul_ps(b, b));
 144                 t = _mm_mul_ps(_mm_mul_ps(k4, t), t);
 145                 a = _mm_sub_ps(a, t);
 146                 _mm_storeu_ps(dst + j, a);
 147             }
 148         }
 149     #elif CV_NEON
 150         float32x4_t v_k = vdupq_n_f32((float)k);
 151
 152         for( ; j <= size.width - 4; j += 4 )
 153         {
 154             float32x4x3_t v_src = vld3q_f32(cov + j * 3);
 155             float32x4_t v_a = v_src.val[0], v_b = v_src.val[1], v_c = v_src.val[2];
 156             float32x4_t v_ac_bb = vmlsq_f32(vmulq_f32(v_a, v_c), v_b, v_b);
 157             float32x4_t v_ac = vaddq_f32(v_a, v_c);
 158             vst1q_f32(dst + j, vmlsq_f32(v_ac_bb, v_k, vmulq_f32(v_ac, v_ac)));
 159         }
 160     #endif
 161
 162         for( ; j < size.width; j++ )
 163         {
 164             float a = cov[j*3];
 165             float b = cov[j*3+1];
 166             float c = cov[j*3+2];
 167             dst[j] = (float)(a*c - b*b - k*(a + c)*(a + c));
 168         }
 169     }
 170 }
 171
 172
 173 static void eigen2x2( const float* cov, float* dst, int n )
 174 {
 175     for( int j = 0; j < n; j++ )
 176     {
 177         double a = cov[j*3];
 178         double b = cov[j*3+1];
 179         double c = cov[j*3+2];
 180
 181         double u = (a + c)*0.5;
 182         double v = std::sqrt((a - c)*(a - c)*0.25 + b*b);
 183         double l1 = u + v;
 184         double l2 = u - v;
 185
 186         double x = b;
 187         double y = l1 - a;
 188         double e = fabs(x);
 189
 190         if( e + fabs(y) < 1e-4 )
 191         {
 192             y = b;
 193             x = l1 - c;
 194             e = fabs(x);
 195             if( e + fabs(y) < 1e-4 )
 196             {
 197                 e = 1./(e + fabs(y) + FLT_EPSILON);
 198                 x *= e, y *= e;
 199             }
 200         }
 201
 202         double d = 1./std::sqrt(x*x + y*y + DBL_EPSILON);
 203         dst[6*j] = (float)l1;
 204         dst[6*j + 2] = (float)(x*d);
 205         dst[6*j + 3] = (float)(y*d);
 206
 207         x = b;
 208         y = l2 - a;
 209         e = fabs(x);
 210
 211         if( e + fabs(y) < 1e-4 )
 212         {
 213             y = b;
 214             x = l2 - c;
 215             e = fabs(x);
 216             if( e + fabs(y) < 1e-4 )
 217             {
 218                 e = 1./(e + fabs(y) + FLT_EPSILON);
 219                 x *= e, y *= e;
 220             }
 221         }
 222
 223         d = 1./std::sqrt(x*x + y*y + DBL_EPSILON);
 224         dst[6*j + 1] = (float)l2;
 225         dst[6*j + 4] = (float)(x*d);
 226         dst[6*j + 5] = (float)(y*d);
 227     }
 228 }
 229
 230 static void calcEigenValsVecs( const Mat& _cov, Mat& _dst )
 231 {
 232     Size size = _cov.size();
 233     if( _cov.isContinuous() && _dst.isContinuous() )
 234     {
 235         size.width *= size.height;
 236         size.height = 1;
 237     }
 238
 239     for( int i = 0; i < size.height; i++ )
 240     {
 241         const float* cov = _cov.ptr<float>(i);
 242         float* dst = _dst.ptr<float>(i);
 243
 244         eigen2x2(cov, dst, size.width);
 245     }
 246 }
 247
 248
 249 enum { MINEIGENVAL=0, HARRIS=1, EIGENVALSVECS=2 };
 250
 251
 252 static void
 253 cornerEigenValsVecs( const Mat& src, Mat& eigenv, int block_size,
 254                      int aperture_size, int op_type, double k=0.,
 255                      int borderType=BORDER_DEFAULT )
 256 {
 257 #ifdef HAVE_TEGRA_OPTIMIZATION
 258     if (tegra::cornerEigenValsVecs(src, eigenv, block_size, aperture_size, op_type, k, borderType))
 259         return;
 260 #endif
 261
 262     int depth = src.depth();
 263     double scale = (double)(1 << ((aperture_size > 0 ? aperture_size : 3) - 1)) * block_size;
 264     if( aperture_size < 0 )
 265         scale *= 2.0;
 266     if( depth == CV_8U )
 267         scale *= 255.0;
 268     scale = 1.0/scale;
 269
 270     CV_Assert( src.type() == CV_8UC1 || src.type() == CV_32FC1 );
 271
 272     Mat Dx, Dy;
 273     if( aperture_size > 0 )
 274     {
 275         Sobel( src, Dx, CV_32F, 1, 0, aperture_size, scale, 0, borderType );
 276         Sobel( src, Dy, CV_32F, 0, 1, aperture_size, scale, 0, borderType );
 277     }
 278     else
 279     {
 280         Scharr( src, Dx, CV_32F, 1, 0, scale, 0, borderType );
 281         Scharr( src, Dy, CV_32F, 0, 1, scale, 0, borderType );
 282     }
 283
 284     Size size = src.size();
 285     Mat cov( size, CV_32FC3 );
 286     int i, j;
 287
 288     for( i = 0; i < size.height; i++ )
 289     {
 290         float* cov_data = cov.ptr<float>(i);
 291         const float* dxdata = Dx.ptr<float>(i);
 292         const float* dydata = Dy.ptr<float>(i);
 293
 294         for( j = 0; j < size.width; j++ )
 295         {
 296             float dx = dxdata[j];
 297             float dy = dydata[j];
 298
 299             cov_data[j*3] = dx*dx;
 300             cov_data[j*3+1] = dx*dy;
 301             cov_data[j*3+2] = dy*dy;
 302         }
 303     }
 304
 305     boxFilter(cov, cov, cov.depth(), Size(block_size, block_size),
 306         Point(-1,-1), false, borderType );
 307
 308     if( op_type == MINEIGENVAL )
 309         calcMinEigenVal( cov, eigenv );
 310     else if( op_type == HARRIS )
 311         calcHarris( cov, eigenv, k );
 312     else if( op_type == EIGENVALSVECS )
 313         calcEigenValsVecs( cov, eigenv );
 314 }
 315
 316 #ifdef HAVE_OPENCL
 317
 318 static bool extractCovData(InputArray _src, UMat & Dx, UMat & Dy, int depth,
 319                            float scale, int aperture_size, int borderType)
 320 {
 321     UMat src = _src.getUMat();
 322
 323     Size wholeSize;
 324     Point ofs;
 325     src.locateROI(wholeSize, ofs);
 326
 327     const int sobel_lsz = 16;
 328     if ((aperture_size == 3 || aperture_size == 5 || aperture_size == 7 || aperture_size == -1) &&
 329         wholeSize.height > sobel_lsz + (aperture_size >> 1) &&
 330         wholeSize.width > sobel_lsz + (aperture_size >> 1))
 331     {
 332         CV_Assert(depth == CV_8U || depth == CV_32F);
 333
 334         Dx.create(src.size(), CV_32FC1);
 335         Dy.create(src.size(), CV_32FC1);
 336
 337         size_t localsize[2] = { sobel_lsz, sobel_lsz };
 338         size_t globalsize[2] = { localsize[0] * (1 + (src.cols - 1) / localsize[0]),
 339                                  localsize[1] * (1 + (src.rows - 1) / localsize[1]) };
 340
 341         int src_offset_x = (int)((src.offset % src.step) / src.elemSize());
 342         int src_offset_y = (int)(src.offset / src.step);
 343
 344         const char * const borderTypes[] = { "BORDER_CONSTANT", "BORDER_REPLICATE", "BORDER_REFLECT",
 345                                              "BORDER_WRAP", "BORDER_REFLECT101" };
 346
 347         ocl::Kernel k(format("sobel%d", aperture_size).c_str(), ocl::imgproc::covardata_oclsrc,
 348                       cv::format("-D BLK_X=%d -D BLK_Y=%d -D %s -D SRCTYPE=%s%s",
 349                                  (int)localsize[0], (int)localsize[1], borderTypes[borderType], ocl::typeToStr(depth),
 350                                  aperture_size < 0 ? " -D SCHARR" : ""));
 351         if (k.empty())
 352             return false;
 353
 354         k.args(ocl::KernelArg::PtrReadOnly(src), (int)src.step, src_offset_x, src_offset_y,
 355                ocl::KernelArg::WriteOnlyNoSize(Dx), ocl::KernelArg::WriteOnly(Dy),
 356                wholeSize.height, wholeSize.width, scale);
 357
 358         return k.run(2, globalsize, localsize, false);
 359     }
 360     else
 361     {
 362         if (aperture_size > 0)
 363         {
 364             Sobel(_src, Dx, CV_32F, 1, 0, aperture_size, scale, 0, borderType);
 365             Sobel(_src, Dy, CV_32F, 0, 1, aperture_size, scale, 0, borderType);
 366         }
 367         else
 368         {
 369             Scharr(_src, Dx, CV_32F, 1, 0, scale, 0, borderType);
 370             Scharr(_src, Dy, CV_32F, 0, 1, scale, 0, borderType);
 371         }
 372     }
 373
 374     return true;
 375 }
 376
 377 static bool ocl_cornerMinEigenValVecs(InputArray _src, OutputArray _dst, int block_size,
 378                                       int aperture_size, double k, int borderType, int op_type)
 379 {
 380     CV_Assert(op_type == HARRIS || op_type == MINEIGENVAL);
 381
 382     if ( !(borderType == BORDER_CONSTANT || borderType == BORDER_REPLICATE ||
 383            borderType == BORDER_REFLECT || borderType == BORDER_REFLECT_101) )
 384         return false;
 385
 386     int type = _src.type(), depth = CV_MAT_DEPTH(type);
 387     if ( !(type == CV_8UC1 || type == CV_32FC1) )
 388         return false;
 389
 390     const char * const borderTypes[] = { "BORDER_CONSTANT", "BORDER_REPLICATE", "BORDER_REFLECT",
 391                                          "BORDER_WRAP", "BORDER_REFLECT101" };
 392     const char * const cornerType[] = { "CORNER_MINEIGENVAL", "CORNER_HARRIS", 0 };
 393
 394
 395     double scale = (double)(1 << ((aperture_size > 0 ? aperture_size : 3) - 1)) * block_size;
 396     if (aperture_size < 0)
 397         scale *= 2.0;
 398     if (depth == CV_8U)
 399         scale *= 255.0;
 400     scale = 1.0 / scale;
 401
 402     UMat Dx, Dy;
 403     if (!extractCovData(_src, Dx, Dy, depth, (float)scale, aperture_size, borderType))
 404         return false;
 405
 406     ocl::Kernel cornelKernel("corner", ocl::imgproc::corner_oclsrc,
 407                              format("-D anX=%d -D anY=%d -D ksX=%d -D ksY=%d -D %s -D %s",
 408                                     block_size / 2, block_size / 2, block_size, block_size,
 409                                     borderTypes[borderType], cornerType[op_type]));
 410     if (cornelKernel.empty())
 411         return false;
 412
 413     _dst.createSameSize(_src, CV_32FC1);
 414     UMat dst = _dst.getUMat();
 415
 416     cornelKernel.args(ocl::KernelArg::ReadOnly(Dx), ocl::KernelArg::ReadOnly(Dy),
 417                       ocl::KernelArg::WriteOnly(dst), (float)k);
 418
 419     size_t blockSizeX = 256, blockSizeY = 1;
 420     size_t gSize = blockSizeX - block_size / 2 * 2;
 421     size_t globalSizeX = (Dx.cols) % gSize == 0 ? Dx.cols / gSize * blockSizeX : (Dx.cols / gSize + 1) * blockSizeX;
 422     size_t rows_per_thread = 2;
 423     size_t globalSizeY = ((Dx.rows + rows_per_thread - 1) / rows_per_thread) % blockSizeY == 0 ?
 424                          ((Dx.rows + rows_per_thread - 1) / rows_per_thread) :
 425                          (((Dx.rows + rows_per_thread - 1) / rows_per_thread) / blockSizeY + 1) * blockSizeY;
 426
 427     size_t globalsize[2] = { globalSizeX, globalSizeY }, localsize[2] = { blockSizeX, blockSizeY };
 428     return cornelKernel.run(2, globalsize, localsize, false);
 429 }
 430
 431 static bool ocl_preCornerDetect( InputArray _src, OutputArray _dst, int ksize, int borderType, int depth )
 432 {
 433     UMat Dx, Dy, D2x, D2y, Dxy;
 434
 435     if (!extractCovData(_src, Dx, Dy, depth, 1, ksize, borderType))
 436         return false;
 437
 438     Sobel( _src, D2x, CV_32F, 2, 0, ksize, 1, 0, borderType );
 439     Sobel( _src, D2y, CV_32F, 0, 2, ksize, 1, 0, borderType );
 440     Sobel( _src, Dxy, CV_32F, 1, 1, ksize, 1, 0, borderType );
 441
 442     _dst.create( _src.size(), CV_32FC1 );
 443     UMat dst = _dst.getUMat();
 444
 445     double factor = 1 << (ksize - 1);
 446     if( depth == CV_8U )
 447         factor *= 255;
 448     factor = 1./(factor * factor * factor);
 449
 450     ocl::Kernel k("preCornerDetect", ocl::imgproc::precornerdetect_oclsrc);
 451     if (k.empty())
 452         return false;
 453
 454     k.args(ocl::KernelArg::ReadOnlyNoSize(Dx), ocl::KernelArg::ReadOnlyNoSize(Dy),
 455            ocl::KernelArg::ReadOnlyNoSize(D2x), ocl::KernelArg::ReadOnlyNoSize(D2y),
 456            ocl::KernelArg::ReadOnlyNoSize(Dxy), ocl::KernelArg::WriteOnly(dst), (float)factor);
 457
 458     size_t globalsize[2] = { dst.cols, dst.rows };
 459     return k.run(2, globalsize, NULL, false);
 460 }
 461
 462 #endif
 463
 464 }
 465
 466 void cv::cornerMinEigenVal( InputArray _src, OutputArray _dst, int blockSize, int ksize, int borderType )
 467 {
 468     CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
 469                ocl_cornerMinEigenValVecs(_src, _dst, blockSize, ksize, 0.0, borderType, MINEIGENVAL))
 470
 471     Mat src = _src.getMat();
 472     _dst.create( src.size(), CV_32FC1 );
 473     Mat dst = _dst.getMat();
 474 #if defined(HAVE_IPP) && (IPP_VERSION_MAJOR >= 8)
 475     typedef IppStatus (CV_STDCALL * ippiMinEigenValGetBufferSize)(IppiSize, int, int, int*);
 476     typedef IppStatus (CV_STDCALL * ippiMinEigenVal)(const void*, int, Ipp32f*, int, IppiSize, IppiKernelType, int, int, Ipp8u*);
 477     IppiKernelType kerType;
 478     int kerSize = ksize;
 479     if (ksize < 0)
 480     {
 481         kerType = ippKernelScharr;
 482         kerSize = 3;
 483     } else
 484     {
 485         kerType = ippKernelSobel;
 486     }
 487     bool isolated = (borderType & BORDER_ISOLATED) != 0;
 488     int borderTypeNI = borderType & ~BORDER_ISOLATED;
 489     if ((borderTypeNI == BORDER_REPLICATE && (!src.isSubmatrix() || isolated)) &&
 490         (kerSize == 3 || kerSize == 5) && (blockSize == 3 || blockSize == 5))
 491     {
 492         ippiMinEigenValGetBufferSize getBufferSizeFunc = 0;
 493         ippiMinEigenVal minEigenValFunc = 0;
 494         float norm_coef = 0.f;
 495
 496         if (src.type() == CV_8UC1)
 497         {
 498             getBufferSizeFunc = (ippiMinEigenValGetBufferSize) ippiMinEigenValGetBufferSize_8u32f_C1R;
 499             minEigenValFunc = (ippiMinEigenVal) ippiMinEigenVal_8u32f_C1R;
 500             norm_coef = 1.f / 255.f;
 501         } else if (src.type() == CV_32FC1)
 502         {
 503             getBufferSizeFunc = (ippiMinEigenValGetBufferSize) ippiMinEigenValGetBufferSize_32f_C1R;
 504             minEigenValFunc = (ippiMinEigenVal) ippiMinEigenVal_32f_C1R;
 505             norm_coef = 255.f;
 506         }
 507         norm_coef = kerType == ippKernelSobel ? norm_coef : norm_coef / 2.45f;
 508
 509         if (getBufferSizeFunc && minEigenValFunc)
 510         {
 511             int bufferSize;
 512             IppiSize srcRoi = { src.cols, src.rows };
 513             IppStatus ok = getBufferSizeFunc(srcRoi, kerSize, blockSize, &bufferSize);
 514             if (ok >= 0)
 515             {
 516                 AutoBuffer<uchar> buffer(bufferSize);
 517                 ok = minEigenValFunc(src.ptr(), (int) src.step, dst.ptr<Ipp32f>(), (int) dst.step, srcRoi, kerType, kerSize, blockSize, buffer);
 518                 CV_SUPPRESS_DEPRECATED_START
 519                 if (ok >= 0) ok = ippiMulC_32f_C1IR(norm_coef, dst.ptr<Ipp32f>(), (int) dst.step, srcRoi);
 520                 CV_SUPPRESS_DEPRECATED_END
 521                 if (ok >= 0)
 522                     return;
 523             }
 524             setIppErrorStatus();
 525         }
 526     }
 527 #endif
 528     cornerEigenValsVecs( src, dst, blockSize, ksize, MINEIGENVAL, 0, borderType );
 529 }
 530
 531 void cv::cornerHarris( InputArray _src, OutputArray _dst, int blockSize, int ksize, double k, int borderType )
 532 {
 533     CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
 534                ocl_cornerMinEigenValVecs(_src, _dst, blockSize, ksize, k, borderType, HARRIS))
 535
 536     Mat src = _src.getMat();
 537     _dst.create( src.size(), CV_32FC1 );
 538     Mat dst = _dst.getMat();
 539
 540 #if IPP_VERSION_X100 >= 801 && 0
 541     int type = src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
 542     int borderTypeNI = borderType & ~BORDER_ISOLATED;
 543     bool isolated = (borderType & BORDER_ISOLATED) != 0;
 544
 545     if ( (ksize == 3 || ksize == 5) && (type == CV_8UC1 || type == CV_32FC1) &&
 546         (borderTypeNI == BORDER_CONSTANT || borderTypeNI == BORDER_REPLICATE) && cn == 1 && (!src.isSubmatrix() || isolated) )
 547     {
 548         IppiSize roisize = { src.cols, src.rows };
 549         IppiMaskSize masksize = ksize == 5 ? ippMskSize5x5 : ippMskSize3x3;
 550         IppDataType datatype = type == CV_8UC1 ? ipp8u : ipp32f;
 551         Ipp32s bufsize = 0;
 552
 553         double scale = (double)(1 << ((ksize > 0 ? ksize : 3) - 1)) * blockSize;
 554         if (ksize < 0)
 555             scale *= 2.0;
 556         if (depth == CV_8U)
 557             scale *= 255.0;
 558         scale = std::pow(scale, -4.0);
 559
 560         if (ippiHarrisCornerGetBufferSize(roisize, masksize, blockSize, datatype, cn, &bufsize) >= 0)
 561         {
 562             Ipp8u * buffer = ippsMalloc_8u(bufsize);
 563             IppiDifferentialKernel filterType = ksize > 0 ? ippFilterSobel : ippFilterScharr;
 564             IppiBorderType borderTypeIpp = borderTypeNI == BORDER_CONSTANT ? ippBorderConst : ippBorderRepl;
 565             IppStatus status = (IppStatus)-1;
 566
 567             if (depth == CV_8U)
 568                 status = ippiHarrisCorner_8u32f_C1R((const Ipp8u *)src.data, (int)src.step, (Ipp32f *)dst.data, (int)dst.step, roisize,
 569                                                     filterType, masksize, blockSize, (Ipp32f)k, (Ipp32f)scale, borderTypeIpp, 0, buffer);
 570             else if (depth == CV_32F)
 571                 status = ippiHarrisCorner_32f_C1R((const Ipp32f *)src.data, (int)src.step, (Ipp32f *)dst.data, (int)dst.step, roisize,
 572                                                   filterType, masksize, blockSize, (Ipp32f)k, (Ipp32f)scale, borderTypeIpp, 0, buffer);
 573             ippsFree(buffer);
 574
 575             if (status >= 0)
 576                 return;
 577         }
 578         setIppErrorStatus();
 579     }
 580 #endif
 581
 582     cornerEigenValsVecs( src, dst, blockSize, ksize, HARRIS, k, borderType );
 583 }
 584
 585
 586 void cv::cornerEigenValsAndVecs( InputArray _src, OutputArray _dst, int blockSize, int ksize, int borderType )
 587 {
 588     Mat src = _src.getMat();
 589     Size dsz = _dst.size();
 590     int dtype = _dst.type();
 591
 592     if( dsz.height != src.rows || dsz.width*CV_MAT_CN(dtype) != src.cols*6 || CV_MAT_DEPTH(dtype) != CV_32F )
 593         _dst.create( src.size(), CV_32FC(6) );
 594     Mat dst = _dst.getMat();
 595     cornerEigenValsVecs( src, dst, blockSize, ksize, EIGENVALSVECS, 0, borderType );
 596 }
 597
 598
 599 void cv::preCornerDetect( InputArray _src, OutputArray _dst, int ksize, int borderType )
 600 {
 601     int type = _src.type();
 602     CV_Assert( type == CV_8UC1 || type == CV_32FC1 );
 603
 604     CV_OCL_RUN( _src.dims() <= 2 && _dst.isUMat(),
 605                 ocl_preCornerDetect(_src, _dst, ksize, borderType, CV_MAT_DEPTH(type)))
 606
 607     Mat Dx, Dy, D2x, D2y, Dxy, src = _src.getMat();
 608     _dst.create( src.size(), CV_32FC1 );
 609     Mat dst = _dst.getMat();
 610
 611     Sobel( src, Dx, CV_32F, 1, 0, ksize, 1, 0, borderType );
 612     Sobel( src, Dy, CV_32F, 0, 1, ksize, 1, 0, borderType );
 613     Sobel( src, D2x, CV_32F, 2, 0, ksize, 1, 0, borderType );
 614     Sobel( src, D2y, CV_32F, 0, 2, ksize, 1, 0, borderType );
 615     Sobel( src, Dxy, CV_32F, 1, 1, ksize, 1, 0, borderType );
 616
 617     double factor = 1 << (ksize - 1);
 618     if( src.depth() == CV_8U )
 619         factor *= 255;
 620     factor = 1./(factor * factor * factor);
 621 #if CV_NEON || CV_SSE2
 622     float factor_f = (float)factor;
 623 #endif
 624
 625 #if CV_SSE2
 626     volatile bool haveSSE2 = cv::checkHardwareSupport(CV_CPU_SSE2);
 627     __m128 v_factor = _mm_set1_ps(factor_f), v_m2 = _mm_set1_ps(-2.0f);
 628 #endif
 629
 630     Size size = src.size();
 631     int i, j;
 632     for( i = 0; i < size.height; i++ )
 633     {
 634         float* dstdata = dst.ptr<float>(i);
 635         const float* dxdata = Dx.ptr<float>(i);
 636         const float* dydata = Dy.ptr<float>(i);
 637         const float* d2xdata = D2x.ptr<float>(i);
 638         const float* d2ydata = D2y.ptr<float>(i);
 639         const float* dxydata = Dxy.ptr<float>(i);
 640
 641         j = 0;
 642
 643 #if CV_SSE2
 644         if (haveSSE2)
 645         {
 646             for( ; j <= size.width - 4; j += 4 )
 647             {
 648                 __m128 v_dx = _mm_loadu_ps((const float *)(dxdata + j));
 649                 __m128 v_dy = _mm_loadu_ps((const float *)(dydata + j));
 650
 651                 __m128 v_s1 = _mm_mul_ps(_mm_mul_ps(v_dx, v_dx), _mm_loadu_ps((const float *)(d2ydata + j)));
 652                 __m128 v_s2 = _mm_mul_ps(_mm_mul_ps(v_dy, v_dy), _mm_loadu_ps((const float *)(d2xdata + j)));
 653                 __m128 v_s3 = _mm_mul_ps(_mm_mul_ps(v_dx, v_dy), _mm_loadu_ps((const float *)(dxydata + j)));
 654                 v_s1 = _mm_mul_ps(v_factor, _mm_add_ps(v_s1, _mm_add_ps(v_s2, _mm_mul_ps(v_s3, v_m2))));
 655                 _mm_storeu_ps(dstdata + j, v_s1);
 656             }
 657         }
 658 #elif CV_NEON
 659         for( ; j <= size.width - 4; j += 4 )
 660         {
 661             float32x4_t v_dx = vld1q_f32(dxdata + j), v_dy = vld1q_f32(dydata + j);
 662             float32x4_t v_s = vmulq_f32(v_dx, vmulq_f32(v_dx, vld1q_f32(d2ydata + j)));
 663             v_s = vmlaq_f32(v_s, vld1q_f32(d2xdata + j), vmulq_f32(v_dy, v_dy));
 664             v_s = vmlaq_f32(v_s, vld1q_f32(dxydata + j), vmulq_n_f32(vmulq_f32(v_dy, v_dx), -2));
 665             vst1q_f32(dstdata + j, vmulq_n_f32(v_s, factor_f));
 666         }
 667 #endif
 668
 669         for( ; j < size.width; j++ )
 670         {
 671             float dx = dxdata[j];
 672             float dy = dydata[j];
 673             dstdata[j] = (float)(factor*(dx*dx*d2ydata[j] + dy*dy*d2xdata[j] - 2*dx*dy*dxydata[j]));
 674         }
 675     }
 676 }
 677
 678 CV_IMPL void
 679 cvCornerMinEigenVal( const CvArr* srcarr, CvArr* dstarr,
 680                      int block_size, int aperture_size )
 681 {
 682     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
 683
 684     CV_Assert( src.size() == dst.size() && dst.type() == CV_32FC1 );
 685     cv::cornerMinEigenVal( src, dst, block_size, aperture_size, cv::BORDER_REPLICATE );
 686 }
 687
 688 CV_IMPL void
 689 cvCornerHarris( const CvArr* srcarr, CvArr* dstarr,
 690                 int block_size, int aperture_size, double k )
 691 {
 692     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
 693
 694     CV_Assert( src.size() == dst.size() && dst.type() == CV_32FC1 );
 695     cv::cornerHarris( src, dst, block_size, aperture_size, k, cv::BORDER_REPLICATE );
 696 }
 697
 698
 699 CV_IMPL void
 700 cvCornerEigenValsAndVecs( const void* srcarr, void* dstarr,
 701                           int block_size, int aperture_size )
 702 {
 703     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
 704
 705     CV_Assert( src.rows == dst.rows && src.cols*6 == dst.cols*dst.channels() && dst.depth() == CV_32F );
 706     cv::cornerEigenValsAndVecs( src, dst, block_size, aperture_size, cv::BORDER_REPLICATE );
 707 }
 708
 709
 710 CV_IMPL void
 711 cvPreCornerDetect( const void* srcarr, void* dstarr, int aperture_size )
 712 {
 713     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
 714
 715     CV_Assert( src.size() == dst.size() && dst.type() == CV_32FC1 );
 716     cv::preCornerDetect( src, dst, aperture_size, cv::BORDER_REPLICATE );
 717 }
 718
 719 /* End of file */