modules/imgproc/src/imgwarp.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 /* ////////////////////////////////////////////////////////////////////
  44 //
  45 //  Geometrical transforms on images and matrices: rotation, zoom etc.
  46 //
  47 // */
  48
  49 #include "precomp.hpp"
  50 #include "opencl_kernels.hpp"
  51
  52 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
  53 static IppStatus sts = ippInit();
  54 #endif
  55
  56 namespace cv
  57 {
  58 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR*100 + IPP_VERSION_MINOR >= 701)
  59     typedef IppStatus (CV_STDCALL* ippiResizeFunc)(const void*, int, const void*, int, IppiPoint, IppiSize, IppiBorderType, void*, void*, Ipp8u*);
  60 #endif
  61
  62 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
  63     typedef IppStatus (CV_STDCALL* ippiSetFunc)(const void*, void *, int, IppiSize);
  64     typedef IppStatus (CV_STDCALL* ippiWarpPerspectiveBackFunc)(const void*, IppiSize, int, IppiRect, void *, int, IppiRect, double [3][3], int);
  65     typedef IppStatus (CV_STDCALL* ippiWarpAffineBackFunc)(const void*, IppiSize, int, IppiRect, void *, int, IppiRect, double [2][3], int);
  66
  67     template <int channels, typename Type>
  68     bool IPPSetSimple(cv::Scalar value, void *dataPointer, int step, IppiSize &size, ippiSetFunc func)
  69     {
  70         Type values[channels];
  71         for( int i = 0; i < channels; i++ )
  72             values[i] = (Type)value[i];
  73         return func(values, dataPointer, step, size) >= 0;
  74     }
  75
  76     bool IPPSet(const cv::Scalar &value, void *dataPointer, int step, IppiSize &size, int channels, int depth)
  77     {
  78         if( channels == 1 )
  79         {
  80             switch( depth )
  81             {
  82             case CV_8U:
  83                 return ippiSet_8u_C1R((Ipp8u)value[0], (Ipp8u *)dataPointer, step, size) >= 0;
  84             case CV_16U:
  85                 return ippiSet_16u_C1R((Ipp16u)value[0], (Ipp16u *)dataPointer, step, size) >= 0;
  86             case CV_32F:
  87                 return ippiSet_32f_C1R((Ipp32f)value[0], (Ipp32f *)dataPointer, step, size) >= 0;
  88             }
  89         }
  90         else
  91         {
  92             if( channels == 3 )
  93             {
  94                 switch( depth )
  95                 {
  96                 case CV_8U:
  97                     return IPPSetSimple<3, Ipp8u>(value, dataPointer, step, size, (ippiSetFunc)ippiSet_8u_C3R);
  98                 case CV_16U:
  99                     return IPPSetSimple<3, Ipp16u>(value, dataPointer, step, size, (ippiSetFunc)ippiSet_16u_C3R);
 100                 case CV_32F:
 101                     return IPPSetSimple<3, Ipp32f>(value, dataPointer, step, size, (ippiSetFunc)ippiSet_32f_C3R);
 102                 }
 103             }
 104             else if( channels == 4 )
 105             {
 106                 switch( depth )
 107                 {
 108                 case CV_8U:
 109                     return IPPSetSimple<4, Ipp8u>(value, dataPointer, step, size, (ippiSetFunc)ippiSet_8u_C4R);
 110                 case CV_16U:
 111                     return IPPSetSimple<4, Ipp16u>(value, dataPointer, step, size, (ippiSetFunc)ippiSet_16u_C4R);
 112                 case CV_32F:
 113                     return IPPSetSimple<4, Ipp32f>(value, dataPointer, step, size, (ippiSetFunc)ippiSet_32f_C4R);
 114                 }
 115             }
 116         }
 117         return false;
 118     }
 119 #endif
 120
 121 /************** interpolation formulas and tables ***************/
 122
 123 const int INTER_RESIZE_COEF_BITS=11;
 124 const int INTER_RESIZE_COEF_SCALE=1 << INTER_RESIZE_COEF_BITS;
 125
 126 const int INTER_REMAP_COEF_BITS=15;
 127 const int INTER_REMAP_COEF_SCALE=1 << INTER_REMAP_COEF_BITS;
 128
 129 static uchar NNDeltaTab_i[INTER_TAB_SIZE2][2];
 130
 131 static float BilinearTab_f[INTER_TAB_SIZE2][2][2];
 132 static short BilinearTab_i[INTER_TAB_SIZE2][2][2];
 133
 134 #if CV_SSE2
 135 static short BilinearTab_iC4_buf[INTER_TAB_SIZE2+2][2][8];
 136 static short (*BilinearTab_iC4)[2][8] = (short (*)[2][8])alignPtr(BilinearTab_iC4_buf, 16);
 137 #endif
 138
 139 static float BicubicTab_f[INTER_TAB_SIZE2][4][4];
 140 static short BicubicTab_i[INTER_TAB_SIZE2][4][4];
 141
 142 static float Lanczos4Tab_f[INTER_TAB_SIZE2][8][8];
 143 static short Lanczos4Tab_i[INTER_TAB_SIZE2][8][8];
 144
 145 static inline void interpolateLinear( float x, float* coeffs )
 146 {
 147     coeffs[0] = 1.f - x;
 148     coeffs[1] = x;
 149 }
 150
 151 static inline void interpolateCubic( float x, float* coeffs )
 152 {
 153     const float A = -0.75f;
 154
 155     coeffs[0] = ((A*(x + 1) - 5*A)*(x + 1) + 8*A)*(x + 1) - 4*A;
 156     coeffs[1] = ((A + 2)*x - (A + 3))*x*x + 1;
 157     coeffs[2] = ((A + 2)*(1 - x) - (A + 3))*(1 - x)*(1 - x) + 1;
 158     coeffs[3] = 1.f - coeffs[0] - coeffs[1] - coeffs[2];
 159 }
 160
 161 static inline void interpolateLanczos4( float x, float* coeffs )
 162 {
 163     static const double s45 = 0.70710678118654752440084436210485;
 164     static const double cs[][2]=
 165     {{1, 0}, {-s45, -s45}, {0, 1}, {s45, -s45}, {-1, 0}, {s45, s45}, {0, -1}, {-s45, s45}};
 166
 167     if( x < FLT_EPSILON )
 168     {
 169         for( int i = 0; i < 8; i++ )
 170             coeffs[i] = 0;
 171         coeffs[3] = 1;
 172         return;
 173     }
 174
 175     float sum = 0;
 176     double y0=-(x+3)*CV_PI*0.25, s0 = sin(y0), c0=cos(y0);
 177     for(int i = 0; i < 8; i++ )
 178     {
 179         double y = -(x+3-i)*CV_PI*0.25;
 180         coeffs[i] = (float)((cs[i][0]*s0 + cs[i][1]*c0)/(y*y));
 181         sum += coeffs[i];
 182     }
 183
 184     sum = 1.f/sum;
 185     for(int i = 0; i < 8; i++ )
 186         coeffs[i] *= sum;
 187 }
 188
 189 static void initInterTab1D(int method, float* tab, int tabsz)
 190 {
 191     float scale = 1.f/tabsz;
 192     if( method == INTER_LINEAR )
 193     {
 194         for( int i = 0; i < tabsz; i++, tab += 2 )
 195             interpolateLinear( i*scale, tab );
 196     }
 197     else if( method == INTER_CUBIC )
 198     {
 199         for( int i = 0; i < tabsz; i++, tab += 4 )
 200             interpolateCubic( i*scale, tab );
 201     }
 202     else if( method == INTER_LANCZOS4 )
 203     {
 204         for( int i = 0; i < tabsz; i++, tab += 8 )
 205             interpolateLanczos4( i*scale, tab );
 206     }
 207     else
 208         CV_Error( CV_StsBadArg, "Unknown interpolation method" );
 209 }
 210
 211
 212 static const void* initInterTab2D( int method, bool fixpt )
 213 {
 214     static bool inittab[INTER_MAX+1] = {false};
 215     float* tab = 0;
 216     short* itab = 0;
 217     int ksize = 0;
 218     if( method == INTER_LINEAR )
 219         tab = BilinearTab_f[0][0], itab = BilinearTab_i[0][0], ksize=2;
 220     else if( method == INTER_CUBIC )
 221         tab = BicubicTab_f[0][0], itab = BicubicTab_i[0][0], ksize=4;
 222     else if( method == INTER_LANCZOS4 )
 223         tab = Lanczos4Tab_f[0][0], itab = Lanczos4Tab_i[0][0], ksize=8;
 224     else
 225         CV_Error( CV_StsBadArg, "Unknown/unsupported interpolation type" );
 226
 227     if( !inittab[method] )
 228     {
 229         AutoBuffer<float> _tab(8*INTER_TAB_SIZE);
 230         int i, j, k1, k2;
 231         initInterTab1D(method, _tab, INTER_TAB_SIZE);
 232         for( i = 0; i < INTER_TAB_SIZE; i++ )
 233             for( j = 0; j < INTER_TAB_SIZE; j++, tab += ksize*ksize, itab += ksize*ksize )
 234             {
 235                 int isum = 0;
 236                 NNDeltaTab_i[i*INTER_TAB_SIZE+j][0] = j < INTER_TAB_SIZE/2;
 237                 NNDeltaTab_i[i*INTER_TAB_SIZE+j][1] = i < INTER_TAB_SIZE/2;
 238
 239                 for( k1 = 0; k1 < ksize; k1++ )
 240                 {
 241                     float vy = _tab[i*ksize + k1];
 242                     for( k2 = 0; k2 < ksize; k2++ )
 243                     {
 244                         float v = vy*_tab[j*ksize + k2];
 245                         tab[k1*ksize + k2] = v;
 246                         isum += itab[k1*ksize + k2] = saturate_cast<short>(v*INTER_REMAP_COEF_SCALE);
 247                     }
 248                 }
 249
 250                 if( isum != INTER_REMAP_COEF_SCALE )
 251                 {
 252                     int diff = isum - INTER_REMAP_COEF_SCALE;
 253                     int ksize2 = ksize/2, Mk1=ksize2, Mk2=ksize2, mk1=ksize2, mk2=ksize2;
 254                     for( k1 = ksize2; k1 < ksize2+2; k1++ )
 255                         for( k2 = ksize2; k2 < ksize2+2; k2++ )
 256                         {
 257                             if( itab[k1*ksize+k2] < itab[mk1*ksize+mk2] )
 258                                 mk1 = k1, mk2 = k2;
 259                             else if( itab[k1*ksize+k2] > itab[Mk1*ksize+Mk2] )
 260                                 Mk1 = k1, Mk2 = k2;
 261                         }
 262                     if( diff < 0 )
 263                         itab[Mk1*ksize + Mk2] = (short)(itab[Mk1*ksize + Mk2] - diff);
 264                     else
 265                         itab[mk1*ksize + mk2] = (short)(itab[mk1*ksize + mk2] - diff);
 266                 }
 267             }
 268         tab -= INTER_TAB_SIZE2*ksize*ksize;
 269         itab -= INTER_TAB_SIZE2*ksize*ksize;
 270 #if CV_SSE2
 271         if( method == INTER_LINEAR )
 272         {
 273             for( i = 0; i < INTER_TAB_SIZE2; i++ )
 274                 for( j = 0; j < 4; j++ )
 275                 {
 276                     BilinearTab_iC4[i][0][j*2] = BilinearTab_i[i][0][0];
 277                     BilinearTab_iC4[i][0][j*2+1] = BilinearTab_i[i][0][1];
 278                     BilinearTab_iC4[i][1][j*2] = BilinearTab_i[i][1][0];
 279                     BilinearTab_iC4[i][1][j*2+1] = BilinearTab_i[i][1][1];
 280                 }
 281         }
 282 #endif
 283         inittab[method] = true;
 284     }
 285     return fixpt ? (const void*)itab : (const void*)tab;
 286 }
 287
 288 #ifndef __MINGW32__
 289 static bool initAllInterTab2D()
 290 {
 291     return  initInterTab2D( INTER_LINEAR, false ) &&
 292             initInterTab2D( INTER_LINEAR, true ) &&
 293             initInterTab2D( INTER_CUBIC, false ) &&
 294             initInterTab2D( INTER_CUBIC, true ) &&
 295             initInterTab2D( INTER_LANCZOS4, false ) &&
 296             initInterTab2D( INTER_LANCZOS4, true );
 297 }
 298
 299 static volatile bool doInitAllInterTab2D = initAllInterTab2D();
 300 #endif
 301
 302 template<typename ST, typename DT> struct Cast
 303 {
 304     typedef ST type1;
 305     typedef DT rtype;
 306
 307     DT operator()(ST val) const { return saturate_cast<DT>(val); }
 308 };
 309
 310 template<typename ST, typename DT, int bits> struct FixedPtCast
 311 {
 312     typedef ST type1;
 313     typedef DT rtype;
 314     enum { SHIFT = bits, DELTA = 1 << (bits-1) };
 315
 316     DT operator()(ST val) const { return saturate_cast<DT>((val + DELTA)>>SHIFT); }
 317 };
 318
 319 /****************************************************************************************\
 320 *                                         Resize                                         *
 321 \****************************************************************************************/
 322
 323 class resizeNNInvoker :
 324     public ParallelLoopBody
 325 {
 326 public:
 327     resizeNNInvoker(const Mat& _src, Mat &_dst, int *_x_ofs, int _pix_size4, double _ify) :
 328         ParallelLoopBody(), src(_src), dst(_dst), x_ofs(_x_ofs), pix_size4(_pix_size4),
 329         ify(_ify)
 330     {
 331     }
 332
 333     virtual void operator() (const Range& range) const
 334     {
 335         Size ssize = src.size(), dsize = dst.size();
 336         int y, x, pix_size = (int)src.elemSize();
 337
 338         for( y = range.start; y < range.end; y++ )
 339         {
 340             uchar* D = dst.data + dst.step*y;
 341             int sy = std::min(cvFloor(y*ify), ssize.height-1);
 342             const uchar* S = src.data + src.step*sy;
 343
 344             switch( pix_size )
 345             {
 346             case 1:
 347                 for( x = 0; x <= dsize.width - 2; x += 2 )
 348                 {
 349                     uchar t0 = S[x_ofs[x]];
 350                     uchar t1 = S[x_ofs[x+1]];
 351                     D[x] = t0;
 352                     D[x+1] = t1;
 353                 }
 354
 355                 for( ; x < dsize.width; x++ )
 356                     D[x] = S[x_ofs[x]];
 357                 break;
 358             case 2:
 359                 for( x = 0; x < dsize.width; x++ )
 360                     *(ushort*)(D + x*2) = *(ushort*)(S + x_ofs[x]);
 361                 break;
 362             case 3:
 363                 for( x = 0; x < dsize.width; x++, D += 3 )
 364                 {
 365                     const uchar* _tS = S + x_ofs[x];
 366                     D[0] = _tS[0]; D[1] = _tS[1]; D[2] = _tS[2];
 367                 }
 368                 break;
 369             case 4:
 370                 for( x = 0; x < dsize.width; x++ )
 371                     *(int*)(D + x*4) = *(int*)(S + x_ofs[x]);
 372                 break;
 373             case 6:
 374                 for( x = 0; x < dsize.width; x++, D += 6 )
 375                 {
 376                     const ushort* _tS = (const ushort*)(S + x_ofs[x]);
 377                     ushort* _tD = (ushort*)D;
 378                     _tD[0] = _tS[0]; _tD[1] = _tS[1]; _tD[2] = _tS[2];
 379                 }
 380                 break;
 381             case 8:
 382                 for( x = 0; x < dsize.width; x++, D += 8 )
 383                 {
 384                     const int* _tS = (const int*)(S + x_ofs[x]);
 385                     int* _tD = (int*)D;
 386                     _tD[0] = _tS[0]; _tD[1] = _tS[1];
 387                 }
 388                 break;
 389             case 12:
 390                 for( x = 0; x < dsize.width; x++, D += 12 )
 391                 {
 392                     const int* _tS = (const int*)(S + x_ofs[x]);
 393                     int* _tD = (int*)D;
 394                     _tD[0] = _tS[0]; _tD[1] = _tS[1]; _tD[2] = _tS[2];
 395                 }
 396                 break;
 397             default:
 398                 for( x = 0; x < dsize.width; x++, D += pix_size )
 399                 {
 400                     const int* _tS = (const int*)(S + x_ofs[x]);
 401                     int* _tD = (int*)D;
 402                     for( int k = 0; k < pix_size4; k++ )
 403                         _tD[k] = _tS[k];
 404                 }
 405             }
 406         }
 407     }
 408
 409 private:
 410     const Mat src;
 411     Mat dst;
 412     int* x_ofs, pix_size4;
 413     double ify;
 414
 415     resizeNNInvoker(const resizeNNInvoker&);
 416     resizeNNInvoker& operator=(const resizeNNInvoker&);
 417 };
 418
 419 static void
 420 resizeNN( const Mat& src, Mat& dst, double fx, double fy )
 421 {
 422     Size ssize = src.size(), dsize = dst.size();
 423     AutoBuffer<int> _x_ofs(dsize.width);
 424     int* x_ofs = _x_ofs;
 425     int pix_size = (int)src.elemSize();
 426     int pix_size4 = (int)(pix_size / sizeof(int));
 427     double ifx = 1./fx, ify = 1./fy;
 428     int x;
 429
 430     for( x = 0; x < dsize.width; x++ )
 431     {
 432         int sx = cvFloor(x*ifx);
 433         x_ofs[x] = std::min(sx, ssize.width-1)*pix_size;
 434     }
 435
 436     Range range(0, dsize.height);
 437     resizeNNInvoker invoker(src, dst, x_ofs, pix_size4, ify);
 438     parallel_for_(range, invoker, dst.total()/(double)(1<<16));
 439 }
 440
 441
 442 struct VResizeNoVec
 443 {
 444     int operator()(const uchar**, uchar*, const uchar*, int ) const { return 0; }
 445 };
 446
 447 struct HResizeNoVec
 448 {
 449     int operator()(const uchar**, uchar**, int, const int*,
 450         const uchar*, int, int, int, int, int) const { return 0; }
 451 };
 452
 453 #if CV_SSE2
 454
 455 struct VResizeLinearVec_32s8u
 456 {
 457     int operator()(const uchar** _src, uchar* dst, const uchar* _beta, int width ) const
 458     {
 459         if( !checkHardwareSupport(CV_CPU_SSE2) )
 460             return 0;
 461
 462         const int** src = (const int**)_src;
 463         const short* beta = (const short*)_beta;
 464         const int *S0 = src[0], *S1 = src[1];
 465         int x = 0;
 466         __m128i b0 = _mm_set1_epi16(beta[0]), b1 = _mm_set1_epi16(beta[1]);
 467         __m128i delta = _mm_set1_epi16(2);
 468
 469         if( (((size_t)S0|(size_t)S1)&15) == 0 )
 470             for( ; x <= width - 16; x += 16 )
 471             {
 472                 __m128i x0, x1, x2, y0, y1, y2;
 473                 x0 = _mm_load_si128((const __m128i*)(S0 + x));
 474                 x1 = _mm_load_si128((const __m128i*)(S0 + x + 4));
 475                 y0 = _mm_load_si128((const __m128i*)(S1 + x));
 476                 y1 = _mm_load_si128((const __m128i*)(S1 + x + 4));
 477                 x0 = _mm_packs_epi32(_mm_srai_epi32(x0, 4), _mm_srai_epi32(x1, 4));
 478                 y0 = _mm_packs_epi32(_mm_srai_epi32(y0, 4), _mm_srai_epi32(y1, 4));
 479
 480                 x1 = _mm_load_si128((const __m128i*)(S0 + x + 8));
 481                 x2 = _mm_load_si128((const __m128i*)(S0 + x + 12));
 482                 y1 = _mm_load_si128((const __m128i*)(S1 + x + 8));
 483                 y2 = _mm_load_si128((const __m128i*)(S1 + x + 12));
 484                 x1 = _mm_packs_epi32(_mm_srai_epi32(x1, 4), _mm_srai_epi32(x2, 4));
 485                 y1 = _mm_packs_epi32(_mm_srai_epi32(y1, 4), _mm_srai_epi32(y2, 4));
 486
 487                 x0 = _mm_adds_epi16(_mm_mulhi_epi16( x0, b0 ), _mm_mulhi_epi16( y0, b1 ));
 488                 x1 = _mm_adds_epi16(_mm_mulhi_epi16( x1, b0 ), _mm_mulhi_epi16( y1, b1 ));
 489
 490                 x0 = _mm_srai_epi16(_mm_adds_epi16(x0, delta), 2);
 491                 x1 = _mm_srai_epi16(_mm_adds_epi16(x1, delta), 2);
 492                 _mm_storeu_si128( (__m128i*)(dst + x), _mm_packus_epi16(x0, x1));
 493             }
 494         else
 495             for( ; x <= width - 16; x += 16 )
 496             {
 497                 __m128i x0, x1, x2, y0, y1, y2;
 498                 x0 = _mm_loadu_si128((const __m128i*)(S0 + x));
 499                 x1 = _mm_loadu_si128((const __m128i*)(S0 + x + 4));
 500                 y0 = _mm_loadu_si128((const __m128i*)(S1 + x));
 501                 y1 = _mm_loadu_si128((const __m128i*)(S1 + x + 4));
 502                 x0 = _mm_packs_epi32(_mm_srai_epi32(x0, 4), _mm_srai_epi32(x1, 4));
 503                 y0 = _mm_packs_epi32(_mm_srai_epi32(y0, 4), _mm_srai_epi32(y1, 4));
 504
 505                 x1 = _mm_loadu_si128((const __m128i*)(S0 + x + 8));
 506                 x2 = _mm_loadu_si128((const __m128i*)(S0 + x + 12));
 507                 y1 = _mm_loadu_si128((const __m128i*)(S1 + x + 8));
 508                 y2 = _mm_loadu_si128((const __m128i*)(S1 + x + 12));
 509                 x1 = _mm_packs_epi32(_mm_srai_epi32(x1, 4), _mm_srai_epi32(x2, 4));
 510                 y1 = _mm_packs_epi32(_mm_srai_epi32(y1, 4), _mm_srai_epi32(y2, 4));
 511
 512                 x0 = _mm_adds_epi16(_mm_mulhi_epi16( x0, b0 ), _mm_mulhi_epi16( y0, b1 ));
 513                 x1 = _mm_adds_epi16(_mm_mulhi_epi16( x1, b0 ), _mm_mulhi_epi16( y1, b1 ));
 514
 515                 x0 = _mm_srai_epi16(_mm_adds_epi16(x0, delta), 2);
 516                 x1 = _mm_srai_epi16(_mm_adds_epi16(x1, delta), 2);
 517                 _mm_storeu_si128( (__m128i*)(dst + x), _mm_packus_epi16(x0, x1));
 518             }
 519
 520         for( ; x < width - 4; x += 4 )
 521         {
 522             __m128i x0, y0;
 523             x0 = _mm_srai_epi32(_mm_loadu_si128((const __m128i*)(S0 + x)), 4);
 524             y0 = _mm_srai_epi32(_mm_loadu_si128((const __m128i*)(S1 + x)), 4);
 525             x0 = _mm_packs_epi32(x0, x0);
 526             y0 = _mm_packs_epi32(y0, y0);
 527             x0 = _mm_adds_epi16(_mm_mulhi_epi16(x0, b0), _mm_mulhi_epi16(y0, b1));
 528             x0 = _mm_srai_epi16(_mm_adds_epi16(x0, delta), 2);
 529             x0 = _mm_packus_epi16(x0, x0);
 530             *(int*)(dst + x) = _mm_cvtsi128_si32(x0);
 531         }
 532
 533         return x;
 534     }
 535 };
 536
 537
 538 template<int shiftval> struct VResizeLinearVec_32f16
 539 {
 540     int operator()(const uchar** _src, uchar* _dst, const uchar* _beta, int width ) const
 541     {
 542         if( !checkHardwareSupport(CV_CPU_SSE2) )
 543             return 0;
 544
 545         const float** src = (const float**)_src;
 546         const float* beta = (const float*)_beta;
 547         const float *S0 = src[0], *S1 = src[1];
 548         ushort* dst = (ushort*)_dst;
 549         int x = 0;
 550
 551         __m128 b0 = _mm_set1_ps(beta[0]), b1 = _mm_set1_ps(beta[1]);
 552         __m128i preshift = _mm_set1_epi32(shiftval);
 553         __m128i postshift = _mm_set1_epi16((short)shiftval);
 554
 555         if( (((size_t)S0|(size_t)S1)&15) == 0 )
 556             for( ; x <= width - 16; x += 16 )
 557             {
 558                 __m128 x0, x1, y0, y1;
 559                 __m128i t0, t1, t2;
 560                 x0 = _mm_load_ps(S0 + x);
 561                 x1 = _mm_load_ps(S0 + x + 4);
 562                 y0 = _mm_load_ps(S1 + x);
 563                 y1 = _mm_load_ps(S1 + x + 4);
 564
 565                 x0 = _mm_add_ps(_mm_mul_ps(x0, b0), _mm_mul_ps(y0, b1));
 566                 x1 = _mm_add_ps(_mm_mul_ps(x1, b0), _mm_mul_ps(y1, b1));
 567                 t0 = _mm_add_epi32(_mm_cvtps_epi32(x0), preshift);
 568                 t2 = _mm_add_epi32(_mm_cvtps_epi32(x1), preshift);
 569                 t0 = _mm_add_epi16(_mm_packs_epi32(t0, t2), postshift);
 570
 571                 x0 = _mm_load_ps(S0 + x + 8);
 572                 x1 = _mm_load_ps(S0 + x + 12);
 573                 y0 = _mm_load_ps(S1 + x + 8);
 574                 y1 = _mm_load_ps(S1 + x + 12);
 575
 576                 x0 = _mm_add_ps(_mm_mul_ps(x0, b0), _mm_mul_ps(y0, b1));
 577                 x1 = _mm_add_ps(_mm_mul_ps(x1, b0), _mm_mul_ps(y1, b1));
 578                 t1 = _mm_add_epi32(_mm_cvtps_epi32(x0), preshift);
 579                 t2 = _mm_add_epi32(_mm_cvtps_epi32(x1), preshift);
 580                 t1 = _mm_add_epi16(_mm_packs_epi32(t1, t2), postshift);
 581
 582                 _mm_storeu_si128( (__m128i*)(dst + x), t0);
 583                 _mm_storeu_si128( (__m128i*)(dst + x + 8), t1);
 584             }
 585         else
 586             for( ; x <= width - 16; x += 16 )
 587             {
 588                 __m128 x0, x1, y0, y1;
 589                 __m128i t0, t1, t2;
 590                 x0 = _mm_loadu_ps(S0 + x);
 591                 x1 = _mm_loadu_ps(S0 + x + 4);
 592                 y0 = _mm_loadu_ps(S1 + x);
 593                 y1 = _mm_loadu_ps(S1 + x + 4);
 594
 595                 x0 = _mm_add_ps(_mm_mul_ps(x0, b0), _mm_mul_ps(y0, b1));
 596                 x1 = _mm_add_ps(_mm_mul_ps(x1, b0), _mm_mul_ps(y1, b1));
 597                 t0 = _mm_add_epi32(_mm_cvtps_epi32(x0), preshift);
 598                 t2 = _mm_add_epi32(_mm_cvtps_epi32(x1), preshift);
 599                 t0 = _mm_add_epi16(_mm_packs_epi32(t0, t2), postshift);
 600
 601                 x0 = _mm_loadu_ps(S0 + x + 8);
 602                 x1 = _mm_loadu_ps(S0 + x + 12);
 603                 y0 = _mm_loadu_ps(S1 + x + 8);
 604                 y1 = _mm_loadu_ps(S1 + x + 12);
 605
 606                 x0 = _mm_add_ps(_mm_mul_ps(x0, b0), _mm_mul_ps(y0, b1));
 607                 x1 = _mm_add_ps(_mm_mul_ps(x1, b0), _mm_mul_ps(y1, b1));
 608                 t1 = _mm_add_epi32(_mm_cvtps_epi32(x0), preshift);
 609                 t2 = _mm_add_epi32(_mm_cvtps_epi32(x1), preshift);
 610                 t1 = _mm_add_epi16(_mm_packs_epi32(t1, t2), postshift);
 611
 612                 _mm_storeu_si128( (__m128i*)(dst + x), t0);
 613                 _mm_storeu_si128( (__m128i*)(dst + x + 8), t1);
 614             }
 615
 616         for( ; x < width - 4; x += 4 )
 617         {
 618             __m128 x0, y0;
 619             __m128i t0;
 620             x0 = _mm_loadu_ps(S0 + x);
 621             y0 = _mm_loadu_ps(S1 + x);
 622
 623             x0 = _mm_add_ps(_mm_mul_ps(x0, b0), _mm_mul_ps(y0, b1));
 624             t0 = _mm_add_epi32(_mm_cvtps_epi32(x0), preshift);
 625             t0 = _mm_add_epi16(_mm_packs_epi32(t0, t0), postshift);
 626             _mm_storel_epi64( (__m128i*)(dst + x), t0);
 627         }
 628
 629         return x;
 630     }
 631 };
 632
 633 typedef VResizeLinearVec_32f16<SHRT_MIN> VResizeLinearVec_32f16u;
 634 typedef VResizeLinearVec_32f16<0> VResizeLinearVec_32f16s;
 635
 636 struct VResizeLinearVec_32f
 637 {
 638     int operator()(const uchar** _src, uchar* _dst, const uchar* _beta, int width ) const
 639     {
 640         if( !checkHardwareSupport(CV_CPU_SSE) )
 641             return 0;
 642
 643         const float** src = (const float**)_src;
 644         const float* beta = (const float*)_beta;
 645         const float *S0 = src[0], *S1 = src[1];
 646         float* dst = (float*)_dst;
 647         int x = 0;
 648
 649         __m128 b0 = _mm_set1_ps(beta[0]), b1 = _mm_set1_ps(beta[1]);
 650
 651         if( (((size_t)S0|(size_t)S1)&15) == 0 )
 652             for( ; x <= width - 8; x += 8 )
 653             {
 654                 __m128 x0, x1, y0, y1;
 655                 x0 = _mm_load_ps(S0 + x);
 656                 x1 = _mm_load_ps(S0 + x + 4);
 657                 y0 = _mm_load_ps(S1 + x);
 658                 y1 = _mm_load_ps(S1 + x + 4);
 659
 660                 x0 = _mm_add_ps(_mm_mul_ps(x0, b0), _mm_mul_ps(y0, b1));
 661                 x1 = _mm_add_ps(_mm_mul_ps(x1, b0), _mm_mul_ps(y1, b1));
 662
 663                 _mm_storeu_ps( dst + x, x0);
 664                 _mm_storeu_ps( dst + x + 4, x1);
 665             }
 666         else
 667             for( ; x <= width - 8; x += 8 )
 668             {
 669                 __m128 x0, x1, y0, y1;
 670                 x0 = _mm_loadu_ps(S0 + x);
 671                 x1 = _mm_loadu_ps(S0 + x + 4);
 672                 y0 = _mm_loadu_ps(S1 + x);
 673                 y1 = _mm_loadu_ps(S1 + x + 4);
 674
 675                 x0 = _mm_add_ps(_mm_mul_ps(x0, b0), _mm_mul_ps(y0, b1));
 676                 x1 = _mm_add_ps(_mm_mul_ps(x1, b0), _mm_mul_ps(y1, b1));
 677
 678                 _mm_storeu_ps( dst + x, x0);
 679                 _mm_storeu_ps( dst + x + 4, x1);
 680             }
 681
 682         return x;
 683     }
 684 };
 685
 686
 687 struct VResizeCubicVec_32s8u
 688 {
 689     int operator()(const uchar** _src, uchar* dst, const uchar* _beta, int width ) const
 690     {
 691         if( !checkHardwareSupport(CV_CPU_SSE2) )
 692             return 0;
 693
 694         const int** src = (const int**)_src;
 695         const short* beta = (const short*)_beta;
 696         const int *S0 = src[0], *S1 = src[1], *S2 = src[2], *S3 = src[3];
 697         int x = 0;
 698         float scale = 1.f/(INTER_RESIZE_COEF_SCALE*INTER_RESIZE_COEF_SCALE);
 699         __m128 b0 = _mm_set1_ps(beta[0]*scale), b1 = _mm_set1_ps(beta[1]*scale),
 700             b2 = _mm_set1_ps(beta[2]*scale), b3 = _mm_set1_ps(beta[3]*scale);
 701
 702         if( (((size_t)S0|(size_t)S1|(size_t)S2|(size_t)S3)&15) == 0 )
 703             for( ; x <= width - 8; x += 8 )
 704             {
 705                 __m128i x0, x1, y0, y1;
 706                 __m128 s0, s1, f0, f1;
 707                 x0 = _mm_load_si128((const __m128i*)(S0 + x));
 708                 x1 = _mm_load_si128((const __m128i*)(S0 + x + 4));
 709                 y0 = _mm_load_si128((const __m128i*)(S1 + x));
 710                 y1 = _mm_load_si128((const __m128i*)(S1 + x + 4));
 711
 712                 s0 = _mm_mul_ps(_mm_cvtepi32_ps(x0), b0);
 713                 s1 = _mm_mul_ps(_mm_cvtepi32_ps(x1), b0);
 714                 f0 = _mm_mul_ps(_mm_cvtepi32_ps(y0), b1);
 715                 f1 = _mm_mul_ps(_mm_cvtepi32_ps(y1), b1);
 716                 s0 = _mm_add_ps(s0, f0);
 717                 s1 = _mm_add_ps(s1, f1);
 718
 719                 x0 = _mm_load_si128((const __m128i*)(S2 + x));
 720                 x1 = _mm_load_si128((const __m128i*)(S2 + x + 4));
 721                 y0 = _mm_load_si128((const __m128i*)(S3 + x));
 722                 y1 = _mm_load_si128((const __m128i*)(S3 + x + 4));
 723
 724                 f0 = _mm_mul_ps(_mm_cvtepi32_ps(x0), b2);
 725                 f1 = _mm_mul_ps(_mm_cvtepi32_ps(x1), b2);
 726                 s0 = _mm_add_ps(s0, f0);
 727                 s1 = _mm_add_ps(s1, f1);
 728                 f0 = _mm_mul_ps(_mm_cvtepi32_ps(y0), b3);
 729                 f1 = _mm_mul_ps(_mm_cvtepi32_ps(y1), b3);
 730                 s0 = _mm_add_ps(s0, f0);
 731                 s1 = _mm_add_ps(s1, f1);
 732
 733                 x0 = _mm_cvtps_epi32(s0);
 734                 x1 = _mm_cvtps_epi32(s1);
 735
 736                 x0 = _mm_packs_epi32(x0, x1);
 737                 _mm_storel_epi64( (__m128i*)(dst + x), _mm_packus_epi16(x0, x0));
 738             }
 739         else
 740             for( ; x <= width - 8; x += 8 )
 741             {
 742                 __m128i x0, x1, y0, y1;
 743                 __m128 s0, s1, f0, f1;
 744                 x0 = _mm_loadu_si128((const __m128i*)(S0 + x));
 745                 x1 = _mm_loadu_si128((const __m128i*)(S0 + x + 4));
 746                 y0 = _mm_loadu_si128((const __m128i*)(S1 + x));
 747                 y1 = _mm_loadu_si128((const __m128i*)(S1 + x + 4));
 748
 749                 s0 = _mm_mul_ps(_mm_cvtepi32_ps(x0), b0);
 750                 s1 = _mm_mul_ps(_mm_cvtepi32_ps(x1), b0);
 751                 f0 = _mm_mul_ps(_mm_cvtepi32_ps(y0), b1);
 752                 f1 = _mm_mul_ps(_mm_cvtepi32_ps(y1), b1);
 753                 s0 = _mm_add_ps(s0, f0);
 754                 s1 = _mm_add_ps(s1, f1);
 755
 756                 x0 = _mm_loadu_si128((const __m128i*)(S2 + x));
 757                 x1 = _mm_loadu_si128((const __m128i*)(S2 + x + 4));
 758                 y0 = _mm_loadu_si128((const __m128i*)(S3 + x));
 759                 y1 = _mm_loadu_si128((const __m128i*)(S3 + x + 4));
 760
 761                 f0 = _mm_mul_ps(_mm_cvtepi32_ps(x0), b2);
 762                 f1 = _mm_mul_ps(_mm_cvtepi32_ps(x1), b2);
 763                 s0 = _mm_add_ps(s0, f0);
 764                 s1 = _mm_add_ps(s1, f1);
 765                 f0 = _mm_mul_ps(_mm_cvtepi32_ps(y0), b3);
 766                 f1 = _mm_mul_ps(_mm_cvtepi32_ps(y1), b3);
 767                 s0 = _mm_add_ps(s0, f0);
 768                 s1 = _mm_add_ps(s1, f1);
 769
 770                 x0 = _mm_cvtps_epi32(s0);
 771                 x1 = _mm_cvtps_epi32(s1);
 772
 773                 x0 = _mm_packs_epi32(x0, x1);
 774                 _mm_storel_epi64( (__m128i*)(dst + x), _mm_packus_epi16(x0, x0));
 775             }
 776
 777         return x;
 778     }
 779 };
 780
 781
 782 template<int shiftval> struct VResizeCubicVec_32f16
 783 {
 784     int operator()(const uchar** _src, uchar* _dst, const uchar* _beta, int width ) const
 785     {
 786         if( !checkHardwareSupport(CV_CPU_SSE2) )
 787             return 0;
 788
 789         const float** src = (const float**)_src;
 790         const float* beta = (const float*)_beta;
 791         const float *S0 = src[0], *S1 = src[1], *S2 = src[2], *S3 = src[3];
 792         ushort* dst = (ushort*)_dst;
 793         int x = 0;
 794         __m128 b0 = _mm_set1_ps(beta[0]), b1 = _mm_set1_ps(beta[1]),
 795             b2 = _mm_set1_ps(beta[2]), b3 = _mm_set1_ps(beta[3]);
 796         __m128i preshift = _mm_set1_epi32(shiftval);
 797         __m128i postshift = _mm_set1_epi16((short)shiftval);
 798
 799         for( ; x <= width - 8; x += 8 )
 800         {
 801             __m128 x0, x1, y0, y1, s0, s1;
 802             __m128i t0, t1;
 803             x0 = _mm_loadu_ps(S0 + x);
 804             x1 = _mm_loadu_ps(S0 + x + 4);
 805             y0 = _mm_loadu_ps(S1 + x);
 806             y1 = _mm_loadu_ps(S1 + x + 4);
 807
 808             s0 = _mm_mul_ps(x0, b0);
 809             s1 = _mm_mul_ps(x1, b0);
 810             y0 = _mm_mul_ps(y0, b1);
 811             y1 = _mm_mul_ps(y1, b1);
 812             s0 = _mm_add_ps(s0, y0);
 813             s1 = _mm_add_ps(s1, y1);
 814
 815             x0 = _mm_loadu_ps(S2 + x);
 816             x1 = _mm_loadu_ps(S2 + x + 4);
 817             y0 = _mm_loadu_ps(S3 + x);
 818             y1 = _mm_loadu_ps(S3 + x + 4);
 819
 820             x0 = _mm_mul_ps(x0, b2);
 821             x1 = _mm_mul_ps(x1, b2);
 822             y0 = _mm_mul_ps(y0, b3);
 823             y1 = _mm_mul_ps(y1, b3);
 824             s0 = _mm_add_ps(s0, x0);
 825             s1 = _mm_add_ps(s1, x1);
 826             s0 = _mm_add_ps(s0, y0);
 827             s1 = _mm_add_ps(s1, y1);
 828
 829             t0 = _mm_add_epi32(_mm_cvtps_epi32(s0), preshift);
 830             t1 = _mm_add_epi32(_mm_cvtps_epi32(s1), preshift);
 831
 832             t0 = _mm_add_epi16(_mm_packs_epi32(t0, t1), postshift);
 833             _mm_storeu_si128( (__m128i*)(dst + x), t0);
 834         }
 835
 836         return x;
 837     }
 838 };
 839
 840 typedef VResizeCubicVec_32f16<SHRT_MIN> VResizeCubicVec_32f16u;
 841 typedef VResizeCubicVec_32f16<0> VResizeCubicVec_32f16s;
 842
 843 struct VResizeCubicVec_32f
 844 {
 845     int operator()(const uchar** _src, uchar* _dst, const uchar* _beta, int width ) const
 846     {
 847         if( !checkHardwareSupport(CV_CPU_SSE) )
 848             return 0;
 849
 850         const float** src = (const float**)_src;
 851         const float* beta = (const float*)_beta;
 852         const float *S0 = src[0], *S1 = src[1], *S2 = src[2], *S3 = src[3];
 853         float* dst = (float*)_dst;
 854         int x = 0;
 855         __m128 b0 = _mm_set1_ps(beta[0]), b1 = _mm_set1_ps(beta[1]),
 856             b2 = _mm_set1_ps(beta[2]), b3 = _mm_set1_ps(beta[3]);
 857
 858         for( ; x <= width - 8; x += 8 )
 859         {
 860             __m128 x0, x1, y0, y1, s0, s1;
 861             x0 = _mm_loadu_ps(S0 + x);
 862             x1 = _mm_loadu_ps(S0 + x + 4);
 863             y0 = _mm_loadu_ps(S1 + x);
 864             y1 = _mm_loadu_ps(S1 + x + 4);
 865
 866             s0 = _mm_mul_ps(x0, b0);
 867             s1 = _mm_mul_ps(x1, b0);
 868             y0 = _mm_mul_ps(y0, b1);
 869             y1 = _mm_mul_ps(y1, b1);
 870             s0 = _mm_add_ps(s0, y0);
 871             s1 = _mm_add_ps(s1, y1);
 872
 873             x0 = _mm_loadu_ps(S2 + x);
 874             x1 = _mm_loadu_ps(S2 + x + 4);
 875             y0 = _mm_loadu_ps(S3 + x);
 876             y1 = _mm_loadu_ps(S3 + x + 4);
 877
 878             x0 = _mm_mul_ps(x0, b2);
 879             x1 = _mm_mul_ps(x1, b2);
 880             y0 = _mm_mul_ps(y0, b3);
 881             y1 = _mm_mul_ps(y1, b3);
 882             s0 = _mm_add_ps(s0, x0);
 883             s1 = _mm_add_ps(s1, x1);
 884             s0 = _mm_add_ps(s0, y0);
 885             s1 = _mm_add_ps(s1, y1);
 886
 887             _mm_storeu_ps( dst + x, s0);
 888             _mm_storeu_ps( dst + x + 4, s1);
 889         }
 890
 891         return x;
 892     }
 893 };
 894
 895 #else
 896
 897 typedef VResizeNoVec VResizeLinearVec_32s8u;
 898 typedef VResizeNoVec VResizeLinearVec_32f16u;
 899 typedef VResizeNoVec VResizeLinearVec_32f16s;
 900 typedef VResizeNoVec VResizeLinearVec_32f;
 901
 902 typedef VResizeNoVec VResizeCubicVec_32s8u;
 903 typedef VResizeNoVec VResizeCubicVec_32f16u;
 904 typedef VResizeNoVec VResizeCubicVec_32f16s;
 905 typedef VResizeNoVec VResizeCubicVec_32f;
 906
 907 #endif
 908
 909 typedef HResizeNoVec HResizeLinearVec_8u32s;
 910 typedef HResizeNoVec HResizeLinearVec_16u32f;
 911 typedef HResizeNoVec HResizeLinearVec_16s32f;
 912 typedef HResizeNoVec HResizeLinearVec_32f;
 913 typedef HResizeNoVec HResizeLinearVec_64f;
 914
 915
 916 template<typename T, typename WT, typename AT, int ONE, class VecOp>
 917 struct HResizeLinear
 918 {
 919     typedef T value_type;
 920     typedef WT buf_type;
 921     typedef AT alpha_type;
 922
 923     void operator()(const T** src, WT** dst, int count,
 924                     const int* xofs, const AT* alpha,
 925                     int swidth, int dwidth, int cn, int xmin, int xmax ) const
 926     {
 927         int dx, k;
 928         VecOp vecOp;
 929
 930         int dx0 = vecOp((const uchar**)src, (uchar**)dst, count,
 931             xofs, (const uchar*)alpha, swidth, dwidth, cn, xmin, xmax );
 932
 933         for( k = 0; k <= count - 2; k++ )
 934         {
 935             const T *S0 = src[k], *S1 = src[k+1];
 936             WT *D0 = dst[k], *D1 = dst[k+1];
 937             for( dx = dx0; dx < xmax; dx++ )
 938             {
 939                 int sx = xofs[dx];
 940                 WT a0 = alpha[dx*2], a1 = alpha[dx*2+1];
 941                 WT t0 = S0[sx]*a0 + S0[sx + cn]*a1;
 942                 WT t1 = S1[sx]*a0 + S1[sx + cn]*a1;
 943                 D0[dx] = t0; D1[dx] = t1;
 944             }
 945
 946             for( ; dx < dwidth; dx++ )
 947             {
 948                 int sx = xofs[dx];
 949                 D0[dx] = WT(S0[sx]*ONE); D1[dx] = WT(S1[sx]*ONE);
 950             }
 951         }
 952
 953         for( ; k < count; k++ )
 954         {
 955             const T *S = src[k];
 956             WT *D = dst[k];
 957             for( dx = 0; dx < xmax; dx++ )
 958             {
 959                 int sx = xofs[dx];
 960                 D[dx] = S[sx]*alpha[dx*2] + S[sx+cn]*alpha[dx*2+1];
 961             }
 962
 963             for( ; dx < dwidth; dx++ )
 964                 D[dx] = WT(S[xofs[dx]]*ONE);
 965         }
 966     }
 967 };
 968
 969
 970 template<typename T, typename WT, typename AT, class CastOp, class VecOp>
 971 struct VResizeLinear
 972 {
 973     typedef T value_type;
 974     typedef WT buf_type;
 975     typedef AT alpha_type;
 976
 977     void operator()(const WT** src, T* dst, const AT* beta, int width ) const
 978     {
 979         WT b0 = beta[0], b1 = beta[1];
 980         const WT *S0 = src[0], *S1 = src[1];
 981         CastOp castOp;
 982         VecOp vecOp;
 983
 984         int x = vecOp((const uchar**)src, (uchar*)dst, (const uchar*)beta, width);
 985         #if CV_ENABLE_UNROLLED
 986         for( ; x <= width - 4; x += 4 )
 987         {
 988             WT t0, t1;
 989             t0 = S0[x]*b0 + S1[x]*b1;
 990             t1 = S0[x+1]*b0 + S1[x+1]*b1;
 991             dst[x] = castOp(t0); dst[x+1] = castOp(t1);
 992             t0 = S0[x+2]*b0 + S1[x+2]*b1;
 993             t1 = S0[x+3]*b0 + S1[x+3]*b1;
 994             dst[x+2] = castOp(t0); dst[x+3] = castOp(t1);
 995         }
 996         #endif
 997         for( ; x < width; x++ )
 998             dst[x] = castOp(S0[x]*b0 + S1[x]*b1);
 999     }
1000 };
1001
1002 template<>
1003 struct VResizeLinear<uchar, int, short, FixedPtCast<int, uchar, INTER_RESIZE_COEF_BITS*2>, VResizeLinearVec_32s8u>
1004 {
1005     typedef uchar value_type;
1006     typedef int buf_type;
1007     typedef short alpha_type;
1008
1009     void operator()(const buf_type** src, value_type* dst, const alpha_type* beta, int width ) const
1010     {
1011         alpha_type b0 = beta[0], b1 = beta[1];
1012         const buf_type *S0 = src[0], *S1 = src[1];
1013         VResizeLinearVec_32s8u vecOp;
1014
1015         int x = vecOp((const uchar**)src, (uchar*)dst, (const uchar*)beta, width);
1016         #if CV_ENABLE_UNROLLED
1017         for( ; x <= width - 4; x += 4 )
1018         {
1019             dst[x+0] = uchar(( ((b0 * (S0[x+0] >> 4)) >> 16) + ((b1 * (S1[x+0] >> 4)) >> 16) + 2)>>2);
1020             dst[x+1] = uchar(( ((b0 * (S0[x+1] >> 4)) >> 16) + ((b1 * (S1[x+1] >> 4)) >> 16) + 2)>>2);
1021             dst[x+2] = uchar(( ((b0 * (S0[x+2] >> 4)) >> 16) + ((b1 * (S1[x+2] >> 4)) >> 16) + 2)>>2);
1022             dst[x+3] = uchar(( ((b0 * (S0[x+3] >> 4)) >> 16) + ((b1 * (S1[x+3] >> 4)) >> 16) + 2)>>2);
1023         }
1024         #endif
1025         for( ; x < width; x++ )
1026             dst[x] = uchar(( ((b0 * (S0[x] >> 4)) >> 16) + ((b1 * (S1[x] >> 4)) >> 16) + 2)>>2);
1027     }
1028 };
1029
1030
1031 template<typename T, typename WT, typename AT>
1032 struct HResizeCubic
1033 {
1034     typedef T value_type;
1035     typedef WT buf_type;
1036     typedef AT alpha_type;
1037
1038     void operator()(const T** src, WT** dst, int count,
1039                     const int* xofs, const AT* alpha,
1040                     int swidth, int dwidth, int cn, int xmin, int xmax ) const
1041     {
1042         for( int k = 0; k < count; k++ )
1043         {
1044             const T *S = src[k];
1045             WT *D = dst[k];
1046             int dx = 0, limit = xmin;
1047             for(;;)
1048             {
1049                 for( ; dx < limit; dx++, alpha += 4 )
1050                 {
1051                     int j, sx = xofs[dx] - cn;
1052                     WT v = 0;
1053                     for( j = 0; j < 4; j++ )
1054                     {
1055                         int sxj = sx + j*cn;
1056                         if( (unsigned)sxj >= (unsigned)swidth )
1057                         {
1058                             while( sxj < 0 )
1059                                 sxj += cn;
1060                             while( sxj >= swidth )
1061                                 sxj -= cn;
1062                         }
1063                         v += S[sxj]*alpha[j];
1064                     }
1065                     D[dx] = v;
1066                 }
1067                 if( limit == dwidth )
1068                     break;
1069                 for( ; dx < xmax; dx++, alpha += 4 )
1070                 {
1071                     int sx = xofs[dx];
1072                     D[dx] = S[sx-cn]*alpha[0] + S[sx]*alpha[1] +
1073                         S[sx+cn]*alpha[2] + S[sx+cn*2]*alpha[3];
1074                 }
1075                 limit = dwidth;
1076             }
1077             alpha -= dwidth*4;
1078         }
1079     }
1080 };
1081
1082
1083 template<typename T, typename WT, typename AT, class CastOp, class VecOp>
1084 struct VResizeCubic
1085 {
1086     typedef T value_type;
1087     typedef WT buf_type;
1088     typedef AT alpha_type;
1089
1090     void operator()(const WT** src, T* dst, const AT* beta, int width ) const
1091     {
1092         WT b0 = beta[0], b1 = beta[1], b2 = beta[2], b3 = beta[3];
1093         const WT *S0 = src[0], *S1 = src[1], *S2 = src[2], *S3 = src[3];
1094         CastOp castOp;
1095         VecOp vecOp;
1096
1097         int x = vecOp((const uchar**)src, (uchar*)dst, (const uchar*)beta, width);
1098         for( ; x < width; x++ )
1099             dst[x] = castOp(S0[x]*b0 + S1[x]*b1 + S2[x]*b2 + S3[x]*b3);
1100     }
1101 };
1102
1103
1104 template<typename T, typename WT, typename AT>
1105 struct HResizeLanczos4
1106 {
1107     typedef T value_type;
1108     typedef WT buf_type;
1109     typedef AT alpha_type;
1110
1111     void operator()(const T** src, WT** dst, int count,
1112                     const int* xofs, const AT* alpha,
1113                     int swidth, int dwidth, int cn, int xmin, int xmax ) const
1114     {
1115         for( int k = 0; k < count; k++ )
1116         {
1117             const T *S = src[k];
1118             WT *D = dst[k];
1119             int dx = 0, limit = xmin;
1120             for(;;)
1121             {
1122                 for( ; dx < limit; dx++, alpha += 8 )
1123                 {
1124                     int j, sx = xofs[dx] - cn*3;
1125                     WT v = 0;
1126                     for( j = 0; j < 8; j++ )
1127                     {
1128                         int sxj = sx + j*cn;
1129                         if( (unsigned)sxj >= (unsigned)swidth )
1130                         {
1131                             while( sxj < 0 )
1132                                 sxj += cn;
1133                             while( sxj >= swidth )
1134                                 sxj -= cn;
1135                         }
1136                         v += S[sxj]*alpha[j];
1137                     }
1138                     D[dx] = v;
1139                 }
1140                 if( limit == dwidth )
1141                     break;
1142                 for( ; dx < xmax; dx++, alpha += 8 )
1143                 {
1144                     int sx = xofs[dx];
1145                     D[dx] = S[sx-cn*3]*alpha[0] + S[sx-cn*2]*alpha[1] +
1146                         S[sx-cn]*alpha[2] + S[sx]*alpha[3] +
1147                         S[sx+cn]*alpha[4] + S[sx+cn*2]*alpha[5] +
1148                         S[sx+cn*3]*alpha[6] + S[sx+cn*4]*alpha[7];
1149                 }
1150                 limit = dwidth;
1151             }
1152             alpha -= dwidth*8;
1153         }
1154     }
1155 };
1156
1157
1158 template<typename T, typename WT, typename AT, class CastOp, class VecOp>
1159 struct VResizeLanczos4
1160 {
1161     typedef T value_type;
1162     typedef WT buf_type;
1163     typedef AT alpha_type;
1164
1165     void operator()(const WT** src, T* dst, const AT* beta, int width ) const
1166     {
1167         CastOp castOp;
1168         VecOp vecOp;
1169         int k, x = vecOp((const uchar**)src, (uchar*)dst, (const uchar*)beta, width);
1170         #if CV_ENABLE_UNROLLED
1171         for( ; x <= width - 4; x += 4 )
1172         {
1173             WT b = beta[0];
1174             const WT* S = src[0];
1175             WT s0 = S[x]*b, s1 = S[x+1]*b, s2 = S[x+2]*b, s3 = S[x+3]*b;
1176
1177             for( k = 1; k < 8; k++ )
1178             {
1179                 b = beta[k]; S = src[k];
1180                 s0 += S[x]*b; s1 += S[x+1]*b;
1181                 s2 += S[x+2]*b; s3 += S[x+3]*b;
1182             }
1183
1184             dst[x] = castOp(s0); dst[x+1] = castOp(s1);
1185             dst[x+2] = castOp(s2); dst[x+3] = castOp(s3);
1186         }
1187         #endif
1188         for( ; x < width; x++ )
1189         {
1190             dst[x] = castOp(src[0][x]*beta[0] + src[1][x]*beta[1] +
1191                 src[2][x]*beta[2] + src[3][x]*beta[3] + src[4][x]*beta[4] +
1192                 src[5][x]*beta[5] + src[6][x]*beta[6] + src[7][x]*beta[7]);
1193         }
1194     }
1195 };
1196
1197
1198 static inline int clip(int x, int a, int b)
1199 {
1200     return x >= a ? (x < b ? x : b-1) : a;
1201 }
1202
1203 static const int MAX_ESIZE=16;
1204
1205 template <typename HResize, typename VResize>
1206 class resizeGeneric_Invoker :
1207     public ParallelLoopBody
1208 {
1209 public:
1210     typedef typename HResize::value_type T;
1211     typedef typename HResize::buf_type WT;
1212     typedef typename HResize::alpha_type AT;
1213
1214     resizeGeneric_Invoker(const Mat& _src, Mat &_dst, const int *_xofs, const int *_yofs,
1215         const AT* _alpha, const AT* __beta, const Size& _ssize, const Size &_dsize,
1216         int _ksize, int _xmin, int _xmax) :
1217         ParallelLoopBody(), src(_src), dst(_dst), xofs(_xofs), yofs(_yofs),
1218         alpha(_alpha), _beta(__beta), ssize(_ssize), dsize(_dsize),
1219         ksize(_ksize), xmin(_xmin), xmax(_xmax)
1220     {
1221     }
1222
1223     virtual void operator() (const Range& range) const
1224     {
1225         int dy, cn = src.channels();
1226         HResize hresize;
1227         VResize vresize;
1228
1229         int bufstep = (int)alignSize(dsize.width, 16);
1230         AutoBuffer<WT> _buffer(bufstep*ksize);
1231         const T* srows[MAX_ESIZE]={0};
1232         WT* rows[MAX_ESIZE]={0};
1233         int prev_sy[MAX_ESIZE];
1234
1235         for(int k = 0; k < ksize; k++ )
1236         {
1237             prev_sy[k] = -1;
1238             rows[k] = (WT*)_buffer + bufstep*k;
1239         }
1240
1241         const AT* beta = _beta + ksize * range.start;
1242
1243         for( dy = range.start; dy < range.end; dy++, beta += ksize )
1244         {
1245             int sy0 = yofs[dy], k0=ksize, k1=0, ksize2 = ksize/2;
1246
1247             for(int k = 0; k < ksize; k++ )
1248             {
1249                 int sy = clip(sy0 - ksize2 + 1 + k, 0, ssize.height);
1250                 for( k1 = std::max(k1, k); k1 < ksize; k1++ )
1251                 {
1252                     if( sy == prev_sy[k1] ) // if the sy-th row has been computed already, reuse it.
1253                     {
1254                         if( k1 > k )
1255                             memcpy( rows[k], rows[k1], bufstep*sizeof(rows[0][0]) );
1256                         break;
1257                     }
1258                 }
1259                 if( k1 == ksize )
1260                     k0 = std::min(k0, k); // remember the first row that needs to be computed
1261                 srows[k] = (T*)(src.data + src.step*sy);
1262                 prev_sy[k] = sy;
1263             }
1264
1265             if( k0 < ksize )
1266                 hresize( (const T**)(srows + k0), (WT**)(rows + k0), ksize - k0, xofs, (const AT*)(alpha),
1267                         ssize.width, dsize.width, cn, xmin, xmax );
1268             vresize( (const WT**)rows, (T*)(dst.data + dst.step*dy), beta, dsize.width );
1269         }
1270     }
1271
1272 private:
1273     Mat src;
1274     Mat dst;
1275     const int* xofs, *yofs;
1276     const AT* alpha, *_beta;
1277     Size ssize, dsize;
1278     int ksize, xmin, xmax;
1279 };
1280
1281 template<class HResize, class VResize>
1282 static void resizeGeneric_( const Mat& src, Mat& dst,
1283                             const int* xofs, const void* _alpha,
1284                             const int* yofs, const void* _beta,
1285                             int xmin, int xmax, int ksize )
1286 {
1287     typedef typename HResize::alpha_type AT;
1288
1289     const AT* beta = (const AT*)_beta;
1290     Size ssize = src.size(), dsize = dst.size();
1291     int cn = src.channels();
1292     ssize.width *= cn;
1293     dsize.width *= cn;
1294     xmin *= cn;
1295     xmax *= cn;
1296     // image resize is a separable operation. In case of not too strong
1297
1298     Range range(0, dsize.height);
1299     resizeGeneric_Invoker<HResize, VResize> invoker(src, dst, xofs, yofs, (const AT*)_alpha, beta,
1300         ssize, dsize, ksize, xmin, xmax);
1301     parallel_for_(range, invoker, dst.total()/(double)(1<<16));
1302 }
1303
1304 template <typename T, typename WT>
1305 struct ResizeAreaFastNoVec
1306 {
1307     ResizeAreaFastNoVec(int, int) { }
1308     ResizeAreaFastNoVec(int, int, int, int) { }
1309     int operator() (const T*, T*, int) const
1310     { return 0; }
1311 };
1312
1313 #if CV_SSE2
1314 class ResizeAreaFastVec_SIMD_8u
1315 {
1316 public:
1317     ResizeAreaFastVec_SIMD_8u(int _cn, int _step) :
1318         cn(_cn), step(_step)
1319     {
1320         use_simd = checkHardwareSupport(CV_CPU_SSE2);
1321     }
1322
1323     int operator() (const uchar* S, uchar* D, int w) const
1324     {
1325         if (!use_simd)
1326             return 0;
1327
1328         int dx = 0;
1329         const uchar* S0 = S;
1330         const uchar* S1 = S0 + step;
1331         __m128i zero = _mm_setzero_si128();
1332         __m128i delta2 = _mm_set1_epi16(2);
1333
1334         if (cn == 1)
1335         {
1336             __m128i masklow = _mm_set1_epi16(0x00ff);
1337             for ( ; dx <= w - 8; dx += 8, S0 += 16, S1 += 16, D += 8)
1338             {
1339                 __m128i r0 = _mm_loadu_si128((const __m128i*)S0);
1340                 __m128i r1 = _mm_loadu_si128((const __m128i*)S1);
1341
1342                 __m128i s0 = _mm_add_epi16(_mm_srli_epi16(r0, 8), _mm_and_si128(r0, masklow));
1343                 __m128i s1 = _mm_add_epi16(_mm_srli_epi16(r1, 8), _mm_and_si128(r1, masklow));
1344                 s0 = _mm_add_epi16(_mm_add_epi16(s0, s1), delta2);
1345                 s0 = _mm_packus_epi16(_mm_srli_epi16(s0, 2), zero);
1346
1347                 _mm_storel_epi64((__m128i*)D, s0);
1348             }
1349         }
1350         else if (cn == 3)
1351             for ( ; dx <= w - 11; dx += 6, S0 += 12, S1 += 12, D += 6)
1352             {
1353                 __m128i r0 = _mm_loadu_si128((const __m128i*)S0);
1354                 __m128i r1 = _mm_loadu_si128((const __m128i*)S1);
1355
1356                 __m128i r0_16l = _mm_unpacklo_epi8(r0, zero);
1357                 __m128i r0_16h = _mm_unpacklo_epi8(_mm_srli_si128(r0, 6), zero);
1358                 __m128i r1_16l = _mm_unpacklo_epi8(r1, zero);
1359                 __m128i r1_16h = _mm_unpacklo_epi8(_mm_srli_si128(r1, 6), zero);
1360
1361                 __m128i s0 = _mm_add_epi16(r0_16l, _mm_srli_si128(r0_16l, 6));
1362                 __m128i s1 = _mm_add_epi16(r1_16l, _mm_srli_si128(r1_16l, 6));
1363                 s0 = _mm_add_epi16(s1, _mm_add_epi16(s0, delta2));
1364                 s0 = _mm_packus_epi16(_mm_srli_epi16(s0, 2), zero);
1365                 _mm_storel_epi64((__m128i*)D, s0);
1366
1367                 s0 = _mm_add_epi16(r0_16h, _mm_srli_si128(r0_16h, 6));
1368                 s1 = _mm_add_epi16(r1_16h, _mm_srli_si128(r1_16h, 6));
1369                 s0 = _mm_add_epi16(s1, _mm_add_epi16(s0, delta2));
1370                 s0 = _mm_packus_epi16(_mm_srli_epi16(s0, 2), zero);
1371                 _mm_storel_epi64((__m128i*)(D+3), s0);
1372             }
1373         else
1374         {
1375             CV_Assert(cn == 4);
1376             int v[] = { 0, 0, -1, -1 };
1377             __m128i mask = _mm_loadu_si128((const __m128i*)v);
1378
1379             for ( ; dx <= w - 8; dx += 8, S0 += 16, S1 += 16, D += 8)
1380             {
1381                 __m128i r0 = _mm_loadu_si128((const __m128i*)S0);
1382                 __m128i r1 = _mm_loadu_si128((const __m128i*)S1);
1383
1384                 __m128i r0_16l = _mm_unpacklo_epi8(r0, zero);
1385                 __m128i r0_16h = _mm_unpackhi_epi8(r0, zero);
1386                 __m128i r1_16l = _mm_unpacklo_epi8(r1, zero);
1387                 __m128i r1_16h = _mm_unpackhi_epi8(r1, zero);
1388
1389                 __m128i s0 = _mm_add_epi16(r0_16l, _mm_srli_si128(r0_16l, 8));
1390                 __m128i s1 = _mm_add_epi16(r1_16l, _mm_srli_si128(r1_16l, 8));
1391                 s0 = _mm_add_epi16(s1, _mm_add_epi16(s0, delta2));
1392                 __m128i res0 = _mm_srli_epi16(s0, 2);
1393
1394                 s0 = _mm_add_epi16(r0_16h, _mm_srli_si128(r0_16h, 8));
1395                 s1 = _mm_add_epi16(r1_16h, _mm_srli_si128(r1_16h, 8));
1396                 s0 = _mm_add_epi16(s1, _mm_add_epi16(s0, delta2));
1397                 __m128i res1 = _mm_srli_epi16(s0, 2);
1398                 s0 = _mm_packus_epi16(_mm_or_si128(_mm_andnot_si128(mask, res0),
1399                                                    _mm_and_si128(mask, _mm_slli_si128(res1, 8))), zero);
1400                 _mm_storel_epi64((__m128i*)(D), s0);
1401             }
1402         }
1403
1404         return dx;
1405     }
1406
1407 private:
1408     int cn;
1409     bool use_simd;
1410     int step;
1411 };
1412
1413 class ResizeAreaFastVec_SIMD_16u
1414 {
1415 public:
1416     ResizeAreaFastVec_SIMD_16u(int _cn, int _step) :
1417         cn(_cn), step(_step)
1418     {
1419         use_simd = checkHardwareSupport(CV_CPU_SSE2);
1420     }
1421
1422     int operator() (const ushort* S, ushort* D, int w) const
1423     {
1424         if (!use_simd)
1425             return 0;
1426
1427         int dx = 0;
1428         const ushort* S0 = (const ushort*)S;
1429         const ushort* S1 = (const ushort*)((const uchar*)(S) + step);
1430         __m128i masklow = _mm_set1_epi32(0x0000ffff);
1431         __m128i zero = _mm_setzero_si128();
1432         __m128i delta2 = _mm_set1_epi32(2);
1433
1434 #define _mm_packus_epi32(a, zero) _mm_packs_epi32(_mm_srai_epi32(_mm_slli_epi32(a, 16), 16), zero)
1435
1436         if (cn == 1)
1437         {
1438             for ( ; dx <= w - 4; dx += 4, S0 += 8, S1 += 8, D += 4)
1439             {
1440                 __m128i r0 = _mm_loadu_si128((const __m128i*)S0);
1441                 __m128i r1 = _mm_loadu_si128((const __m128i*)S1);
1442
1443                 __m128i s0 = _mm_add_epi32(_mm_srli_epi32(r0, 16), _mm_and_si128(r0, masklow));
1444                 __m128i s1 = _mm_add_epi32(_mm_srli_epi32(r1, 16), _mm_and_si128(r1, masklow));
1445                 s0 = _mm_add_epi32(_mm_add_epi32(s0, s1), delta2);
1446                 s0 = _mm_srli_epi32(s0, 2);
1447                 s0 = _mm_packus_epi32(s0, zero);
1448
1449                 _mm_storel_epi64((__m128i*)D, s0);
1450             }
1451         }
1452         else if (cn == 3)
1453             for ( ; dx <= w - 4; dx += 3, S0 += 6, S1 += 6, D += 3)
1454             {
1455                 __m128i r0 = _mm_loadu_si128((const __m128i*)S0);
1456                 __m128i r1 = _mm_loadu_si128((const __m128i*)S1);
1457
1458                 __m128i r0_16l = _mm_unpacklo_epi16(r0, zero);
1459                 __m128i r0_16h = _mm_unpacklo_epi16(_mm_srli_si128(r0, 6), zero);
1460                 __m128i r1_16l = _mm_unpacklo_epi16(r1, zero);
1461                 __m128i r1_16h = _mm_unpacklo_epi16(_mm_srli_si128(r1, 6), zero);
1462
1463                 __m128i s0 = _mm_add_epi32(r0_16l, r0_16h);
1464                 __m128i s1 = _mm_add_epi32(r1_16l, r1_16h);
1465                 s0 = _mm_add_epi32(delta2, _mm_add_epi32(s0, s1));
1466                 s0 = _mm_packus_epi32(_mm_srli_epi32(s0, 2), zero);
1467                 _mm_storel_epi64((__m128i*)D, s0);
1468             }
1469         else
1470         {
1471             CV_Assert(cn == 4);
1472             for ( ; dx <= w - 4; dx += 4, S0 += 8, S1 += 8, D += 4)
1473             {
1474                 __m128i r0 = _mm_loadu_si128((const __m128i*)S0);
1475                 __m128i r1 = _mm_loadu_si128((const __m128i*)S1);
1476
1477                 __m128i r0_32l = _mm_unpacklo_epi16(r0, zero);
1478                 __m128i r0_32h = _mm_unpackhi_epi16(r0, zero);
1479                 __m128i r1_32l = _mm_unpacklo_epi16(r1, zero);
1480                 __m128i r1_32h = _mm_unpackhi_epi16(r1, zero);
1481
1482                 __m128i s0 = _mm_add_epi32(r0_32l, r0_32h);
1483                 __m128i s1 = _mm_add_epi32(r1_32l, r1_32h);
1484                 s0 = _mm_add_epi32(s1, _mm_add_epi32(s0, delta2));
1485                 s0 = _mm_packus_epi32(_mm_srli_epi32(s0, 2), zero);
1486                 _mm_storel_epi64((__m128i*)D, s0);
1487             }
1488         }
1489
1490 #undef _mm_packus_epi32
1491
1492         return dx;
1493     }
1494
1495 private:
1496     int cn;
1497     int step;
1498     bool use_simd;
1499 };
1500
1501 #else
1502 typedef ResizeAreaFastNoVec<uchar, uchar> ResizeAreaFastVec_SIMD_8u;
1503 typedef ResizeAreaFastNoVec<ushort, ushort> ResizeAreaFastVec_SIMD_16u;
1504 #endif
1505
1506 template<typename T, typename SIMDVecOp>
1507 struct ResizeAreaFastVec
1508 {
1509     ResizeAreaFastVec(int _scale_x, int _scale_y, int _cn, int _step) :
1510         scale_x(_scale_x), scale_y(_scale_y), cn(_cn), step(_step), vecOp(_cn, _step)
1511     {
1512         fast_mode = scale_x == 2 && scale_y == 2 && (cn == 1 || cn == 3 || cn == 4);
1513     }
1514
1515     int operator() (const T* S, T* D, int w) const
1516     {
1517         if (!fast_mode)
1518             return 0;
1519
1520         const T* nextS = (const T*)((const uchar*)S + step);
1521         int dx = vecOp(S, D, w);
1522
1523         if (cn == 1)
1524             for( ; dx < w; ++dx )
1525             {
1526                 int index = dx*2;
1527                 D[dx] = (T)((S[index] + S[index+1] + nextS[index] + nextS[index+1] + 2) >> 2);
1528             }
1529         else if (cn == 3)
1530             for( ; dx < w; dx += 3 )
1531             {
1532                 int index = dx*2;
1533                 D[dx] = (T)((S[index] + S[index+3] + nextS[index] + nextS[index+3] + 2) >> 2);
1534                 D[dx+1] = (T)((S[index+1] + S[index+4] + nextS[index+1] + nextS[index+4] + 2) >> 2);
1535                 D[dx+2] = (T)((S[index+2] + S[index+5] + nextS[index+2] + nextS[index+5] + 2) >> 2);
1536             }
1537         else
1538             {
1539                 CV_Assert(cn == 4);
1540                 for( ; dx < w; dx += 4 )
1541                 {
1542                     int index = dx*2;
1543                     D[dx] = (T)((S[index] + S[index+4] + nextS[index] + nextS[index+4] + 2) >> 2);
1544                     D[dx+1] = (T)((S[index+1] + S[index+5] + nextS[index+1] + nextS[index+5] + 2) >> 2);
1545                     D[dx+2] = (T)((S[index+2] + S[index+6] + nextS[index+2] + nextS[index+6] + 2) >> 2);
1546                     D[dx+3] = (T)((S[index+3] + S[index+7] + nextS[index+3] + nextS[index+7] + 2) >> 2);
1547                 }
1548             }
1549
1550         return dx;
1551     }
1552
1553 private:
1554     int scale_x, scale_y;
1555     int cn;
1556     bool fast_mode;
1557     int step;
1558     SIMDVecOp vecOp;
1559 };
1560
1561 template <typename T, typename WT, typename VecOp>
1562 class resizeAreaFast_Invoker :
1563     public ParallelLoopBody
1564 {
1565 public:
1566     resizeAreaFast_Invoker(const Mat &_src, Mat &_dst,
1567         int _scale_x, int _scale_y, const int* _ofs, const int* _xofs) :
1568         ParallelLoopBody(), src(_src), dst(_dst), scale_x(_scale_x),
1569         scale_y(_scale_y), ofs(_ofs), xofs(_xofs)
1570     {
1571     }
1572
1573     virtual void operator() (const Range& range) const
1574     {
1575         Size ssize = src.size(), dsize = dst.size();
1576         int cn = src.channels();
1577         int area = scale_x*scale_y;
1578         float scale = 1.f/(area);
1579         int dwidth1 = (ssize.width/scale_x)*cn;
1580         dsize.width *= cn;
1581         ssize.width *= cn;
1582         int dy, dx, k = 0;
1583
1584         VecOp vop(scale_x, scale_y, src.channels(), (int)src.step/*, area_ofs*/);
1585
1586         for( dy = range.start; dy < range.end; dy++ )
1587         {
1588             T* D = (T*)(dst.data + dst.step*dy);
1589             int sy0 = dy*scale_y;
1590             int w = sy0 + scale_y <= ssize.height ? dwidth1 : 0;
1591
1592             if( sy0 >= ssize.height )
1593             {
1594                 for( dx = 0; dx < dsize.width; dx++ )
1595                     D[dx] = 0;
1596                 continue;
1597             }
1598
1599             dx = vop((const T*)(src.data + src.step * sy0), D, w);
1600             for( ; dx < w; dx++ )
1601             {
1602                 const T* S = (const T*)(src.data + src.step * sy0) + xofs[dx];
1603                 WT sum = 0;
1604                 k = 0;
1605                 #if CV_ENABLE_UNROLLED
1606                 for( ; k <= area - 4; k += 4 )
1607                     sum += S[ofs[k]] + S[ofs[k+1]] + S[ofs[k+2]] + S[ofs[k+3]];
1608                 #endif
1609                 for( ; k < area; k++ )
1610                     sum += S[ofs[k]];
1611
1612                 D[dx] = saturate_cast<T>(sum * scale);
1613             }
1614
1615             for( ; dx < dsize.width; dx++ )
1616             {
1617                 WT sum = 0;
1618                 int count = 0, sx0 = xofs[dx];
1619                 if( sx0 >= ssize.width )
1620                     D[dx] = 0;
1621
1622                 for( int sy = 0; sy < scale_y; sy++ )
1623                 {
1624                     if( sy0 + sy >= ssize.height )
1625                         break;
1626                     const T* S = (const T*)(src.data + src.step*(sy0 + sy)) + sx0;
1627                     for( int sx = 0; sx < scale_x*cn; sx += cn )
1628                     {
1629                         if( sx0 + sx >= ssize.width )
1630                             break;
1631                         sum += S[sx];
1632                         count++;
1633                     }
1634                 }
1635
1636                 D[dx] = saturate_cast<T>((float)sum/count);
1637             }
1638         }
1639     }
1640
1641 private:
1642     Mat src;
1643     Mat dst;
1644     int scale_x, scale_y;
1645     const int *ofs, *xofs;
1646 };
1647
1648 template<typename T, typename WT, typename VecOp>
1649 static void resizeAreaFast_( const Mat& src, Mat& dst, const int* ofs, const int* xofs,
1650                              int scale_x, int scale_y )
1651 {
1652     Range range(0, dst.rows);
1653     resizeAreaFast_Invoker<T, WT, VecOp> invoker(src, dst, scale_x,
1654         scale_y, ofs, xofs);
1655     parallel_for_(range, invoker, dst.total()/(double)(1<<16));
1656 }
1657
1658 struct DecimateAlpha
1659 {
1660     int si, di;
1661     float alpha;
1662 };
1663
1664
1665 template<typename T, typename WT> class ResizeArea_Invoker :
1666     public ParallelLoopBody
1667 {
1668 public:
1669     ResizeArea_Invoker( const Mat& _src, Mat& _dst,
1670                         const DecimateAlpha* _xtab, int _xtab_size,
1671                         const DecimateAlpha* _ytab, int _ytab_size,
1672                         const int* _tabofs )
1673     {
1674         src = &_src;
1675         dst = &_dst;
1676         xtab0 = _xtab;
1677         xtab_size0 = _xtab_size;
1678         ytab = _ytab;
1679         ytab_size = _ytab_size;
1680         tabofs = _tabofs;
1681     }
1682
1683     virtual void operator() (const Range& range) const
1684     {
1685         Size dsize = dst->size();
1686         int cn = dst->channels();
1687         dsize.width *= cn;
1688         AutoBuffer<WT> _buffer(dsize.width*2);
1689         const DecimateAlpha* xtab = xtab0;
1690         int xtab_size = xtab_size0;
1691         WT *buf = _buffer, *sum = buf + dsize.width;
1692         int j_start = tabofs[range.start], j_end = tabofs[range.end], j, k, dx, prev_dy = ytab[j_start].di;
1693
1694         for( dx = 0; dx < dsize.width; dx++ )
1695             sum[dx] = (WT)0;
1696
1697         for( j = j_start; j < j_end; j++ )
1698         {
1699             WT beta = ytab[j].alpha;
1700             int dy = ytab[j].di;
1701             int sy = ytab[j].si;
1702
1703             {
1704                 const T* S = (const T*)(src->data + src->step*sy);
1705                 for( dx = 0; dx < dsize.width; dx++ )
1706                     buf[dx] = (WT)0;
1707
1708                 if( cn == 1 )
1709                     for( k = 0; k < xtab_size; k++ )
1710                     {
1711                         int dxn = xtab[k].di;
1712                         WT alpha = xtab[k].alpha;
1713                         buf[dxn] += S[xtab[k].si]*alpha;
1714                     }
1715                 else if( cn == 2 )
1716                     for( k = 0; k < xtab_size; k++ )
1717                     {
1718                         int sxn = xtab[k].si;
1719                         int dxn = xtab[k].di;
1720                         WT alpha = xtab[k].alpha;
1721                         WT t0 = buf[dxn] + S[sxn]*alpha;
1722                         WT t1 = buf[dxn+1] + S[sxn+1]*alpha;
1723                         buf[dxn] = t0; buf[dxn+1] = t1;
1724                     }
1725                 else if( cn == 3 )
1726                     for( k = 0; k < xtab_size; k++ )
1727                     {
1728                         int sxn = xtab[k].si;
1729                         int dxn = xtab[k].di;
1730                         WT alpha = xtab[k].alpha;
1731                         WT t0 = buf[dxn] + S[sxn]*alpha;
1732                         WT t1 = buf[dxn+1] + S[sxn+1]*alpha;
1733                         WT t2 = buf[dxn+2] + S[sxn+2]*alpha;
1734                         buf[dxn] = t0; buf[dxn+1] = t1; buf[dxn+2] = t2;
1735                     }
1736                 else if( cn == 4 )
1737                 {
1738                     for( k = 0; k < xtab_size; k++ )
1739                     {
1740                         int sxn = xtab[k].si;
1741                         int dxn = xtab[k].di;
1742                         WT alpha = xtab[k].alpha;
1743                         WT t0 = buf[dxn] + S[sxn]*alpha;
1744                         WT t1 = buf[dxn+1] + S[sxn+1]*alpha;
1745                         buf[dxn] = t0; buf[dxn+1] = t1;
1746                         t0 = buf[dxn+2] + S[sxn+2]*alpha;
1747                         t1 = buf[dxn+3] + S[sxn+3]*alpha;
1748                         buf[dxn+2] = t0; buf[dxn+3] = t1;
1749                     }
1750                 }
1751                 else
1752                 {
1753                     for( k = 0; k < xtab_size; k++ )
1754                     {
1755                         int sxn = xtab[k].si;
1756                         int dxn = xtab[k].di;
1757                         WT alpha = xtab[k].alpha;
1758                         for( int c = 0; c < cn; c++ )
1759                             buf[dxn + c] += S[sxn + c]*alpha;
1760                     }
1761                 }
1762             }
1763
1764             if( dy != prev_dy )
1765             {
1766                 T* D = (T*)(dst->data + dst->step*prev_dy);
1767
1768                 for( dx = 0; dx < dsize.width; dx++ )
1769                 {
1770                     D[dx] = saturate_cast<T>(sum[dx]);
1771                     sum[dx] = beta*buf[dx];
1772                 }
1773                 prev_dy = dy;
1774             }
1775             else
1776             {
1777                 for( dx = 0; dx < dsize.width; dx++ )
1778                     sum[dx] += beta*buf[dx];
1779             }
1780         }
1781
1782         {
1783         T* D = (T*)(dst->data + dst->step*prev_dy);
1784         for( dx = 0; dx < dsize.width; dx++ )
1785             D[dx] = saturate_cast<T>(sum[dx]);
1786         }
1787     }
1788
1789 private:
1790     const Mat* src;
1791     Mat* dst;
1792     const DecimateAlpha* xtab0;
1793     const DecimateAlpha* ytab;
1794     int xtab_size0, ytab_size;
1795     const int* tabofs;
1796 };
1797
1798
1799 template <typename T, typename WT>
1800 static void resizeArea_( const Mat& src, Mat& dst,
1801                          const DecimateAlpha* xtab, int xtab_size,
1802                          const DecimateAlpha* ytab, int ytab_size,
1803                          const int* tabofs )
1804 {
1805     parallel_for_(Range(0, dst.rows),
1806                  ResizeArea_Invoker<T, WT>(src, dst, xtab, xtab_size, ytab, ytab_size, tabofs),
1807                  dst.total()/((double)(1 << 16)));
1808 }
1809
1810
1811 typedef void (*ResizeFunc)( const Mat& src, Mat& dst,
1812                             const int* xofs, const void* alpha,
1813                             const int* yofs, const void* beta,
1814                             int xmin, int xmax, int ksize );
1815
1816 typedef void (*ResizeAreaFastFunc)( const Mat& src, Mat& dst,
1817                                     const int* ofs, const int *xofs,
1818                                     int scale_x, int scale_y );
1819
1820 typedef void (*ResizeAreaFunc)( const Mat& src, Mat& dst,
1821                                 const DecimateAlpha* xtab, int xtab_size,
1822                                 const DecimateAlpha* ytab, int ytab_size,
1823                                 const int* yofs);
1824
1825
1826 static int computeResizeAreaTab( int ssize, int dsize, int cn, double scale, DecimateAlpha* tab )
1827 {
1828     int k = 0;
1829     for(int dx = 0; dx < dsize; dx++ )
1830     {
1831         double fsx1 = dx * scale;
1832         double fsx2 = fsx1 + scale;
1833         double cellWidth = std::min(scale, ssize - fsx1);
1834
1835         int sx1 = cvCeil(fsx1), sx2 = cvFloor(fsx2);
1836
1837         sx2 = std::min(sx2, ssize - 1);
1838         sx1 = std::min(sx1, sx2);
1839
1840         if( sx1 - fsx1 > 1e-3 )
1841         {
1842             assert( k < ssize*2 );
1843             tab[k].di = dx * cn;
1844             tab[k].si = (sx1 - 1) * cn;
1845             tab[k++].alpha = (float)((sx1 - fsx1) / cellWidth);
1846         }
1847
1848         for(int sx = sx1; sx < sx2; sx++ )
1849         {
1850             assert( k < ssize*2 );
1851             tab[k].di = dx * cn;
1852             tab[k].si = sx * cn;
1853             tab[k++].alpha = float(1.0 / cellWidth);
1854         }
1855
1856         if( fsx2 - sx2 > 1e-3 )
1857         {
1858             assert( k < ssize*2 );
1859             tab[k].di = dx * cn;
1860             tab[k].si = sx2 * cn;
1861             tab[k++].alpha = (float)(std::min(std::min(fsx2 - sx2, 1.), cellWidth) / cellWidth);
1862         }
1863     }
1864     return k;
1865 }
1866
1867 #define CHECK_FUNC(FUNC) if( FUNC==0 ) { *ok = false; return;}
1868 #define CHECK_STATUS(STATUS) if( STATUS!=ippStsNoErr ) { *ok = false; return;}
1869
1870 #define SET_IPP_RESIZE_LINEAR_FUNC_PTR(TYPE, CN) \
1871     func = (ippiResizeFunc)ippiResizeLinear_##TYPE##_##CN##R; CHECK_FUNC(func);\
1872     status = ippiResizeGetSize_##TYPE(srcSize, dstSize, (IppiInterpolationType)mode, 0, &specSize, &initSize); CHECK_STATUS(status)\
1873     specBuf.allocate(specSize);\
1874     pSpec = (uchar*)specBuf;\
1875     status = ippiResizeLinearInit_##TYPE(srcSize, dstSize, (IppiResizeSpec_32f*)pSpec); CHECK_STATUS(status);
1876
1877 #define SET_IPP_RESIZE_LINEAR_FUNC_64_PTR(TYPE, CN) \
1878     func = (ippiResizeFunc)ippiResizeLinear_##TYPE##_##CN##R; CHECK_FUNC(func);\
1879     status = ippiResizeGetSize_##TYPE(srcSize, dstSize, (IppiInterpolationType)mode, 0, &specSize, &initSize); CHECK_STATUS(status)\
1880     specBuf.allocate(specSize);\
1881     pSpec = (uchar*)specBuf;\
1882     status = ippiResizeLinearInit_##TYPE(srcSize, dstSize, (IppiResizeSpec_64f*)pSpec); CHECK_STATUS(status);
1883
1884 #define SET_IPP_RESIZE_CUBIC_FUNC_PTR(TYPE, CN) \
1885     func = (ippiResizeFunc)ippiResizeCubic_##TYPE##_##CN##R; CHECK_FUNC(func);\
1886     status = ippiResizeGetSize_##TYPE(srcSize, dstSize, (IppiInterpolationType)mode, 0, &specSize, &initSize); CHECK_STATUS(status)\
1887     specBuf.allocate(specSize);\
1888     pSpec = (uchar*)specBuf;\
1889     status = ippiResizeCubicInit_##TYPE(srcSize, dstSize,  0.f, 0.75f, (IppiResizeSpec_32f*)pSpec, pInit); CHECK_STATUS(status);
1890
1891 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR*100 + IPP_VERSION_MINOR >= 701)
1892 class IPPresizeInvoker :
1893     public ParallelLoopBody
1894 {
1895 public:
1896     IPPresizeInvoker(Mat &_src, Mat &_dst, double _inv_scale_x, double _inv_scale_y, int _mode, bool *_ok) :
1897       ParallelLoopBody(), src(_src), dst(_dst), inv_scale_x(_inv_scale_x), inv_scale_y(_inv_scale_y), mode(_mode), ok(_ok)
1898       {
1899           IppStatus status = ippStsNotSupportedModeErr;
1900           *ok = true;
1901           IppiSize srcSize, dstSize;
1902           int type = src.type();
1903           int specSize = 0, initSize = 0;
1904           srcSize.width  = src.cols;
1905           srcSize.height = src.rows;
1906           dstSize.width  = dst.cols;
1907           dstSize.height = dst.rows;
1908
1909           if (mode == (int)ippLinear)
1910           {
1911               switch (type)
1912               {
1913               case CV_8UC1:  SET_IPP_RESIZE_LINEAR_FUNC_PTR(8u,C1);  break;
1914               case CV_8UC3:  SET_IPP_RESIZE_LINEAR_FUNC_PTR(8u,C3);  break;
1915               case CV_8UC4:  SET_IPP_RESIZE_LINEAR_FUNC_PTR(8u,C4);  break;
1916               case CV_16UC1: SET_IPP_RESIZE_LINEAR_FUNC_PTR(16u,C1); break;
1917               case CV_16UC3: SET_IPP_RESIZE_LINEAR_FUNC_PTR(16u,C3); break;
1918               case CV_16UC4: SET_IPP_RESIZE_LINEAR_FUNC_PTR(16u,C4); break;
1919               case CV_16SC1: SET_IPP_RESIZE_LINEAR_FUNC_PTR(16s,C1); break;
1920               case CV_16SC3: SET_IPP_RESIZE_LINEAR_FUNC_PTR(16s,C3); break;
1921               case CV_16SC4: SET_IPP_RESIZE_LINEAR_FUNC_PTR(16s,C4); break;
1922               case CV_32FC1: SET_IPP_RESIZE_LINEAR_FUNC_PTR(32f,C1); break;
1923               case CV_32FC3: SET_IPP_RESIZE_LINEAR_FUNC_PTR(32f,C3); break;
1924               case CV_32FC4: SET_IPP_RESIZE_LINEAR_FUNC_PTR(32f,C4); break;
1925               case CV_64FC1: SET_IPP_RESIZE_LINEAR_FUNC_64_PTR(64f,C1); break;
1926               case CV_64FC3: SET_IPP_RESIZE_LINEAR_FUNC_64_PTR(64f,C3); break;
1927               case CV_64FC4: SET_IPP_RESIZE_LINEAR_FUNC_64_PTR(64f,C4); break;
1928               default: { *ok = false; return;} break;
1929               }
1930           }
1931           else if (mode == (int)ippCubic)
1932           {
1933               AutoBuffer<uchar> buf(initSize);
1934               uchar* pInit = (uchar*)buf;
1935               switch (type)
1936               {
1937               case CV_8UC1:  SET_IPP_RESIZE_CUBIC_FUNC_PTR(8u,C1);  break;
1938               case CV_8UC3:  SET_IPP_RESIZE_CUBIC_FUNC_PTR(8u,C3);  break;
1939               case CV_8UC4:  SET_IPP_RESIZE_CUBIC_FUNC_PTR(8u,C4);  break;
1940               case CV_16UC1: SET_IPP_RESIZE_CUBIC_FUNC_PTR(16u,C1); break;
1941               case CV_16UC3: SET_IPP_RESIZE_CUBIC_FUNC_PTR(16u,C3); break;
1942               case CV_16UC4: SET_IPP_RESIZE_CUBIC_FUNC_PTR(16u,C4); break;
1943               case CV_16SC1: SET_IPP_RESIZE_CUBIC_FUNC_PTR(16s,C1); break;
1944               case CV_16SC3: SET_IPP_RESIZE_CUBIC_FUNC_PTR(16s,C3); break;
1945               case CV_16SC4: SET_IPP_RESIZE_CUBIC_FUNC_PTR(16s,C4); break;
1946               case CV_32FC1: SET_IPP_RESIZE_CUBIC_FUNC_PTR(32f,C1); break;
1947               case CV_32FC3: SET_IPP_RESIZE_CUBIC_FUNC_PTR(32f,C3); break;
1948               case CV_32FC4: SET_IPP_RESIZE_CUBIC_FUNC_PTR(32f,C4); break;
1949               default: { *ok = false; return;} break;
1950               }
1951           }
1952       }
1953
1954       ~IPPresizeInvoker()
1955       {
1956       }
1957
1958       virtual void operator() (const Range& range) const
1959       {
1960           if (*ok == false) return;
1961
1962           int cn = src.channels();
1963           int dsty = CV_IMIN(cvRound(range.start * inv_scale_y), dst.rows);
1964           int dstwidth  = CV_IMIN(cvRound(src.cols * inv_scale_x), dst.cols);
1965           int dstheight = CV_IMIN(cvRound(range.end * inv_scale_y), dst.rows);
1966
1967           IppiPoint dstOffset = { 0, dsty }, srcOffset = {0, 0};
1968           IppiSize  dstSize   = { dstwidth, dstheight - dsty };
1969           int bufsize = 0, itemSize = 0;
1970
1971           IppStatus status = ippStsNotSupportedModeErr;
1972
1973           switch (src.depth())
1974           {
1975           case CV_8U:
1976               itemSize = 1;
1977               status = ippiResizeGetBufferSize_8u((IppiResizeSpec_32f*)pSpec, dstSize, cn, &bufsize);
1978               CHECK_STATUS(status);
1979               status = ippiResizeGetSrcOffset_8u((IppiResizeSpec_32f*)pSpec, dstOffset, &srcOffset);
1980               break;
1981           case CV_16U:
1982               itemSize = 2;
1983               status = ippiResizeGetBufferSize_16u((IppiResizeSpec_32f*)pSpec, dstSize, cn, &bufsize);
1984               CHECK_STATUS(status);
1985               status = ippiResizeGetSrcOffset_16u((IppiResizeSpec_32f*)pSpec, dstOffset, &srcOffset);
1986               break;
1987           case CV_16S:
1988               itemSize = 2;
1989               status = ippiResizeGetBufferSize_16s((IppiResizeSpec_32f*)pSpec, dstSize, cn, &bufsize);
1990               CHECK_STATUS(status);
1991               status = ippiResizeGetSrcOffset_16s((IppiResizeSpec_32f*)pSpec, dstOffset, &srcOffset);
1992               break;
1993           case CV_32F:
1994               itemSize = 4;
1995               status = ippiResizeGetBufferSize_32f((IppiResizeSpec_32f*)pSpec, dstSize, cn, &bufsize);
1996               CHECK_STATUS(status);
1997               status = ippiResizeGetSrcOffset_32f((IppiResizeSpec_32f*)pSpec, dstOffset, &srcOffset);
1998               break;
1999           case CV_64F:
2000               itemSize = 8;
2001               status = ippiResizeGetBufferSize_64f((IppiResizeSpec_64f*)pSpec, dstSize, cn, &bufsize);
2002               CHECK_STATUS(status);
2003               status = ippiResizeGetSrcOffset_64f((IppiResizeSpec_64f*)pSpec, dstOffset, &srcOffset);
2004               break;
2005           }
2006
2007           CHECK_STATUS(status);
2008
2009           Ipp8u* pSrc = (Ipp8u*)src.data + (int)src.step[0] * srcOffset.y + srcOffset.x * cn * itemSize;
2010           Ipp8u* pDst = (Ipp8u*)dst.data + (int)dst.step[0] * dstOffset.y + dstOffset.x * cn * itemSize;
2011
2012           AutoBuffer<uchar> buf(bufsize + 64);
2013           uchar* bufptr = alignPtr((uchar*)buf, 32);
2014
2015           CHECK_FUNC(func);
2016
2017           if( func( pSrc, (int)src.step[0], pDst, (int)dst.step[0], dstOffset, dstSize, ippBorderRepl, 0, pSpec, bufptr ) < 0 )
2018               *ok = false;
2019       }
2020 private:
2021     Mat &src;
2022     Mat &dst;
2023     double inv_scale_x;
2024     double inv_scale_y;
2025     void *pSpec;
2026     AutoBuffer<uchar>   specBuf;
2027     int mode;
2028     ippiResizeFunc func;
2029     bool *ok;
2030     const IPPresizeInvoker& operator= (const IPPresizeInvoker&);
2031 };
2032 #endif
2033
2034 #ifdef HAVE_OPENCL
2035
2036 static void ocl_computeResizeAreaTabs(int ssize, int dsize, double scale, int * const map_tab,
2037                                           float * const alpha_tab, int * const ofs_tab)
2038 {
2039     int k = 0, dx = 0;
2040     for ( ; dx < dsize; dx++)
2041     {
2042         ofs_tab[dx] = k;
2043
2044         double fsx1 = dx * scale;
2045         double fsx2 = fsx1 + scale;
2046         double cellWidth = std::min(scale, ssize - fsx1);
2047
2048         int sx1 = cvCeil(fsx1), sx2 = cvFloor(fsx2);
2049
2050         sx2 = std::min(sx2, ssize - 1);
2051         sx1 = std::min(sx1, sx2);
2052
2053         if (sx1 - fsx1 > 1e-3)
2054         {
2055             map_tab[k] = sx1 - 1;
2056             alpha_tab[k++] = (float)((sx1 - fsx1) / cellWidth);
2057         }
2058
2059         for (int sx = sx1; sx < sx2; sx++)
2060         {
2061             map_tab[k] = sx;
2062             alpha_tab[k++] = float(1.0 / cellWidth);
2063         }
2064
2065         if (fsx2 - sx2 > 1e-3)
2066         {
2067             map_tab[k] = sx2;
2068             alpha_tab[k++] = (float)(std::min(std::min(fsx2 - sx2, 1.), cellWidth) / cellWidth);
2069         }
2070     }
2071     ofs_tab[dx] = k;
2072 }
2073
2074 static void ocl_computeResizeAreaFastTabs(int * dmap_tab, int * smap_tab, int scale, int dcols, int scol)
2075 {
2076     for (int i = 0; i < dcols; ++i)
2077         dmap_tab[i] = scale * i;
2078
2079     for (int i = 0, size = dcols * scale; i < size; ++i)
2080         smap_tab[i] = std::min(scol - 1, i);
2081 }
2082
2083 static bool ocl_resize( InputArray _src, OutputArray _dst, Size dsize,
2084                         double fx, double fy, int interpolation)
2085 {
2086     int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
2087
2088     double inv_fx = 1. / fx, inv_fy = 1. / fy;
2089     float inv_fxf = (float)inv_fx, inv_fyf = (float)inv_fy;
2090
2091     if( !(cn <= 4 &&
2092            (interpolation == INTER_NEAREST || interpolation == INTER_LINEAR ||
2093             (interpolation == INTER_AREA && inv_fx >= 1 && inv_fy >= 1) )) )
2094         return false;
2095
2096     UMat src = _src.getUMat();
2097     _dst.create(dsize, type);
2098     UMat dst = _dst.getUMat();
2099
2100     ocl::Kernel k;
2101     size_t globalsize[] = { dst.cols, dst.rows };
2102
2103     if (interpolation == INTER_LINEAR)
2104     {
2105         int wdepth = std::max(depth, CV_32S);
2106         int wtype = CV_MAKETYPE(wdepth, cn);
2107         char buf[2][32];
2108         k.create("resizeLN", ocl::imgproc::resize_oclsrc,
2109                  format("-D INTER_LINEAR -D depth=%d -D PIXTYPE=%s -D PIXTYPE1=%s "
2110                         "-D WORKTYPE=%s -D convertToWT=%s -D convertToDT=%s -D cn=%d",
2111                         depth, ocl::typeToStr(type), ocl::typeToStr(depth), ocl::typeToStr(wtype),
2112                         ocl::convertTypeStr(depth, wdepth, cn, buf[0]),
2113                         ocl::convertTypeStr(wdepth, depth, cn, buf[1]),
2114                         cn));
2115     }
2116     else if (interpolation == INTER_NEAREST)
2117     {
2118         k.create("resizeNN", ocl::imgproc::resize_oclsrc,
2119                  format("-D INTER_NEAREST -D PIXTYPE=%s -D PIXTYPE1=%s -D cn=%d",
2120                         ocl::memopTypeToStr(type), ocl::memopTypeToStr(depth), cn));
2121     }
2122     else if (interpolation == INTER_AREA)
2123     {
2124         int iscale_x = saturate_cast<int>(inv_fx);
2125         int iscale_y = saturate_cast<int>(inv_fy);
2126         bool is_area_fast = std::abs(inv_fx - iscale_x) < DBL_EPSILON &&
2127                         std::abs(inv_fy - iscale_y) < DBL_EPSILON;
2128         int wdepth = std::max(depth, is_area_fast ? CV_32S : CV_32F);
2129         int wtype = CV_MAKE_TYPE(wdepth, cn);
2130
2131         char cvt[2][40];
2132         String buildOption = format("-D INTER_AREA -D PIXTYPE=%s -D PIXTYPE1=%s -D WTV=%s -D convertToWTV=%s -D cn=%d",
2133                                     ocl::typeToStr(type), ocl::typeToStr(depth), ocl::typeToStr(wtype),
2134                                     ocl::convertTypeStr(depth, wdepth, cn, cvt[0]), cn);
2135
2136         UMat alphaOcl, tabofsOcl, mapOcl;
2137         UMat dmap, smap;
2138
2139         if (is_area_fast)
2140         {
2141             int wdepth2 = std::max(CV_32F, depth), wtype2 = CV_MAKE_TYPE(wdepth2, cn);
2142             buildOption = buildOption + format(" -D convertToPIXTYPE=%s -D WT2V=%s -D convertToWT2V=%s -D INTER_AREA_FAST"
2143                                                " -D XSCALE=%d -D YSCALE=%d -D SCALE=%ff",
2144                                                ocl::convertTypeStr(wdepth2, depth, cn, cvt[0]),
2145                                                ocl::typeToStr(wtype2), ocl::convertTypeStr(wdepth, wdepth2, cn, cvt[1]),
2146                                   iscale_x, iscale_y, 1.0f / (iscale_x * iscale_y));
2147
2148             k.create("resizeAREA_FAST", ocl::imgproc::resize_oclsrc, buildOption);
2149             if (k.empty())
2150                 return false;
2151
2152             int smap_tab_size = dst.cols * iscale_x + dst.rows * iscale_y;
2153             AutoBuffer<int> dmap_tab(dst.cols + dst.rows), smap_tab(smap_tab_size);
2154             int * dxmap_tab = dmap_tab, * dymap_tab = dxmap_tab + dst.cols;
2155             int * sxmap_tab = smap_tab, * symap_tab = smap_tab + dst.cols * iscale_y;
2156
2157             ocl_computeResizeAreaFastTabs(dxmap_tab, sxmap_tab, iscale_x, dst.cols, src.cols);
2158             ocl_computeResizeAreaFastTabs(dymap_tab, symap_tab, iscale_y, dst.rows, src.rows);
2159
2160             Mat(1, dst.cols + dst.rows, CV_32SC1, (void *)dmap_tab).copyTo(dmap);
2161             Mat(1, smap_tab_size, CV_32SC1, (void *)smap_tab).copyTo(smap);
2162         }
2163         else
2164         {
2165             buildOption = buildOption + format(" -D convertToPIXTYPE=%s", ocl::convertTypeStr(wdepth, depth, cn, cvt[0]));
2166             k.create("resizeAREA", ocl::imgproc::resize_oclsrc, buildOption);
2167             if (k.empty())
2168                 return false;
2169
2170             Size ssize = src.size();
2171             int xytab_size = (ssize.width + ssize.height) << 1;
2172             int tabofs_size = dsize.height + dsize.width + 2;
2173
2174             AutoBuffer<int> _xymap_tab(xytab_size), _xyofs_tab(tabofs_size);
2175             AutoBuffer<float> _xyalpha_tab(xytab_size);
2176             int * xmap_tab = _xymap_tab, * ymap_tab = _xymap_tab + (ssize.width << 1);
2177             float * xalpha_tab = _xyalpha_tab, * yalpha_tab = _xyalpha_tab + (ssize.width << 1);
2178             int * xofs_tab = _xyofs_tab, * yofs_tab = _xyofs_tab + dsize.width + 1;
2179
2180             ocl_computeResizeAreaTabs(ssize.width, dsize.width, inv_fx, xmap_tab, xalpha_tab, xofs_tab);
2181             ocl_computeResizeAreaTabs(ssize.height, dsize.height, inv_fy, ymap_tab, yalpha_tab, yofs_tab);
2182
2183             // loading precomputed arrays to GPU
2184             Mat(1, xytab_size, CV_32FC1, (void *)_xyalpha_tab).copyTo(alphaOcl);
2185             Mat(1, xytab_size, CV_32SC1, (void *)_xymap_tab).copyTo(mapOcl);
2186             Mat(1, tabofs_size, CV_32SC1, (void *)_xyofs_tab).copyTo(tabofsOcl);
2187         }
2188
2189         ocl::KernelArg srcarg = ocl::KernelArg::ReadOnly(src), dstarg = ocl::KernelArg::WriteOnly(dst);
2190
2191         if (is_area_fast)
2192             k.args(srcarg, dstarg, ocl::KernelArg::PtrReadOnly(dmap), ocl::KernelArg::PtrReadOnly(smap));
2193         else
2194             k.args(srcarg, dstarg, inv_fxf, inv_fyf, ocl::KernelArg::PtrReadOnly(tabofsOcl),
2195                    ocl::KernelArg::PtrReadOnly(mapOcl), ocl::KernelArg::PtrReadOnly(alphaOcl));
2196
2197         return k.run(2, globalsize, NULL, false);
2198     }
2199
2200     if( k.empty() )
2201         return false;
2202     k.args(ocl::KernelArg::ReadOnly(src), ocl::KernelArg::WriteOnly(dst),
2203            (float)inv_fx, (float)inv_fy);
2204
2205     return k.run(2, globalsize, 0, false);
2206 }
2207
2208 #endif
2209
2210 }
2211
2212 //////////////////////////////////////////////////////////////////////////////////////////
2213
2214 void cv::resize( InputArray _src, OutputArray _dst, Size dsize,
2215                  double inv_scale_x, double inv_scale_y, int interpolation )
2216 {
2217     static ResizeFunc linear_tab[] =
2218     {
2219         resizeGeneric_<
2220             HResizeLinear<uchar, int, short,
2221                 INTER_RESIZE_COEF_SCALE,
2222                 HResizeLinearVec_8u32s>,
2223             VResizeLinear<uchar, int, short,
2224                 FixedPtCast<int, uchar, INTER_RESIZE_COEF_BITS*2>,
2225                 VResizeLinearVec_32s8u> >,
2226         0,
2227         resizeGeneric_<
2228             HResizeLinear<ushort, float, float, 1,
2229                 HResizeLinearVec_16u32f>,
2230             VResizeLinear<ushort, float, float, Cast<float, ushort>,
2231                 VResizeLinearVec_32f16u> >,
2232         resizeGeneric_<
2233             HResizeLinear<short, float, float, 1,
2234                 HResizeLinearVec_16s32f>,
2235             VResizeLinear<short, float, float, Cast<float, short>,
2236                 VResizeLinearVec_32f16s> >,
2237         0,
2238         resizeGeneric_<
2239             HResizeLinear<float, float, float, 1,
2240                 HResizeLinearVec_32f>,
2241             VResizeLinear<float, float, float, Cast<float, float>,
2242                 VResizeLinearVec_32f> >,
2243         resizeGeneric_<
2244             HResizeLinear<double, double, float, 1,
2245                 HResizeNoVec>,
2246             VResizeLinear<double, double, float, Cast<double, double>,
2247                 VResizeNoVec> >,
2248         0
2249     };
2250
2251     static ResizeFunc cubic_tab[] =
2252     {
2253         resizeGeneric_<
2254             HResizeCubic<uchar, int, short>,
2255             VResizeCubic<uchar, int, short,
2256                 FixedPtCast<int, uchar, INTER_RESIZE_COEF_BITS*2>,
2257                 VResizeCubicVec_32s8u> >,
2258         0,
2259         resizeGeneric_<
2260             HResizeCubic<ushort, float, float>,
2261             VResizeCubic<ushort, float, float, Cast<float, ushort>,
2262             VResizeCubicVec_32f16u> >,
2263         resizeGeneric_<
2264             HResizeCubic<short, float, float>,
2265             VResizeCubic<short, float, float, Cast<float, short>,
2266             VResizeCubicVec_32f16s> >,
2267         0,
2268         resizeGeneric_<
2269             HResizeCubic<float, float, float>,
2270             VResizeCubic<float, float, float, Cast<float, float>,
2271             VResizeCubicVec_32f> >,
2272         resizeGeneric_<
2273             HResizeCubic<double, double, float>,
2274             VResizeCubic<double, double, float, Cast<double, double>,
2275             VResizeNoVec> >,
2276         0
2277     };
2278
2279     static ResizeFunc lanczos4_tab[] =
2280     {
2281         resizeGeneric_<HResizeLanczos4<uchar, int, short>,
2282             VResizeLanczos4<uchar, int, short,
2283             FixedPtCast<int, uchar, INTER_RESIZE_COEF_BITS*2>,
2284             VResizeNoVec> >,
2285         0,
2286         resizeGeneric_<HResizeLanczos4<ushort, float, float>,
2287             VResizeLanczos4<ushort, float, float, Cast<float, ushort>,
2288             VResizeNoVec> >,
2289         resizeGeneric_<HResizeLanczos4<short, float, float>,
2290             VResizeLanczos4<short, float, float, Cast<float, short>,
2291             VResizeNoVec> >,
2292         0,
2293         resizeGeneric_<HResizeLanczos4<float, float, float>,
2294             VResizeLanczos4<float, float, float, Cast<float, float>,
2295             VResizeNoVec> >,
2296         resizeGeneric_<HResizeLanczos4<double, double, float>,
2297             VResizeLanczos4<double, double, float, Cast<double, double>,
2298             VResizeNoVec> >,
2299         0
2300     };
2301
2302     static ResizeAreaFastFunc areafast_tab[] =
2303     {
2304         resizeAreaFast_<uchar, int, ResizeAreaFastVec<uchar, ResizeAreaFastVec_SIMD_8u> >,
2305         0,
2306         resizeAreaFast_<ushort, float, ResizeAreaFastVec<ushort, ResizeAreaFastVec_SIMD_16u> >,
2307         resizeAreaFast_<short, float, ResizeAreaFastVec<short, ResizeAreaFastNoVec<short, float> > >,
2308         0,
2309         resizeAreaFast_<float, float, ResizeAreaFastNoVec<float, float> >,
2310         resizeAreaFast_<double, double, ResizeAreaFastNoVec<double, double> >,
2311         0
2312     };
2313
2314     static ResizeAreaFunc area_tab[] =
2315     {
2316         resizeArea_<uchar, float>, 0, resizeArea_<ushort, float>,
2317         resizeArea_<short, float>, 0, resizeArea_<float, float>,
2318         resizeArea_<double, double>, 0
2319     };
2320
2321     Size ssize = _src.size();
2322
2323     CV_Assert( ssize.area() > 0 );
2324     CV_Assert( dsize.area() > 0 || (inv_scale_x > 0 && inv_scale_y > 0) );
2325     if( dsize.area() == 0 )
2326     {
2327         dsize = Size(saturate_cast<int>(ssize.width*inv_scale_x),
2328                      saturate_cast<int>(ssize.height*inv_scale_y));
2329         CV_Assert( dsize.area() > 0 );
2330     }
2331     else
2332     {
2333         inv_scale_x = (double)dsize.width/ssize.width;
2334         inv_scale_y = (double)dsize.height/ssize.height;
2335     }
2336
2337     CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
2338                ocl_resize(_src, _dst, dsize, inv_scale_x, inv_scale_y, interpolation))
2339
2340     Mat src = _src.getMat();
2341     _dst.create(dsize, src.type());
2342     Mat dst = _dst.getMat();
2343
2344 #ifdef HAVE_TEGRA_OPTIMIZATION
2345     if (tegra::resize(src, dst, (float)inv_scale_x, (float)inv_scale_y, interpolation))
2346         return;
2347 #endif
2348
2349     int depth = src.depth(), cn = src.channels();
2350     double scale_x = 1./inv_scale_x, scale_y = 1./inv_scale_y;
2351     int k, sx, sy, dx, dy;
2352
2353 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR*100 + IPP_VERSION_MINOR >= 701)
2354 #define IPP_RESIZE_EPS    1.e-10
2355
2356     double ex = fabs((double)dsize.width/src.cols  - inv_scale_x)/inv_scale_x;
2357     double ey = fabs((double)dsize.height/src.rows - inv_scale_y)/inv_scale_y;
2358
2359     if ((ex < IPP_RESIZE_EPS && ey < IPP_RESIZE_EPS && depth != CV_64F) ||
2360         (ex == 0 && ey == 0 && depth == CV_64F))
2361     {
2362         int mode = 0;
2363         if (interpolation == INTER_LINEAR && src.rows >= 2 && src.cols >= 2)
2364         {
2365             mode = ippLinear;
2366         }
2367         else if (interpolation == INTER_CUBIC && src.rows >= 4 && src.cols >= 4)
2368         {
2369             mode = ippCubic;
2370         }
2371         if( mode != 0 && (cn == 1 || cn ==3 || cn == 4) &&
2372             (depth == CV_8U || depth == CV_16U || depth == CV_16S || depth == CV_32F ||
2373             (depth == CV_64F && mode == ippLinear)))
2374         {
2375             bool ok = true;
2376             Range range(0, src.rows);
2377             IPPresizeInvoker invoker(src, dst, inv_scale_x, inv_scale_y, mode, &ok);
2378             parallel_for_(range, invoker, dst.total()/(double)(1<<16));
2379             if( ok )
2380                 return;
2381         }
2382     }
2383 #undef IPP_RESIZE_EPS
2384 #endif
2385
2386     if( interpolation == INTER_NEAREST )
2387     {
2388         resizeNN( src, dst, inv_scale_x, inv_scale_y );
2389         return;
2390     }
2391
2392     {
2393         int iscale_x = saturate_cast<int>(scale_x);
2394         int iscale_y = saturate_cast<int>(scale_y);
2395
2396         bool is_area_fast = std::abs(scale_x - iscale_x) < DBL_EPSILON &&
2397                 std::abs(scale_y - iscale_y) < DBL_EPSILON;
2398
2399         // in case of scale_x && scale_y is equal to 2
2400         // INTER_AREA (fast) also is equal to INTER_LINEAR
2401         if( interpolation == INTER_LINEAR && is_area_fast && iscale_x == 2 && iscale_y == 2 )
2402             interpolation = INTER_AREA;
2403
2404         // true "area" interpolation is only implemented for the case (scale_x <= 1 && scale_y <= 1).
2405         // In other cases it is emulated using some variant of bilinear interpolation
2406         if( interpolation == INTER_AREA && scale_x >= 1 && scale_y >= 1 )
2407         {
2408             if( is_area_fast )
2409             {
2410                 int area = iscale_x*iscale_y;
2411                 size_t srcstep = src.step / src.elemSize1();
2412                 AutoBuffer<int> _ofs(area + dsize.width*cn);
2413                 int* ofs = _ofs;
2414                 int* xofs = ofs + area;
2415                 ResizeAreaFastFunc func = areafast_tab[depth];
2416                 CV_Assert( func != 0 );
2417
2418                 for( sy = 0, k = 0; sy < iscale_y; sy++ )
2419                     for( sx = 0; sx < iscale_x; sx++ )
2420                         ofs[k++] = (int)(sy*srcstep + sx*cn);
2421
2422                 for( dx = 0; dx < dsize.width; dx++ )
2423                 {
2424                     int j = dx * cn;
2425                     sx = iscale_x * j;
2426                     for( k = 0; k < cn; k++ )
2427                         xofs[j + k] = sx + k;
2428                 }
2429
2430                 func( src, dst, ofs, xofs, iscale_x, iscale_y );
2431                 return;
2432             }
2433
2434             ResizeAreaFunc func = area_tab[depth];
2435             CV_Assert( func != 0 && cn <= 4 );
2436
2437             AutoBuffer<DecimateAlpha> _xytab((ssize.width + ssize.height)*2);
2438             DecimateAlpha* xtab = _xytab, *ytab = xtab + ssize.width*2;
2439
2440             int xtab_size = computeResizeAreaTab(ssize.width, dsize.width, cn, scale_x, xtab);
2441             int ytab_size = computeResizeAreaTab(ssize.height, dsize.height, 1, scale_y, ytab);
2442
2443             AutoBuffer<int> _tabofs(dsize.height + 1);
2444             int* tabofs = _tabofs;
2445             for( k = 0, dy = 0; k < ytab_size; k++ )
2446             {
2447                 if( k == 0 || ytab[k].di != ytab[k-1].di )
2448                 {
2449                     assert( ytab[k].di == dy );
2450                     tabofs[dy++] = k;
2451                 }
2452             }
2453             tabofs[dy] = ytab_size;
2454
2455             func( src, dst, xtab, xtab_size, ytab, ytab_size, tabofs );
2456             return;
2457         }
2458     }
2459
2460     int xmin = 0, xmax = dsize.width, width = dsize.width*cn;
2461     bool area_mode = interpolation == INTER_AREA;
2462     bool fixpt = depth == CV_8U;
2463     float fx, fy;
2464     ResizeFunc func=0;
2465     int ksize=0, ksize2;
2466     if( interpolation == INTER_CUBIC )
2467         ksize = 4, func = cubic_tab[depth];
2468     else if( interpolation == INTER_LANCZOS4 )
2469         ksize = 8, func = lanczos4_tab[depth];
2470     else if( interpolation == INTER_LINEAR || interpolation == INTER_AREA )
2471         ksize = 2, func = linear_tab[depth];
2472     else
2473         CV_Error( CV_StsBadArg, "Unknown interpolation method" );
2474     ksize2 = ksize/2;
2475
2476     CV_Assert( func != 0 );
2477
2478     AutoBuffer<uchar> _buffer((width + dsize.height)*(sizeof(int) + sizeof(float)*ksize));
2479     int* xofs = (int*)(uchar*)_buffer;
2480     int* yofs = xofs + width;
2481     float* alpha = (float*)(yofs + dsize.height);
2482     short* ialpha = (short*)alpha;
2483     float* beta = alpha + width*ksize;
2484     short* ibeta = ialpha + width*ksize;
2485     float cbuf[MAX_ESIZE];
2486
2487     for( dx = 0; dx < dsize.width; dx++ )
2488     {
2489         if( !area_mode )
2490         {
2491             fx = (float)((dx+0.5)*scale_x - 0.5);
2492             sx = cvFloor(fx);
2493             fx -= sx;
2494         }
2495         else
2496         {
2497             sx = cvFloor(dx*scale_x);
2498             fx = (float)((dx+1) - (sx+1)*inv_scale_x);
2499             fx = fx <= 0 ? 0.f : fx - cvFloor(fx);
2500         }
2501
2502         if( sx < ksize2-1 )
2503         {
2504             xmin = dx+1;
2505             if( sx < 0 && (interpolation != INTER_CUBIC && interpolation != INTER_LANCZOS4))
2506                 fx = 0, sx = 0;
2507         }
2508
2509         if( sx + ksize2 >= ssize.width )
2510         {
2511             xmax = std::min( xmax, dx );
2512             if( sx >= ssize.width-1 && (interpolation != INTER_CUBIC && interpolation != INTER_LANCZOS4))
2513                 fx = 0, sx = ssize.width-1;
2514         }
2515
2516         for( k = 0, sx *= cn; k < cn; k++ )
2517             xofs[dx*cn + k] = sx + k;
2518
2519         if( interpolation == INTER_CUBIC )
2520             interpolateCubic( fx, cbuf );
2521         else if( interpolation == INTER_LANCZOS4 )
2522             interpolateLanczos4( fx, cbuf );
2523         else
2524         {
2525             cbuf[0] = 1.f - fx;
2526             cbuf[1] = fx;
2527         }
2528         if( fixpt )
2529         {
2530             for( k = 0; k < ksize; k++ )
2531                 ialpha[dx*cn*ksize + k] = saturate_cast<short>(cbuf[k]*INTER_RESIZE_COEF_SCALE);
2532             for( ; k < cn*ksize; k++ )
2533                 ialpha[dx*cn*ksize + k] = ialpha[dx*cn*ksize + k - ksize];
2534         }
2535         else
2536         {
2537             for( k = 0; k < ksize; k++ )
2538                 alpha[dx*cn*ksize + k] = cbuf[k];
2539             for( ; k < cn*ksize; k++ )
2540                 alpha[dx*cn*ksize + k] = alpha[dx*cn*ksize + k - ksize];
2541         }
2542     }
2543
2544     for( dy = 0; dy < dsize.height; dy++ )
2545     {
2546         if( !area_mode )
2547         {
2548             fy = (float)((dy+0.5)*scale_y - 0.5);
2549             sy = cvFloor(fy);
2550             fy -= sy;
2551         }
2552         else
2553         {
2554             sy = cvFloor(dy*scale_y);
2555             fy = (float)((dy+1) - (sy+1)*inv_scale_y);
2556             fy = fy <= 0 ? 0.f : fy - cvFloor(fy);
2557         }
2558
2559         yofs[dy] = sy;
2560         if( interpolation == INTER_CUBIC )
2561             interpolateCubic( fy, cbuf );
2562         else if( interpolation == INTER_LANCZOS4 )
2563             interpolateLanczos4( fy, cbuf );
2564         else
2565         {
2566             cbuf[0] = 1.f - fy;
2567             cbuf[1] = fy;
2568         }
2569
2570         if( fixpt )
2571         {
2572             for( k = 0; k < ksize; k++ )
2573                 ibeta[dy*ksize + k] = saturate_cast<short>(cbuf[k]*INTER_RESIZE_COEF_SCALE);
2574         }
2575         else
2576         {
2577             for( k = 0; k < ksize; k++ )
2578                 beta[dy*ksize + k] = cbuf[k];
2579         }
2580     }
2581
2582     func( src, dst, xofs, fixpt ? (void*)ialpha : (void*)alpha, yofs,
2583           fixpt ? (void*)ibeta : (void*)beta, xmin, xmax, ksize );
2584 }
2585
2586
2587 /****************************************************************************************\
2588 *                       General warping (affine, perspective, remap)                     *
2589 \****************************************************************************************/
2590
2591 namespace cv
2592 {
2593
2594 template<typename T>
2595 static void remapNearest( const Mat& _src, Mat& _dst, const Mat& _xy,
2596                           int borderType, const Scalar& _borderValue )
2597 {
2598     Size ssize = _src.size(), dsize = _dst.size();
2599     int cn = _src.channels();
2600     const T* S0 = (const T*)_src.data;
2601     size_t sstep = _src.step/sizeof(S0[0]);
2602     Scalar_<T> cval(saturate_cast<T>(_borderValue[0]),
2603         saturate_cast<T>(_borderValue[1]),
2604         saturate_cast<T>(_borderValue[2]),
2605         saturate_cast<T>(_borderValue[3]));
2606     int dx, dy;
2607
2608     unsigned width1 = ssize.width, height1 = ssize.height;
2609
2610     if( _dst.isContinuous() && _xy.isContinuous() )
2611     {
2612         dsize.width *= dsize.height;
2613         dsize.height = 1;
2614     }
2615
2616     for( dy = 0; dy < dsize.height; dy++ )
2617     {
2618         T* D = (T*)(_dst.data + _dst.step*dy);
2619         const short* XY = (const short*)(_xy.data + _xy.step*dy);
2620
2621         if( cn == 1 )
2622         {
2623             for( dx = 0; dx < dsize.width; dx++ )
2624             {
2625                 int sx = XY[dx*2], sy = XY[dx*2+1];
2626                 if( (unsigned)sx < width1 && (unsigned)sy < height1 )
2627                     D[dx] = S0[sy*sstep + sx];
2628                 else
2629                 {
2630                     if( borderType == BORDER_REPLICATE )
2631                     {
2632                         sx = clip(sx, 0, ssize.width);
2633                         sy = clip(sy, 0, ssize.height);
2634                         D[dx] = S0[sy*sstep + sx];
2635                     }
2636                     else if( borderType == BORDER_CONSTANT )
2637                         D[dx] = cval[0];
2638                     else if( borderType != BORDER_TRANSPARENT )
2639                     {
2640                         sx = borderInterpolate(sx, ssize.width, borderType);
2641                         sy = borderInterpolate(sy, ssize.height, borderType);
2642                         D[dx] = S0[sy*sstep + sx];
2643                     }
2644                 }
2645             }
2646         }
2647         else
2648         {
2649             for( dx = 0; dx < dsize.width; dx++, D += cn )
2650             {
2651                 int sx = XY[dx*2], sy = XY[dx*2+1], k;
2652                 const T *S;
2653                 if( (unsigned)sx < width1 && (unsigned)sy < height1 )
2654                 {
2655                     if( cn == 3 )
2656                     {
2657                         S = S0 + sy*sstep + sx*3;
2658                         D[0] = S[0], D[1] = S[1], D[2] = S[2];
2659                     }
2660                     else if( cn == 4 )
2661                     {
2662                         S = S0 + sy*sstep + sx*4;
2663                         D[0] = S[0], D[1] = S[1], D[2] = S[2], D[3] = S[3];
2664                     }
2665                     else
2666                     {
2667                         S = S0 + sy*sstep + sx*cn;
2668                         for( k = 0; k < cn; k++ )
2669                             D[k] = S[k];
2670                     }
2671                 }
2672                 else if( borderType != BORDER_TRANSPARENT )
2673                 {
2674                     if( borderType == BORDER_REPLICATE )
2675                     {
2676                         sx = clip(sx, 0, ssize.width);
2677                         sy = clip(sy, 0, ssize.height);
2678                         S = S0 + sy*sstep + sx*cn;
2679                     }
2680                     else if( borderType == BORDER_CONSTANT )
2681                         S = &cval[0];
2682                     else
2683                     {
2684                         sx = borderInterpolate(sx, ssize.width, borderType);
2685                         sy = borderInterpolate(sy, ssize.height, borderType);
2686                         S = S0 + sy*sstep + sx*cn;
2687                     }
2688                     for( k = 0; k < cn; k++ )
2689                         D[k] = S[k];
2690                 }
2691             }
2692         }
2693     }
2694 }
2695
2696
2697 struct RemapNoVec
2698 {
2699     int operator()( const Mat&, void*, const short*, const ushort*,
2700                     const void*, int ) const { return 0; }
2701 };
2702
2703 #if CV_SSE2
2704
2705 struct RemapVec_8u
2706 {
2707     int operator()( const Mat& _src, void* _dst, const short* XY,
2708                     const ushort* FXY, const void* _wtab, int width ) const
2709     {
2710         int cn = _src.channels(), x = 0, sstep = (int)_src.step;
2711
2712         if( (cn != 1 && cn != 3 && cn != 4) || !checkHardwareSupport(CV_CPU_SSE2) ||
2713             sstep > 0x8000 )
2714             return 0;
2715
2716         const uchar *S0 = _src.data, *S1 = _src.data + _src.step;
2717         const short* wtab = cn == 1 ? (const short*)_wtab : &BilinearTab_iC4[0][0][0];
2718         uchar* D = (uchar*)_dst;
2719         __m128i delta = _mm_set1_epi32(INTER_REMAP_COEF_SCALE/2);
2720         __m128i xy2ofs = _mm_set1_epi32(cn + (sstep << 16));
2721         __m128i z = _mm_setzero_si128();
2722         int CV_DECL_ALIGNED(16) iofs0[4], iofs1[4];
2723
2724         if( cn == 1 )
2725         {
2726             for( ; x <= width - 8; x += 8 )
2727             {
2728                 __m128i xy0 = _mm_loadu_si128( (const __m128i*)(XY + x*2));
2729                 __m128i xy1 = _mm_loadu_si128( (const __m128i*)(XY + x*2 + 8));
2730                 __m128i v0, v1, v2, v3, a0, a1, b0, b1;
2731                 unsigned i0, i1;
2732
2733                 xy0 = _mm_madd_epi16( xy0, xy2ofs );
2734                 xy1 = _mm_madd_epi16( xy1, xy2ofs );
2735                 _mm_store_si128( (__m128i*)iofs0, xy0 );
2736                 _mm_store_si128( (__m128i*)iofs1, xy1 );
2737
2738                 i0 = *(ushort*)(S0 + iofs0[0]) + (*(ushort*)(S0 + iofs0[1]) << 16);
2739                 i1 = *(ushort*)(S0 + iofs0[2]) + (*(ushort*)(S0 + iofs0[3]) << 16);
2740                 v0 = _mm_unpacklo_epi32(_mm_cvtsi32_si128(i0), _mm_cvtsi32_si128(i1));
2741                 i0 = *(ushort*)(S1 + iofs0[0]) + (*(ushort*)(S1 + iofs0[1]) << 16);
2742                 i1 = *(ushort*)(S1 + iofs0[2]) + (*(ushort*)(S1 + iofs0[3]) << 16);
2743                 v1 = _mm_unpacklo_epi32(_mm_cvtsi32_si128(i0), _mm_cvtsi32_si128(i1));
2744                 v0 = _mm_unpacklo_epi8(v0, z);
2745                 v1 = _mm_unpacklo_epi8(v1, z);
2746
2747                 a0 = _mm_unpacklo_epi32(_mm_loadl_epi64((__m128i*)(wtab+FXY[x]*4)),
2748                                         _mm_loadl_epi64((__m128i*)(wtab+FXY[x+1]*4)));
2749                 a1 = _mm_unpacklo_epi32(_mm_loadl_epi64((__m128i*)(wtab+FXY[x+2]*4)),
2750                                         _mm_loadl_epi64((__m128i*)(wtab+FXY[x+3]*4)));
2751                 b0 = _mm_unpacklo_epi64(a0, a1);
2752                 b1 = _mm_unpackhi_epi64(a0, a1);
2753                 v0 = _mm_madd_epi16(v0, b0);
2754                 v1 = _mm_madd_epi16(v1, b1);
2755                 v0 = _mm_add_epi32(_mm_add_epi32(v0, v1), delta);
2756
2757                 i0 = *(ushort*)(S0 + iofs1[0]) + (*(ushort*)(S0 + iofs1[1]) << 16);
2758                 i1 = *(ushort*)(S0 + iofs1[2]) + (*(ushort*)(S0 + iofs1[3]) << 16);
2759                 v2 = _mm_unpacklo_epi32(_mm_cvtsi32_si128(i0), _mm_cvtsi32_si128(i1));
2760                 i0 = *(ushort*)(S1 + iofs1[0]) + (*(ushort*)(S1 + iofs1[1]) << 16);
2761                 i1 = *(ushort*)(S1 + iofs1[2]) + (*(ushort*)(S1 + iofs1[3]) << 16);
2762                 v3 = _mm_unpacklo_epi32(_mm_cvtsi32_si128(i0), _mm_cvtsi32_si128(i1));
2763                 v2 = _mm_unpacklo_epi8(v2, z);
2764                 v3 = _mm_unpacklo_epi8(v3, z);
2765
2766                 a0 = _mm_unpacklo_epi32(_mm_loadl_epi64((__m128i*)(wtab+FXY[x+4]*4)),
2767                                         _mm_loadl_epi64((__m128i*)(wtab+FXY[x+5]*4)));
2768                 a1 = _mm_unpacklo_epi32(_mm_loadl_epi64((__m128i*)(wtab+FXY[x+6]*4)),
2769                                         _mm_loadl_epi64((__m128i*)(wtab+FXY[x+7]*4)));
2770                 b0 = _mm_unpacklo_epi64(a0, a1);
2771                 b1 = _mm_unpackhi_epi64(a0, a1);
2772                 v2 = _mm_madd_epi16(v2, b0);
2773                 v3 = _mm_madd_epi16(v3, b1);
2774                 v2 = _mm_add_epi32(_mm_add_epi32(v2, v3), delta);
2775
2776                 v0 = _mm_srai_epi32(v0, INTER_REMAP_COEF_BITS);
2777                 v2 = _mm_srai_epi32(v2, INTER_REMAP_COEF_BITS);
2778                 v0 = _mm_packus_epi16(_mm_packs_epi32(v0, v2), z);
2779                 _mm_storel_epi64( (__m128i*)(D + x), v0 );
2780             }
2781         }
2782         else if( cn == 3 )
2783         {
2784             for( ; x <= width - 5; x += 4, D += 12 )
2785             {
2786                 __m128i xy0 = _mm_loadu_si128( (const __m128i*)(XY + x*2));
2787                 __m128i u0, v0, u1, v1;
2788
2789                 xy0 = _mm_madd_epi16( xy0, xy2ofs );
2790                 _mm_store_si128( (__m128i*)iofs0, xy0 );
2791                 const __m128i *w0, *w1;
2792                 w0 = (const __m128i*)(wtab + FXY[x]*16);
2793                 w1 = (const __m128i*)(wtab + FXY[x+1]*16);
2794
2795                 u0 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S0 + iofs0[0])),
2796                                        _mm_cvtsi32_si128(*(int*)(S0 + iofs0[0] + 3)));
2797                 v0 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S1 + iofs0[0])),
2798                                        _mm_cvtsi32_si128(*(int*)(S1 + iofs0[0] + 3)));
2799                 u1 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S0 + iofs0[1])),
2800                                        _mm_cvtsi32_si128(*(int*)(S0 + iofs0[1] + 3)));
2801                 v1 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S1 + iofs0[1])),
2802                                        _mm_cvtsi32_si128(*(int*)(S1 + iofs0[1] + 3)));
2803                 u0 = _mm_unpacklo_epi8(u0, z);
2804                 v0 = _mm_unpacklo_epi8(v0, z);
2805                 u1 = _mm_unpacklo_epi8(u1, z);
2806                 v1 = _mm_unpacklo_epi8(v1, z);
2807                 u0 = _mm_add_epi32(_mm_madd_epi16(u0, w0[0]), _mm_madd_epi16(v0, w0[1]));
2808                 u1 = _mm_add_epi32(_mm_madd_epi16(u1, w1[0]), _mm_madd_epi16(v1, w1[1]));
2809                 u0 = _mm_srai_epi32(_mm_add_epi32(u0, delta), INTER_REMAP_COEF_BITS);
2810                 u1 = _mm_srai_epi32(_mm_add_epi32(u1, delta), INTER_REMAP_COEF_BITS);
2811                 u0 = _mm_slli_si128(u0, 4);
2812                 u0 = _mm_packs_epi32(u0, u1);
2813                 u0 = _mm_packus_epi16(u0, u0);
2814                 _mm_storel_epi64((__m128i*)D, _mm_srli_si128(u0,1));
2815
2816                 w0 = (const __m128i*)(wtab + FXY[x+2]*16);
2817                 w1 = (const __m128i*)(wtab + FXY[x+3]*16);
2818
2819                 u0 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S0 + iofs0[2])),
2820                                        _mm_cvtsi32_si128(*(int*)(S0 + iofs0[2] + 3)));
2821                 v0 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S1 + iofs0[2])),
2822                                        _mm_cvtsi32_si128(*(int*)(S1 + iofs0[2] + 3)));
2823                 u1 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S0 + iofs0[3])),
2824                                        _mm_cvtsi32_si128(*(int*)(S0 + iofs0[3] + 3)));
2825                 v1 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S1 + iofs0[3])),
2826                                        _mm_cvtsi32_si128(*(int*)(S1 + iofs0[3] + 3)));
2827                 u0 = _mm_unpacklo_epi8(u0, z);
2828                 v0 = _mm_unpacklo_epi8(v0, z);
2829                 u1 = _mm_unpacklo_epi8(u1, z);
2830                 v1 = _mm_unpacklo_epi8(v1, z);
2831                 u0 = _mm_add_epi32(_mm_madd_epi16(u0, w0[0]), _mm_madd_epi16(v0, w0[1]));
2832                 u1 = _mm_add_epi32(_mm_madd_epi16(u1, w1[0]), _mm_madd_epi16(v1, w1[1]));
2833                 u0 = _mm_srai_epi32(_mm_add_epi32(u0, delta), INTER_REMAP_COEF_BITS);
2834                 u1 = _mm_srai_epi32(_mm_add_epi32(u1, delta), INTER_REMAP_COEF_BITS);
2835                 u0 = _mm_slli_si128(u0, 4);
2836                 u0 = _mm_packs_epi32(u0, u1);
2837                 u0 = _mm_packus_epi16(u0, u0);
2838                 _mm_storel_epi64((__m128i*)(D + 6), _mm_srli_si128(u0,1));
2839             }
2840         }
2841         else if( cn == 4 )
2842         {
2843             for( ; x <= width - 4; x += 4, D += 16 )
2844             {
2845                 __m128i xy0 = _mm_loadu_si128( (const __m128i*)(XY + x*2));
2846                 __m128i u0, v0, u1, v1;
2847
2848                 xy0 = _mm_madd_epi16( xy0, xy2ofs );
2849                 _mm_store_si128( (__m128i*)iofs0, xy0 );
2850                 const __m128i *w0, *w1;
2851                 w0 = (const __m128i*)(wtab + FXY[x]*16);
2852                 w1 = (const __m128i*)(wtab + FXY[x+1]*16);
2853
2854                 u0 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S0 + iofs0[0])),
2855                                        _mm_cvtsi32_si128(*(int*)(S0 + iofs0[0] + 4)));
2856                 v0 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S1 + iofs0[0])),
2857                                        _mm_cvtsi32_si128(*(int*)(S1 + iofs0[0] + 4)));
2858                 u1 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S0 + iofs0[1])),
2859                                        _mm_cvtsi32_si128(*(int*)(S0 + iofs0[1] + 4)));
2860                 v1 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S1 + iofs0[1])),
2861                                        _mm_cvtsi32_si128(*(int*)(S1 + iofs0[1] + 4)));
2862                 u0 = _mm_unpacklo_epi8(u0, z);
2863                 v0 = _mm_unpacklo_epi8(v0, z);
2864                 u1 = _mm_unpacklo_epi8(u1, z);
2865                 v1 = _mm_unpacklo_epi8(v1, z);
2866                 u0 = _mm_add_epi32(_mm_madd_epi16(u0, w0[0]), _mm_madd_epi16(v0, w0[1]));
2867                 u1 = _mm_add_epi32(_mm_madd_epi16(u1, w1[0]), _mm_madd_epi16(v1, w1[1]));
2868                 u0 = _mm_srai_epi32(_mm_add_epi32(u0, delta), INTER_REMAP_COEF_BITS);
2869                 u1 = _mm_srai_epi32(_mm_add_epi32(u1, delta), INTER_REMAP_COEF_BITS);
2870                 u0 = _mm_packs_epi32(u0, u1);
2871                 u0 = _mm_packus_epi16(u0, u0);
2872                 _mm_storel_epi64((__m128i*)D, u0);
2873
2874                 w0 = (const __m128i*)(wtab + FXY[x+2]*16);
2875                 w1 = (const __m128i*)(wtab + FXY[x+3]*16);
2876
2877                 u0 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S0 + iofs0[2])),
2878                                        _mm_cvtsi32_si128(*(int*)(S0 + iofs0[2] + 4)));
2879                 v0 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S1 + iofs0[2])),
2880                                        _mm_cvtsi32_si128(*(int*)(S1 + iofs0[2] + 4)));
2881                 u1 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S0 + iofs0[3])),
2882                                        _mm_cvtsi32_si128(*(int*)(S0 + iofs0[3] + 4)));
2883                 v1 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S1 + iofs0[3])),
2884                                        _mm_cvtsi32_si128(*(int*)(S1 + iofs0[3] + 4)));
2885                 u0 = _mm_unpacklo_epi8(u0, z);
2886                 v0 = _mm_unpacklo_epi8(v0, z);
2887                 u1 = _mm_unpacklo_epi8(u1, z);
2888                 v1 = _mm_unpacklo_epi8(v1, z);
2889                 u0 = _mm_add_epi32(_mm_madd_epi16(u0, w0[0]), _mm_madd_epi16(v0, w0[1]));
2890                 u1 = _mm_add_epi32(_mm_madd_epi16(u1, w1[0]), _mm_madd_epi16(v1, w1[1]));
2891                 u0 = _mm_srai_epi32(_mm_add_epi32(u0, delta), INTER_REMAP_COEF_BITS);
2892                 u1 = _mm_srai_epi32(_mm_add_epi32(u1, delta), INTER_REMAP_COEF_BITS);
2893                 u0 = _mm_packs_epi32(u0, u1);
2894                 u0 = _mm_packus_epi16(u0, u0);
2895                 _mm_storel_epi64((__m128i*)(D + 8), u0);
2896             }
2897         }
2898
2899         return x;
2900     }
2901 };
2902
2903 #else
2904
2905 typedef RemapNoVec RemapVec_8u;
2906
2907 #endif
2908
2909
2910 template<class CastOp, class VecOp, typename AT>
2911 static void remapBilinear( const Mat& _src, Mat& _dst, const Mat& _xy,
2912                            const Mat& _fxy, const void* _wtab,
2913                            int borderType, const Scalar& _borderValue )
2914 {
2915     typedef typename CastOp::rtype T;
2916     typedef typename CastOp::type1 WT;
2917     Size ssize = _src.size(), dsize = _dst.size();
2918     int cn = _src.channels();
2919     const AT* wtab = (const AT*)_wtab;
2920     const T* S0 = (const T*)_src.data;
2921     size_t sstep = _src.step/sizeof(S0[0]);
2922     Scalar_<T> cval(saturate_cast<T>(_borderValue[0]),
2923         saturate_cast<T>(_borderValue[1]),
2924         saturate_cast<T>(_borderValue[2]),
2925         saturate_cast<T>(_borderValue[3]));
2926     int dx, dy;
2927     CastOp castOp;
2928     VecOp vecOp;
2929
2930     unsigned width1 = std::max(ssize.width-1, 0), height1 = std::max(ssize.height-1, 0);
2931     CV_Assert( cn <= 4 && ssize.area() > 0 );
2932 #if CV_SSE2
2933     if( _src.type() == CV_8UC3 )
2934         width1 = std::max(ssize.width-2, 0);
2935 #endif
2936
2937     for( dy = 0; dy < dsize.height; dy++ )
2938     {
2939         T* D = (T*)(_dst.data + _dst.step*dy);
2940         const short* XY = (const short*)(_xy.data + _xy.step*dy);
2941         const ushort* FXY = (const ushort*)(_fxy.data + _fxy.step*dy);
2942         int X0 = 0;
2943         bool prevInlier = false;
2944
2945         for( dx = 0; dx <= dsize.width; dx++ )
2946         {
2947             bool curInlier = dx < dsize.width ?
2948                 (unsigned)XY[dx*2] < width1 &&
2949                 (unsigned)XY[dx*2+1] < height1 : !prevInlier;
2950             if( curInlier == prevInlier )
2951                 continue;
2952
2953             int X1 = dx;
2954             dx = X0;
2955             X0 = X1;
2956             prevInlier = curInlier;
2957
2958             if( !curInlier )
2959             {
2960                 int len = vecOp( _src, D, XY + dx*2, FXY + dx, wtab, X1 - dx );
2961                 D += len*cn;
2962                 dx += len;
2963
2964                 if( cn == 1 )
2965                 {
2966                     for( ; dx < X1; dx++, D++ )
2967                     {
2968                         int sx = XY[dx*2], sy = XY[dx*2+1];
2969                         const AT* w = wtab + FXY[dx]*4;
2970                         const T* S = S0 + sy*sstep + sx;
2971                         *D = castOp(WT(S[0]*w[0] + S[1]*w[1] + S[sstep]*w[2] + S[sstep+1]*w[3]));
2972                     }
2973                 }
2974                 else if( cn == 2 )
2975                     for( ; dx < X1; dx++, D += 2 )
2976                     {
2977                         int sx = XY[dx*2], sy = XY[dx*2+1];
2978                         const AT* w = wtab + FXY[dx]*4;
2979                         const T* S = S0 + sy*sstep + sx*2;
2980                         WT t0 = S[0]*w[0] + S[2]*w[1] + S[sstep]*w[2] + S[sstep+2]*w[3];
2981                         WT t1 = S[1]*w[0] + S[3]*w[1] + S[sstep+1]*w[2] + S[sstep+3]*w[3];
2982                         D[0] = castOp(t0); D[1] = castOp(t1);
2983                     }
2984                 else if( cn == 3 )
2985                     for( ; dx < X1; dx++, D += 3 )
2986                     {
2987                         int sx = XY[dx*2], sy = XY[dx*2+1];
2988                         const AT* w = wtab + FXY[dx]*4;
2989                         const T* S = S0 + sy*sstep + sx*3;
2990                         WT t0 = S[0]*w[0] + S[3]*w[1] + S[sstep]*w[2] + S[sstep+3]*w[3];
2991                         WT t1 = S[1]*w[0] + S[4]*w[1] + S[sstep+1]*w[2] + S[sstep+4]*w[3];
2992                         WT t2 = S[2]*w[0] + S[5]*w[1] + S[sstep+2]*w[2] + S[sstep+5]*w[3];
2993                         D[0] = castOp(t0); D[1] = castOp(t1); D[2] = castOp(t2);
2994                     }
2995                 else
2996                     for( ; dx < X1; dx++, D += 4 )
2997                     {
2998                         int sx = XY[dx*2], sy = XY[dx*2+1];
2999                         const AT* w = wtab + FXY[dx]*4;
3000                         const T* S = S0 + sy*sstep + sx*4;
3001                         WT t0 = S[0]*w[0] + S[4]*w[1] + S[sstep]*w[2] + S[sstep+4]*w[3];
3002                         WT t1 = S[1]*w[0] + S[5]*w[1] + S[sstep+1]*w[2] + S[sstep+5]*w[3];
3003                         D[0] = castOp(t0); D[1] = castOp(t1);
3004                         t0 = S[2]*w[0] + S[6]*w[1] + S[sstep+2]*w[2] + S[sstep+6]*w[3];
3005                         t1 = S[3]*w[0] + S[7]*w[1] + S[sstep+3]*w[2] + S[sstep+7]*w[3];
3006                         D[2] = castOp(t0); D[3] = castOp(t1);
3007                     }
3008             }
3009             else
3010             {
3011                 if( borderType == BORDER_TRANSPARENT && cn != 3 )
3012                 {
3013                     D += (X1 - dx)*cn;
3014                     dx = X1;
3015                     continue;
3016                 }
3017
3018                 if( cn == 1 )
3019                     for( ; dx < X1; dx++, D++ )
3020                     {
3021                         int sx = XY[dx*2], sy = XY[dx*2+1];
3022                         if( borderType == BORDER_CONSTANT &&
3023                             (sx >= ssize.width || sx+1 < 0 ||
3024                              sy >= ssize.height || sy+1 < 0) )
3025                         {
3026                             D[0] = cval[0];
3027                         }
3028                         else
3029                         {
3030                             int sx0, sx1, sy0, sy1;
3031                             T v0, v1, v2, v3;
3032                             const AT* w = wtab + FXY[dx]*4;
3033                             if( borderType == BORDER_REPLICATE )
3034                             {
3035                                 sx0 = clip(sx, 0, ssize.width);
3036                                 sx1 = clip(sx+1, 0, ssize.width);
3037                                 sy0 = clip(sy, 0, ssize.height);
3038                                 sy1 = clip(sy+1, 0, ssize.height);
3039                                 v0 = S0[sy0*sstep + sx0];
3040                                 v1 = S0[sy0*sstep + sx1];
3041                                 v2 = S0[sy1*sstep + sx0];
3042                                 v3 = S0[sy1*sstep + sx1];
3043                             }
3044                             else
3045                             {
3046                                 sx0 = borderInterpolate(sx, ssize.width, borderType);
3047                                 sx1 = borderInterpolate(sx+1, ssize.width, borderType);
3048                                 sy0 = borderInterpolate(sy, ssize.height, borderType);
3049                                 sy1 = borderInterpolate(sy+1, ssize.height, borderType);
3050                                 v0 = sx0 >= 0 && sy0 >= 0 ? S0[sy0*sstep + sx0] : cval[0];
3051                                 v1 = sx1 >= 0 && sy0 >= 0 ? S0[sy0*sstep + sx1] : cval[0];
3052                                 v2 = sx0 >= 0 && sy1 >= 0 ? S0[sy1*sstep + sx0] : cval[0];
3053                                 v3 = sx1 >= 0 && sy1 >= 0 ? S0[sy1*sstep + sx1] : cval[0];
3054                             }
3055                             D[0] = castOp(WT(v0*w[0] + v1*w[1] + v2*w[2] + v3*w[3]));
3056                         }
3057                     }
3058                 else
3059                     for( ; dx < X1; dx++, D += cn )
3060                     {
3061                         int sx = XY[dx*2], sy = XY[dx*2+1], k;
3062                         if( borderType == BORDER_CONSTANT &&
3063                             (sx >= ssize.width || sx+1 < 0 ||
3064                              sy >= ssize.height || sy+1 < 0) )
3065                         {
3066                             for( k = 0; k < cn; k++ )
3067                                 D[k] = cval[k];
3068                         }
3069                         else
3070                         {
3071                             int sx0, sx1, sy0, sy1;
3072                             const T *v0, *v1, *v2, *v3;
3073                             const AT* w = wtab + FXY[dx]*4;
3074                             if( borderType == BORDER_REPLICATE )
3075                             {
3076                                 sx0 = clip(sx, 0, ssize.width);
3077                                 sx1 = clip(sx+1, 0, ssize.width);
3078                                 sy0 = clip(sy, 0, ssize.height);
3079                                 sy1 = clip(sy+1, 0, ssize.height);
3080                                 v0 = S0 + sy0*sstep + sx0*cn;
3081                                 v1 = S0 + sy0*sstep + sx1*cn;
3082                                 v2 = S0 + sy1*sstep + sx0*cn;
3083                                 v3 = S0 + sy1*sstep + sx1*cn;
3084                             }
3085                             else if( borderType == BORDER_TRANSPARENT &&
3086                                 ((unsigned)sx >= (unsigned)(ssize.width-1) ||
3087                                 (unsigned)sy >= (unsigned)(ssize.height-1)))
3088                                 continue;
3089                             else
3090                             {
3091                                 sx0 = borderInterpolate(sx, ssize.width, borderType);
3092                                 sx1 = borderInterpolate(sx+1, ssize.width, borderType);
3093                                 sy0 = borderInterpolate(sy, ssize.height, borderType);
3094                                 sy1 = borderInterpolate(sy+1, ssize.height, borderType);
3095                                 v0 = sx0 >= 0 && sy0 >= 0 ? S0 + sy0*sstep + sx0*cn : &cval[0];
3096                                 v1 = sx1 >= 0 && sy0 >= 0 ? S0 + sy0*sstep + sx1*cn : &cval[0];
3097                                 v2 = sx0 >= 0 && sy1 >= 0 ? S0 + sy1*sstep + sx0*cn : &cval[0];
3098                                 v3 = sx1 >= 0 && sy1 >= 0 ? S0 + sy1*sstep + sx1*cn : &cval[0];
3099                             }
3100                             for( k = 0; k < cn; k++ )
3101                                 D[k] = castOp(WT(v0[k]*w[0] + v1[k]*w[1] + v2[k]*w[2] + v3[k]*w[3]));
3102                         }
3103                     }
3104             }
3105         }
3106     }
3107 }
3108
3109
3110 template<class CastOp, typename AT, int ONE>
3111 static void remapBicubic( const Mat& _src, Mat& _dst, const Mat& _xy,
3112                           const Mat& _fxy, const void* _wtab,
3113                           int borderType, const Scalar& _borderValue )
3114 {
3115     typedef typename CastOp::rtype T;
3116     typedef typename CastOp::type1 WT;
3117     Size ssize = _src.size(), dsize = _dst.size();
3118     int cn = _src.channels();
3119     const AT* wtab = (const AT*)_wtab;
3120     const T* S0 = (const T*)_src.data;
3121     size_t sstep = _src.step/sizeof(S0[0]);
3122     Scalar_<T> cval(saturate_cast<T>(_borderValue[0]),
3123         saturate_cast<T>(_borderValue[1]),
3124         saturate_cast<T>(_borderValue[2]),
3125         saturate_cast<T>(_borderValue[3]));
3126     int dx, dy;
3127     CastOp castOp;
3128     int borderType1 = borderType != BORDER_TRANSPARENT ? borderType : BORDER_REFLECT_101;
3129
3130     unsigned width1 = std::max(ssize.width-3, 0), height1 = std::max(ssize.height-3, 0);
3131
3132     if( _dst.isContinuous() && _xy.isContinuous() && _fxy.isContinuous() )
3133     {
3134         dsize.width *= dsize.height;
3135         dsize.height = 1;
3136     }
3137
3138     for( dy = 0; dy < dsize.height; dy++ )
3139     {
3140         T* D = (T*)(_dst.data + _dst.step*dy);
3141         const short* XY = (const short*)(_xy.data + _xy.step*dy);
3142         const ushort* FXY = (const ushort*)(_fxy.data + _fxy.step*dy);
3143
3144         for( dx = 0; dx < dsize.width; dx++, D += cn )
3145         {
3146             int sx = XY[dx*2]-1, sy = XY[dx*2+1]-1;
3147             const AT* w = wtab + FXY[dx]*16;
3148             int i, k;
3149             if( (unsigned)sx < width1 && (unsigned)sy < height1 )
3150             {
3151                 const T* S = S0 + sy*sstep + sx*cn;
3152                 for( k = 0; k < cn; k++ )
3153                 {
3154                     WT sum = S[0]*w[0] + S[cn]*w[1] + S[cn*2]*w[2] + S[cn*3]*w[3];
3155                     S += sstep;
3156                     sum += S[0]*w[4] + S[cn]*w[5] + S[cn*2]*w[6] + S[cn*3]*w[7];
3157                     S += sstep;
3158                     sum += S[0]*w[8] + S[cn]*w[9] + S[cn*2]*w[10] + S[cn*3]*w[11];
3159                     S += sstep;
3160                     sum += S[0]*w[12] + S[cn]*w[13] + S[cn*2]*w[14] + S[cn*3]*w[15];
3161                     S += 1 - sstep*3;
3162                     D[k] = castOp(sum);
3163                 }
3164             }
3165             else
3166             {
3167                 int x[4], y[4];
3168                 if( borderType == BORDER_TRANSPARENT &&
3169                     ((unsigned)(sx+1) >= (unsigned)ssize.width ||
3170                     (unsigned)(sy+1) >= (unsigned)ssize.height) )
3171                     continue;
3172
3173                 if( borderType1 == BORDER_CONSTANT &&
3174                     (sx >= ssize.width || sx+4 <= 0 ||
3175                     sy >= ssize.height || sy+4 <= 0))
3176                 {
3177                     for( k = 0; k < cn; k++ )
3178                         D[k] = cval[k];
3179                     continue;
3180                 }
3181
3182                 for( i = 0; i < 4; i++ )
3183                 {
3184                     x[i] = borderInterpolate(sx + i, ssize.width, borderType1)*cn;
3185                     y[i] = borderInterpolate(sy + i, ssize.height, borderType1);
3186                 }
3187
3188                 for( k = 0; k < cn; k++, S0++, w -= 16 )
3189                 {
3190                     WT cv = cval[k], sum = cv*ONE;
3191                     for( i = 0; i < 4; i++, w += 4 )
3192                     {
3193                         int yi = y[i];
3194                         const T* S = S0 + yi*sstep;
3195                         if( yi < 0 )
3196                             continue;
3197                         if( x[0] >= 0 )
3198                             sum += (S[x[0]] - cv)*w[0];
3199                         if( x[1] >= 0 )
3200                             sum += (S[x[1]] - cv)*w[1];
3201                         if( x[2] >= 0 )
3202                             sum += (S[x[2]] - cv)*w[2];
3203                         if( x[3] >= 0 )
3204                             sum += (S[x[3]] - cv)*w[3];
3205                     }
3206                     D[k] = castOp(sum);
3207                 }
3208                 S0 -= cn;
3209             }
3210         }
3211     }
3212 }
3213
3214
3215 template<class CastOp, typename AT, int ONE>
3216 static void remapLanczos4( const Mat& _src, Mat& _dst, const Mat& _xy,
3217                            const Mat& _fxy, const void* _wtab,
3218                            int borderType, const Scalar& _borderValue )
3219 {
3220     typedef typename CastOp::rtype T;
3221     typedef typename CastOp::type1 WT;
3222     Size ssize = _src.size(), dsize = _dst.size();
3223     int cn = _src.channels();
3224     const AT* wtab = (const AT*)_wtab;
3225     const T* S0 = (const T*)_src.data;
3226     size_t sstep = _src.step/sizeof(S0[0]);
3227     Scalar_<T> cval(saturate_cast<T>(_borderValue[0]),
3228         saturate_cast<T>(_borderValue[1]),
3229         saturate_cast<T>(_borderValue[2]),
3230         saturate_cast<T>(_borderValue[3]));
3231     int dx, dy;
3232     CastOp castOp;
3233     int borderType1 = borderType != BORDER_TRANSPARENT ? borderType : BORDER_REFLECT_101;
3234
3235     unsigned width1 = std::max(ssize.width-7, 0), height1 = std::max(ssize.height-7, 0);
3236
3237     if( _dst.isContinuous() && _xy.isContinuous() && _fxy.isContinuous() )
3238     {
3239         dsize.width *= dsize.height;
3240         dsize.height = 1;
3241     }
3242
3243     for( dy = 0; dy < dsize.height; dy++ )
3244     {
3245         T* D = (T*)(_dst.data + _dst.step*dy);
3246         const short* XY = (const short*)(_xy.data + _xy.step*dy);
3247         const ushort* FXY = (const ushort*)(_fxy.data + _fxy.step*dy);
3248
3249         for( dx = 0; dx < dsize.width; dx++, D += cn )
3250         {
3251             int sx = XY[dx*2]-3, sy = XY[dx*2+1]-3;
3252             const AT* w = wtab + FXY[dx]*64;
3253             const T* S = S0 + sy*sstep + sx*cn;
3254             int i, k;
3255             if( (unsigned)sx < width1 && (unsigned)sy < height1 )
3256             {
3257                 for( k = 0; k < cn; k++ )
3258                 {
3259                     WT sum = 0;
3260                     for( int r = 0; r < 8; r++, S += sstep, w += 8 )
3261                         sum += S[0]*w[0] + S[cn]*w[1] + S[cn*2]*w[2] + S[cn*3]*w[3] +
3262                             S[cn*4]*w[4] + S[cn*5]*w[5] + S[cn*6]*w[6] + S[cn*7]*w[7];
3263                     w -= 64;
3264                     S -= sstep*8 - 1;
3265                     D[k] = castOp(sum);
3266                 }
3267             }
3268             else
3269             {
3270                 int x[8], y[8];
3271                 if( borderType == BORDER_TRANSPARENT &&
3272                     ((unsigned)(sx+3) >= (unsigned)ssize.width ||
3273                     (unsigned)(sy+3) >= (unsigned)ssize.height) )
3274                     continue;
3275
3276                 if( borderType1 == BORDER_CONSTANT &&
3277                     (sx >= ssize.width || sx+8 <= 0 ||
3278                     sy >= ssize.height || sy+8 <= 0))
3279                 {
3280                     for( k = 0; k < cn; k++ )
3281                         D[k] = cval[k];
3282                     continue;
3283                 }
3284
3285                 for( i = 0; i < 8; i++ )
3286                 {
3287                     x[i] = borderInterpolate(sx + i, ssize.width, borderType1)*cn;
3288                     y[i] = borderInterpolate(sy + i, ssize.height, borderType1);
3289                 }
3290
3291                 for( k = 0; k < cn; k++, S0++, w -= 64 )
3292                 {
3293                     WT cv = cval[k], sum = cv*ONE;
3294                     for( i = 0; i < 8; i++, w += 8 )
3295                     {
3296                         int yi = y[i];
3297                         const T* S1 = S0 + yi*sstep;
3298                         if( yi < 0 )
3299                             continue;
3300                         if( x[0] >= 0 )
3301                             sum += (S1[x[0]] - cv)*w[0];
3302                         if( x[1] >= 0 )
3303                             sum += (S1[x[1]] - cv)*w[1];
3304                         if( x[2] >= 0 )
3305                             sum += (S1[x[2]] - cv)*w[2];
3306                         if( x[3] >= 0 )
3307                             sum += (S1[x[3]] - cv)*w[3];
3308                         if( x[4] >= 0 )
3309                             sum += (S1[x[4]] - cv)*w[4];
3310                         if( x[5] >= 0 )
3311                             sum += (S1[x[5]] - cv)*w[5];
3312                         if( x[6] >= 0 )
3313                             sum += (S1[x[6]] - cv)*w[6];
3314                         if( x[7] >= 0 )
3315                             sum += (S1[x[7]] - cv)*w[7];
3316                     }
3317                     D[k] = castOp(sum);
3318                 }
3319                 S0 -= cn;
3320             }
3321         }
3322     }
3323 }
3324
3325
3326 typedef void (*RemapNNFunc)(const Mat& _src, Mat& _dst, const Mat& _xy,
3327                             int borderType, const Scalar& _borderValue );
3328
3329 typedef void (*RemapFunc)(const Mat& _src, Mat& _dst, const Mat& _xy,
3330                           const Mat& _fxy, const void* _wtab,
3331                           int borderType, const Scalar& _borderValue);
3332
3333 class RemapInvoker :
3334     public ParallelLoopBody
3335 {
3336 public:
3337     RemapInvoker(const Mat& _src, Mat& _dst, const Mat *_m1,
3338                  const Mat *_m2, int _borderType, const Scalar &_borderValue,
3339                  int _planar_input, RemapNNFunc _nnfunc, RemapFunc _ifunc, const void *_ctab) :
3340         ParallelLoopBody(), src(&_src), dst(&_dst), m1(_m1), m2(_m2),
3341         borderType(_borderType), borderValue(_borderValue),
3342         planar_input(_planar_input), nnfunc(_nnfunc), ifunc(_ifunc), ctab(_ctab)
3343     {
3344     }
3345
3346     virtual void operator() (const Range& range) const
3347     {
3348         int x, y, x1, y1;
3349         const int buf_size = 1 << 14;
3350         int brows0 = std::min(128, dst->rows), map_depth = m1->depth();
3351         int bcols0 = std::min(buf_size/brows0, dst->cols);
3352         brows0 = std::min(buf_size/bcols0, dst->rows);
3353     #if CV_SSE2
3354         bool useSIMD = checkHardwareSupport(CV_CPU_SSE2);
3355     #endif
3356
3357         Mat _bufxy(brows0, bcols0, CV_16SC2), _bufa;
3358         if( !nnfunc )
3359             _bufa.create(brows0, bcols0, CV_16UC1);
3360
3361         for( y = range.start; y < range.end; y += brows0 )
3362         {
3363             for( x = 0; x < dst->cols; x += bcols0 )
3364             {
3365                 int brows = std::min(brows0, range.end - y);
3366                 int bcols = std::min(bcols0, dst->cols - x);
3367                 Mat dpart(*dst, Rect(x, y, bcols, brows));
3368                 Mat bufxy(_bufxy, Rect(0, 0, bcols, brows));
3369
3370                 if( nnfunc )
3371                 {
3372                     if( m1->type() == CV_16SC2 && !m2->data ) // the data is already in the right format
3373                         bufxy = (*m1)(Rect(x, y, bcols, brows));
3374                     else if( map_depth != CV_32F )
3375                     {
3376                         for( y1 = 0; y1 < brows; y1++ )
3377                         {
3378                             short* XY = (short*)(bufxy.data + bufxy.step*y1);
3379                             const short* sXY = (const short*)(m1->data + m1->step*(y+y1)) + x*2;
3380                             const ushort* sA = (const ushort*)(m2->data + m2->step*(y+y1)) + x;
3381
3382                             for( x1 = 0; x1 < bcols; x1++ )
3383                             {
3384                                 int a = sA[x1] & (INTER_TAB_SIZE2-1);
3385                                 XY[x1*2] = sXY[x1*2] + NNDeltaTab_i[a][0];
3386                                 XY[x1*2+1] = sXY[x1*2+1] + NNDeltaTab_i[a][1];
3387                             }
3388                         }
3389                     }
3390                     else if( !planar_input )
3391                         (*m1)(Rect(x, y, bcols, brows)).convertTo(bufxy, bufxy.depth());
3392                     else
3393                     {
3394                         for( y1 = 0; y1 < brows; y1++ )
3395                         {
3396                             short* XY = (short*)(bufxy.data + bufxy.step*y1);
3397                             const float* sX = (const float*)(m1->data + m1->step*(y+y1)) + x;
3398                             const float* sY = (const float*)(m2->data + m2->step*(y+y1)) + x;
3399                             x1 = 0;
3400
3401                         #if CV_SSE2
3402                             if( useSIMD )
3403                             {
3404                                 for( ; x1 <= bcols - 8; x1 += 8 )
3405                                 {
3406                                     __m128 fx0 = _mm_loadu_ps(sX + x1);
3407                                     __m128 fx1 = _mm_loadu_ps(sX + x1 + 4);
3408                                     __m128 fy0 = _mm_loadu_ps(sY + x1);
3409                                     __m128 fy1 = _mm_loadu_ps(sY + x1 + 4);
3410                                     __m128i ix0 = _mm_cvtps_epi32(fx0);
3411                                     __m128i ix1 = _mm_cvtps_epi32(fx1);
3412                                     __m128i iy0 = _mm_cvtps_epi32(fy0);
3413                                     __m128i iy1 = _mm_cvtps_epi32(fy1);
3414                                     ix0 = _mm_packs_epi32(ix0, ix1);
3415                                     iy0 = _mm_packs_epi32(iy0, iy1);
3416                                     ix1 = _mm_unpacklo_epi16(ix0, iy0);
3417                                     iy1 = _mm_unpackhi_epi16(ix0, iy0);
3418                                     _mm_storeu_si128((__m128i*)(XY + x1*2), ix1);
3419                                     _mm_storeu_si128((__m128i*)(XY + x1*2 + 8), iy1);
3420                                 }
3421                             }
3422                         #endif
3423
3424                             for( ; x1 < bcols; x1++ )
3425                             {
3426                                 XY[x1*2] = saturate_cast<short>(sX[x1]);
3427                                 XY[x1*2+1] = saturate_cast<short>(sY[x1]);
3428                             }
3429                         }
3430                     }
3431                     nnfunc( *src, dpart, bufxy, borderType, borderValue );
3432                     continue;
3433                 }
3434
3435                 Mat bufa(_bufa, Rect(0, 0, bcols, brows));
3436                 for( y1 = 0; y1 < brows; y1++ )
3437                 {
3438                     short* XY = (short*)(bufxy.data + bufxy.step*y1);
3439                     ushort* A = (ushort*)(bufa.data + bufa.step*y1);
3440
3441                     if( m1->type() == CV_16SC2 && (m2->type() == CV_16UC1 || m2->type() == CV_16SC1) )
3442                     {
3443                         bufxy = (*m1)(Rect(x, y, bcols, brows));
3444
3445                         const ushort* sA = (const ushort*)(m2->data + m2->step*(y+y1)) + x;
3446                         for( x1 = 0; x1 < bcols; x1++ )
3447                             A[x1] = (ushort)(sA[x1] & (INTER_TAB_SIZE2-1));
3448                     }
3449                     else if( planar_input )
3450                     {
3451                         const float* sX = (const float*)(m1->data + m1->step*(y+y1)) + x;
3452                         const float* sY = (const float*)(m2->data + m2->step*(y+y1)) + x;
3453
3454                         x1 = 0;
3455                     #if CV_SSE2
3456                         if( useSIMD )
3457                         {
3458                             __m128 scale = _mm_set1_ps((float)INTER_TAB_SIZE);
3459                             __m128i mask = _mm_set1_epi32(INTER_TAB_SIZE-1);
3460                             for( ; x1 <= bcols - 8; x1 += 8 )
3461                             {
3462                                 __m128 fx0 = _mm_loadu_ps(sX + x1);
3463                                 __m128 fx1 = _mm_loadu_ps(sX + x1 + 4);
3464                                 __m128 fy0 = _mm_loadu_ps(sY + x1);
3465                                 __m128 fy1 = _mm_loadu_ps(sY + x1 + 4);
3466                                 __m128i ix0 = _mm_cvtps_epi32(_mm_mul_ps(fx0, scale));
3467                                 __m128i ix1 = _mm_cvtps_epi32(_mm_mul_ps(fx1, scale));
3468                                 __m128i iy0 = _mm_cvtps_epi32(_mm_mul_ps(fy0, scale));
3469                                 __m128i iy1 = _mm_cvtps_epi32(_mm_mul_ps(fy1, scale));
3470                                 __m128i mx0 = _mm_and_si128(ix0, mask);
3471                                 __m128i mx1 = _mm_and_si128(ix1, mask);
3472                                 __m128i my0 = _mm_and_si128(iy0, mask);
3473                                 __m128i my1 = _mm_and_si128(iy1, mask);
3474                                 mx0 = _mm_packs_epi32(mx0, mx1);
3475                                 my0 = _mm_packs_epi32(my0, my1);
3476                                 my0 = _mm_slli_epi16(my0, INTER_BITS);
3477                                 mx0 = _mm_or_si128(mx0, my0);
3478                                 _mm_storeu_si128((__m128i*)(A + x1), mx0);
3479                                 ix0 = _mm_srai_epi32(ix0, INTER_BITS);
3480                                 ix1 = _mm_srai_epi32(ix1, INTER_BITS);
3481                                 iy0 = _mm_srai_epi32(iy0, INTER_BITS);
3482                                 iy1 = _mm_srai_epi32(iy1, INTER_BITS);
3483                                 ix0 = _mm_packs_epi32(ix0, ix1);
3484                                 iy0 = _mm_packs_epi32(iy0, iy1);
3485                                 ix1 = _mm_unpacklo_epi16(ix0, iy0);
3486                                 iy1 = _mm_unpackhi_epi16(ix0, iy0);
3487                                 _mm_storeu_si128((__m128i*)(XY + x1*2), ix1);
3488                                 _mm_storeu_si128((__m128i*)(XY + x1*2 + 8), iy1);
3489                             }
3490                         }
3491                     #endif
3492
3493                         for( ; x1 < bcols; x1++ )
3494                         {
3495                             int sx = cvRound(sX[x1]*INTER_TAB_SIZE);
3496                             int sy = cvRound(sY[x1]*INTER_TAB_SIZE);
3497                             int v = (sy & (INTER_TAB_SIZE-1))*INTER_TAB_SIZE + (sx & (INTER_TAB_SIZE-1));
3498                             XY[x1*2] = saturate_cast<short>(sx >> INTER_BITS);
3499                             XY[x1*2+1] = saturate_cast<short>(sy >> INTER_BITS);
3500                             A[x1] = (ushort)v;
3501                         }
3502                     }
3503                     else
3504                     {
3505                         const float* sXY = (const float*)(m1->data + m1->step*(y+y1)) + x*2;
3506
3507                         for( x1 = 0; x1 < bcols; x1++ )
3508                         {
3509                             int sx = cvRound(sXY[x1*2]*INTER_TAB_SIZE);
3510                             int sy = cvRound(sXY[x1*2+1]*INTER_TAB_SIZE);
3511                             int v = (sy & (INTER_TAB_SIZE-1))*INTER_TAB_SIZE + (sx & (INTER_TAB_SIZE-1));
3512                             XY[x1*2] = saturate_cast<short>(sx >> INTER_BITS);
3513                             XY[x1*2+1] = saturate_cast<short>(sy >> INTER_BITS);
3514                             A[x1] = (ushort)v;
3515                         }
3516                     }
3517                 }
3518                 ifunc(*src, dpart, bufxy, bufa, ctab, borderType, borderValue);
3519             }
3520         }
3521     }
3522
3523 private:
3524     const Mat* src;
3525     Mat* dst;
3526     const Mat *m1, *m2;
3527     int borderType;
3528     Scalar borderValue;
3529     int planar_input;
3530     RemapNNFunc nnfunc;
3531     RemapFunc ifunc;
3532     const void *ctab;
3533 };
3534
3535 #ifdef HAVE_OPENCL
3536
3537 static bool ocl_remap(InputArray _src, OutputArray _dst, InputArray _map1, InputArray _map2,
3538                       int interpolation, int borderType, const Scalar& borderValue)
3539 {
3540     int cn = _src.channels(), type = _src.type(), depth = _src.depth();
3541
3542     if (borderType == BORDER_TRANSPARENT || cn == 3 || !(interpolation == INTER_LINEAR || interpolation == INTER_NEAREST)
3543             || _map1.type() == CV_16SC1 || _map2.type() == CV_16SC1)
3544         return false;
3545
3546     UMat src = _src.getUMat(), map1 = _map1.getUMat(), map2 = _map2.getUMat();
3547
3548     if( (map1.type() == CV_16SC2 && (map2.type() == CV_16UC1 || map2.empty())) ||
3549         (map2.type() == CV_16SC2 && (map1.type() == CV_16UC1 || map1.empty())) )
3550     {
3551         if (map1.type() != CV_16SC2)
3552             std::swap(map1, map2);
3553     }
3554     else
3555         CV_Assert( map1.type() == CV_32FC2 || (map1.type() == CV_32FC1 && map2.type() == CV_32FC1) );
3556
3557     _dst.create(map1.size(), type);
3558     UMat dst = _dst.getUMat();
3559
3560     String kernelName = "remap";
3561     if (map1.type() == CV_32FC2 && map2.empty())
3562         kernelName += "_32FC2";
3563     else if (map1.type() == CV_16SC2)
3564     {
3565         kernelName += "_16SC2";
3566         if (!map2.empty())
3567             kernelName += "_16UC1";
3568     }
3569     else if (map1.type() == CV_32FC1 && map2.type() == CV_32FC1)
3570         kernelName += "_2_32FC1";
3571     else
3572         CV_Error(Error::StsBadArg, "Unsupported map types");
3573
3574     static const char * const interMap[] = { "INTER_NEAREST", "INTER_LINEAR", "INTER_CUBIC", "INTER_LINEAR", "INTER_LANCZOS" };
3575     static const char * const borderMap[] = { "BORDER_CONSTANT", "BORDER_REPLICATE", "BORDER_REFLECT", "BORDER_WRAP",
3576                            "BORDER_REFLECT_101", "BORDER_TRANSPARENT" };
3577     String buildOptions = format("-D %s -D %s -D T=%s", interMap[interpolation], borderMap[borderType], ocl::typeToStr(type));
3578
3579     if (interpolation != INTER_NEAREST)
3580     {
3581         char cvt[3][40];
3582         int wdepth = std::max(CV_32F, dst.depth());
3583         buildOptions = buildOptions
3584                       + format(" -D WT=%s -D convertToT=%s -D convertToWT=%s"
3585                                " -D convertToWT2=%s -D WT2=%s",
3586                                ocl::typeToStr(CV_MAKE_TYPE(wdepth, cn)),
3587                                ocl::convertTypeStr(wdepth, depth, cn, cvt[0]),
3588                                ocl::convertTypeStr(depth, wdepth, cn, cvt[1]),
3589                                ocl::convertTypeStr(CV_32S, wdepth, 2, cvt[2]),
3590                                ocl::typeToStr(CV_MAKE_TYPE(wdepth, 2)));
3591     }
3592
3593     ocl::Kernel k(kernelName.c_str(), ocl::imgproc::remap_oclsrc, buildOptions);
3594
3595     Mat scalar(1, 1, type, borderValue);
3596     ocl::KernelArg srcarg = ocl::KernelArg::ReadOnly(src), dstarg = ocl::KernelArg::WriteOnly(dst),
3597             map1arg = ocl::KernelArg::ReadOnlyNoSize(map1),
3598             scalararg = ocl::KernelArg::Constant((void*)scalar.data, scalar.elemSize());
3599
3600     if (map2.empty())
3601         k.args(srcarg, dstarg, map1arg, scalararg);
3602     else
3603         k.args(srcarg, dstarg, map1arg, ocl::KernelArg::ReadOnlyNoSize(map2), scalararg);
3604
3605     size_t globalThreads[2] = { dst.cols, dst.rows };
3606     return k.run(2, globalThreads, NULL, false);
3607 }
3608
3609 #endif
3610
3611 }
3612
3613 void cv::remap( InputArray _src, OutputArray _dst,
3614                 InputArray _map1, InputArray _map2,
3615                 int interpolation, int borderType, const Scalar& borderValue )
3616 {
3617     static RemapNNFunc nn_tab[] =
3618     {
3619         remapNearest<uchar>, remapNearest<schar>, remapNearest<ushort>, remapNearest<short>,
3620         remapNearest<int>, remapNearest<float>, remapNearest<double>, 0
3621     };
3622
3623     static RemapFunc linear_tab[] =
3624     {
3625         remapBilinear<FixedPtCast<int, uchar, INTER_REMAP_COEF_BITS>, RemapVec_8u, short>, 0,
3626         remapBilinear<Cast<float, ushort>, RemapNoVec, float>,
3627         remapBilinear<Cast<float, short>, RemapNoVec, float>, 0,
3628         remapBilinear<Cast<float, float>, RemapNoVec, float>,
3629         remapBilinear<Cast<double, double>, RemapNoVec, float>, 0
3630     };
3631
3632     static RemapFunc cubic_tab[] =
3633     {
3634         remapBicubic<FixedPtCast<int, uchar, INTER_REMAP_COEF_BITS>, short, INTER_REMAP_COEF_SCALE>, 0,
3635         remapBicubic<Cast<float, ushort>, float, 1>,
3636         remapBicubic<Cast<float, short>, float, 1>, 0,
3637         remapBicubic<Cast<float, float>, float, 1>,
3638         remapBicubic<Cast<double, double>, float, 1>, 0
3639     };
3640
3641     static RemapFunc lanczos4_tab[] =
3642     {
3643         remapLanczos4<FixedPtCast<int, uchar, INTER_REMAP_COEF_BITS>, short, INTER_REMAP_COEF_SCALE>, 0,
3644         remapLanczos4<Cast<float, ushort>, float, 1>,
3645         remapLanczos4<Cast<float, short>, float, 1>, 0,
3646         remapLanczos4<Cast<float, float>, float, 1>,
3647         remapLanczos4<Cast<double, double>, float, 1>, 0
3648     };
3649
3650     CV_Assert( _map1.size().area() > 0 );
3651     CV_Assert( _map2.empty() || (_map2.size() == _map1.size()));
3652
3653     CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
3654                ocl_remap(_src, _dst, _map1, _map2, interpolation, borderType, borderValue))
3655
3656     Mat src = _src.getMat(), map1 = _map1.getMat(), map2 = _map2.getMat();
3657     _dst.create( map1.size(), src.type() );
3658     Mat dst = _dst.getMat();
3659     if( dst.data == src.data )
3660         src = src.clone();
3661
3662     int depth = src.depth();
3663     RemapNNFunc nnfunc = 0;
3664     RemapFunc ifunc = 0;
3665     const void* ctab = 0;
3666     bool fixpt = depth == CV_8U;
3667     bool planar_input = false;
3668
3669     if( interpolation == INTER_NEAREST )
3670     {
3671         nnfunc = nn_tab[depth];
3672         CV_Assert( nnfunc != 0 );
3673     }
3674     else
3675     {
3676         if( interpolation == INTER_AREA )
3677             interpolation = INTER_LINEAR;
3678
3679         if( interpolation == INTER_LINEAR )
3680             ifunc = linear_tab[depth];
3681         else if( interpolation == INTER_CUBIC )
3682             ifunc = cubic_tab[depth];
3683         else if( interpolation == INTER_LANCZOS4 )
3684             ifunc = lanczos4_tab[depth];
3685         else
3686             CV_Error( CV_StsBadArg, "Unknown interpolation method" );
3687         CV_Assert( ifunc != 0 );
3688         ctab = initInterTab2D( interpolation, fixpt );
3689     }
3690
3691     const Mat *m1 = &map1, *m2 = &map2;
3692
3693     if( (map1.type() == CV_16SC2 && (map2.type() == CV_16UC1 || map2.type() == CV_16SC1 || !map2.data)) ||
3694         (map2.type() == CV_16SC2 && (map1.type() == CV_16UC1 || map1.type() == CV_16SC1 || !map1.data)) )
3695     {
3696         if( map1.type() != CV_16SC2 )
3697             std::swap(m1, m2);
3698     }
3699     else
3700     {
3701         CV_Assert( ((map1.type() == CV_32FC2 || map1.type() == CV_16SC2) && !map2.data) ||
3702             (map1.type() == CV_32FC1 && map2.type() == CV_32FC1) );
3703         planar_input = map1.channels() == 1;
3704     }
3705
3706     RemapInvoker invoker(src, dst, m1, m2,
3707                          borderType, borderValue, planar_input, nnfunc, ifunc,
3708                          ctab);
3709     parallel_for_(Range(0, dst.rows), invoker, dst.total()/(double)(1<<16));
3710 }
3711
3712
3713 void cv::convertMaps( InputArray _map1, InputArray _map2,
3714                       OutputArray _dstmap1, OutputArray _dstmap2,
3715                       int dstm1type, bool nninterpolate )
3716 {
3717     Mat map1 = _map1.getMat(), map2 = _map2.getMat(), dstmap1, dstmap2;
3718     Size size = map1.size();
3719     const Mat *m1 = &map1, *m2 = &map2;
3720     int m1type = m1->type(), m2type = m2->type();
3721
3722     CV_Assert( (m1type == CV_16SC2 && (nninterpolate || m2type == CV_16UC1 || m2type == CV_16SC1)) ||
3723                (m2type == CV_16SC2 && (nninterpolate || m1type == CV_16UC1 || m1type == CV_16SC1)) ||
3724                (m1type == CV_32FC1 && m2type == CV_32FC1) ||
3725                (m1type == CV_32FC2 && !m2->data) );
3726
3727     if( m2type == CV_16SC2 )
3728     {
3729         std::swap( m1, m2 );
3730         std::swap( m1type, m2type );
3731     }
3732
3733     if( dstm1type <= 0 )
3734         dstm1type = m1type == CV_16SC2 ? CV_32FC2 : CV_16SC2;
3735     CV_Assert( dstm1type == CV_16SC2 || dstm1type == CV_32FC1 || dstm1type == CV_32FC2 );
3736     _dstmap1.create( size, dstm1type );
3737     dstmap1 = _dstmap1.getMat();
3738
3739     if( !nninterpolate && dstm1type != CV_32FC2 )
3740     {
3741         _dstmap2.create( size, dstm1type == CV_16SC2 ? CV_16UC1 : CV_32FC1 );
3742         dstmap2 = _dstmap2.getMat();
3743     }
3744     else
3745         _dstmap2.release();
3746
3747     if( m1type == dstm1type || (nninterpolate &&
3748         ((m1type == CV_16SC2 && dstm1type == CV_32FC2) ||
3749         (m1type == CV_32FC2 && dstm1type == CV_16SC2))) )
3750     {
3751         m1->convertTo( dstmap1, dstmap1.type() );
3752         if( dstmap2.data && dstmap2.type() == m2->type() )
3753             m2->copyTo( dstmap2 );
3754         return;
3755     }
3756
3757     if( m1type == CV_32FC1 && dstm1type == CV_32FC2 )
3758     {
3759         Mat vdata[] = { *m1, *m2 };
3760         merge( vdata, 2, dstmap1 );
3761         return;
3762     }
3763
3764     if( m1type == CV_32FC2 && dstm1type == CV_32FC1 )
3765     {
3766         Mat mv[] = { dstmap1, dstmap2 };
3767         split( *m1, mv );
3768         return;
3769     }
3770
3771     if( m1->isContinuous() && (!m2->data || m2->isContinuous()) &&
3772         dstmap1.isContinuous() && (!dstmap2.data || dstmap2.isContinuous()) )
3773     {
3774         size.width *= size.height;
3775         size.height = 1;
3776     }
3777
3778     const float scale = 1.f/INTER_TAB_SIZE;
3779     int x, y;
3780     for( y = 0; y < size.height; y++ )
3781     {
3782         const float* src1f = (const float*)(m1->data + m1->step*y);
3783         const float* src2f = (const float*)(m2->data + m2->step*y);
3784         const short* src1 = (const short*)src1f;
3785         const ushort* src2 = (const ushort*)src2f;
3786
3787         float* dst1f = (float*)(dstmap1.data + dstmap1.step*y);
3788         float* dst2f = (float*)(dstmap2.data + dstmap2.step*y);
3789         short* dst1 = (short*)dst1f;
3790         ushort* dst2 = (ushort*)dst2f;
3791
3792         if( m1type == CV_32FC1 && dstm1type == CV_16SC2 )
3793         {
3794             if( nninterpolate )
3795                 for( x = 0; x < size.width; x++ )
3796                 {
3797                     dst1[x*2] = saturate_cast<short>(src1f[x]);
3798                     dst1[x*2+1] = saturate_cast<short>(src2f[x]);
3799                 }
3800             else
3801                 for( x = 0; x < size.width; x++ )
3802                 {
3803                     int ix = saturate_cast<int>(src1f[x]*INTER_TAB_SIZE);
3804                     int iy = saturate_cast<int>(src2f[x]*INTER_TAB_SIZE);
3805                     dst1[x*2] = saturate_cast<short>(ix >> INTER_BITS);
3806                     dst1[x*2+1] = saturate_cast<short>(iy >> INTER_BITS);
3807                     dst2[x] = (ushort)((iy & (INTER_TAB_SIZE-1))*INTER_TAB_SIZE + (ix & (INTER_TAB_SIZE-1)));
3808                 }
3809         }
3810         else if( m1type == CV_32FC2 && dstm1type == CV_16SC2 )
3811         {
3812             if( nninterpolate )
3813                 for( x = 0; x < size.width; x++ )
3814                 {
3815                     dst1[x*2] = saturate_cast<short>(src1f[x*2]);
3816                     dst1[x*2+1] = saturate_cast<short>(src1f[x*2+1]);
3817                 }
3818             else
3819                 for( x = 0; x < size.width; x++ )
3820                 {
3821                     int ix = saturate_cast<int>(src1f[x*2]*INTER_TAB_SIZE);
3822                     int iy = saturate_cast<int>(src1f[x*2+1]*INTER_TAB_SIZE);
3823                     dst1[x*2] = saturate_cast<short>(ix >> INTER_BITS);
3824                     dst1[x*2+1] = saturate_cast<short>(iy >> INTER_BITS);
3825                     dst2[x] = (ushort)((iy & (INTER_TAB_SIZE-1))*INTER_TAB_SIZE + (ix & (INTER_TAB_SIZE-1)));
3826                 }
3827         }
3828         else if( m1type == CV_16SC2 && dstm1type == CV_32FC1 )
3829         {
3830             for( x = 0; x < size.width; x++ )
3831             {
3832                 int fxy = src2 ? src2[x] & (INTER_TAB_SIZE2-1) : 0;
3833                 dst1f[x] = src1[x*2] + (fxy & (INTER_TAB_SIZE-1))*scale;
3834                 dst2f[x] = src1[x*2+1] + (fxy >> INTER_BITS)*scale;
3835             }
3836         }
3837         else if( m1type == CV_16SC2 && dstm1type == CV_32FC2 )
3838         {
3839             for( x = 0; x < size.width; x++ )
3840             {
3841                 int fxy = src2 ? src2[x] & (INTER_TAB_SIZE2-1): 0;
3842                 dst1f[x*2] = src1[x*2] + (fxy & (INTER_TAB_SIZE-1))*scale;
3843                 dst1f[x*2+1] = src1[x*2+1] + (fxy >> INTER_BITS)*scale;
3844             }
3845         }
3846         else
3847             CV_Error( CV_StsNotImplemented, "Unsupported combination of input/output matrices" );
3848     }
3849 }
3850
3851
3852 namespace cv
3853 {
3854
3855 class warpAffineInvoker :
3856     public ParallelLoopBody
3857 {
3858 public:
3859     warpAffineInvoker(const Mat &_src, Mat &_dst, int _interpolation, int _borderType,
3860                       const Scalar &_borderValue, int *_adelta, int *_bdelta, double *_M) :
3861         ParallelLoopBody(), src(_src), dst(_dst), interpolation(_interpolation),
3862         borderType(_borderType), borderValue(_borderValue), adelta(_adelta), bdelta(_bdelta),
3863         M(_M)
3864     {
3865     }
3866
3867     virtual void operator() (const Range& range) const
3868     {
3869         const int BLOCK_SZ = 64;
3870         short XY[BLOCK_SZ*BLOCK_SZ*2], A[BLOCK_SZ*BLOCK_SZ];
3871         const int AB_BITS = MAX(10, (int)INTER_BITS);
3872         const int AB_SCALE = 1 << AB_BITS;
3873         int round_delta = interpolation == INTER_NEAREST ? AB_SCALE/2 : AB_SCALE/INTER_TAB_SIZE/2, x, y, x1, y1;
3874     #if CV_SSE2
3875         bool useSIMD = checkHardwareSupport(CV_CPU_SSE2);
3876     #endif
3877
3878         int bh0 = std::min(BLOCK_SZ/2, dst.rows);
3879         int bw0 = std::min(BLOCK_SZ*BLOCK_SZ/bh0, dst.cols);
3880         bh0 = std::min(BLOCK_SZ*BLOCK_SZ/bw0, dst.rows);
3881
3882         for( y = range.start; y < range.end; y += bh0 )
3883         {
3884             for( x = 0; x < dst.cols; x += bw0 )
3885             {
3886                 int bw = std::min( bw0, dst.cols - x);
3887                 int bh = std::min( bh0, range.end - y);
3888
3889                 Mat _XY(bh, bw, CV_16SC2, XY), matA;
3890                 Mat dpart(dst, Rect(x, y, bw, bh));
3891
3892                 for( y1 = 0; y1 < bh; y1++ )
3893                 {
3894                     short* xy = XY + y1*bw*2;
3895                     int X0 = saturate_cast<int>((M[1]*(y + y1) + M[2])*AB_SCALE) + round_delta;
3896                     int Y0 = saturate_cast<int>((M[4]*(y + y1) + M[5])*AB_SCALE) + round_delta;
3897
3898                     if( interpolation == INTER_NEAREST )
3899                         for( x1 = 0; x1 < bw; x1++ )
3900                         {
3901                             int X = (X0 + adelta[x+x1]) >> AB_BITS;
3902                             int Y = (Y0 + bdelta[x+x1]) >> AB_BITS;
3903                             xy[x1*2] = saturate_cast<short>(X);
3904                             xy[x1*2+1] = saturate_cast<short>(Y);
3905                         }
3906                     else
3907                     {
3908                         short* alpha = A + y1*bw;
3909                         x1 = 0;
3910                     #if CV_SSE2
3911                         if( useSIMD )
3912                         {
3913                             __m128i fxy_mask = _mm_set1_epi32(INTER_TAB_SIZE - 1);
3914                             __m128i XX = _mm_set1_epi32(X0), YY = _mm_set1_epi32(Y0);
3915                             for( ; x1 <= bw - 8; x1 += 8 )
3916                             {
3917                                 __m128i tx0, tx1, ty0, ty1;
3918                                 tx0 = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(adelta + x + x1)), XX);
3919                                 ty0 = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(bdelta + x + x1)), YY);
3920                                 tx1 = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(adelta + x + x1 + 4)), XX);
3921                                 ty1 = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(bdelta + x + x1 + 4)), YY);
3922
3923                                 tx0 = _mm_srai_epi32(tx0, AB_BITS - INTER_BITS);
3924                                 ty0 = _mm_srai_epi32(ty0, AB_BITS - INTER_BITS);
3925                                 tx1 = _mm_srai_epi32(tx1, AB_BITS - INTER_BITS);
3926                                 ty1 = _mm_srai_epi32(ty1, AB_BITS - INTER_BITS);
3927
3928                                 __m128i fx_ = _mm_packs_epi32(_mm_and_si128(tx0, fxy_mask),
3929                                                             _mm_and_si128(tx1, fxy_mask));
3930                                 __m128i fy_ = _mm_packs_epi32(_mm_and_si128(ty0, fxy_mask),
3931                                                             _mm_and_si128(ty1, fxy_mask));
3932                                 tx0 = _mm_packs_epi32(_mm_srai_epi32(tx0, INTER_BITS),
3933                                                             _mm_srai_epi32(tx1, INTER_BITS));
3934                                 ty0 = _mm_packs_epi32(_mm_srai_epi32(ty0, INTER_BITS),
3935                                                     _mm_srai_epi32(ty1, INTER_BITS));
3936                                 fx_ = _mm_adds_epi16(fx_, _mm_slli_epi16(fy_, INTER_BITS));
3937
3938                                 _mm_storeu_si128((__m128i*)(xy + x1*2), _mm_unpacklo_epi16(tx0, ty0));
3939                                 _mm_storeu_si128((__m128i*)(xy + x1*2 + 8), _mm_unpackhi_epi16(tx0, ty0));
3940                                 _mm_storeu_si128((__m128i*)(alpha + x1), fx_);
3941                             }
3942                         }
3943                     #endif
3944                         for( ; x1 < bw; x1++ )
3945                         {
3946                             int X = (X0 + adelta[x+x1]) >> (AB_BITS - INTER_BITS);
3947                             int Y = (Y0 + bdelta[x+x1]) >> (AB_BITS - INTER_BITS);
3948                             xy[x1*2] = saturate_cast<short>(X >> INTER_BITS);
3949                             xy[x1*2+1] = saturate_cast<short>(Y >> INTER_BITS);
3950                             alpha[x1] = (short)((Y & (INTER_TAB_SIZE-1))*INTER_TAB_SIZE +
3951                                     (X & (INTER_TAB_SIZE-1)));
3952                         }
3953                     }
3954                 }
3955
3956                 if( interpolation == INTER_NEAREST )
3957                     remap( src, dpart, _XY, Mat(), interpolation, borderType, borderValue );
3958                 else
3959                 {
3960                     Mat _matA(bh, bw, CV_16U, A);
3961                     remap( src, dpart, _XY, _matA, interpolation, borderType, borderValue );
3962                 }
3963             }
3964         }
3965     }
3966
3967 private:
3968     Mat src;
3969     Mat dst;
3970     int interpolation, borderType;
3971     Scalar borderValue;
3972     int *adelta, *bdelta;
3973     double *M;
3974 };
3975
3976 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
3977 class IPPwarpAffineInvoker :
3978     public ParallelLoopBody
3979 {
3980 public:
3981     IPPwarpAffineInvoker(Mat &_src, Mat &_dst, double (&_coeffs)[2][3], int &_interpolation, int &_borderType, const Scalar &_borderValue, ippiWarpAffineBackFunc _func, bool *_ok) :
3982       ParallelLoopBody(), src(_src), dst(_dst), mode(_interpolation), coeffs(_coeffs), borderType(_borderType), borderValue(_borderValue), func(_func), ok(_ok)
3983       {
3984           *ok = true;
3985       }
3986
3987       virtual void operator() (const Range& range) const
3988       {
3989           IppiSize srcsize = { src.cols, src.rows };
3990           IppiRect srcroi = { 0, 0, src.cols, src.rows };
3991           IppiRect dstroi = { 0, range.start, dst.cols, range.end - range.start };
3992           int cnn = src.channels();
3993           if( borderType == BORDER_CONSTANT )
3994           {
3995               IppiSize setSize = { dst.cols, range.end - range.start };
3996               void *dataPointer = dst.data + dst.step[0] * range.start;
3997               if( !IPPSet( borderValue, dataPointer, (int)dst.step[0], setSize, cnn, src.depth() ) )
3998               {
3999                   *ok = false;
4000                   return;
4001               }
4002           }
4003           if( func( src.data, srcsize, (int)src.step[0], srcroi, dst.data, (int)dst.step[0], dstroi, coeffs, mode ) < 0) ////Aug 2013: problem in IPP 7.1, 8.0 : sometimes function return ippStsCoeffErr
4004               *ok = false;
4005       }
4006 private:
4007     Mat &src;
4008     Mat &dst;
4009     double (&coeffs)[2][3];
4010     int mode;
4011     int borderType;
4012     Scalar borderValue;
4013     ippiWarpAffineBackFunc func;
4014     bool *ok;
4015     const IPPwarpAffineInvoker& operator= (const IPPwarpAffineInvoker&);
4016 };
4017 #endif
4018
4019 #ifdef HAVE_OPENCL
4020
4021 enum { OCL_OP_PERSPECTIVE = 1, OCL_OP_AFFINE = 0 };
4022
4023 static bool ocl_warpTransform(InputArray _src, OutputArray _dst, InputArray _M0,
4024                               Size dsize, int flags, int borderType, const Scalar& borderValue,
4025                               int op_type)
4026 {
4027     CV_Assert(op_type == OCL_OP_AFFINE || op_type == OCL_OP_PERSPECTIVE);
4028
4029     int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
4030     double doubleSupport = ocl::Device::getDefault().doubleFPConfig() > 0;
4031
4032     int interpolation = flags & INTER_MAX;
4033     if( interpolation == INTER_AREA )
4034         interpolation = INTER_LINEAR;
4035
4036     if ( !(borderType == cv::BORDER_CONSTANT &&
4037            (interpolation == cv::INTER_NEAREST || interpolation == cv::INTER_LINEAR || interpolation == cv::INTER_CUBIC)) ||
4038          (!doubleSupport && depth == CV_64F) || cn > 4)
4039         return false;
4040
4041     const char * const interpolationMap[3] = { "NEAREST", "LINEAR", "CUBIC" };
4042     ocl::ProgramSource program = op_type == OCL_OP_AFFINE ?
4043                 ocl::imgproc::warp_affine_oclsrc : ocl::imgproc::warp_perspective_oclsrc;
4044     const char * const kernelName = op_type == OCL_OP_AFFINE ? "warpAffine" : "warpPerspective";
4045
4046     int scalarcn = cn == 3 ? 4 : cn;
4047     int wdepth = interpolation == INTER_NEAREST ? depth : std::max(CV_32S, depth);
4048     int sctype = CV_MAKETYPE(wdepth, scalarcn);
4049
4050     ocl::Kernel k;
4051     String opts;
4052     if (interpolation == INTER_NEAREST)
4053     {
4054         opts = format("-D INTER_NEAREST -D T=%s%s -D T1=%s -D ST=%s -D cn=%d", ocl::typeToStr(type),
4055                       doubleSupport ? " -D DOUBLE_SUPPORT" : "",
4056                       ocl::typeToStr(CV_MAT_DEPTH(type)),
4057                       ocl::typeToStr(sctype),
4058                       cn);
4059     }
4060     else
4061     {
4062         char cvt[2][50];
4063         opts = format("-D INTER_%s -D T=%s -D T1=%s -D ST=%s -D WT=%s -D depth=%d -D convertToWT=%s -D convertToT=%s%s -D cn=%d",
4064                       interpolationMap[interpolation], ocl::typeToStr(type),
4065                       ocl::typeToStr(CV_MAT_DEPTH(type)),
4066                       ocl::typeToStr(sctype),
4067                       ocl::typeToStr(CV_MAKE_TYPE(wdepth, cn)), depth,
4068                       ocl::convertTypeStr(depth, wdepth, cn, cvt[0]),
4069                       ocl::convertTypeStr(wdepth, depth, cn, cvt[1]),
4070                       doubleSupport ? " -D DOUBLE_SUPPORT" : "", cn);
4071     }
4072
4073     k.create(kernelName, program, opts);
4074     if (k.empty())
4075         return false;
4076
4077     double borderBuf[] = {0, 0, 0, 0};
4078     scalarToRawData(borderValue, borderBuf, sctype);
4079
4080     UMat src = _src.getUMat(), M0;
4081     _dst.create( dsize.area() == 0 ? src.size() : dsize, src.type() );
4082     UMat dst = _dst.getUMat();
4083
4084     double M[9];
4085     int matRows = (op_type == OCL_OP_AFFINE ? 2 : 3);
4086     Mat matM(matRows, 3, CV_64F, M), M1 = _M0.getMat();
4087     CV_Assert( (M1.type() == CV_32F || M1.type() == CV_64F) &&
4088                M1.rows == matRows && M1.cols == 3 );
4089     M1.convertTo(matM, matM.type());
4090
4091     if( !(flags & WARP_INVERSE_MAP) )
4092     {
4093         if (op_type == OCL_OP_PERSPECTIVE)
4094             invert(matM, matM);
4095         else
4096         {
4097             double D = M[0]*M[4] - M[1]*M[3];
4098             D = D != 0 ? 1./D : 0;
4099             double A11 = M[4]*D, A22=M[0]*D;
4100             M[0] = A11; M[1] *= -D;
4101             M[3] *= -D; M[4] = A22;
4102             double b1 = -M[0]*M[2] - M[1]*M[5];
4103             double b2 = -M[3]*M[2] - M[4]*M[5];
4104             M[2] = b1; M[5] = b2;
4105         }
4106     }
4107     matM.convertTo(M0, doubleSupport ? CV_64F : CV_32F);
4108
4109     k.args(ocl::KernelArg::ReadOnly(src), ocl::KernelArg::WriteOnly(dst), ocl::KernelArg::PtrReadOnly(M0),
4110            ocl::KernelArg(0, 0, 0, borderBuf, CV_ELEM_SIZE(sctype)));
4111
4112     size_t globalThreads[2] = { dst.cols, dst.rows };
4113     return k.run(2, globalThreads, NULL, false);
4114 }
4115
4116 #endif
4117
4118 }
4119
4120
4121 void cv::warpAffine( InputArray _src, OutputArray _dst,
4122                      InputArray _M0, Size dsize,
4123                      int flags, int borderType, const Scalar& borderValue )
4124 {
4125     CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
4126                ocl_warpTransform(_src, _dst, _M0, dsize, flags, borderType,
4127                                  borderValue, OCL_OP_AFFINE))
4128
4129     Mat src = _src.getMat(), M0 = _M0.getMat();
4130     _dst.create( dsize.area() == 0 ? src.size() : dsize, src.type() );
4131     Mat dst = _dst.getMat();
4132     CV_Assert( src.cols > 0 && src.rows > 0 );
4133     if( dst.data == src.data )
4134         src = src.clone();
4135
4136     double M[6];
4137     Mat matM(2, 3, CV_64F, M);
4138     int interpolation = flags & INTER_MAX;
4139     if( interpolation == INTER_AREA )
4140         interpolation = INTER_LINEAR;
4141
4142     CV_Assert( (M0.type() == CV_32F || M0.type() == CV_64F) && M0.rows == 2 && M0.cols == 3 );
4143     M0.convertTo(matM, matM.type());
4144
4145 #ifdef HAVE_TEGRA_OPTIMIZATION
4146     if( tegra::warpAffine(src, dst, M, flags, borderType, borderValue) )
4147         return;
4148 #endif
4149
4150     if( !(flags & WARP_INVERSE_MAP) )
4151     {
4152         double D = M[0]*M[4] - M[1]*M[3];
4153         D = D != 0 ? 1./D : 0;
4154         double A11 = M[4]*D, A22=M[0]*D;
4155         M[0] = A11; M[1] *= -D;
4156         M[3] *= -D; M[4] = A22;
4157         double b1 = -M[0]*M[2] - M[1]*M[5];
4158         double b2 = -M[3]*M[2] - M[4]*M[5];
4159         M[2] = b1; M[5] = b2;
4160     }
4161
4162     int x;
4163     AutoBuffer<int> _abdelta(dst.cols*2);
4164     int* adelta = &_abdelta[0], *bdelta = adelta + dst.cols;
4165     const int AB_BITS = MAX(10, (int)INTER_BITS);
4166     const int AB_SCALE = 1 << AB_BITS;
4167 /*
4168 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
4169     int depth = src.depth();
4170     int channels = src.channels();
4171     if( ( depth == CV_8U || depth == CV_16U || depth == CV_32F ) &&
4172         ( channels == 1 || channels == 3 || channels == 4 ) &&
4173         ( borderType == cv::BORDER_TRANSPARENT || ( borderType == cv::BORDER_CONSTANT ) ) )
4174     {
4175         int type = src.type();
4176         ippiWarpAffineBackFunc ippFunc =
4177             type == CV_8UC1 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_8u_C1R :
4178             type == CV_8UC3 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_8u_C3R :
4179             type == CV_8UC4 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_8u_C4R :
4180             type == CV_16UC1 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_16u_C1R :
4181             type == CV_16UC3 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_16u_C3R :
4182             type == CV_16UC4 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_16u_C4R :
4183             type == CV_32FC1 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_32f_C1R :
4184             type == CV_32FC3 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_32f_C3R :
4185             type == CV_32FC4 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_32f_C4R :
4186             0;
4187         int mode =
4188             flags == INTER_LINEAR ? IPPI_INTER_LINEAR :
4189             flags == INTER_NEAREST ? IPPI_INTER_NN :
4190             flags == INTER_CUBIC ? IPPI_INTER_CUBIC :
4191             0;
4192         if( mode && ippFunc )
4193         {
4194             double coeffs[2][3];
4195             for( int i = 0; i < 2; i++ )
4196             {
4197                 for( int j = 0; j < 3; j++ )
4198                 {
4199                     coeffs[i][j] = matM.at<double>(i, j);
4200                 }
4201             }
4202             bool ok;
4203             Range range(0, dst.rows);
4204             IPPwarpAffineInvoker invoker(src, dst, coeffs, mode, borderType, borderValue, ippFunc, &ok);
4205             parallel_for_(range, invoker, dst.total()/(double)(1<<16));
4206             if( ok )
4207                 return;
4208         }
4209     }
4210 #endif
4211 */
4212     for( x = 0; x < dst.cols; x++ )
4213     {
4214         adelta[x] = saturate_cast<int>(M[0]*x*AB_SCALE);
4215         bdelta[x] = saturate_cast<int>(M[3]*x*AB_SCALE);
4216     }
4217
4218     Range range(0, dst.rows);
4219     warpAffineInvoker invoker(src, dst, interpolation, borderType,
4220                               borderValue, adelta, bdelta, M);
4221     parallel_for_(range, invoker, dst.total()/(double)(1<<16));
4222 }
4223
4224
4225 namespace cv
4226 {
4227
4228 class warpPerspectiveInvoker :
4229     public ParallelLoopBody
4230 {
4231 public:
4232
4233     warpPerspectiveInvoker(const Mat &_src, Mat &_dst, double *_M, int _interpolation,
4234                            int _borderType, const Scalar &_borderValue) :
4235         ParallelLoopBody(), src(_src), dst(_dst), M(_M), interpolation(_interpolation),
4236         borderType(_borderType), borderValue(_borderValue)
4237     {
4238     }
4239
4240     virtual void operator() (const Range& range) const
4241     {
4242         const int BLOCK_SZ = 32;
4243         short XY[BLOCK_SZ*BLOCK_SZ*2], A[BLOCK_SZ*BLOCK_SZ];
4244         int x, y, x1, y1, width = dst.cols, height = dst.rows;
4245
4246         int bh0 = std::min(BLOCK_SZ/2, height);
4247         int bw0 = std::min(BLOCK_SZ*BLOCK_SZ/bh0, width);
4248         bh0 = std::min(BLOCK_SZ*BLOCK_SZ/bw0, height);
4249
4250         for( y = range.start; y < range.end; y += bh0 )
4251         {
4252             for( x = 0; x < width; x += bw0 )
4253             {
4254                 int bw = std::min( bw0, width - x);
4255                 int bh = std::min( bh0, range.end - y); // height
4256
4257                 Mat _XY(bh, bw, CV_16SC2, XY), matA;
4258                 Mat dpart(dst, Rect(x, y, bw, bh));
4259
4260                 for( y1 = 0; y1 < bh; y1++ )
4261                 {
4262                     short* xy = XY + y1*bw*2;
4263                     double X0 = M[0]*x + M[1]*(y + y1) + M[2];
4264                     double Y0 = M[3]*x + M[4]*(y + y1) + M[5];
4265                     double W0 = M[6]*x + M[7]*(y + y1) + M[8];
4266
4267                     if( interpolation == INTER_NEAREST )
4268                         for( x1 = 0; x1 < bw; x1++ )
4269                         {
4270                             double W = W0 + M[6]*x1;
4271                             W = W ? 1./W : 0;
4272                             double fX = std::max((double)INT_MIN, std::min((double)INT_MAX, (X0 + M[0]*x1)*W));
4273                             double fY = std::max((double)INT_MIN, std::min((double)INT_MAX, (Y0 + M[3]*x1)*W));
4274                             int X = saturate_cast<int>(fX);
4275                             int Y = saturate_cast<int>(fY);
4276
4277                             xy[x1*2] = saturate_cast<short>(X);
4278                             xy[x1*2+1] = saturate_cast<short>(Y);
4279                         }
4280                     else
4281                     {
4282                         short* alpha = A + y1*bw;
4283                         for( x1 = 0; x1 < bw; x1++ )
4284                         {
4285                             double W = W0 + M[6]*x1;
4286                             W = W ? INTER_TAB_SIZE/W : 0;
4287                             double fX = std::max((double)INT_MIN, std::min((double)INT_MAX, (X0 + M[0]*x1)*W));
4288                             double fY = std::max((double)INT_MIN, std::min((double)INT_MAX, (Y0 + M[3]*x1)*W));
4289                             int X = saturate_cast<int>(fX);
4290                             int Y = saturate_cast<int>(fY);
4291
4292                             xy[x1*2] = saturate_cast<short>(X >> INTER_BITS);
4293                             xy[x1*2+1] = saturate_cast<short>(Y >> INTER_BITS);
4294                             alpha[x1] = (short)((Y & (INTER_TAB_SIZE-1))*INTER_TAB_SIZE +
4295                                                 (X & (INTER_TAB_SIZE-1)));
4296                         }
4297                     }
4298                 }
4299
4300                 if( interpolation == INTER_NEAREST )
4301                     remap( src, dpart, _XY, Mat(), interpolation, borderType, borderValue );
4302                 else
4303                 {
4304                     Mat _matA(bh, bw, CV_16U, A);
4305                     remap( src, dpart, _XY, _matA, interpolation, borderType, borderValue );
4306                 }
4307             }
4308         }
4309     }
4310
4311 private:
4312     Mat src;
4313     Mat dst;
4314     double* M;
4315     int interpolation, borderType;
4316     Scalar borderValue;
4317 };
4318
4319 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
4320 class IPPwarpPerspectiveInvoker :
4321     public ParallelLoopBody
4322 {
4323 public:
4324     IPPwarpPerspectiveInvoker(Mat &_src, Mat &_dst, double (&_coeffs)[3][3], int &_interpolation, int &_borderType, const Scalar &_borderValue, ippiWarpPerspectiveBackFunc _func, bool *_ok) :
4325       ParallelLoopBody(), src(_src), dst(_dst), mode(_interpolation), coeffs(_coeffs), borderType(_borderType), borderValue(_borderValue), func(_func), ok(_ok)
4326       {
4327           *ok = true;
4328       }
4329
4330       virtual void operator() (const Range& range) const
4331       {
4332           IppiSize srcsize = {src.cols, src.rows};
4333           IppiRect srcroi = {0, 0, src.cols, src.rows};
4334           IppiRect dstroi = {0, range.start, dst.cols, range.end - range.start};
4335           int cnn = src.channels();
4336
4337           if( borderType == BORDER_CONSTANT )
4338           {
4339               IppiSize setSize = {dst.cols, range.end - range.start};
4340               void *dataPointer = dst.data + dst.step[0] * range.start;
4341               if( !IPPSet( borderValue, dataPointer, (int)dst.step[0], setSize, cnn, src.depth() ) )
4342               {
4343                   *ok = false;
4344                   return;
4345               }
4346           }
4347           if( func(src.data, srcsize, (int)src.step[0], srcroi, dst.data, (int)dst.step[0], dstroi, coeffs, mode) < 0)
4348               *ok = false;
4349       }
4350 private:
4351     Mat &src;
4352     Mat &dst;
4353     double (&coeffs)[3][3];
4354     int mode;
4355     int borderType;
4356     const Scalar borderValue;
4357     ippiWarpPerspectiveBackFunc func;
4358     bool *ok;
4359     const IPPwarpPerspectiveInvoker& operator= (const IPPwarpPerspectiveInvoker&);
4360 };
4361 #endif
4362
4363 }
4364
4365 void cv::warpPerspective( InputArray _src, OutputArray _dst, InputArray _M0,
4366                           Size dsize, int flags, int borderType, const Scalar& borderValue )
4367 {
4368     CV_Assert( _src.total() > 0 );
4369
4370     CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
4371                ocl_warpTransform(_src, _dst, _M0, dsize, flags, borderType, borderValue,
4372                               OCL_OP_PERSPECTIVE))
4373
4374     Mat src = _src.getMat(), M0 = _M0.getMat();
4375     _dst.create( dsize.area() == 0 ? src.size() : dsize, src.type() );
4376     Mat dst = _dst.getMat();
4377
4378     if( dst.data == src.data )
4379         src = src.clone();
4380
4381     double M[9];
4382     Mat matM(3, 3, CV_64F, M);
4383     int interpolation = flags & INTER_MAX;
4384     if( interpolation == INTER_AREA )
4385         interpolation = INTER_LINEAR;
4386
4387     CV_Assert( (M0.type() == CV_32F || M0.type() == CV_64F) && M0.rows == 3 && M0.cols == 3 );
4388     M0.convertTo(matM, matM.type());
4389
4390 #ifdef HAVE_TEGRA_OPTIMIZATION
4391     if( tegra::warpPerspective(src, dst, M, flags, borderType, borderValue) )
4392         return;
4393 #endif
4394
4395     if( !(flags & WARP_INVERSE_MAP) )
4396          invert(matM, matM);
4397 /*
4398 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
4399     int depth = src.depth();
4400     int channels = src.channels();
4401     if( ( depth == CV_8U || depth == CV_16U || depth == CV_32F ) &&
4402         ( channels == 1 || channels == 3 || channels == 4 ) &&
4403         ( borderType == cv::BORDER_TRANSPARENT || borderType == cv::BORDER_CONSTANT ) )
4404     {
4405         int type = src.type();
4406         ippiWarpPerspectiveBackFunc ippFunc =
4407             type == CV_8UC1 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_8u_C1R :
4408             type == CV_8UC3 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_8u_C3R :
4409             type == CV_8UC4 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_8u_C4R :
4410             type == CV_16UC1 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_16u_C1R :
4411             type == CV_16UC3 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_16u_C3R :
4412             type == CV_16UC4 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_16u_C4R :
4413             type == CV_32FC1 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_32f_C1R :
4414             type == CV_32FC3 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_32f_C3R :
4415             type == CV_32FC4 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_32f_C4R :
4416             0;
4417         int mode =
4418             flags == INTER_LINEAR ? IPPI_INTER_LINEAR :
4419             flags == INTER_NEAREST ? IPPI_INTER_NN :
4420             flags == INTER_CUBIC ? IPPI_INTER_CUBIC :
4421             0;
4422         if( mode && ippFunc )
4423         {
4424             double coeffs[3][3];
4425             for( int i = 0; i < 3; i++ )
4426             {
4427                 for( int j = 0; j < 3; j++ )
4428                 {
4429                     coeffs[i][j] = matM.at<double>(i, j);
4430                 }
4431             }
4432             bool ok;
4433             Range range(0, dst.rows);
4434             IPPwarpPerspectiveInvoker invoker(src, dst, coeffs, mode, borderType, borderValue, ippFunc, &ok);
4435             parallel_for_(range, invoker, dst.total()/(double)(1<<16));
4436             if( ok )
4437                 return;
4438         }
4439     }
4440 #endif
4441 */
4442     Range range(0, dst.rows);
4443     warpPerspectiveInvoker invoker(src, dst, M, interpolation, borderType, borderValue);
4444     parallel_for_(range, invoker, dst.total()/(double)(1<<16));
4445 }
4446
4447
4448 cv::Mat cv::getRotationMatrix2D( Point2f center, double angle, double scale )
4449 {
4450     angle *= CV_PI/180;
4451     double alpha = cos(angle)*scale;
4452     double beta = sin(angle)*scale;
4453
4454     Mat M(2, 3, CV_64F);
4455     double* m = (double*)M.data;
4456
4457     m[0] = alpha;
4458     m[1] = beta;
4459     m[2] = (1-alpha)*center.x - beta*center.y;
4460     m[3] = -beta;
4461     m[4] = alpha;
4462     m[5] = beta*center.x + (1-alpha)*center.y;
4463
4464     return M;
4465 }
4466
4467 /* Calculates coefficients of perspective transformation
4468  * which maps (xi,yi) to (ui,vi), (i=1,2,3,4):
4469  *
4470  *      c00*xi + c01*yi + c02
4471  * ui = ---------------------
4472  *      c20*xi + c21*yi + c22
4473  *
4474  *      c10*xi + c11*yi + c12
4475  * vi = ---------------------
4476  *      c20*xi + c21*yi + c22
4477  *
4478  * Coefficients are calculated by solving linear system:
4479  * / x0 y0  1  0  0  0 -x0*u0 -y0*u0 \ /c00\ /u0\
4480  * | x1 y1  1  0  0  0 -x1*u1 -y1*u1 | |c01| |u1|
4481  * | x2 y2  1  0  0  0 -x2*u2 -y2*u2 | |c02| |u2|
4482  * | x3 y3  1  0  0  0 -x3*u3 -y3*u3 |.|c10|=|u3|,
4483  * |  0  0  0 x0 y0  1 -x0*v0 -y0*v0 | |c11| |v0|
4484  * |  0  0  0 x1 y1  1 -x1*v1 -y1*v1 | |c12| |v1|
4485  * |  0  0  0 x2 y2  1 -x2*v2 -y2*v2 | |c20| |v2|
4486  * \  0  0  0 x3 y3  1 -x3*v3 -y3*v3 / \c21/ \v3/
4487  *
4488  * where:
4489  *   cij - matrix coefficients, c22 = 1
4490  */
4491 cv::Mat cv::getPerspectiveTransform( const Point2f src[], const Point2f dst[] )
4492 {
4493     Mat M(3, 3, CV_64F), X(8, 1, CV_64F, M.data);
4494     double a[8][8], b[8];
4495     Mat A(8, 8, CV_64F, a), B(8, 1, CV_64F, b);
4496
4497     for( int i = 0; i < 4; ++i )
4498     {
4499         a[i][0] = a[i+4][3] = src[i].x;
4500         a[i][1] = a[i+4][4] = src[i].y;
4501         a[i][2] = a[i+4][5] = 1;
4502         a[i][3] = a[i][4] = a[i][5] =
4503         a[i+4][0] = a[i+4][1] = a[i+4][2] = 0;
4504         a[i][6] = -src[i].x*dst[i].x;
4505         a[i][7] = -src[i].y*dst[i].x;
4506         a[i+4][6] = -src[i].x*dst[i].y;
4507         a[i+4][7] = -src[i].y*dst[i].y;
4508         b[i] = dst[i].x;
4509         b[i+4] = dst[i].y;
4510     }
4511
4512     solve( A, B, X, DECOMP_SVD );
4513     ((double*)M.data)[8] = 1.;
4514
4515     return M;
4516 }
4517
4518 /* Calculates coefficients of affine transformation
4519  * which maps (xi,yi) to (ui,vi), (i=1,2,3):
4520  *
4521  * ui = c00*xi + c01*yi + c02
4522  *
4523  * vi = c10*xi + c11*yi + c12
4524  *
4525  * Coefficients are calculated by solving linear system:
4526  * / x0 y0  1  0  0  0 \ /c00\ /u0\
4527  * | x1 y1  1  0  0  0 | |c01| |u1|
4528  * | x2 y2  1  0  0  0 | |c02| |u2|
4529  * |  0  0  0 x0 y0  1 | |c10| |v0|
4530  * |  0  0  0 x1 y1  1 | |c11| |v1|
4531  * \  0  0  0 x2 y2  1 / |c12| |v2|
4532  *
4533  * where:
4534  *   cij - matrix coefficients
4535  */
4536
4537 cv::Mat cv::getAffineTransform( const Point2f src[], const Point2f dst[] )
4538 {
4539     Mat M(2, 3, CV_64F), X(6, 1, CV_64F, M.data);
4540     double a[6*6], b[6];
4541     Mat A(6, 6, CV_64F, a), B(6, 1, CV_64F, b);
4542
4543     for( int i = 0; i < 3; i++ )
4544     {
4545         int j = i*12;
4546         int k = i*12+6;
4547         a[j] = a[k+3] = src[i].x;
4548         a[j+1] = a[k+4] = src[i].y;
4549         a[j+2] = a[k+5] = 1;
4550         a[j+3] = a[j+4] = a[j+5] = 0;
4551         a[k] = a[k+1] = a[k+2] = 0;
4552         b[i*2] = dst[i].x;
4553         b[i*2+1] = dst[i].y;
4554     }
4555
4556     solve( A, B, X );
4557     return M;
4558 }
4559
4560 void cv::invertAffineTransform(InputArray _matM, OutputArray __iM)
4561 {
4562     Mat matM = _matM.getMat();
4563     CV_Assert(matM.rows == 2 && matM.cols == 3);
4564     __iM.create(2, 3, matM.type());
4565     Mat _iM = __iM.getMat();
4566
4567     if( matM.type() == CV_32F )
4568     {
4569         const float* M = (const float*)matM.data;
4570         float* iM = (float*)_iM.data;
4571         int step = (int)(matM.step/sizeof(M[0])), istep = (int)(_iM.step/sizeof(iM[0]));
4572
4573         double D = M[0]*M[step+1] - M[1]*M[step];
4574         D = D != 0 ? 1./D : 0;
4575         double A11 = M[step+1]*D, A22 = M[0]*D, A12 = -M[1]*D, A21 = -M[step]*D;
4576         double b1 = -A11*M[2] - A12*M[step+2];
4577         double b2 = -A21*M[2] - A22*M[step+2];
4578
4579         iM[0] = (float)A11; iM[1] = (float)A12; iM[2] = (float)b1;
4580         iM[istep] = (float)A21; iM[istep+1] = (float)A22; iM[istep+2] = (float)b2;
4581     }
4582     else if( matM.type() == CV_64F )
4583     {
4584         const double* M = (const double*)matM.data;
4585         double* iM = (double*)_iM.data;
4586         int step = (int)(matM.step/sizeof(M[0])), istep = (int)(_iM.step/sizeof(iM[0]));
4587
4588         double D = M[0]*M[step+1] - M[1]*M[step];
4589         D = D != 0 ? 1./D : 0;
4590         double A11 = M[step+1]*D, A22 = M[0]*D, A12 = -M[1]*D, A21 = -M[step]*D;
4591         double b1 = -A11*M[2] - A12*M[step+2];
4592         double b2 = -A21*M[2] - A22*M[step+2];
4593
4594         iM[0] = A11; iM[1] = A12; iM[2] = b1;
4595         iM[istep] = A21; iM[istep+1] = A22; iM[istep+2] = b2;
4596     }
4597     else
4598         CV_Error( CV_StsUnsupportedFormat, "" );
4599 }
4600
4601 cv::Mat cv::getPerspectiveTransform(InputArray _src, InputArray _dst)
4602 {
4603     Mat src = _src.getMat(), dst = _dst.getMat();
4604     CV_Assert(src.checkVector(2, CV_32F) == 4 && dst.checkVector(2, CV_32F) == 4);
4605     return getPerspectiveTransform((const Point2f*)src.data, (const Point2f*)dst.data);
4606 }
4607
4608 cv::Mat cv::getAffineTransform(InputArray _src, InputArray _dst)
4609 {
4610     Mat src = _src.getMat(), dst = _dst.getMat();
4611     CV_Assert(src.checkVector(2, CV_32F) == 3 && dst.checkVector(2, CV_32F) == 3);
4612     return getAffineTransform((const Point2f*)src.data, (const Point2f*)dst.data);
4613 }
4614
4615 CV_IMPL void
4616 cvResize( const CvArr* srcarr, CvArr* dstarr, int method )
4617 {
4618     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
4619     CV_Assert( src.type() == dst.type() );
4620     cv::resize( src, dst, dst.size(), (double)dst.cols/src.cols,
4621         (double)dst.rows/src.rows, method );
4622 }
4623
4624
4625 CV_IMPL void
4626 cvWarpAffine( const CvArr* srcarr, CvArr* dstarr, const CvMat* marr,
4627               int flags, CvScalar fillval )
4628 {
4629     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
4630     cv::Mat matrix = cv::cvarrToMat(marr);
4631     CV_Assert( src.type() == dst.type() );
4632     cv::warpAffine( src, dst, matrix, dst.size(), flags,
4633         (flags & CV_WARP_FILL_OUTLIERS) ? cv::BORDER_CONSTANT : cv::BORDER_TRANSPARENT,
4634         fillval );
4635 }
4636
4637 CV_IMPL void
4638 cvWarpPerspective( const CvArr* srcarr, CvArr* dstarr, const CvMat* marr,
4639                    int flags, CvScalar fillval )
4640 {
4641     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
4642     cv::Mat matrix = cv::cvarrToMat(marr);
4643     CV_Assert( src.type() == dst.type() );
4644     cv::warpPerspective( src, dst, matrix, dst.size(), flags,
4645         (flags & CV_WARP_FILL_OUTLIERS) ? cv::BORDER_CONSTANT : cv::BORDER_TRANSPARENT,
4646         fillval );
4647 }
4648
4649 CV_IMPL void
4650 cvRemap( const CvArr* srcarr, CvArr* dstarr,
4651          const CvArr* _mapx, const CvArr* _mapy,
4652          int flags, CvScalar fillval )
4653 {
4654     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr), dst0 = dst;
4655     cv::Mat mapx = cv::cvarrToMat(_mapx), mapy = cv::cvarrToMat(_mapy);
4656     CV_Assert( src.type() == dst.type() && dst.size() == mapx.size() );
4657     cv::remap( src, dst, mapx, mapy, flags & cv::INTER_MAX,
4658         (flags & CV_WARP_FILL_OUTLIERS) ? cv::BORDER_CONSTANT : cv::BORDER_TRANSPARENT,
4659         fillval );
4660     CV_Assert( dst0.data == dst.data );
4661 }
4662
4663
4664 CV_IMPL CvMat*
4665 cv2DRotationMatrix( CvPoint2D32f center, double angle,
4666                     double scale, CvMat* matrix )
4667 {
4668     cv::Mat M0 = cv::cvarrToMat(matrix), M = cv::getRotationMatrix2D(center, angle, scale);
4669     CV_Assert( M.size() == M0.size() );
4670     M.convertTo(M0, M0.type());
4671     return matrix;
4672 }
4673
4674
4675 CV_IMPL CvMat*
4676 cvGetPerspectiveTransform( const CvPoint2D32f* src,
4677                           const CvPoint2D32f* dst,
4678                           CvMat* matrix )
4679 {
4680     cv::Mat M0 = cv::cvarrToMat(matrix),
4681         M = cv::getPerspectiveTransform((const cv::Point2f*)src, (const cv::Point2f*)dst);
4682     CV_Assert( M.size() == M0.size() );
4683     M.convertTo(M0, M0.type());
4684     return matrix;
4685 }
4686
4687
4688 CV_IMPL CvMat*
4689 cvGetAffineTransform( const CvPoint2D32f* src,
4690                           const CvPoint2D32f* dst,
4691                           CvMat* matrix )
4692 {
4693     cv::Mat M0 = cv::cvarrToMat(matrix),
4694         M = cv::getAffineTransform((const cv::Point2f*)src, (const cv::Point2f*)dst);
4695     CV_Assert( M.size() == M0.size() );
4696     M.convertTo(M0, M0.type());
4697     return matrix;
4698 }
4699
4700
4701 CV_IMPL void
4702 cvConvertMaps( const CvArr* arr1, const CvArr* arr2, CvArr* dstarr1, CvArr* dstarr2 )
4703 {
4704     cv::Mat map1 = cv::cvarrToMat(arr1), map2;
4705     cv::Mat dstmap1 = cv::cvarrToMat(dstarr1), dstmap2;
4706
4707     if( arr2 )
4708         map2 = cv::cvarrToMat(arr2);
4709     if( dstarr2 )
4710     {
4711         dstmap2 = cv::cvarrToMat(dstarr2);
4712         if( dstmap2.type() == CV_16SC1 )
4713             dstmap2 = cv::Mat(dstmap2.size(), CV_16UC1, dstmap2.data, dstmap2.step);
4714     }
4715
4716     cv::convertMaps( map1, map2, dstmap1, dstmap2, dstmap1.type(), false );
4717 }
4718
4719 /****************************************************************************************\
4720 *                                   Log-Polar Transform                                  *
4721 \****************************************************************************************/
4722
4723 /* now it is done via Remap; more correct implementation should use
4724    some super-sampling technique outside of the "fovea" circle */
4725 CV_IMPL void
4726 cvLogPolar( const CvArr* srcarr, CvArr* dstarr,
4727             CvPoint2D32f center, double M, int flags )
4728 {
4729     cv::Ptr<CvMat> mapx, mapy;
4730
4731     CvMat srcstub, *src = cvGetMat(srcarr, &srcstub);
4732     CvMat dststub, *dst = cvGetMat(dstarr, &dststub);
4733     CvSize ssize, dsize;
4734
4735     if( !CV_ARE_TYPES_EQ( src, dst ))
4736         CV_Error( CV_StsUnmatchedFormats, "" );
4737
4738     if( M <= 0 )
4739         CV_Error( CV_StsOutOfRange, "M should be >0" );
4740
4741     ssize = cvGetMatSize(src);
4742     dsize = cvGetMatSize(dst);
4743
4744     mapx.reset(cvCreateMat( dsize.height, dsize.width, CV_32F ));
4745     mapy.reset(cvCreateMat( dsize.height, dsize.width, CV_32F ));
4746
4747     if( !(flags & CV_WARP_INVERSE_MAP) )
4748     {
4749         int phi, rho;
4750         cv::AutoBuffer<double> _exp_tab(dsize.width);
4751         double* exp_tab = _exp_tab;
4752
4753         for( rho = 0; rho < dst->width; rho++ )
4754             exp_tab[rho] = std::exp(rho/M);
4755
4756         for( phi = 0; phi < dsize.height; phi++ )
4757         {
4758             double cp = cos(phi*2*CV_PI/dsize.height);
4759             double sp = sin(phi*2*CV_PI/dsize.height);
4760             float* mx = (float*)(mapx->data.ptr + phi*mapx->step);
4761             float* my = (float*)(mapy->data.ptr + phi*mapy->step);
4762
4763             for( rho = 0; rho < dsize.width; rho++ )
4764             {
4765                 double r = exp_tab[rho];
4766                 double x = r*cp + center.x;
4767                 double y = r*sp + center.y;
4768
4769                 mx[rho] = (float)x;
4770                 my[rho] = (float)y;
4771             }
4772         }
4773     }
4774     else
4775     {
4776         int x, y;
4777         CvMat bufx, bufy, bufp, bufa;
4778         double ascale = ssize.height/(2*CV_PI);
4779         cv::AutoBuffer<float> _buf(4*dsize.width);
4780         float* buf = _buf;
4781
4782         bufx = cvMat( 1, dsize.width, CV_32F, buf );
4783         bufy = cvMat( 1, dsize.width, CV_32F, buf + dsize.width );
4784         bufp = cvMat( 1, dsize.width, CV_32F, buf + dsize.width*2 );
4785         bufa = cvMat( 1, dsize.width, CV_32F, buf + dsize.width*3 );
4786
4787         for( x = 0; x < dsize.width; x++ )
4788             bufx.data.fl[x] = (float)x - center.x;
4789
4790         for( y = 0; y < dsize.height; y++ )
4791         {
4792             float* mx = (float*)(mapx->data.ptr + y*mapx->step);
4793             float* my = (float*)(mapy->data.ptr + y*mapy->step);
4794
4795             for( x = 0; x < dsize.width; x++ )
4796                 bufy.data.fl[x] = (float)y - center.y;
4797
4798 #if 1
4799             cvCartToPolar( &bufx, &bufy, &bufp, &bufa );
4800
4801             for( x = 0; x < dsize.width; x++ )
4802                 bufp.data.fl[x] += 1.f;
4803
4804             cvLog( &bufp, &bufp );
4805
4806             for( x = 0; x < dsize.width; x++ )
4807             {
4808                 double rho = bufp.data.fl[x]*M;
4809                 double phi = bufa.data.fl[x]*ascale;
4810
4811                 mx[x] = (float)rho;
4812                 my[x] = (float)phi;
4813             }
4814 #else
4815             for( x = 0; x < dsize.width; x++ )
4816             {
4817                 double xx = bufx.data.fl[x];
4818                 double yy = bufy.data.fl[x];
4819
4820                 double p = log(std::sqrt(xx*xx + yy*yy) + 1.)*M;
4821                 double a = atan2(yy,xx);
4822                 if( a < 0 )
4823                     a = 2*CV_PI + a;
4824                 a *= ascale;
4825
4826                 mx[x] = (float)p;
4827                 my[x] = (float)a;
4828             }
4829 #endif
4830         }
4831     }
4832
4833     cvRemap( src, dst, mapx, mapy, flags, cvScalarAll(0) );
4834 }
4835
4836 void cv::logPolar( InputArray _src, OutputArray _dst,
4837                    Point2f center, double M, int flags )
4838 {
4839     Mat src = _src.getMat();
4840     _dst.create( src.size(), src.type() );
4841     CvMat c_src = src, c_dst = _dst.getMat();
4842     cvLogPolar( &c_src, &c_dst, center, M, flags );
4843 }
4844
4845 /****************************************************************************************
4846                                    Linear-Polar Transform
4847   J.L. Blanco, Apr 2009
4848  ****************************************************************************************/
4849 CV_IMPL
4850 void cvLinearPolar( const CvArr* srcarr, CvArr* dstarr,
4851             CvPoint2D32f center, double maxRadius, int flags )
4852 {
4853     cv::Ptr<CvMat> mapx, mapy;
4854
4855     CvMat srcstub, *src = (CvMat*)srcarr;
4856     CvMat dststub, *dst = (CvMat*)dstarr;
4857     CvSize ssize, dsize;
4858
4859     src = cvGetMat( srcarr, &srcstub,0,0 );
4860     dst = cvGetMat( dstarr, &dststub,0,0 );
4861
4862     if( !CV_ARE_TYPES_EQ( src, dst ))
4863         CV_Error( CV_StsUnmatchedFormats, "" );
4864
4865     ssize.width = src->cols;
4866     ssize.height = src->rows;
4867     dsize.width = dst->cols;
4868     dsize.height = dst->rows;
4869
4870     mapx.reset(cvCreateMat( dsize.height, dsize.width, CV_32F ));
4871     mapy.reset(cvCreateMat( dsize.height, dsize.width, CV_32F ));
4872
4873     if( !(flags & CV_WARP_INVERSE_MAP) )
4874     {
4875         int phi, rho;
4876
4877         for( phi = 0; phi < dsize.height; phi++ )
4878         {
4879             double cp = cos(phi*2*CV_PI/dsize.height);
4880             double sp = sin(phi*2*CV_PI/dsize.height);
4881             float* mx = (float*)(mapx->data.ptr + phi*mapx->step);
4882             float* my = (float*)(mapy->data.ptr + phi*mapy->step);
4883
4884             for( rho = 0; rho < dsize.width; rho++ )
4885             {
4886                 double r = maxRadius*(rho+1)/dsize.width;
4887                 double x = r*cp + center.x;
4888                 double y = r*sp + center.y;
4889
4890                 mx[rho] = (float)x;
4891                 my[rho] = (float)y;
4892             }
4893         }
4894     }
4895     else
4896     {
4897         int x, y;
4898         CvMat bufx, bufy, bufp, bufa;
4899         const double ascale = ssize.height/(2*CV_PI);
4900         const double pscale = ssize.width/maxRadius;
4901
4902         cv::AutoBuffer<float> _buf(4*dsize.width);
4903         float* buf = _buf;
4904
4905         bufx = cvMat( 1, dsize.width, CV_32F, buf );
4906         bufy = cvMat( 1, dsize.width, CV_32F, buf + dsize.width );
4907         bufp = cvMat( 1, dsize.width, CV_32F, buf + dsize.width*2 );
4908         bufa = cvMat( 1, dsize.width, CV_32F, buf + dsize.width*3 );
4909
4910         for( x = 0; x < dsize.width; x++ )
4911             bufx.data.fl[x] = (float)x - center.x;
4912
4913         for( y = 0; y < dsize.height; y++ )
4914         {
4915             float* mx = (float*)(mapx->data.ptr + y*mapx->step);
4916             float* my = (float*)(mapy->data.ptr + y*mapy->step);
4917
4918             for( x = 0; x < dsize.width; x++ )
4919                 bufy.data.fl[x] = (float)y - center.y;
4920
4921             cvCartToPolar( &bufx, &bufy, &bufp, &bufa, 0 );
4922
4923             for( x = 0; x < dsize.width; x++ )
4924                 bufp.data.fl[x] += 1.f;
4925
4926             for( x = 0; x < dsize.width; x++ )
4927             {
4928                 double rho = bufp.data.fl[x]*pscale;
4929                 double phi = bufa.data.fl[x]*ascale;
4930                 mx[x] = (float)rho;
4931                 my[x] = (float)phi;
4932             }
4933         }
4934     }
4935
4936     cvRemap( src, dst, mapx, mapy, flags, cvScalarAll(0) );
4937 }
4938
4939 void cv::linearPolar( InputArray _src, OutputArray _dst,
4940                       Point2f center, double maxRadius, int flags )
4941 {
4942     Mat src = _src.getMat();
4943     _dst.create( src.size(), src.type() );
4944     CvMat c_src = src, c_dst = _dst.getMat();
4945     cvLinearPolar( &c_src, &c_dst, center, maxRadius, flags );
4946 }
4947
4948 /* End of file. */