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 == 7 && IPP_VERSION_MINOR >= 1) || IPP_VERSION_MAJOR > 7)
  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 #if defined (HAVE_IPP) && ((IPP_VERSION_MAJOR == 7 && IPP_VERSION_MINOR >= 1) || IPP_VERSION_MAJOR > 7)
1868 class IPPresizeInvoker :
1869     public ParallelLoopBody
1870 {
1871 public:
1872     IPPresizeInvoker(Mat &_src, Mat &_dst, double _inv_scale_x, double _inv_scale_y, int _mode, bool *_ok) :
1873       ParallelLoopBody(), src(_src), dst(_dst), inv_scale_x(_inv_scale_x), inv_scale_y(_inv_scale_y), mode(_mode), ok(_ok)
1874       {
1875           IppStatus status = ippStsNotSupportedModeErr;
1876           *ok = true;
1877           IppiSize srcSize, dstSize;
1878           int type = src.type();
1879           int specSize = 0, initSize = 0;
1880           srcSize.width  = src.cols;
1881           srcSize.height = src.rows;
1882           dstSize.width  = dst.cols;
1883           dstSize.height = dst.rows;
1884
1885           if (mode == (int)ippLinear)
1886           {
1887               func =
1888                   type == CV_8UC1  ? (ippiResizeFunc)ippiResizeLinear_8u_C1R :
1889                   type == CV_8UC3  ? (ippiResizeFunc)ippiResizeLinear_8u_C3R :
1890                   type == CV_8UC4  ? (ippiResizeFunc)ippiResizeLinear_8u_C4R :
1891                   type == CV_16UC1 ? (ippiResizeFunc)ippiResizeLinear_16u_C1R :
1892                   type == CV_16UC3 ? (ippiResizeFunc)ippiResizeLinear_16u_C3R :
1893                   type == CV_16UC4 ? (ippiResizeFunc)ippiResizeLinear_16u_C4R :
1894                   type == CV_16SC1 ? (ippiResizeFunc)ippiResizeLinear_16s_C1R :
1895                   type == CV_16SC3 ? (ippiResizeFunc)ippiResizeLinear_16s_C3R :
1896                   type == CV_16SC4 ? (ippiResizeFunc)ippiResizeLinear_16s_C4R :
1897                   type == CV_32FC1 ? (ippiResizeFunc)ippiResizeLinear_32f_C1R :
1898                   type == CV_32FC3 ? (ippiResizeFunc)ippiResizeLinear_32f_C3R :
1899                   type == CV_32FC4 ? (ippiResizeFunc)ippiResizeLinear_32f_C4R :
1900                   type == CV_64FC1 ? (ippiResizeFunc)ippiResizeLinear_64f_C1R :
1901                   type == CV_64FC3 ? (ippiResizeFunc)ippiResizeLinear_64f_C3R :
1902                   type == CV_64FC4 ? (ippiResizeFunc)ippiResizeLinear_64f_C4R :
1903                   0;
1904           }
1905           else if (mode == (int)ippCubic)
1906           {
1907               func =
1908                   type == CV_8UC1  ? (ippiResizeFunc)ippiResizeCubic_8u_C1R :
1909                   type == CV_8UC3  ? (ippiResizeFunc)ippiResizeCubic_8u_C3R :
1910                   type == CV_8UC4  ? (ippiResizeFunc)ippiResizeCubic_8u_C4R :
1911                   type == CV_16UC1 ? (ippiResizeFunc)ippiResizeCubic_16u_C1R :
1912                   type == CV_16UC3 ? (ippiResizeFunc)ippiResizeCubic_16u_C3R :
1913                   type == CV_16UC4 ? (ippiResizeFunc)ippiResizeCubic_16u_C4R :
1914                   type == CV_16SC1 ? (ippiResizeFunc)ippiResizeCubic_16s_C1R :
1915                   type == CV_16SC3 ? (ippiResizeFunc)ippiResizeCubic_16s_C3R :
1916                   type == CV_16SC4 ? (ippiResizeFunc)ippiResizeCubic_16s_C4R :
1917                   type == CV_32FC1 ? (ippiResizeFunc)ippiResizeCubic_32f_C1R :
1918                   type == CV_32FC3 ? (ippiResizeFunc)ippiResizeCubic_32f_C3R :
1919                   type == CV_32FC4 ? (ippiResizeFunc)ippiResizeCubic_32f_C4R :
1920                   0;
1921           }
1922
1923           switch (src.depth())
1924           {
1925           case CV_8U:
1926               status = ippiResizeGetSize_8u(srcSize, dstSize, (IppiInterpolationType)mode, 0, &specSize, &initSize);
1927               break;
1928           case CV_16U:
1929               status = ippiResizeGetSize_16u(srcSize, dstSize, (IppiInterpolationType)mode, 0, &specSize, &initSize);
1930               break;
1931           case CV_16S:
1932               status = ippiResizeGetSize_16s(srcSize, dstSize, (IppiInterpolationType)mode, 0, &specSize, &initSize);
1933               break;
1934           case CV_32F:
1935               status = ippiResizeGetSize_32f(srcSize, dstSize, (IppiInterpolationType)mode, 0, &specSize, &initSize);
1936               break;
1937           case CV_64F:
1938               status = ippiResizeGetSize_64f(srcSize, dstSize, (IppiInterpolationType)mode, 0, &specSize, &initSize);
1939               break;
1940           }
1941           if (status != ippStsNoErr)
1942           {
1943               *ok = false;
1944               return;
1945           }
1946
1947           specBuf.allocate(specSize);
1948           pSpec = (uchar*)specBuf;
1949
1950           status = ippStsNotSupportedModeErr;
1951           if (mode == (int)ippLinear)
1952           {
1953               switch (src.depth())
1954               {
1955               case CV_8U:
1956                   status = ippiResizeLinearInit_8u(srcSize, dstSize, (IppiResizeSpec_32f*)pSpec);
1957                   break;
1958               case CV_16U:
1959                   status = ippiResizeLinearInit_16u(srcSize, dstSize, (IppiResizeSpec_32f*)pSpec);
1960                   break;
1961               case CV_16S:
1962                   status = ippiResizeLinearInit_16s(srcSize, dstSize, (IppiResizeSpec_32f*)pSpec);
1963                   break;
1964               case CV_32F:
1965                   status = ippiResizeLinearInit_32f(srcSize, dstSize, (IppiResizeSpec_32f*)pSpec);
1966                   break;
1967               case CV_64F:
1968                   status = ippiResizeLinearInit_64f(srcSize, dstSize, (IppiResizeSpec_64f*)pSpec);
1969                   break;
1970               }
1971               if (status != ippStsNoErr)
1972               {
1973                   *ok = false;
1974                   return;
1975               }
1976           }
1977           else if (mode == (int)ippCubic)
1978           {
1979               AutoBuffer<uchar> buf(initSize);
1980               uchar* pInit = (uchar*)buf;
1981
1982               switch (src.depth())
1983               {
1984               case CV_8U:
1985                   status = ippiResizeCubicInit_8u(srcSize, dstSize,  0.f, 0.75f, (IppiResizeSpec_32f*)pSpec, pInit);
1986                   break;
1987               case CV_16U:
1988                   status = ippiResizeCubicInit_16u(srcSize, dstSize, 0.f, 0.75f, (IppiResizeSpec_32f*)pSpec, pInit);
1989                   break;
1990               case CV_16S:
1991                   status = ippiResizeCubicInit_16s(srcSize, dstSize, 0.f, 0.75f, (IppiResizeSpec_32f*)pSpec, pInit);
1992                   break;
1993               case CV_32F:
1994                   status = ippiResizeCubicInit_32f(srcSize, dstSize, 0.f, 0.75f, (IppiResizeSpec_32f*)pSpec, pInit);
1995                   break;
1996               }
1997               if (status != ippStsNoErr) *ok = false;
1998           }
1999       }
2000
2001       ~IPPresizeInvoker()
2002       {
2003       }
2004
2005       virtual void operator() (const Range& range) const
2006       {
2007           if (*ok == false) return;
2008
2009           int cn = src.channels();
2010           int dsty = CV_IMIN(cvRound(range.start * inv_scale_y), dst.rows);
2011           int dstwidth  = CV_IMIN(cvRound(src.cols * inv_scale_x), dst.cols);
2012           int dstheight = CV_IMIN(cvRound(range.end * inv_scale_y), dst.rows);
2013
2014           IppiPoint dstOffset = { 0, dsty }, srcOffset = {0, 0};
2015           IppiSize  dstSize   = { dstwidth, dstheight - dsty };
2016           int bufsize = 0, itemSize = 0;
2017
2018           IppStatus status = ippStsNotSupportedModeErr;
2019
2020           switch (src.depth())
2021           {
2022           case CV_8U:
2023               itemSize = 1;
2024               status = ippiResizeGetBufferSize_8u((IppiResizeSpec_32f*)pSpec, dstSize, cn, &bufsize);
2025               if (status == ippStsNoErr)
2026                   status = ippiResizeGetSrcOffset_8u((IppiResizeSpec_32f*)pSpec, dstOffset, &srcOffset);
2027               break;
2028           case CV_16U:
2029               itemSize = 2;
2030               status = ippiResizeGetBufferSize_16u((IppiResizeSpec_32f*)pSpec, dstSize, cn, &bufsize);
2031               if (status == ippStsNoErr)
2032                   status = ippiResizeGetSrcOffset_16u((IppiResizeSpec_32f*)pSpec, dstOffset, &srcOffset);
2033               break;
2034           case CV_16S:
2035               itemSize = 2;
2036               status = ippiResizeGetBufferSize_16s((IppiResizeSpec_32f*)pSpec, dstSize, cn, &bufsize);
2037               if (status == ippStsNoErr)
2038                   status = ippiResizeGetSrcOffset_16s((IppiResizeSpec_32f*)pSpec, dstOffset, &srcOffset);
2039               break;
2040           case CV_32F:
2041               itemSize = 4;
2042               status = ippiResizeGetBufferSize_32f((IppiResizeSpec_32f*)pSpec, dstSize, cn, &bufsize);
2043               if (status == ippStsNoErr)
2044                   status = ippiResizeGetSrcOffset_32f((IppiResizeSpec_32f*)pSpec, dstOffset, &srcOffset);
2045               break;
2046           case CV_64F:
2047               itemSize = 4;
2048               status = ippiResizeGetBufferSize_64f((IppiResizeSpec_64f*)pSpec, dstSize, cn, &bufsize);
2049               if (status == ippStsNoErr)
2050                   status = ippiResizeGetSrcOffset_64f((IppiResizeSpec_64f*)pSpec, dstOffset, &srcOffset);
2051               break;
2052           }
2053           if (status != ippStsNoErr)
2054           {
2055               *ok = false;
2056               return;
2057           }
2058
2059           Ipp8u* pSrc = (Ipp8u*)src.data + (int)src.step[0] * srcOffset.y + srcOffset.x * cn * itemSize;
2060           Ipp8u* pDst = (Ipp8u*)dst.data + (int)dst.step[0] * dstOffset.y + dstOffset.x * cn * itemSize;
2061
2062           AutoBuffer<uchar> buf(bufsize + 64);
2063           uchar* bufptr = alignPtr((uchar*)buf, 32);
2064           if( func( pSrc, (int)src.step[0], pDst, (int)dst.step[0], dstOffset, dstSize, ippBorderRepl, 0, pSpec, bufptr ) <= 0 )
2065               *ok = false;
2066       }
2067 private:
2068     Mat &src;
2069     Mat &dst;
2070     double inv_scale_x;
2071     double inv_scale_y;
2072     void *pSpec;
2073     AutoBuffer<uchar>   specBuf;
2074     int mode;
2075     ippiResizeFunc func;
2076     bool *ok;
2077     const IPPresizeInvoker& operator= (const IPPresizeInvoker&);
2078 };
2079 #endif
2080
2081 #ifdef HAVE_OPENCL
2082
2083 static void ocl_computeResizeAreaTabs(int ssize, int dsize, double scale, int * const map_tab,
2084                                           float * const alpha_tab, int * const ofs_tab)
2085 {
2086     int k = 0, dx = 0;
2087     for ( ; dx < dsize; dx++)
2088     {
2089         ofs_tab[dx] = k;
2090
2091         double fsx1 = dx * scale;
2092         double fsx2 = fsx1 + scale;
2093         double cellWidth = std::min(scale, ssize - fsx1);
2094
2095         int sx1 = cvCeil(fsx1), sx2 = cvFloor(fsx2);
2096
2097         sx2 = std::min(sx2, ssize - 1);
2098         sx1 = std::min(sx1, sx2);
2099
2100         if (sx1 - fsx1 > 1e-3)
2101         {
2102             map_tab[k] = sx1 - 1;
2103             alpha_tab[k++] = (float)((sx1 - fsx1) / cellWidth);
2104         }
2105
2106         for (int sx = sx1; sx < sx2; sx++)
2107         {
2108             map_tab[k] = sx;
2109             alpha_tab[k++] = float(1.0 / cellWidth);
2110         }
2111
2112         if (fsx2 - sx2 > 1e-3)
2113         {
2114             map_tab[k] = sx2;
2115             alpha_tab[k++] = (float)(std::min(std::min(fsx2 - sx2, 1.), cellWidth) / cellWidth);
2116         }
2117     }
2118     ofs_tab[dx] = k;
2119 }
2120
2121 static void ocl_computeResizeAreaFastTabs(int * dmap_tab, int * smap_tab, int scale, int dcols, int scol)
2122 {
2123     for (int i = 0; i < dcols; ++i)
2124         dmap_tab[i] = scale * i;
2125
2126     for (int i = 0, size = dcols * scale; i < size; ++i)
2127         smap_tab[i] = std::min(scol - 1, i);
2128 }
2129
2130 static bool ocl_resize( InputArray _src, OutputArray _dst, Size dsize,
2131                         double fx, double fy, int interpolation)
2132 {
2133     int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
2134
2135     double inv_fx = 1. / fx, inv_fy = 1. / fy;
2136     float inv_fxf = (float)inv_fx, inv_fyf = (float)inv_fy;
2137
2138     if( !(cn <= 4 &&
2139            (interpolation == INTER_NEAREST || interpolation == INTER_LINEAR ||
2140             (interpolation == INTER_AREA && inv_fx >= 1 && inv_fy >= 1) )) )
2141         return false;
2142
2143     UMat src = _src.getUMat();
2144     _dst.create(dsize, type);
2145     UMat dst = _dst.getUMat();
2146
2147     ocl::Kernel k;
2148     size_t globalsize[] = { dst.cols, dst.rows };
2149
2150     if (interpolation == INTER_LINEAR)
2151     {
2152         int wdepth = std::max(depth, CV_32S);
2153         int wtype = CV_MAKETYPE(wdepth, cn);
2154         char buf[2][32];
2155         k.create("resizeLN", ocl::imgproc::resize_oclsrc,
2156                  format("-D INTER_LINEAR -D depth=%d -D PIXTYPE=%s -D PIXTYPE1=%s "
2157                         "-D WORKTYPE=%s -D convertToWT=%s -D convertToDT=%s -D cn=%d",
2158                         depth, ocl::typeToStr(type), ocl::typeToStr(depth), ocl::typeToStr(wtype),
2159                         ocl::convertTypeStr(depth, wdepth, cn, buf[0]),
2160                         ocl::convertTypeStr(wdepth, depth, cn, buf[1]),
2161                         cn));
2162     }
2163     else if (interpolation == INTER_NEAREST)
2164     {
2165         k.create("resizeNN", ocl::imgproc::resize_oclsrc,
2166                  format("-D INTER_NEAREST -D PIXTYPE=%s -D PIXTYPE1=%s -D cn=%d",
2167                         ocl::memopTypeToStr(type), ocl::memopTypeToStr(depth), cn));
2168     }
2169     else if (interpolation == INTER_AREA)
2170     {
2171         int iscale_x = saturate_cast<int>(inv_fx);
2172         int iscale_y = saturate_cast<int>(inv_fy);
2173         bool is_area_fast = std::abs(inv_fx - iscale_x) < DBL_EPSILON &&
2174                         std::abs(inv_fy - iscale_y) < DBL_EPSILON;
2175         int wdepth = std::max(depth, is_area_fast ? CV_32S : CV_32F);
2176         int wtype = CV_MAKE_TYPE(wdepth, cn);
2177
2178         char cvt[2][40];
2179         String buildOption = format("-D INTER_AREA -D PIXTYPE=%s -D PIXTYPE1=%s -D WTV=%s -D convertToWTV=%s -D cn=%d",
2180                                     ocl::typeToStr(type), ocl::typeToStr(depth), ocl::typeToStr(wtype),
2181                                     ocl::convertTypeStr(depth, wdepth, cn, cvt[0]), cn);
2182
2183         UMat alphaOcl, tabofsOcl, mapOcl;
2184         UMat dmap, smap;
2185
2186         if (is_area_fast)
2187         {
2188             int wdepth2 = std::max(CV_32F, depth), wtype2 = CV_MAKE_TYPE(wdepth2, cn);
2189             buildOption = buildOption + format(" -D convertToPIXTYPE=%s -D WT2V=%s -D convertToWT2V=%s -D INTER_AREA_FAST"
2190                                                " -D XSCALE=%d -D YSCALE=%d -D SCALE=%ff",
2191                                                ocl::convertTypeStr(wdepth2, depth, cn, cvt[0]),
2192                                                ocl::typeToStr(wtype2), ocl::convertTypeStr(wdepth, wdepth2, cn, cvt[1]),
2193                                   iscale_x, iscale_y, 1.0f / (iscale_x * iscale_y));
2194
2195             k.create("resizeAREA_FAST", ocl::imgproc::resize_oclsrc, buildOption);
2196             if (k.empty())
2197                 return false;
2198
2199             int smap_tab_size = dst.cols * iscale_x + dst.rows * iscale_y;
2200             AutoBuffer<int> dmap_tab(dst.cols + dst.rows), smap_tab(smap_tab_size);
2201             int * dxmap_tab = dmap_tab, * dymap_tab = dxmap_tab + dst.cols;
2202             int * sxmap_tab = smap_tab, * symap_tab = smap_tab + dst.cols * iscale_y;
2203
2204             ocl_computeResizeAreaFastTabs(dxmap_tab, sxmap_tab, iscale_x, dst.cols, src.cols);
2205             ocl_computeResizeAreaFastTabs(dymap_tab, symap_tab, iscale_y, dst.rows, src.rows);
2206
2207             Mat(1, dst.cols + dst.rows, CV_32SC1, (void *)dmap_tab).copyTo(dmap);
2208             Mat(1, smap_tab_size, CV_32SC1, (void *)smap_tab).copyTo(smap);
2209         }
2210         else
2211         {
2212             buildOption = buildOption + format(" -D convertToPIXTYPE=%s", ocl::convertTypeStr(wdepth, depth, cn, cvt[0]));
2213             k.create("resizeAREA", ocl::imgproc::resize_oclsrc, buildOption);
2214             if (k.empty())
2215                 return false;
2216
2217             Size ssize = src.size();
2218             int xytab_size = (ssize.width + ssize.height) << 1;
2219             int tabofs_size = dsize.height + dsize.width + 2;
2220
2221             AutoBuffer<int> _xymap_tab(xytab_size), _xyofs_tab(tabofs_size);
2222             AutoBuffer<float> _xyalpha_tab(xytab_size);
2223             int * xmap_tab = _xymap_tab, * ymap_tab = _xymap_tab + (ssize.width << 1);
2224             float * xalpha_tab = _xyalpha_tab, * yalpha_tab = _xyalpha_tab + (ssize.width << 1);
2225             int * xofs_tab = _xyofs_tab, * yofs_tab = _xyofs_tab + dsize.width + 1;
2226
2227             ocl_computeResizeAreaTabs(ssize.width, dsize.width, inv_fx, xmap_tab, xalpha_tab, xofs_tab);
2228             ocl_computeResizeAreaTabs(ssize.height, dsize.height, inv_fy, ymap_tab, yalpha_tab, yofs_tab);
2229
2230             // loading precomputed arrays to GPU
2231             Mat(1, xytab_size, CV_32FC1, (void *)_xyalpha_tab).copyTo(alphaOcl);
2232             Mat(1, xytab_size, CV_32SC1, (void *)_xymap_tab).copyTo(mapOcl);
2233             Mat(1, tabofs_size, CV_32SC1, (void *)_xyofs_tab).copyTo(tabofsOcl);
2234         }
2235
2236         ocl::KernelArg srcarg = ocl::KernelArg::ReadOnly(src), dstarg = ocl::KernelArg::WriteOnly(dst);
2237
2238         if (is_area_fast)
2239             k.args(srcarg, dstarg, ocl::KernelArg::PtrReadOnly(dmap), ocl::KernelArg::PtrReadOnly(smap));
2240         else
2241             k.args(srcarg, dstarg, inv_fxf, inv_fyf, ocl::KernelArg::PtrReadOnly(tabofsOcl),
2242                    ocl::KernelArg::PtrReadOnly(mapOcl), ocl::KernelArg::PtrReadOnly(alphaOcl));
2243
2244         return k.run(2, globalsize, NULL, false);
2245     }
2246
2247     if( k.empty() )
2248         return false;
2249     k.args(ocl::KernelArg::ReadOnly(src), ocl::KernelArg::WriteOnly(dst),
2250            (float)inv_fx, (float)inv_fy);
2251
2252     return k.run(2, globalsize, 0, false);
2253 }
2254
2255 #endif
2256
2257 }
2258
2259 //////////////////////////////////////////////////////////////////////////////////////////
2260
2261 void cv::resize( InputArray _src, OutputArray _dst, Size dsize,
2262                  double inv_scale_x, double inv_scale_y, int interpolation )
2263 {
2264     static ResizeFunc linear_tab[] =
2265     {
2266         resizeGeneric_<
2267             HResizeLinear<uchar, int, short,
2268                 INTER_RESIZE_COEF_SCALE,
2269                 HResizeLinearVec_8u32s>,
2270             VResizeLinear<uchar, int, short,
2271                 FixedPtCast<int, uchar, INTER_RESIZE_COEF_BITS*2>,
2272                 VResizeLinearVec_32s8u> >,
2273         0,
2274         resizeGeneric_<
2275             HResizeLinear<ushort, float, float, 1,
2276                 HResizeLinearVec_16u32f>,
2277             VResizeLinear<ushort, float, float, Cast<float, ushort>,
2278                 VResizeLinearVec_32f16u> >,
2279         resizeGeneric_<
2280             HResizeLinear<short, float, float, 1,
2281                 HResizeLinearVec_16s32f>,
2282             VResizeLinear<short, float, float, Cast<float, short>,
2283                 VResizeLinearVec_32f16s> >,
2284         0,
2285         resizeGeneric_<
2286             HResizeLinear<float, float, float, 1,
2287                 HResizeLinearVec_32f>,
2288             VResizeLinear<float, float, float, Cast<float, float>,
2289                 VResizeLinearVec_32f> >,
2290         resizeGeneric_<
2291             HResizeLinear<double, double, float, 1,
2292                 HResizeNoVec>,
2293             VResizeLinear<double, double, float, Cast<double, double>,
2294                 VResizeNoVec> >,
2295         0
2296     };
2297
2298     static ResizeFunc cubic_tab[] =
2299     {
2300         resizeGeneric_<
2301             HResizeCubic<uchar, int, short>,
2302             VResizeCubic<uchar, int, short,
2303                 FixedPtCast<int, uchar, INTER_RESIZE_COEF_BITS*2>,
2304                 VResizeCubicVec_32s8u> >,
2305         0,
2306         resizeGeneric_<
2307             HResizeCubic<ushort, float, float>,
2308             VResizeCubic<ushort, float, float, Cast<float, ushort>,
2309             VResizeCubicVec_32f16u> >,
2310         resizeGeneric_<
2311             HResizeCubic<short, float, float>,
2312             VResizeCubic<short, float, float, Cast<float, short>,
2313             VResizeCubicVec_32f16s> >,
2314         0,
2315         resizeGeneric_<
2316             HResizeCubic<float, float, float>,
2317             VResizeCubic<float, float, float, Cast<float, float>,
2318             VResizeCubicVec_32f> >,
2319         resizeGeneric_<
2320             HResizeCubic<double, double, float>,
2321             VResizeCubic<double, double, float, Cast<double, double>,
2322             VResizeNoVec> >,
2323         0
2324     };
2325
2326     static ResizeFunc lanczos4_tab[] =
2327     {
2328         resizeGeneric_<HResizeLanczos4<uchar, int, short>,
2329             VResizeLanczos4<uchar, int, short,
2330             FixedPtCast<int, uchar, INTER_RESIZE_COEF_BITS*2>,
2331             VResizeNoVec> >,
2332         0,
2333         resizeGeneric_<HResizeLanczos4<ushort, float, float>,
2334             VResizeLanczos4<ushort, float, float, Cast<float, ushort>,
2335             VResizeNoVec> >,
2336         resizeGeneric_<HResizeLanczos4<short, float, float>,
2337             VResizeLanczos4<short, float, float, Cast<float, short>,
2338             VResizeNoVec> >,
2339         0,
2340         resizeGeneric_<HResizeLanczos4<float, float, float>,
2341             VResizeLanczos4<float, float, float, Cast<float, float>,
2342             VResizeNoVec> >,
2343         resizeGeneric_<HResizeLanczos4<double, double, float>,
2344             VResizeLanczos4<double, double, float, Cast<double, double>,
2345             VResizeNoVec> >,
2346         0
2347     };
2348
2349     static ResizeAreaFastFunc areafast_tab[] =
2350     {
2351         resizeAreaFast_<uchar, int, ResizeAreaFastVec<uchar, ResizeAreaFastVec_SIMD_8u> >,
2352         0,
2353         resizeAreaFast_<ushort, float, ResizeAreaFastVec<ushort, ResizeAreaFastVec_SIMD_16u> >,
2354         resizeAreaFast_<short, float, ResizeAreaFastVec<short, ResizeAreaFastNoVec<short, float> > >,
2355         0,
2356         resizeAreaFast_<float, float, ResizeAreaFastNoVec<float, float> >,
2357         resizeAreaFast_<double, double, ResizeAreaFastNoVec<double, double> >,
2358         0
2359     };
2360
2361     static ResizeAreaFunc area_tab[] =
2362     {
2363         resizeArea_<uchar, float>, 0, resizeArea_<ushort, float>,
2364         resizeArea_<short, float>, 0, resizeArea_<float, float>,
2365         resizeArea_<double, double>, 0
2366     };
2367
2368     Size ssize = _src.size();
2369
2370     CV_Assert( ssize.area() > 0 );
2371     CV_Assert( dsize.area() > 0 || (inv_scale_x > 0 && inv_scale_y > 0) );
2372     if( dsize.area() == 0 )
2373     {
2374         dsize = Size(saturate_cast<int>(ssize.width*inv_scale_x),
2375                      saturate_cast<int>(ssize.height*inv_scale_y));
2376         CV_Assert( dsize.area() > 0 );
2377     }
2378     else
2379     {
2380         inv_scale_x = (double)dsize.width/ssize.width;
2381         inv_scale_y = (double)dsize.height/ssize.height;
2382     }
2383
2384     CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
2385                ocl_resize(_src, _dst, dsize, inv_scale_x, inv_scale_y, interpolation))
2386
2387     Mat src = _src.getMat();
2388     _dst.create(dsize, src.type());
2389     Mat dst = _dst.getMat();
2390
2391 #ifdef HAVE_TEGRA_OPTIMIZATION
2392     if (tegra::resize(src, dst, (float)inv_scale_x, (float)inv_scale_y, interpolation))
2393         return;
2394 #endif
2395
2396     int depth = src.depth(), cn = src.channels();
2397     double scale_x = 1./inv_scale_x, scale_y = 1./inv_scale_y;
2398     int k, sx, sy, dx, dy;
2399
2400 #if defined (HAVE_IPP) && ((IPP_VERSION_MAJOR == 7 && IPP_VERSION_MINOR >= 1) || IPP_VERSION_MAJOR > 7)
2401 #define IPP_RESIZE_EPS    1.e-10
2402
2403     double ex = fabs((double)dsize.width/src.cols  - inv_scale_x)/inv_scale_x;
2404     double ey = fabs((double)dsize.height/src.rows - inv_scale_y)/inv_scale_y;
2405
2406     if ((ex < IPP_RESIZE_EPS && ey < IPP_RESIZE_EPS && depth != CV_64F) ||
2407         (ex == 0 && ey == 0 && depth == CV_64F))
2408     {
2409         int mode = 0;
2410         if (interpolation == INTER_LINEAR && src.rows >= 2 && src.cols >= 2)
2411         {
2412             mode = ippLinear;
2413         }
2414         else if (interpolation == INTER_CUBIC && src.rows >= 4 && src.cols >= 4)
2415         {
2416             mode = ippCubic;
2417         }
2418         if( mode != 0 && (cn == 1 || cn ==3 || cn == 4) &&
2419             (depth == CV_8U || depth == CV_16U || depth == CV_16S || depth == CV_32F ||
2420             (depth == CV_64F && mode == ippLinear)))
2421         {
2422             bool ok = true;
2423             Range range(0, src.rows);
2424             IPPresizeInvoker invoker(src, dst, inv_scale_x, inv_scale_y, mode, &ok);
2425             parallel_for_(range, invoker, dst.total()/(double)(1<<16));
2426             if( ok )
2427                 return;
2428         }
2429     }
2430 #undef IPP_RESIZE_EPS
2431 #endif
2432
2433     if( interpolation == INTER_NEAREST )
2434     {
2435         resizeNN( src, dst, inv_scale_x, inv_scale_y );
2436         return;
2437     }
2438
2439     {
2440         int iscale_x = saturate_cast<int>(scale_x);
2441         int iscale_y = saturate_cast<int>(scale_y);
2442
2443         bool is_area_fast = std::abs(scale_x - iscale_x) < DBL_EPSILON &&
2444                 std::abs(scale_y - iscale_y) < DBL_EPSILON;
2445
2446         // in case of scale_x && scale_y is equal to 2
2447         // INTER_AREA (fast) also is equal to INTER_LINEAR
2448         if( interpolation == INTER_LINEAR && is_area_fast && iscale_x == 2 && iscale_y == 2 )
2449             interpolation = INTER_AREA;
2450
2451         // true "area" interpolation is only implemented for the case (scale_x <= 1 && scale_y <= 1).
2452         // In other cases it is emulated using some variant of bilinear interpolation
2453         if( interpolation == INTER_AREA && scale_x >= 1 && scale_y >= 1 )
2454         {
2455             if( is_area_fast )
2456             {
2457                 int area = iscale_x*iscale_y;
2458                 size_t srcstep = src.step / src.elemSize1();
2459                 AutoBuffer<int> _ofs(area + dsize.width*cn);
2460                 int* ofs = _ofs;
2461                 int* xofs = ofs + area;
2462                 ResizeAreaFastFunc func = areafast_tab[depth];
2463                 CV_Assert( func != 0 );
2464
2465                 for( sy = 0, k = 0; sy < iscale_y; sy++ )
2466                     for( sx = 0; sx < iscale_x; sx++ )
2467                         ofs[k++] = (int)(sy*srcstep + sx*cn);
2468
2469                 for( dx = 0; dx < dsize.width; dx++ )
2470                 {
2471                     int j = dx * cn;
2472                     sx = iscale_x * j;
2473                     for( k = 0; k < cn; k++ )
2474                         xofs[j + k] = sx + k;
2475                 }
2476
2477                 func( src, dst, ofs, xofs, iscale_x, iscale_y );
2478                 return;
2479             }
2480
2481             ResizeAreaFunc func = area_tab[depth];
2482             CV_Assert( func != 0 && cn <= 4 );
2483
2484             AutoBuffer<DecimateAlpha> _xytab((ssize.width + ssize.height)*2);
2485             DecimateAlpha* xtab = _xytab, *ytab = xtab + ssize.width*2;
2486
2487             int xtab_size = computeResizeAreaTab(ssize.width, dsize.width, cn, scale_x, xtab);
2488             int ytab_size = computeResizeAreaTab(ssize.height, dsize.height, 1, scale_y, ytab);
2489
2490             AutoBuffer<int> _tabofs(dsize.height + 1);
2491             int* tabofs = _tabofs;
2492             for( k = 0, dy = 0; k < ytab_size; k++ )
2493             {
2494                 if( k == 0 || ytab[k].di != ytab[k-1].di )
2495                 {
2496                     assert( ytab[k].di == dy );
2497                     tabofs[dy++] = k;
2498                 }
2499             }
2500             tabofs[dy] = ytab_size;
2501
2502             func( src, dst, xtab, xtab_size, ytab, ytab_size, tabofs );
2503             return;
2504         }
2505     }
2506
2507     int xmin = 0, xmax = dsize.width, width = dsize.width*cn;
2508     bool area_mode = interpolation == INTER_AREA;
2509     bool fixpt = depth == CV_8U;
2510     float fx, fy;
2511     ResizeFunc func=0;
2512     int ksize=0, ksize2;
2513     if( interpolation == INTER_CUBIC )
2514         ksize = 4, func = cubic_tab[depth];
2515     else if( interpolation == INTER_LANCZOS4 )
2516         ksize = 8, func = lanczos4_tab[depth];
2517     else if( interpolation == INTER_LINEAR || interpolation == INTER_AREA )
2518         ksize = 2, func = linear_tab[depth];
2519     else
2520         CV_Error( CV_StsBadArg, "Unknown interpolation method" );
2521     ksize2 = ksize/2;
2522
2523     CV_Assert( func != 0 );
2524
2525     AutoBuffer<uchar> _buffer((width + dsize.height)*(sizeof(int) + sizeof(float)*ksize));
2526     int* xofs = (int*)(uchar*)_buffer;
2527     int* yofs = xofs + width;
2528     float* alpha = (float*)(yofs + dsize.height);
2529     short* ialpha = (short*)alpha;
2530     float* beta = alpha + width*ksize;
2531     short* ibeta = ialpha + width*ksize;
2532     float cbuf[MAX_ESIZE];
2533
2534     for( dx = 0; dx < dsize.width; dx++ )
2535     {
2536         if( !area_mode )
2537         {
2538             fx = (float)((dx+0.5)*scale_x - 0.5);
2539             sx = cvFloor(fx);
2540             fx -= sx;
2541         }
2542         else
2543         {
2544             sx = cvFloor(dx*scale_x);
2545             fx = (float)((dx+1) - (sx+1)*inv_scale_x);
2546             fx = fx <= 0 ? 0.f : fx - cvFloor(fx);
2547         }
2548
2549         if( sx < ksize2-1 )
2550         {
2551             xmin = dx+1;
2552             if( sx < 0 && (interpolation != INTER_CUBIC && interpolation != INTER_LANCZOS4))
2553                 fx = 0, sx = 0;
2554         }
2555
2556         if( sx + ksize2 >= ssize.width )
2557         {
2558             xmax = std::min( xmax, dx );
2559             if( sx >= ssize.width-1 && (interpolation != INTER_CUBIC && interpolation != INTER_LANCZOS4))
2560                 fx = 0, sx = ssize.width-1;
2561         }
2562
2563         for( k = 0, sx *= cn; k < cn; k++ )
2564             xofs[dx*cn + k] = sx + k;
2565
2566         if( interpolation == INTER_CUBIC )
2567             interpolateCubic( fx, cbuf );
2568         else if( interpolation == INTER_LANCZOS4 )
2569             interpolateLanczos4( fx, cbuf );
2570         else
2571         {
2572             cbuf[0] = 1.f - fx;
2573             cbuf[1] = fx;
2574         }
2575         if( fixpt )
2576         {
2577             for( k = 0; k < ksize; k++ )
2578                 ialpha[dx*cn*ksize + k] = saturate_cast<short>(cbuf[k]*INTER_RESIZE_COEF_SCALE);
2579             for( ; k < cn*ksize; k++ )
2580                 ialpha[dx*cn*ksize + k] = ialpha[dx*cn*ksize + k - ksize];
2581         }
2582         else
2583         {
2584             for( k = 0; k < ksize; k++ )
2585                 alpha[dx*cn*ksize + k] = cbuf[k];
2586             for( ; k < cn*ksize; k++ )
2587                 alpha[dx*cn*ksize + k] = alpha[dx*cn*ksize + k - ksize];
2588         }
2589     }
2590
2591     for( dy = 0; dy < dsize.height; dy++ )
2592     {
2593         if( !area_mode )
2594         {
2595             fy = (float)((dy+0.5)*scale_y - 0.5);
2596             sy = cvFloor(fy);
2597             fy -= sy;
2598         }
2599         else
2600         {
2601             sy = cvFloor(dy*scale_y);
2602             fy = (float)((dy+1) - (sy+1)*inv_scale_y);
2603             fy = fy <= 0 ? 0.f : fy - cvFloor(fy);
2604         }
2605
2606         yofs[dy] = sy;
2607         if( interpolation == INTER_CUBIC )
2608             interpolateCubic( fy, cbuf );
2609         else if( interpolation == INTER_LANCZOS4 )
2610             interpolateLanczos4( fy, cbuf );
2611         else
2612         {
2613             cbuf[0] = 1.f - fy;
2614             cbuf[1] = fy;
2615         }
2616
2617         if( fixpt )
2618         {
2619             for( k = 0; k < ksize; k++ )
2620                 ibeta[dy*ksize + k] = saturate_cast<short>(cbuf[k]*INTER_RESIZE_COEF_SCALE);
2621         }
2622         else
2623         {
2624             for( k = 0; k < ksize; k++ )
2625                 beta[dy*ksize + k] = cbuf[k];
2626         }
2627     }
2628
2629     func( src, dst, xofs, fixpt ? (void*)ialpha : (void*)alpha, yofs,
2630           fixpt ? (void*)ibeta : (void*)beta, xmin, xmax, ksize );
2631 }
2632
2633
2634 /****************************************************************************************\
2635 *                       General warping (affine, perspective, remap)                     *
2636 \****************************************************************************************/
2637
2638 namespace cv
2639 {
2640
2641 template<typename T>
2642 static void remapNearest( const Mat& _src, Mat& _dst, const Mat& _xy,
2643                           int borderType, const Scalar& _borderValue )
2644 {
2645     Size ssize = _src.size(), dsize = _dst.size();
2646     int cn = _src.channels();
2647     const T* S0 = (const T*)_src.data;
2648     size_t sstep = _src.step/sizeof(S0[0]);
2649     Scalar_<T> cval(saturate_cast<T>(_borderValue[0]),
2650         saturate_cast<T>(_borderValue[1]),
2651         saturate_cast<T>(_borderValue[2]),
2652         saturate_cast<T>(_borderValue[3]));
2653     int dx, dy;
2654
2655     unsigned width1 = ssize.width, height1 = ssize.height;
2656
2657     if( _dst.isContinuous() && _xy.isContinuous() )
2658     {
2659         dsize.width *= dsize.height;
2660         dsize.height = 1;
2661     }
2662
2663     for( dy = 0; dy < dsize.height; dy++ )
2664     {
2665         T* D = (T*)(_dst.data + _dst.step*dy);
2666         const short* XY = (const short*)(_xy.data + _xy.step*dy);
2667
2668         if( cn == 1 )
2669         {
2670             for( dx = 0; dx < dsize.width; dx++ )
2671             {
2672                 int sx = XY[dx*2], sy = XY[dx*2+1];
2673                 if( (unsigned)sx < width1 && (unsigned)sy < height1 )
2674                     D[dx] = S0[sy*sstep + sx];
2675                 else
2676                 {
2677                     if( borderType == BORDER_REPLICATE )
2678                     {
2679                         sx = clip(sx, 0, ssize.width);
2680                         sy = clip(sy, 0, ssize.height);
2681                         D[dx] = S0[sy*sstep + sx];
2682                     }
2683                     else if( borderType == BORDER_CONSTANT )
2684                         D[dx] = cval[0];
2685                     else if( borderType != BORDER_TRANSPARENT )
2686                     {
2687                         sx = borderInterpolate(sx, ssize.width, borderType);
2688                         sy = borderInterpolate(sy, ssize.height, borderType);
2689                         D[dx] = S0[sy*sstep + sx];
2690                     }
2691                 }
2692             }
2693         }
2694         else
2695         {
2696             for( dx = 0; dx < dsize.width; dx++, D += cn )
2697             {
2698                 int sx = XY[dx*2], sy = XY[dx*2+1], k;
2699                 const T *S;
2700                 if( (unsigned)sx < width1 && (unsigned)sy < height1 )
2701                 {
2702                     if( cn == 3 )
2703                     {
2704                         S = S0 + sy*sstep + sx*3;
2705                         D[0] = S[0], D[1] = S[1], D[2] = S[2];
2706                     }
2707                     else if( cn == 4 )
2708                     {
2709                         S = S0 + sy*sstep + sx*4;
2710                         D[0] = S[0], D[1] = S[1], D[2] = S[2], D[3] = S[3];
2711                     }
2712                     else
2713                     {
2714                         S = S0 + sy*sstep + sx*cn;
2715                         for( k = 0; k < cn; k++ )
2716                             D[k] = S[k];
2717                     }
2718                 }
2719                 else if( borderType != BORDER_TRANSPARENT )
2720                 {
2721                     if( borderType == BORDER_REPLICATE )
2722                     {
2723                         sx = clip(sx, 0, ssize.width);
2724                         sy = clip(sy, 0, ssize.height);
2725                         S = S0 + sy*sstep + sx*cn;
2726                     }
2727                     else if( borderType == BORDER_CONSTANT )
2728                         S = &cval[0];
2729                     else
2730                     {
2731                         sx = borderInterpolate(sx, ssize.width, borderType);
2732                         sy = borderInterpolate(sy, ssize.height, borderType);
2733                         S = S0 + sy*sstep + sx*cn;
2734                     }
2735                     for( k = 0; k < cn; k++ )
2736                         D[k] = S[k];
2737                 }
2738             }
2739         }
2740     }
2741 }
2742
2743
2744 struct RemapNoVec
2745 {
2746     int operator()( const Mat&, void*, const short*, const ushort*,
2747                     const void*, int ) const { return 0; }
2748 };
2749
2750 #if CV_SSE2
2751
2752 struct RemapVec_8u
2753 {
2754     int operator()( const Mat& _src, void* _dst, const short* XY,
2755                     const ushort* FXY, const void* _wtab, int width ) const
2756     {
2757         int cn = _src.channels(), x = 0, sstep = (int)_src.step;
2758
2759         if( (cn != 1 && cn != 3 && cn != 4) || !checkHardwareSupport(CV_CPU_SSE2)||
2760             sstep > 0x8000 )
2761             return 0;
2762
2763         const uchar *S0 = _src.data, *S1 = _src.data + _src.step;
2764         const short* wtab = cn == 1 ? (const short*)_wtab : &BilinearTab_iC4[0][0][0];
2765         uchar* D = (uchar*)_dst;
2766         __m128i delta = _mm_set1_epi32(INTER_REMAP_COEF_SCALE/2);
2767         __m128i xy2ofs = _mm_set1_epi32(cn + (sstep << 16));
2768         __m128i z = _mm_setzero_si128();
2769         int CV_DECL_ALIGNED(16) iofs0[4], iofs1[4];
2770
2771         if( cn == 1 )
2772         {
2773             for( ; x <= width - 8; x += 8 )
2774             {
2775                 __m128i xy0 = _mm_loadu_si128( (const __m128i*)(XY + x*2));
2776                 __m128i xy1 = _mm_loadu_si128( (const __m128i*)(XY + x*2 + 8));
2777                 __m128i v0, v1, v2, v3, a0, a1, b0, b1;
2778                 unsigned i0, i1;
2779
2780                 xy0 = _mm_madd_epi16( xy0, xy2ofs );
2781                 xy1 = _mm_madd_epi16( xy1, xy2ofs );
2782                 _mm_store_si128( (__m128i*)iofs0, xy0 );
2783                 _mm_store_si128( (__m128i*)iofs1, xy1 );
2784
2785                 i0 = *(ushort*)(S0 + iofs0[0]) + (*(ushort*)(S0 + iofs0[1]) << 16);
2786                 i1 = *(ushort*)(S0 + iofs0[2]) + (*(ushort*)(S0 + iofs0[3]) << 16);
2787                 v0 = _mm_unpacklo_epi32(_mm_cvtsi32_si128(i0), _mm_cvtsi32_si128(i1));
2788                 i0 = *(ushort*)(S1 + iofs0[0]) + (*(ushort*)(S1 + iofs0[1]) << 16);
2789                 i1 = *(ushort*)(S1 + iofs0[2]) + (*(ushort*)(S1 + iofs0[3]) << 16);
2790                 v1 = _mm_unpacklo_epi32(_mm_cvtsi32_si128(i0), _mm_cvtsi32_si128(i1));
2791                 v0 = _mm_unpacklo_epi8(v0, z);
2792                 v1 = _mm_unpacklo_epi8(v1, z);
2793
2794                 a0 = _mm_unpacklo_epi32(_mm_loadl_epi64((__m128i*)(wtab+FXY[x]*4)),
2795                                         _mm_loadl_epi64((__m128i*)(wtab+FXY[x+1]*4)));
2796                 a1 = _mm_unpacklo_epi32(_mm_loadl_epi64((__m128i*)(wtab+FXY[x+2]*4)),
2797                                         _mm_loadl_epi64((__m128i*)(wtab+FXY[x+3]*4)));
2798                 b0 = _mm_unpacklo_epi64(a0, a1);
2799                 b1 = _mm_unpackhi_epi64(a0, a1);
2800                 v0 = _mm_madd_epi16(v0, b0);
2801                 v1 = _mm_madd_epi16(v1, b1);
2802                 v0 = _mm_add_epi32(_mm_add_epi32(v0, v1), delta);
2803
2804                 i0 = *(ushort*)(S0 + iofs1[0]) + (*(ushort*)(S0 + iofs1[1]) << 16);
2805                 i1 = *(ushort*)(S0 + iofs1[2]) + (*(ushort*)(S0 + iofs1[3]) << 16);
2806                 v2 = _mm_unpacklo_epi32(_mm_cvtsi32_si128(i0), _mm_cvtsi32_si128(i1));
2807                 i0 = *(ushort*)(S1 + iofs1[0]) + (*(ushort*)(S1 + iofs1[1]) << 16);
2808                 i1 = *(ushort*)(S1 + iofs1[2]) + (*(ushort*)(S1 + iofs1[3]) << 16);
2809                 v3 = _mm_unpacklo_epi32(_mm_cvtsi32_si128(i0), _mm_cvtsi32_si128(i1));
2810                 v2 = _mm_unpacklo_epi8(v2, z);
2811                 v3 = _mm_unpacklo_epi8(v3, z);
2812
2813                 a0 = _mm_unpacklo_epi32(_mm_loadl_epi64((__m128i*)(wtab+FXY[x+4]*4)),
2814                                         _mm_loadl_epi64((__m128i*)(wtab+FXY[x+5]*4)));
2815                 a1 = _mm_unpacklo_epi32(_mm_loadl_epi64((__m128i*)(wtab+FXY[x+6]*4)),
2816                                         _mm_loadl_epi64((__m128i*)(wtab+FXY[x+7]*4)));
2817                 b0 = _mm_unpacklo_epi64(a0, a1);
2818                 b1 = _mm_unpackhi_epi64(a0, a1);
2819                 v2 = _mm_madd_epi16(v2, b0);
2820                 v3 = _mm_madd_epi16(v3, b1);
2821                 v2 = _mm_add_epi32(_mm_add_epi32(v2, v3), delta);
2822
2823                 v0 = _mm_srai_epi32(v0, INTER_REMAP_COEF_BITS);
2824                 v2 = _mm_srai_epi32(v2, INTER_REMAP_COEF_BITS);
2825                 v0 = _mm_packus_epi16(_mm_packs_epi32(v0, v2), z);
2826                 _mm_storel_epi64( (__m128i*)(D + x), v0 );
2827             }
2828         }
2829         else if( cn == 3 )
2830         {
2831             for( ; x <= width - 5; x += 4, D += 12 )
2832             {
2833                 __m128i xy0 = _mm_loadu_si128( (const __m128i*)(XY + x*2));
2834                 __m128i u0, v0, u1, v1;
2835
2836                 xy0 = _mm_madd_epi16( xy0, xy2ofs );
2837                 _mm_store_si128( (__m128i*)iofs0, xy0 );
2838                 const __m128i *w0, *w1;
2839                 w0 = (const __m128i*)(wtab + FXY[x]*16);
2840                 w1 = (const __m128i*)(wtab + FXY[x+1]*16);
2841
2842                 u0 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S0 + iofs0[0])),
2843                                        _mm_cvtsi32_si128(*(int*)(S0 + iofs0[0] + 3)));
2844                 v0 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S1 + iofs0[0])),
2845                                        _mm_cvtsi32_si128(*(int*)(S1 + iofs0[0] + 3)));
2846                 u1 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S0 + iofs0[1])),
2847                                        _mm_cvtsi32_si128(*(int*)(S0 + iofs0[1] + 3)));
2848                 v1 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S1 + iofs0[1])),
2849                                        _mm_cvtsi32_si128(*(int*)(S1 + iofs0[1] + 3)));
2850                 u0 = _mm_unpacklo_epi8(u0, z);
2851                 v0 = _mm_unpacklo_epi8(v0, z);
2852                 u1 = _mm_unpacklo_epi8(u1, z);
2853                 v1 = _mm_unpacklo_epi8(v1, z);
2854                 u0 = _mm_add_epi32(_mm_madd_epi16(u0, w0[0]), _mm_madd_epi16(v0, w0[1]));
2855                 u1 = _mm_add_epi32(_mm_madd_epi16(u1, w1[0]), _mm_madd_epi16(v1, w1[1]));
2856                 u0 = _mm_srai_epi32(_mm_add_epi32(u0, delta), INTER_REMAP_COEF_BITS);
2857                 u1 = _mm_srai_epi32(_mm_add_epi32(u1, delta), INTER_REMAP_COEF_BITS);
2858                 u0 = _mm_slli_si128(u0, 4);
2859                 u0 = _mm_packs_epi32(u0, u1);
2860                 u0 = _mm_packus_epi16(u0, u0);
2861                 _mm_storel_epi64((__m128i*)D, _mm_srli_si128(u0,1));
2862
2863                 w0 = (const __m128i*)(wtab + FXY[x+2]*16);
2864                 w1 = (const __m128i*)(wtab + FXY[x+3]*16);
2865
2866                 u0 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S0 + iofs0[2])),
2867                                        _mm_cvtsi32_si128(*(int*)(S0 + iofs0[2] + 3)));
2868                 v0 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S1 + iofs0[2])),
2869                                        _mm_cvtsi32_si128(*(int*)(S1 + iofs0[2] + 3)));
2870                 u1 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S0 + iofs0[3])),
2871                                        _mm_cvtsi32_si128(*(int*)(S0 + iofs0[3] + 3)));
2872                 v1 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S1 + iofs0[3])),
2873                                        _mm_cvtsi32_si128(*(int*)(S1 + iofs0[3] + 3)));
2874                 u0 = _mm_unpacklo_epi8(u0, z);
2875                 v0 = _mm_unpacklo_epi8(v0, z);
2876                 u1 = _mm_unpacklo_epi8(u1, z);
2877                 v1 = _mm_unpacklo_epi8(v1, z);
2878                 u0 = _mm_add_epi32(_mm_madd_epi16(u0, w0[0]), _mm_madd_epi16(v0, w0[1]));
2879                 u1 = _mm_add_epi32(_mm_madd_epi16(u1, w1[0]), _mm_madd_epi16(v1, w1[1]));
2880                 u0 = _mm_srai_epi32(_mm_add_epi32(u0, delta), INTER_REMAP_COEF_BITS);
2881                 u1 = _mm_srai_epi32(_mm_add_epi32(u1, delta), INTER_REMAP_COEF_BITS);
2882                 u0 = _mm_slli_si128(u0, 4);
2883                 u0 = _mm_packs_epi32(u0, u1);
2884                 u0 = _mm_packus_epi16(u0, u0);
2885                 _mm_storel_epi64((__m128i*)(D + 6), _mm_srli_si128(u0,1));
2886             }
2887         }
2888         else if( cn == 4 )
2889         {
2890             for( ; x <= width - 4; x += 4, D += 16 )
2891             {
2892                 __m128i xy0 = _mm_loadu_si128( (const __m128i*)(XY + x*2));
2893                 __m128i u0, v0, u1, v1;
2894
2895                 xy0 = _mm_madd_epi16( xy0, xy2ofs );
2896                 _mm_store_si128( (__m128i*)iofs0, xy0 );
2897                 const __m128i *w0, *w1;
2898                 w0 = (const __m128i*)(wtab + FXY[x]*16);
2899                 w1 = (const __m128i*)(wtab + FXY[x+1]*16);
2900
2901                 u0 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S0 + iofs0[0])),
2902                                        _mm_cvtsi32_si128(*(int*)(S0 + iofs0[0] + 4)));
2903                 v0 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S1 + iofs0[0])),
2904                                        _mm_cvtsi32_si128(*(int*)(S1 + iofs0[0] + 4)));
2905                 u1 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S0 + iofs0[1])),
2906                                        _mm_cvtsi32_si128(*(int*)(S0 + iofs0[1] + 4)));
2907                 v1 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S1 + iofs0[1])),
2908                                        _mm_cvtsi32_si128(*(int*)(S1 + iofs0[1] + 4)));
2909                 u0 = _mm_unpacklo_epi8(u0, z);
2910                 v0 = _mm_unpacklo_epi8(v0, z);
2911                 u1 = _mm_unpacklo_epi8(u1, z);
2912                 v1 = _mm_unpacklo_epi8(v1, z);
2913                 u0 = _mm_add_epi32(_mm_madd_epi16(u0, w0[0]), _mm_madd_epi16(v0, w0[1]));
2914                 u1 = _mm_add_epi32(_mm_madd_epi16(u1, w1[0]), _mm_madd_epi16(v1, w1[1]));
2915                 u0 = _mm_srai_epi32(_mm_add_epi32(u0, delta), INTER_REMAP_COEF_BITS);
2916                 u1 = _mm_srai_epi32(_mm_add_epi32(u1, delta), INTER_REMAP_COEF_BITS);
2917                 u0 = _mm_packs_epi32(u0, u1);
2918                 u0 = _mm_packus_epi16(u0, u0);
2919                 _mm_storel_epi64((__m128i*)D, u0);
2920
2921                 w0 = (const __m128i*)(wtab + FXY[x+2]*16);
2922                 w1 = (const __m128i*)(wtab + FXY[x+3]*16);
2923
2924                 u0 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S0 + iofs0[2])),
2925                                        _mm_cvtsi32_si128(*(int*)(S0 + iofs0[2] + 4)));
2926                 v0 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S1 + iofs0[2])),
2927                                        _mm_cvtsi32_si128(*(int*)(S1 + iofs0[2] + 4)));
2928                 u1 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S0 + iofs0[3])),
2929                                        _mm_cvtsi32_si128(*(int*)(S0 + iofs0[3] + 4)));
2930                 v1 = _mm_unpacklo_epi8(_mm_cvtsi32_si128(*(int*)(S1 + iofs0[3])),
2931                                        _mm_cvtsi32_si128(*(int*)(S1 + iofs0[3] + 4)));
2932                 u0 = _mm_unpacklo_epi8(u0, z);
2933                 v0 = _mm_unpacklo_epi8(v0, z);
2934                 u1 = _mm_unpacklo_epi8(u1, z);
2935                 v1 = _mm_unpacklo_epi8(v1, z);
2936                 u0 = _mm_add_epi32(_mm_madd_epi16(u0, w0[0]), _mm_madd_epi16(v0, w0[1]));
2937                 u1 = _mm_add_epi32(_mm_madd_epi16(u1, w1[0]), _mm_madd_epi16(v1, w1[1]));
2938                 u0 = _mm_srai_epi32(_mm_add_epi32(u0, delta), INTER_REMAP_COEF_BITS);
2939                 u1 = _mm_srai_epi32(_mm_add_epi32(u1, delta), INTER_REMAP_COEF_BITS);
2940                 u0 = _mm_packs_epi32(u0, u1);
2941                 u0 = _mm_packus_epi16(u0, u0);
2942                 _mm_storel_epi64((__m128i*)(D + 8), u0);
2943             }
2944         }
2945
2946         return x;
2947     }
2948 };
2949
2950 #else
2951
2952 typedef RemapNoVec RemapVec_8u;
2953
2954 #endif
2955
2956
2957 template<class CastOp, class VecOp, typename AT>
2958 static void remapBilinear( const Mat& _src, Mat& _dst, const Mat& _xy,
2959                            const Mat& _fxy, const void* _wtab,
2960                            int borderType, const Scalar& _borderValue )
2961 {
2962     typedef typename CastOp::rtype T;
2963     typedef typename CastOp::type1 WT;
2964     Size ssize = _src.size(), dsize = _dst.size();
2965     int cn = _src.channels();
2966     const AT* wtab = (const AT*)_wtab;
2967     const T* S0 = (const T*)_src.data;
2968     size_t sstep = _src.step/sizeof(S0[0]);
2969     Scalar_<T> cval(saturate_cast<T>(_borderValue[0]),
2970         saturate_cast<T>(_borderValue[1]),
2971         saturate_cast<T>(_borderValue[2]),
2972         saturate_cast<T>(_borderValue[3]));
2973     int dx, dy;
2974     CastOp castOp;
2975     VecOp vecOp;
2976
2977     unsigned width1 = std::max(ssize.width-1, 0), height1 = std::max(ssize.height-1, 0);
2978     CV_Assert( cn <= 4 && ssize.area() > 0 );
2979 #if CV_SSE2
2980     if( _src.type() == CV_8UC3 )
2981         width1 = std::max(ssize.width-2, 0);
2982 #endif
2983
2984     for( dy = 0; dy < dsize.height; dy++ )
2985     {
2986         T* D = (T*)(_dst.data + _dst.step*dy);
2987         const short* XY = (const short*)(_xy.data + _xy.step*dy);
2988         const ushort* FXY = (const ushort*)(_fxy.data + _fxy.step*dy);
2989         int X0 = 0;
2990         bool prevInlier = false;
2991
2992         for( dx = 0; dx <= dsize.width; dx++ )
2993         {
2994             bool curInlier = dx < dsize.width ?
2995                 (unsigned)XY[dx*2] < width1 &&
2996                 (unsigned)XY[dx*2+1] < height1 : !prevInlier;
2997             if( curInlier == prevInlier )
2998                 continue;
2999
3000             int X1 = dx;
3001             dx = X0;
3002             X0 = X1;
3003             prevInlier = curInlier;
3004
3005             if( !curInlier )
3006             {
3007                 int len = vecOp( _src, D, XY + dx*2, FXY + dx, wtab, X1 - dx );
3008                 D += len*cn;
3009                 dx += len;
3010
3011                 if( cn == 1 )
3012                 {
3013                     for( ; dx < X1; dx++, D++ )
3014                     {
3015                         int sx = XY[dx*2], sy = XY[dx*2+1];
3016                         const AT* w = wtab + FXY[dx]*4;
3017                         const T* S = S0 + sy*sstep + sx;
3018                         *D = castOp(WT(S[0]*w[0] + S[1]*w[1] + S[sstep]*w[2] + S[sstep+1]*w[3]));
3019                     }
3020                 }
3021                 else if( cn == 2 )
3022                     for( ; dx < X1; dx++, D += 2 )
3023                     {
3024                         int sx = XY[dx*2], sy = XY[dx*2+1];
3025                         const AT* w = wtab + FXY[dx]*4;
3026                         const T* S = S0 + sy*sstep + sx*2;
3027                         WT t0 = S[0]*w[0] + S[2]*w[1] + S[sstep]*w[2] + S[sstep+2]*w[3];
3028                         WT t1 = S[1]*w[0] + S[3]*w[1] + S[sstep+1]*w[2] + S[sstep+3]*w[3];
3029                         D[0] = castOp(t0); D[1] = castOp(t1);
3030                     }
3031                 else if( cn == 3 )
3032                     for( ; dx < X1; dx++, D += 3 )
3033                     {
3034                         int sx = XY[dx*2], sy = XY[dx*2+1];
3035                         const AT* w = wtab + FXY[dx]*4;
3036                         const T* S = S0 + sy*sstep + sx*3;
3037                         WT t0 = S[0]*w[0] + S[3]*w[1] + S[sstep]*w[2] + S[sstep+3]*w[3];
3038                         WT t1 = S[1]*w[0] + S[4]*w[1] + S[sstep+1]*w[2] + S[sstep+4]*w[3];
3039                         WT t2 = S[2]*w[0] + S[5]*w[1] + S[sstep+2]*w[2] + S[sstep+5]*w[3];
3040                         D[0] = castOp(t0); D[1] = castOp(t1); D[2] = castOp(t2);
3041                     }
3042                 else
3043                     for( ; dx < X1; dx++, D += 4 )
3044                     {
3045                         int sx = XY[dx*2], sy = XY[dx*2+1];
3046                         const AT* w = wtab + FXY[dx]*4;
3047                         const T* S = S0 + sy*sstep + sx*4;
3048                         WT t0 = S[0]*w[0] + S[4]*w[1] + S[sstep]*w[2] + S[sstep+4]*w[3];
3049                         WT t1 = S[1]*w[0] + S[5]*w[1] + S[sstep+1]*w[2] + S[sstep+5]*w[3];
3050                         D[0] = castOp(t0); D[1] = castOp(t1);
3051                         t0 = S[2]*w[0] + S[6]*w[1] + S[sstep+2]*w[2] + S[sstep+6]*w[3];
3052                         t1 = S[3]*w[0] + S[7]*w[1] + S[sstep+3]*w[2] + S[sstep+7]*w[3];
3053                         D[2] = castOp(t0); D[3] = castOp(t1);
3054                     }
3055             }
3056             else
3057             {
3058                 if( borderType == BORDER_TRANSPARENT && cn != 3 )
3059                 {
3060                     D += (X1 - dx)*cn;
3061                     dx = X1;
3062                     continue;
3063                 }
3064
3065                 if( cn == 1 )
3066                     for( ; dx < X1; dx++, D++ )
3067                     {
3068                         int sx = XY[dx*2], sy = XY[dx*2+1];
3069                         if( borderType == BORDER_CONSTANT &&
3070                             (sx >= ssize.width || sx+1 < 0 ||
3071                              sy >= ssize.height || sy+1 < 0) )
3072                         {
3073                             D[0] = cval[0];
3074                         }
3075                         else
3076                         {
3077                             int sx0, sx1, sy0, sy1;
3078                             T v0, v1, v2, v3;
3079                             const AT* w = wtab + FXY[dx]*4;
3080                             if( borderType == BORDER_REPLICATE )
3081                             {
3082                                 sx0 = clip(sx, 0, ssize.width);
3083                                 sx1 = clip(sx+1, 0, ssize.width);
3084                                 sy0 = clip(sy, 0, ssize.height);
3085                                 sy1 = clip(sy+1, 0, ssize.height);
3086                                 v0 = S0[sy0*sstep + sx0];
3087                                 v1 = S0[sy0*sstep + sx1];
3088                                 v2 = S0[sy1*sstep + sx0];
3089                                 v3 = S0[sy1*sstep + sx1];
3090                             }
3091                             else
3092                             {
3093                                 sx0 = borderInterpolate(sx, ssize.width, borderType);
3094                                 sx1 = borderInterpolate(sx+1, ssize.width, borderType);
3095                                 sy0 = borderInterpolate(sy, ssize.height, borderType);
3096                                 sy1 = borderInterpolate(sy+1, ssize.height, borderType);
3097                                 v0 = sx0 >= 0 && sy0 >= 0 ? S0[sy0*sstep + sx0] : cval[0];
3098                                 v1 = sx1 >= 0 && sy0 >= 0 ? S0[sy0*sstep + sx1] : cval[0];
3099                                 v2 = sx0 >= 0 && sy1 >= 0 ? S0[sy1*sstep + sx0] : cval[0];
3100                                 v3 = sx1 >= 0 && sy1 >= 0 ? S0[sy1*sstep + sx1] : cval[0];
3101                             }
3102                             D[0] = castOp(WT(v0*w[0] + v1*w[1] + v2*w[2] + v3*w[3]));
3103                         }
3104                     }
3105                 else
3106                     for( ; dx < X1; dx++, D += cn )
3107                     {
3108                         int sx = XY[dx*2], sy = XY[dx*2+1], k;
3109                         if( borderType == BORDER_CONSTANT &&
3110                             (sx >= ssize.width || sx+1 < 0 ||
3111                              sy >= ssize.height || sy+1 < 0) )
3112                         {
3113                             for( k = 0; k < cn; k++ )
3114                                 D[k] = cval[k];
3115                         }
3116                         else
3117                         {
3118                             int sx0, sx1, sy0, sy1;
3119                             const T *v0, *v1, *v2, *v3;
3120                             const AT* w = wtab + FXY[dx]*4;
3121                             if( borderType == BORDER_REPLICATE )
3122                             {
3123                                 sx0 = clip(sx, 0, ssize.width);
3124                                 sx1 = clip(sx+1, 0, ssize.width);
3125                                 sy0 = clip(sy, 0, ssize.height);
3126                                 sy1 = clip(sy+1, 0, ssize.height);
3127                                 v0 = S0 + sy0*sstep + sx0*cn;
3128                                 v1 = S0 + sy0*sstep + sx1*cn;
3129                                 v2 = S0 + sy1*sstep + sx0*cn;
3130                                 v3 = S0 + sy1*sstep + sx1*cn;
3131                             }
3132                             else if( borderType == BORDER_TRANSPARENT &&
3133                                 ((unsigned)sx >= (unsigned)(ssize.width-1) ||
3134                                 (unsigned)sy >= (unsigned)(ssize.height-1)))
3135                                 continue;
3136                             else
3137                             {
3138                                 sx0 = borderInterpolate(sx, ssize.width, borderType);
3139                                 sx1 = borderInterpolate(sx+1, ssize.width, borderType);
3140                                 sy0 = borderInterpolate(sy, ssize.height, borderType);
3141                                 sy1 = borderInterpolate(sy+1, ssize.height, borderType);
3142                                 v0 = sx0 >= 0 && sy0 >= 0 ? S0 + sy0*sstep + sx0*cn : &cval[0];
3143                                 v1 = sx1 >= 0 && sy0 >= 0 ? S0 + sy0*sstep + sx1*cn : &cval[0];
3144                                 v2 = sx0 >= 0 && sy1 >= 0 ? S0 + sy1*sstep + sx0*cn : &cval[0];
3145                                 v3 = sx1 >= 0 && sy1 >= 0 ? S0 + sy1*sstep + sx1*cn : &cval[0];
3146                             }
3147                             for( k = 0; k < cn; k++ )
3148                                 D[k] = castOp(WT(v0[k]*w[0] + v1[k]*w[1] + v2[k]*w[2] + v3[k]*w[3]));
3149                         }
3150                     }
3151             }
3152         }
3153     }
3154 }
3155
3156
3157 template<class CastOp, typename AT, int ONE>
3158 static void remapBicubic( const Mat& _src, Mat& _dst, const Mat& _xy,
3159                           const Mat& _fxy, const void* _wtab,
3160                           int borderType, const Scalar& _borderValue )
3161 {
3162     typedef typename CastOp::rtype T;
3163     typedef typename CastOp::type1 WT;
3164     Size ssize = _src.size(), dsize = _dst.size();
3165     int cn = _src.channels();
3166     const AT* wtab = (const AT*)_wtab;
3167     const T* S0 = (const T*)_src.data;
3168     size_t sstep = _src.step/sizeof(S0[0]);
3169     Scalar_<T> cval(saturate_cast<T>(_borderValue[0]),
3170         saturate_cast<T>(_borderValue[1]),
3171         saturate_cast<T>(_borderValue[2]),
3172         saturate_cast<T>(_borderValue[3]));
3173     int dx, dy;
3174     CastOp castOp;
3175     int borderType1 = borderType != BORDER_TRANSPARENT ? borderType : BORDER_REFLECT_101;
3176
3177     unsigned width1 = std::max(ssize.width-3, 0), height1 = std::max(ssize.height-3, 0);
3178
3179     if( _dst.isContinuous() && _xy.isContinuous() && _fxy.isContinuous() )
3180     {
3181         dsize.width *= dsize.height;
3182         dsize.height = 1;
3183     }
3184
3185     for( dy = 0; dy < dsize.height; dy++ )
3186     {
3187         T* D = (T*)(_dst.data + _dst.step*dy);
3188         const short* XY = (const short*)(_xy.data + _xy.step*dy);
3189         const ushort* FXY = (const ushort*)(_fxy.data + _fxy.step*dy);
3190
3191         for( dx = 0; dx < dsize.width; dx++, D += cn )
3192         {
3193             int sx = XY[dx*2]-1, sy = XY[dx*2+1]-1;
3194             const AT* w = wtab + FXY[dx]*16;
3195             int i, k;
3196             if( (unsigned)sx < width1 && (unsigned)sy < height1 )
3197             {
3198                 const T* S = S0 + sy*sstep + sx*cn;
3199                 for( k = 0; k < cn; k++ )
3200                 {
3201                     WT sum = S[0]*w[0] + S[cn]*w[1] + S[cn*2]*w[2] + S[cn*3]*w[3];
3202                     S += sstep;
3203                     sum += S[0]*w[4] + S[cn]*w[5] + S[cn*2]*w[6] + S[cn*3]*w[7];
3204                     S += sstep;
3205                     sum += S[0]*w[8] + S[cn]*w[9] + S[cn*2]*w[10] + S[cn*3]*w[11];
3206                     S += sstep;
3207                     sum += S[0]*w[12] + S[cn]*w[13] + S[cn*2]*w[14] + S[cn*3]*w[15];
3208                     S += 1 - sstep*3;
3209                     D[k] = castOp(sum);
3210                 }
3211             }
3212             else
3213             {
3214                 int x[4], y[4];
3215                 if( borderType == BORDER_TRANSPARENT &&
3216                     ((unsigned)(sx+1) >= (unsigned)ssize.width ||
3217                     (unsigned)(sy+1) >= (unsigned)ssize.height) )
3218                     continue;
3219
3220                 if( borderType1 == BORDER_CONSTANT &&
3221                     (sx >= ssize.width || sx+4 <= 0 ||
3222                     sy >= ssize.height || sy+4 <= 0))
3223                 {
3224                     for( k = 0; k < cn; k++ )
3225                         D[k] = cval[k];
3226                     continue;
3227                 }
3228
3229                 for( i = 0; i < 4; i++ )
3230                 {
3231                     x[i] = borderInterpolate(sx + i, ssize.width, borderType1)*cn;
3232                     y[i] = borderInterpolate(sy + i, ssize.height, borderType1);
3233                 }
3234
3235                 for( k = 0; k < cn; k++, S0++, w -= 16 )
3236                 {
3237                     WT cv = cval[k], sum = cv*ONE;
3238                     for( i = 0; i < 4; i++, w += 4 )
3239                     {
3240                         int yi = y[i];
3241                         const T* S = S0 + yi*sstep;
3242                         if( yi < 0 )
3243                             continue;
3244                         if( x[0] >= 0 )
3245                             sum += (S[x[0]] - cv)*w[0];
3246                         if( x[1] >= 0 )
3247                             sum += (S[x[1]] - cv)*w[1];
3248                         if( x[2] >= 0 )
3249                             sum += (S[x[2]] - cv)*w[2];
3250                         if( x[3] >= 0 )
3251                             sum += (S[x[3]] - cv)*w[3];
3252                     }
3253                     D[k] = castOp(sum);
3254                 }
3255                 S0 -= cn;
3256             }
3257         }
3258     }
3259 }
3260
3261
3262 template<class CastOp, typename AT, int ONE>
3263 static void remapLanczos4( const Mat& _src, Mat& _dst, const Mat& _xy,
3264                            const Mat& _fxy, const void* _wtab,
3265                            int borderType, const Scalar& _borderValue )
3266 {
3267     typedef typename CastOp::rtype T;
3268     typedef typename CastOp::type1 WT;
3269     Size ssize = _src.size(), dsize = _dst.size();
3270     int cn = _src.channels();
3271     const AT* wtab = (const AT*)_wtab;
3272     const T* S0 = (const T*)_src.data;
3273     size_t sstep = _src.step/sizeof(S0[0]);
3274     Scalar_<T> cval(saturate_cast<T>(_borderValue[0]),
3275         saturate_cast<T>(_borderValue[1]),
3276         saturate_cast<T>(_borderValue[2]),
3277         saturate_cast<T>(_borderValue[3]));
3278     int dx, dy;
3279     CastOp castOp;
3280     int borderType1 = borderType != BORDER_TRANSPARENT ? borderType : BORDER_REFLECT_101;
3281
3282     unsigned width1 = std::max(ssize.width-7, 0), height1 = std::max(ssize.height-7, 0);
3283
3284     if( _dst.isContinuous() && _xy.isContinuous() && _fxy.isContinuous() )
3285     {
3286         dsize.width *= dsize.height;
3287         dsize.height = 1;
3288     }
3289
3290     for( dy = 0; dy < dsize.height; dy++ )
3291     {
3292         T* D = (T*)(_dst.data + _dst.step*dy);
3293         const short* XY = (const short*)(_xy.data + _xy.step*dy);
3294         const ushort* FXY = (const ushort*)(_fxy.data + _fxy.step*dy);
3295
3296         for( dx = 0; dx < dsize.width; dx++, D += cn )
3297         {
3298             int sx = XY[dx*2]-3, sy = XY[dx*2+1]-3;
3299             const AT* w = wtab + FXY[dx]*64;
3300             const T* S = S0 + sy*sstep + sx*cn;
3301             int i, k;
3302             if( (unsigned)sx < width1 && (unsigned)sy < height1 )
3303             {
3304                 for( k = 0; k < cn; k++ )
3305                 {
3306                     WT sum = 0;
3307                     for( int r = 0; r < 8; r++, S += sstep, w += 8 )
3308                         sum += S[0]*w[0] + S[cn]*w[1] + S[cn*2]*w[2] + S[cn*3]*w[3] +
3309                             S[cn*4]*w[4] + S[cn*5]*w[5] + S[cn*6]*w[6] + S[cn*7]*w[7];
3310                     w -= 64;
3311                     S -= sstep*8 - 1;
3312                     D[k] = castOp(sum);
3313                 }
3314             }
3315             else
3316             {
3317                 int x[8], y[8];
3318                 if( borderType == BORDER_TRANSPARENT &&
3319                     ((unsigned)(sx+3) >= (unsigned)ssize.width ||
3320                     (unsigned)(sy+3) >= (unsigned)ssize.height) )
3321                     continue;
3322
3323                 if( borderType1 == BORDER_CONSTANT &&
3324                     (sx >= ssize.width || sx+8 <= 0 ||
3325                     sy >= ssize.height || sy+8 <= 0))
3326                 {
3327                     for( k = 0; k < cn; k++ )
3328                         D[k] = cval[k];
3329                     continue;
3330                 }
3331
3332                 for( i = 0; i < 8; i++ )
3333                 {
3334                     x[i] = borderInterpolate(sx + i, ssize.width, borderType1)*cn;
3335                     y[i] = borderInterpolate(sy + i, ssize.height, borderType1);
3336                 }
3337
3338                 for( k = 0; k < cn; k++, S0++, w -= 64 )
3339                 {
3340                     WT cv = cval[k], sum = cv*ONE;
3341                     for( i = 0; i < 8; i++, w += 8 )
3342                     {
3343                         int yi = y[i];
3344                         const T* S1 = S0 + yi*sstep;
3345                         if( yi < 0 )
3346                             continue;
3347                         if( x[0] >= 0 )
3348                             sum += (S1[x[0]] - cv)*w[0];
3349                         if( x[1] >= 0 )
3350                             sum += (S1[x[1]] - cv)*w[1];
3351                         if( x[2] >= 0 )
3352                             sum += (S1[x[2]] - cv)*w[2];
3353                         if( x[3] >= 0 )
3354                             sum += (S1[x[3]] - cv)*w[3];
3355                         if( x[4] >= 0 )
3356                             sum += (S1[x[4]] - cv)*w[4];
3357                         if( x[5] >= 0 )
3358                             sum += (S1[x[5]] - cv)*w[5];
3359                         if( x[6] >= 0 )
3360                             sum += (S1[x[6]] - cv)*w[6];
3361                         if( x[7] >= 0 )
3362                             sum += (S1[x[7]] - cv)*w[7];
3363                     }
3364                     D[k] = castOp(sum);
3365                 }
3366                 S0 -= cn;
3367             }
3368         }
3369     }
3370 }
3371
3372
3373 typedef void (*RemapNNFunc)(const Mat& _src, Mat& _dst, const Mat& _xy,
3374                             int borderType, const Scalar& _borderValue );
3375
3376 typedef void (*RemapFunc)(const Mat& _src, Mat& _dst, const Mat& _xy,
3377                           const Mat& _fxy, const void* _wtab,
3378                           int borderType, const Scalar& _borderValue);
3379
3380 class RemapInvoker :
3381     public ParallelLoopBody
3382 {
3383 public:
3384     RemapInvoker(const Mat& _src, Mat& _dst, const Mat *_m1,
3385                  const Mat *_m2, int _borderType, const Scalar &_borderValue,
3386                  int _planar_input, RemapNNFunc _nnfunc, RemapFunc _ifunc, const void *_ctab) :
3387         ParallelLoopBody(), src(&_src), dst(&_dst), m1(_m1), m2(_m2),
3388         borderType(_borderType), borderValue(_borderValue),
3389         planar_input(_planar_input), nnfunc(_nnfunc), ifunc(_ifunc), ctab(_ctab)
3390     {
3391     }
3392
3393     virtual void operator() (const Range& range) const
3394     {
3395         int x, y, x1, y1;
3396         const int buf_size = 1 << 14;
3397         int brows0 = std::min(128, dst->rows), map_depth = m1->depth();
3398         int bcols0 = std::min(buf_size/brows0, dst->cols);
3399         brows0 = std::min(buf_size/bcols0, dst->rows);
3400     #if CV_SSE2
3401         bool useSIMD = checkHardwareSupport(CV_CPU_SSE2);
3402     #endif
3403
3404         Mat _bufxy(brows0, bcols0, CV_16SC2), _bufa;
3405         if( !nnfunc )
3406             _bufa.create(brows0, bcols0, CV_16UC1);
3407
3408         for( y = range.start; y < range.end; y += brows0 )
3409         {
3410             for( x = 0; x < dst->cols; x += bcols0 )
3411             {
3412                 int brows = std::min(brows0, range.end - y);
3413                 int bcols = std::min(bcols0, dst->cols - x);
3414                 Mat dpart(*dst, Rect(x, y, bcols, brows));
3415                 Mat bufxy(_bufxy, Rect(0, 0, bcols, brows));
3416
3417                 if( nnfunc )
3418                 {
3419                     if( m1->type() == CV_16SC2 && !m2->data ) // the data is already in the right format
3420                         bufxy = (*m1)(Rect(x, y, bcols, brows));
3421                     else if( map_depth != CV_32F )
3422                     {
3423                         for( y1 = 0; y1 < brows; y1++ )
3424                         {
3425                             short* XY = (short*)(bufxy.data + bufxy.step*y1);
3426                             const short* sXY = (const short*)(m1->data + m1->step*(y+y1)) + x*2;
3427                             const ushort* sA = (const ushort*)(m2->data + m2->step*(y+y1)) + x;
3428
3429                             for( x1 = 0; x1 < bcols; x1++ )
3430                             {
3431                                 int a = sA[x1] & (INTER_TAB_SIZE2-1);
3432                                 XY[x1*2] = sXY[x1*2] + NNDeltaTab_i[a][0];
3433                                 XY[x1*2+1] = sXY[x1*2+1] + NNDeltaTab_i[a][1];
3434                             }
3435                         }
3436                     }
3437                     else if( !planar_input )
3438                         (*m1)(Rect(x, y, bcols, brows)).convertTo(bufxy, bufxy.depth());
3439                     else
3440                     {
3441                         for( y1 = 0; y1 < brows; y1++ )
3442                         {
3443                             short* XY = (short*)(bufxy.data + bufxy.step*y1);
3444                             const float* sX = (const float*)(m1->data + m1->step*(y+y1)) + x;
3445                             const float* sY = (const float*)(m2->data + m2->step*(y+y1)) + x;
3446                             x1 = 0;
3447
3448                         #if CV_SSE2
3449                             if( useSIMD )
3450                             {
3451                                 for( ; x1 <= bcols - 8; x1 += 8 )
3452                                 {
3453                                     __m128 fx0 = _mm_loadu_ps(sX + x1);
3454                                     __m128 fx1 = _mm_loadu_ps(sX + x1 + 4);
3455                                     __m128 fy0 = _mm_loadu_ps(sY + x1);
3456                                     __m128 fy1 = _mm_loadu_ps(sY + x1 + 4);
3457                                     __m128i ix0 = _mm_cvtps_epi32(fx0);
3458                                     __m128i ix1 = _mm_cvtps_epi32(fx1);
3459                                     __m128i iy0 = _mm_cvtps_epi32(fy0);
3460                                     __m128i iy1 = _mm_cvtps_epi32(fy1);
3461                                     ix0 = _mm_packs_epi32(ix0, ix1);
3462                                     iy0 = _mm_packs_epi32(iy0, iy1);
3463                                     ix1 = _mm_unpacklo_epi16(ix0, iy0);
3464                                     iy1 = _mm_unpackhi_epi16(ix0, iy0);
3465                                     _mm_storeu_si128((__m128i*)(XY + x1*2), ix1);
3466                                     _mm_storeu_si128((__m128i*)(XY + x1*2 + 8), iy1);
3467                                 }
3468                             }
3469                         #endif
3470
3471                             for( ; x1 < bcols; x1++ )
3472                             {
3473                                 XY[x1*2] = saturate_cast<short>(sX[x1]);
3474                                 XY[x1*2+1] = saturate_cast<short>(sY[x1]);
3475                             }
3476                         }
3477                     }
3478                     nnfunc( *src, dpart, bufxy, borderType, borderValue );
3479                     continue;
3480                 }
3481
3482                 Mat bufa(_bufa, Rect(0, 0, bcols, brows));
3483                 for( y1 = 0; y1 < brows; y1++ )
3484                 {
3485                     short* XY = (short*)(bufxy.data + bufxy.step*y1);
3486                     ushort* A = (ushort*)(bufa.data + bufa.step*y1);
3487
3488                     if( m1->type() == CV_16SC2 && (m2->type() == CV_16UC1 || m2->type() == CV_16SC1) )
3489                     {
3490                         bufxy = (*m1)(Rect(x, y, bcols, brows));
3491
3492                         const ushort* sA = (const ushort*)(m2->data + m2->step*(y+y1)) + x;
3493                         for( x1 = 0; x1 < bcols; x1++ )
3494                             A[x1] = (ushort)(sA[x1] & (INTER_TAB_SIZE2-1));
3495                     }
3496                     else if( planar_input )
3497                     {
3498                         const float* sX = (const float*)(m1->data + m1->step*(y+y1)) + x;
3499                         const float* sY = (const float*)(m2->data + m2->step*(y+y1)) + x;
3500
3501                         x1 = 0;
3502                     #if CV_SSE2
3503                         if( useSIMD )
3504                         {
3505                             __m128 scale = _mm_set1_ps((float)INTER_TAB_SIZE);
3506                             __m128i mask = _mm_set1_epi32(INTER_TAB_SIZE-1);
3507                             for( ; x1 <= bcols - 8; x1 += 8 )
3508                             {
3509                                 __m128 fx0 = _mm_loadu_ps(sX + x1);
3510                                 __m128 fx1 = _mm_loadu_ps(sX + x1 + 4);
3511                                 __m128 fy0 = _mm_loadu_ps(sY + x1);
3512                                 __m128 fy1 = _mm_loadu_ps(sY + x1 + 4);
3513                                 __m128i ix0 = _mm_cvtps_epi32(_mm_mul_ps(fx0, scale));
3514                                 __m128i ix1 = _mm_cvtps_epi32(_mm_mul_ps(fx1, scale));
3515                                 __m128i iy0 = _mm_cvtps_epi32(_mm_mul_ps(fy0, scale));
3516                                 __m128i iy1 = _mm_cvtps_epi32(_mm_mul_ps(fy1, scale));
3517                                 __m128i mx0 = _mm_and_si128(ix0, mask);
3518                                 __m128i mx1 = _mm_and_si128(ix1, mask);
3519                                 __m128i my0 = _mm_and_si128(iy0, mask);
3520                                 __m128i my1 = _mm_and_si128(iy1, mask);
3521                                 mx0 = _mm_packs_epi32(mx0, mx1);
3522                                 my0 = _mm_packs_epi32(my0, my1);
3523                                 my0 = _mm_slli_epi16(my0, INTER_BITS);
3524                                 mx0 = _mm_or_si128(mx0, my0);
3525                                 _mm_storeu_si128((__m128i*)(A + x1), mx0);
3526                                 ix0 = _mm_srai_epi32(ix0, INTER_BITS);
3527                                 ix1 = _mm_srai_epi32(ix1, INTER_BITS);
3528                                 iy0 = _mm_srai_epi32(iy0, INTER_BITS);
3529                                 iy1 = _mm_srai_epi32(iy1, INTER_BITS);
3530                                 ix0 = _mm_packs_epi32(ix0, ix1);
3531                                 iy0 = _mm_packs_epi32(iy0, iy1);
3532                                 ix1 = _mm_unpacklo_epi16(ix0, iy0);
3533                                 iy1 = _mm_unpackhi_epi16(ix0, iy0);
3534                                 _mm_storeu_si128((__m128i*)(XY + x1*2), ix1);
3535                                 _mm_storeu_si128((__m128i*)(XY + x1*2 + 8), iy1);
3536                             }
3537                         }
3538                     #endif
3539
3540                         for( ; x1 < bcols; x1++ )
3541                         {
3542                             int sx = cvRound(sX[x1]*INTER_TAB_SIZE);
3543                             int sy = cvRound(sY[x1]*INTER_TAB_SIZE);
3544                             int v = (sy & (INTER_TAB_SIZE-1))*INTER_TAB_SIZE + (sx & (INTER_TAB_SIZE-1));
3545                             XY[x1*2] = saturate_cast<short>(sx >> INTER_BITS);
3546                             XY[x1*2+1] = saturate_cast<short>(sy >> INTER_BITS);
3547                             A[x1] = (ushort)v;
3548                         }
3549                     }
3550                     else
3551                     {
3552                         const float* sXY = (const float*)(m1->data + m1->step*(y+y1)) + x*2;
3553
3554                         for( x1 = 0; x1 < bcols; x1++ )
3555                         {
3556                             int sx = cvRound(sXY[x1*2]*INTER_TAB_SIZE);
3557                             int sy = cvRound(sXY[x1*2+1]*INTER_TAB_SIZE);
3558                             int v = (sy & (INTER_TAB_SIZE-1))*INTER_TAB_SIZE + (sx & (INTER_TAB_SIZE-1));
3559                             XY[x1*2] = saturate_cast<short>(sx >> INTER_BITS);
3560                             XY[x1*2+1] = saturate_cast<short>(sy >> INTER_BITS);
3561                             A[x1] = (ushort)v;
3562                         }
3563                     }
3564                 }
3565                 ifunc(*src, dpart, bufxy, bufa, ctab, borderType, borderValue);
3566             }
3567         }
3568     }
3569
3570 private:
3571     const Mat* src;
3572     Mat* dst;
3573     const Mat *m1, *m2;
3574     int borderType;
3575     Scalar borderValue;
3576     int planar_input;
3577     RemapNNFunc nnfunc;
3578     RemapFunc ifunc;
3579     const void *ctab;
3580 };
3581
3582 #ifdef HAVE_OPENCL
3583
3584 static bool ocl_remap(InputArray _src, OutputArray _dst, InputArray _map1, InputArray _map2,
3585                       int interpolation, int borderType, const Scalar& borderValue)
3586 {
3587     int cn = _src.channels(), type = _src.type(), depth = _src.depth();
3588
3589     if (borderType == BORDER_TRANSPARENT || cn == 3 || !(interpolation == INTER_LINEAR || interpolation == INTER_NEAREST)
3590             || _map1.type() == CV_16SC1 || _map2.type() == CV_16SC1)
3591         return false;
3592
3593     UMat src = _src.getUMat(), map1 = _map1.getUMat(), map2 = _map2.getUMat();
3594
3595     if( (map1.type() == CV_16SC2 && (map2.type() == CV_16UC1 || map2.empty())) ||
3596         (map2.type() == CV_16SC2 && (map1.type() == CV_16UC1 || map1.empty())) )
3597     {
3598         if (map1.type() != CV_16SC2)
3599             std::swap(map1, map2);
3600     }
3601     else
3602         CV_Assert( map1.type() == CV_32FC2 || (map1.type() == CV_32FC1 && map2.type() == CV_32FC1) );
3603
3604     _dst.create(map1.size(), type);
3605     UMat dst = _dst.getUMat();
3606
3607     String kernelName = "remap";
3608     if (map1.type() == CV_32FC2 && map2.empty())
3609         kernelName += "_32FC2";
3610     else if (map1.type() == CV_16SC2)
3611     {
3612         kernelName += "_16SC2";
3613         if (!map2.empty())
3614             kernelName += "_16UC1";
3615     }
3616     else if (map1.type() == CV_32FC1 && map2.type() == CV_32FC1)
3617         kernelName += "_2_32FC1";
3618     else
3619         CV_Error(Error::StsBadArg, "Unsupported map types");
3620
3621     static const char * const interMap[] = { "INTER_NEAREST", "INTER_LINEAR", "INTER_CUBIC", "INTER_LINEAR", "INTER_LANCZOS" };
3622     static const char * const borderMap[] = { "BORDER_CONSTANT", "BORDER_REPLICATE", "BORDER_REFLECT", "BORDER_WRAP",
3623                            "BORDER_REFLECT_101", "BORDER_TRANSPARENT" };
3624     String buildOptions = format("-D %s -D %s -D T=%s", interMap[interpolation], borderMap[borderType], ocl::typeToStr(type));
3625
3626     if (interpolation != INTER_NEAREST)
3627     {
3628         char cvt[3][40];
3629         int wdepth = std::max(CV_32F, dst.depth());
3630         buildOptions = buildOptions
3631                       + format(" -D WT=%s -D convertToT=%s -D convertToWT=%s"
3632                                " -D convertToWT2=%s -D WT2=%s",
3633                                ocl::typeToStr(CV_MAKE_TYPE(wdepth, cn)),
3634                                ocl::convertTypeStr(wdepth, depth, cn, cvt[0]),
3635                                ocl::convertTypeStr(depth, wdepth, cn, cvt[1]),
3636                                ocl::convertTypeStr(CV_32S, wdepth, 2, cvt[2]),
3637                                ocl::typeToStr(CV_MAKE_TYPE(wdepth, 2)));
3638     }
3639
3640     ocl::Kernel k(kernelName.c_str(), ocl::imgproc::remap_oclsrc, buildOptions);
3641
3642     Mat scalar(1, 1, type, borderValue);
3643     ocl::KernelArg srcarg = ocl::KernelArg::ReadOnly(src), dstarg = ocl::KernelArg::WriteOnly(dst),
3644             map1arg = ocl::KernelArg::ReadOnlyNoSize(map1),
3645             scalararg = ocl::KernelArg::Constant((void*)scalar.data, scalar.elemSize());
3646
3647     if (map2.empty())
3648         k.args(srcarg, dstarg, map1arg, scalararg);
3649     else
3650         k.args(srcarg, dstarg, map1arg, ocl::KernelArg::ReadOnlyNoSize(map2), scalararg);
3651
3652     size_t globalThreads[2] = { dst.cols, dst.rows };
3653     return k.run(2, globalThreads, NULL, false);
3654 }
3655
3656 #endif
3657
3658 }
3659
3660 void cv::remap( InputArray _src, OutputArray _dst,
3661                 InputArray _map1, InputArray _map2,
3662                 int interpolation, int borderType, const Scalar& borderValue )
3663 {
3664     static RemapNNFunc nn_tab[] =
3665     {
3666         remapNearest<uchar>, remapNearest<schar>, remapNearest<ushort>, remapNearest<short>,
3667         remapNearest<int>, remapNearest<float>, remapNearest<double>, 0
3668     };
3669
3670     static RemapFunc linear_tab[] =
3671     {
3672         remapBilinear<FixedPtCast<int, uchar, INTER_REMAP_COEF_BITS>, RemapVec_8u, short>, 0,
3673         remapBilinear<Cast<float, ushort>, RemapNoVec, float>,
3674         remapBilinear<Cast<float, short>, RemapNoVec, float>, 0,
3675         remapBilinear<Cast<float, float>, RemapNoVec, float>,
3676         remapBilinear<Cast<double, double>, RemapNoVec, float>, 0
3677     };
3678
3679     static RemapFunc cubic_tab[] =
3680     {
3681         remapBicubic<FixedPtCast<int, uchar, INTER_REMAP_COEF_BITS>, short, INTER_REMAP_COEF_SCALE>, 0,
3682         remapBicubic<Cast<float, ushort>, float, 1>,
3683         remapBicubic<Cast<float, short>, float, 1>, 0,
3684         remapBicubic<Cast<float, float>, float, 1>,
3685         remapBicubic<Cast<double, double>, float, 1>, 0
3686     };
3687
3688     static RemapFunc lanczos4_tab[] =
3689     {
3690         remapLanczos4<FixedPtCast<int, uchar, INTER_REMAP_COEF_BITS>, short, INTER_REMAP_COEF_SCALE>, 0,
3691         remapLanczos4<Cast<float, ushort>, float, 1>,
3692         remapLanczos4<Cast<float, short>, float, 1>, 0,
3693         remapLanczos4<Cast<float, float>, float, 1>,
3694         remapLanczos4<Cast<double, double>, float, 1>, 0
3695     };
3696
3697     CV_Assert( _map1.size().area() > 0 );
3698     CV_Assert( _map2.empty() || (_map2.size() == _map1.size()));
3699
3700     CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
3701                ocl_remap(_src, _dst, _map1, _map2, interpolation, borderType, borderValue))
3702
3703     Mat src = _src.getMat(), map1 = _map1.getMat(), map2 = _map2.getMat();
3704     _dst.create( map1.size(), src.type() );
3705     Mat dst = _dst.getMat();
3706     if( dst.data == src.data )
3707         src = src.clone();
3708
3709     int depth = src.depth();
3710     RemapNNFunc nnfunc = 0;
3711     RemapFunc ifunc = 0;
3712     const void* ctab = 0;
3713     bool fixpt = depth == CV_8U;
3714     bool planar_input = false;
3715
3716     if( interpolation == INTER_NEAREST )
3717     {
3718         nnfunc = nn_tab[depth];
3719         CV_Assert( nnfunc != 0 );
3720     }
3721     else
3722     {
3723         if( interpolation == INTER_AREA )
3724             interpolation = INTER_LINEAR;
3725
3726         if( interpolation == INTER_LINEAR )
3727             ifunc = linear_tab[depth];
3728         else if( interpolation == INTER_CUBIC )
3729             ifunc = cubic_tab[depth];
3730         else if( interpolation == INTER_LANCZOS4 )
3731             ifunc = lanczos4_tab[depth];
3732         else
3733             CV_Error( CV_StsBadArg, "Unknown interpolation method" );
3734         CV_Assert( ifunc != 0 );
3735         ctab = initInterTab2D( interpolation, fixpt );
3736     }
3737
3738     const Mat *m1 = &map1, *m2 = &map2;
3739
3740     if( (map1.type() == CV_16SC2 && (map2.type() == CV_16UC1 || map2.type() == CV_16SC1 || !map2.data)) ||
3741         (map2.type() == CV_16SC2 && (map1.type() == CV_16UC1 || map1.type() == CV_16SC1 || !map1.data)) )
3742     {
3743         if( map1.type() != CV_16SC2 )
3744             std::swap(m1, m2);
3745     }
3746     else
3747     {
3748         CV_Assert( ((map1.type() == CV_32FC2 || map1.type() == CV_16SC2) && !map2.data) ||
3749             (map1.type() == CV_32FC1 && map2.type() == CV_32FC1) );
3750         planar_input = map1.channels() == 1;
3751     }
3752
3753     RemapInvoker invoker(src, dst, m1, m2,
3754                          borderType, borderValue, planar_input, nnfunc, ifunc,
3755                          ctab);
3756     parallel_for_(Range(0, dst.rows), invoker, dst.total()/(double)(1<<16));
3757 }
3758
3759
3760 void cv::convertMaps( InputArray _map1, InputArray _map2,
3761                       OutputArray _dstmap1, OutputArray _dstmap2,
3762                       int dstm1type, bool nninterpolate )
3763 {
3764     Mat map1 = _map1.getMat(), map2 = _map2.getMat(), dstmap1, dstmap2;
3765     Size size = map1.size();
3766     const Mat *m1 = &map1, *m2 = &map2;
3767     int m1type = m1->type(), m2type = m2->type();
3768
3769     CV_Assert( (m1type == CV_16SC2 && (nninterpolate || m2type == CV_16UC1 || m2type == CV_16SC1)) ||
3770                (m2type == CV_16SC2 && (nninterpolate || m1type == CV_16UC1 || m1type == CV_16SC1)) ||
3771                (m1type == CV_32FC1 && m2type == CV_32FC1) ||
3772                (m1type == CV_32FC2 && !m2->data) );
3773
3774     if( m2type == CV_16SC2 )
3775     {
3776         std::swap( m1, m2 );
3777         std::swap( m1type, m2type );
3778     }
3779
3780     if( dstm1type <= 0 )
3781         dstm1type = m1type == CV_16SC2 ? CV_32FC2 : CV_16SC2;
3782     CV_Assert( dstm1type == CV_16SC2 || dstm1type == CV_32FC1 || dstm1type == CV_32FC2 );
3783     _dstmap1.create( size, dstm1type );
3784     dstmap1 = _dstmap1.getMat();
3785
3786     if( !nninterpolate && dstm1type != CV_32FC2 )
3787     {
3788         _dstmap2.create( size, dstm1type == CV_16SC2 ? CV_16UC1 : CV_32FC1 );
3789         dstmap2 = _dstmap2.getMat();
3790     }
3791     else
3792         _dstmap2.release();
3793
3794     if( m1type == dstm1type || (nninterpolate &&
3795         ((m1type == CV_16SC2 && dstm1type == CV_32FC2) ||
3796         (m1type == CV_32FC2 && dstm1type == CV_16SC2))) )
3797     {
3798         m1->convertTo( dstmap1, dstmap1.type() );
3799         if( dstmap2.data && dstmap2.type() == m2->type() )
3800             m2->copyTo( dstmap2 );
3801         return;
3802     }
3803
3804     if( m1type == CV_32FC1 && dstm1type == CV_32FC2 )
3805     {
3806         Mat vdata[] = { *m1, *m2 };
3807         merge( vdata, 2, dstmap1 );
3808         return;
3809     }
3810
3811     if( m1type == CV_32FC2 && dstm1type == CV_32FC1 )
3812     {
3813         Mat mv[] = { dstmap1, dstmap2 };
3814         split( *m1, mv );
3815         return;
3816     }
3817
3818     if( m1->isContinuous() && (!m2->data || m2->isContinuous()) &&
3819         dstmap1.isContinuous() && (!dstmap2.data || dstmap2.isContinuous()) )
3820     {
3821         size.width *= size.height;
3822         size.height = 1;
3823     }
3824
3825     const float scale = 1.f/INTER_TAB_SIZE;
3826     int x, y;
3827     for( y = 0; y < size.height; y++ )
3828     {
3829         const float* src1f = (const float*)(m1->data + m1->step*y);
3830         const float* src2f = (const float*)(m2->data + m2->step*y);
3831         const short* src1 = (const short*)src1f;
3832         const ushort* src2 = (const ushort*)src2f;
3833
3834         float* dst1f = (float*)(dstmap1.data + dstmap1.step*y);
3835         float* dst2f = (float*)(dstmap2.data + dstmap2.step*y);
3836         short* dst1 = (short*)dst1f;
3837         ushort* dst2 = (ushort*)dst2f;
3838
3839         if( m1type == CV_32FC1 && dstm1type == CV_16SC2 )
3840         {
3841             if( nninterpolate )
3842                 for( x = 0; x < size.width; x++ )
3843                 {
3844                     dst1[x*2] = saturate_cast<short>(src1f[x]);
3845                     dst1[x*2+1] = saturate_cast<short>(src2f[x]);
3846                 }
3847             else
3848                 for( x = 0; x < size.width; x++ )
3849                 {
3850                     int ix = saturate_cast<int>(src1f[x]*INTER_TAB_SIZE);
3851                     int iy = saturate_cast<int>(src2f[x]*INTER_TAB_SIZE);
3852                     dst1[x*2] = saturate_cast<short>(ix >> INTER_BITS);
3853                     dst1[x*2+1] = saturate_cast<short>(iy >> INTER_BITS);
3854                     dst2[x] = (ushort)((iy & (INTER_TAB_SIZE-1))*INTER_TAB_SIZE + (ix & (INTER_TAB_SIZE-1)));
3855                 }
3856         }
3857         else if( m1type == CV_32FC2 && dstm1type == CV_16SC2 )
3858         {
3859             if( nninterpolate )
3860                 for( x = 0; x < size.width; x++ )
3861                 {
3862                     dst1[x*2] = saturate_cast<short>(src1f[x*2]);
3863                     dst1[x*2+1] = saturate_cast<short>(src1f[x*2+1]);
3864                 }
3865             else
3866                 for( x = 0; x < size.width; x++ )
3867                 {
3868                     int ix = saturate_cast<int>(src1f[x*2]*INTER_TAB_SIZE);
3869                     int iy = saturate_cast<int>(src1f[x*2+1]*INTER_TAB_SIZE);
3870                     dst1[x*2] = saturate_cast<short>(ix >> INTER_BITS);
3871                     dst1[x*2+1] = saturate_cast<short>(iy >> INTER_BITS);
3872                     dst2[x] = (ushort)((iy & (INTER_TAB_SIZE-1))*INTER_TAB_SIZE + (ix & (INTER_TAB_SIZE-1)));
3873                 }
3874         }
3875         else if( m1type == CV_16SC2 && dstm1type == CV_32FC1 )
3876         {
3877             for( x = 0; x < size.width; x++ )
3878             {
3879                 int fxy = src2 ? src2[x] & (INTER_TAB_SIZE2-1) : 0;
3880                 dst1f[x] = src1[x*2] + (fxy & (INTER_TAB_SIZE-1))*scale;
3881                 dst2f[x] = src1[x*2+1] + (fxy >> INTER_BITS)*scale;
3882             }
3883         }
3884         else if( m1type == CV_16SC2 && dstm1type == CV_32FC2 )
3885         {
3886             for( x = 0; x < size.width; x++ )
3887             {
3888                 int fxy = src2 ? src2[x] & (INTER_TAB_SIZE2-1): 0;
3889                 dst1f[x*2] = src1[x*2] + (fxy & (INTER_TAB_SIZE-1))*scale;
3890                 dst1f[x*2+1] = src1[x*2+1] + (fxy >> INTER_BITS)*scale;
3891             }
3892         }
3893         else
3894             CV_Error( CV_StsNotImplemented, "Unsupported combination of input/output matrices" );
3895     }
3896 }
3897
3898
3899 namespace cv
3900 {
3901
3902 class warpAffineInvoker :
3903     public ParallelLoopBody
3904 {
3905 public:
3906     warpAffineInvoker(const Mat &_src, Mat &_dst, int _interpolation, int _borderType,
3907                       const Scalar &_borderValue, int *_adelta, int *_bdelta, double *_M) :
3908         ParallelLoopBody(), src(_src), dst(_dst), interpolation(_interpolation),
3909         borderType(_borderType), borderValue(_borderValue), adelta(_adelta), bdelta(_bdelta),
3910         M(_M)
3911     {
3912     }
3913
3914     virtual void operator() (const Range& range) const
3915     {
3916         const int BLOCK_SZ = 64;
3917         short XY[BLOCK_SZ*BLOCK_SZ*2], A[BLOCK_SZ*BLOCK_SZ];
3918         const int AB_BITS = MAX(10, (int)INTER_BITS);
3919         const int AB_SCALE = 1 << AB_BITS;
3920         int round_delta = interpolation == INTER_NEAREST ? AB_SCALE/2 : AB_SCALE/INTER_TAB_SIZE/2, x, y, x1, y1;
3921     #if CV_SSE2
3922         bool useSIMD = checkHardwareSupport(CV_CPU_SSE2);
3923     #endif
3924
3925         int bh0 = std::min(BLOCK_SZ/2, dst.rows);
3926         int bw0 = std::min(BLOCK_SZ*BLOCK_SZ/bh0, dst.cols);
3927         bh0 = std::min(BLOCK_SZ*BLOCK_SZ/bw0, dst.rows);
3928
3929         for( y = range.start; y < range.end; y += bh0 )
3930         {
3931             for( x = 0; x < dst.cols; x += bw0 )
3932             {
3933                 int bw = std::min( bw0, dst.cols - x);
3934                 int bh = std::min( bh0, range.end - y);
3935
3936                 Mat _XY(bh, bw, CV_16SC2, XY), matA;
3937                 Mat dpart(dst, Rect(x, y, bw, bh));
3938
3939                 for( y1 = 0; y1 < bh; y1++ )
3940                 {
3941                     short* xy = XY + y1*bw*2;
3942                     int X0 = saturate_cast<int>((M[1]*(y + y1) + M[2])*AB_SCALE) + round_delta;
3943                     int Y0 = saturate_cast<int>((M[4]*(y + y1) + M[5])*AB_SCALE) + round_delta;
3944
3945                     if( interpolation == INTER_NEAREST )
3946                         for( x1 = 0; x1 < bw; x1++ )
3947                         {
3948                             int X = (X0 + adelta[x+x1]) >> AB_BITS;
3949                             int Y = (Y0 + bdelta[x+x1]) >> AB_BITS;
3950                             xy[x1*2] = saturate_cast<short>(X);
3951                             xy[x1*2+1] = saturate_cast<short>(Y);
3952                         }
3953                     else
3954                     {
3955                         short* alpha = A + y1*bw;
3956                         x1 = 0;
3957                     #if CV_SSE2
3958                         if( useSIMD )
3959                         {
3960                             __m128i fxy_mask = _mm_set1_epi32(INTER_TAB_SIZE - 1);
3961                             __m128i XX = _mm_set1_epi32(X0), YY = _mm_set1_epi32(Y0);
3962                             for( ; x1 <= bw - 8; x1 += 8 )
3963                             {
3964                                 __m128i tx0, tx1, ty0, ty1;
3965                                 tx0 = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(adelta + x + x1)), XX);
3966                                 ty0 = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(bdelta + x + x1)), YY);
3967                                 tx1 = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(adelta + x + x1 + 4)), XX);
3968                                 ty1 = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(bdelta + x + x1 + 4)), YY);
3969
3970                                 tx0 = _mm_srai_epi32(tx0, AB_BITS - INTER_BITS);
3971                                 ty0 = _mm_srai_epi32(ty0, AB_BITS - INTER_BITS);
3972                                 tx1 = _mm_srai_epi32(tx1, AB_BITS - INTER_BITS);
3973                                 ty1 = _mm_srai_epi32(ty1, AB_BITS - INTER_BITS);
3974
3975                                 __m128i fx_ = _mm_packs_epi32(_mm_and_si128(tx0, fxy_mask),
3976                                                             _mm_and_si128(tx1, fxy_mask));
3977                                 __m128i fy_ = _mm_packs_epi32(_mm_and_si128(ty0, fxy_mask),
3978                                                             _mm_and_si128(ty1, fxy_mask));
3979                                 tx0 = _mm_packs_epi32(_mm_srai_epi32(tx0, INTER_BITS),
3980                                                             _mm_srai_epi32(tx1, INTER_BITS));
3981                                 ty0 = _mm_packs_epi32(_mm_srai_epi32(ty0, INTER_BITS),
3982                                                     _mm_srai_epi32(ty1, INTER_BITS));
3983                                 fx_ = _mm_adds_epi16(fx_, _mm_slli_epi16(fy_, INTER_BITS));
3984
3985                                 _mm_storeu_si128((__m128i*)(xy + x1*2), _mm_unpacklo_epi16(tx0, ty0));
3986                                 _mm_storeu_si128((__m128i*)(xy + x1*2 + 8), _mm_unpackhi_epi16(tx0, ty0));
3987                                 _mm_storeu_si128((__m128i*)(alpha + x1), fx_);
3988                             }
3989                         }
3990                     #endif
3991                         for( ; x1 < bw; x1++ )
3992                         {
3993                             int X = (X0 + adelta[x+x1]) >> (AB_BITS - INTER_BITS);
3994                             int Y = (Y0 + bdelta[x+x1]) >> (AB_BITS - INTER_BITS);
3995                             xy[x1*2] = saturate_cast<short>(X >> INTER_BITS);
3996                             xy[x1*2+1] = saturate_cast<short>(Y >> INTER_BITS);
3997                             alpha[x1] = (short)((Y & (INTER_TAB_SIZE-1))*INTER_TAB_SIZE +
3998                                     (X & (INTER_TAB_SIZE-1)));
3999                         }
4000                     }
4001                 }
4002
4003                 if( interpolation == INTER_NEAREST )
4004                     remap( src, dpart, _XY, Mat(), interpolation, borderType, borderValue );
4005                 else
4006                 {
4007                     Mat _matA(bh, bw, CV_16U, A);
4008                     remap( src, dpart, _XY, _matA, interpolation, borderType, borderValue );
4009                 }
4010             }
4011         }
4012     }
4013
4014 private:
4015     Mat src;
4016     Mat dst;
4017     int interpolation, borderType;
4018     Scalar borderValue;
4019     int *adelta, *bdelta;
4020     double *M;
4021 };
4022
4023 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
4024 class IPPwarpAffineInvoker :
4025     public ParallelLoopBody
4026 {
4027 public:
4028     IPPwarpAffineInvoker(Mat &_src, Mat &_dst, double (&_coeffs)[2][3], int &_interpolation, int &_borderType, const Scalar &_borderValue, ippiWarpAffineBackFunc _func, bool *_ok) :
4029       ParallelLoopBody(), src(_src), dst(_dst), mode(_interpolation), coeffs(_coeffs), borderType(_borderType), borderValue(_borderValue), func(_func), ok(_ok)
4030       {
4031           *ok = true;
4032       }
4033
4034       virtual void operator() (const Range& range) const
4035       {
4036           IppiSize srcsize = { src.cols, src.rows };
4037           IppiRect srcroi = { 0, 0, src.cols, src.rows };
4038           IppiRect dstroi = { 0, range.start, dst.cols, range.end - range.start };
4039           int cnn = src.channels();
4040           if( borderType == BORDER_CONSTANT )
4041           {
4042               IppiSize setSize = { dst.cols, range.end - range.start };
4043               void *dataPointer = dst.data + dst.step[0] * range.start;
4044               if( !IPPSet( borderValue, dataPointer, (int)dst.step[0], setSize, cnn, src.depth() ) )
4045               {
4046                   *ok = false;
4047                   return;
4048               }
4049           }
4050           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
4051               *ok = false;
4052       }
4053 private:
4054     Mat &src;
4055     Mat &dst;
4056     double (&coeffs)[2][3];
4057     int mode;
4058     int borderType;
4059     Scalar borderValue;
4060     ippiWarpAffineBackFunc func;
4061     bool *ok;
4062     const IPPwarpAffineInvoker& operator= (const IPPwarpAffineInvoker&);
4063 };
4064 #endif
4065
4066 #ifdef HAVE_OPENCL
4067
4068 enum { OCL_OP_PERSPECTIVE = 1, OCL_OP_AFFINE = 0 };
4069
4070 static bool ocl_warpTransform(InputArray _src, OutputArray _dst, InputArray _M0,
4071                               Size dsize, int flags, int borderType, const Scalar& borderValue,
4072                               int op_type)
4073 {
4074     CV_Assert(op_type == OCL_OP_AFFINE || op_type == OCL_OP_PERSPECTIVE);
4075
4076     int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
4077     double doubleSupport = ocl::Device::getDefault().doubleFPConfig() > 0;
4078
4079     int interpolation = flags & INTER_MAX;
4080     if( interpolation == INTER_AREA )
4081         interpolation = INTER_LINEAR;
4082
4083     if ( !(borderType == cv::BORDER_CONSTANT &&
4084            (interpolation == cv::INTER_NEAREST || interpolation == cv::INTER_LINEAR || interpolation == cv::INTER_CUBIC)) ||
4085          (!doubleSupport && depth == CV_64F) || cn > 4)
4086         return false;
4087
4088     const char * const interpolationMap[3] = { "NEAREST", "LINEAR", "CUBIC" };
4089     ocl::ProgramSource program = op_type == OCL_OP_AFFINE ?
4090                 ocl::imgproc::warp_affine_oclsrc : ocl::imgproc::warp_perspective_oclsrc;
4091     const char * const kernelName = op_type == OCL_OP_AFFINE ? "warpAffine" : "warpPerspective";
4092
4093     int scalarcn = cn == 3 ? 4 : cn;
4094     int wdepth = interpolation == INTER_NEAREST ? depth : std::max(CV_32S, depth);
4095     int sctype = CV_MAKETYPE(wdepth, scalarcn);
4096
4097     ocl::Kernel k;
4098     String opts;
4099     if (interpolation == INTER_NEAREST)
4100     {
4101         opts = format("-D INTER_NEAREST -D T=%s%s -D T1=%s -D ST=%s -D cn=%d", ocl::typeToStr(type),
4102                       doubleSupport ? " -D DOUBLE_SUPPORT" : "",
4103                       ocl::typeToStr(CV_MAT_DEPTH(type)),
4104                       ocl::typeToStr(sctype),
4105                       cn);
4106     }
4107     else
4108     {
4109         char cvt[2][50];
4110         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",
4111                       interpolationMap[interpolation], ocl::typeToStr(type),
4112                       ocl::typeToStr(CV_MAT_DEPTH(type)),
4113                       ocl::typeToStr(sctype),
4114                       ocl::typeToStr(CV_MAKE_TYPE(wdepth, cn)), depth,
4115                       ocl::convertTypeStr(depth, wdepth, cn, cvt[0]),
4116                       ocl::convertTypeStr(wdepth, depth, cn, cvt[1]),
4117                       doubleSupport ? " -D DOUBLE_SUPPORT" : "", cn);
4118     }
4119
4120     k.create(kernelName, program, opts);
4121     if (k.empty())
4122         return false;
4123
4124     double borderBuf[] = {0, 0, 0, 0};
4125     scalarToRawData(borderValue, borderBuf, sctype);
4126
4127     UMat src = _src.getUMat(), M0;
4128     _dst.create( dsize.area() == 0 ? src.size() : dsize, src.type() );
4129     UMat dst = _dst.getUMat();
4130
4131     double M[9];
4132     int matRows = (op_type == OCL_OP_AFFINE ? 2 : 3);
4133     Mat matM(matRows, 3, CV_64F, M), M1 = _M0.getMat();
4134     CV_Assert( (M1.type() == CV_32F || M1.type() == CV_64F) &&
4135                M1.rows == matRows && M1.cols == 3 );
4136     M1.convertTo(matM, matM.type());
4137
4138     if( !(flags & WARP_INVERSE_MAP) )
4139     {
4140         if (op_type == OCL_OP_PERSPECTIVE)
4141             invert(matM, matM);
4142         else
4143         {
4144             double D = M[0]*M[4] - M[1]*M[3];
4145             D = D != 0 ? 1./D : 0;
4146             double A11 = M[4]*D, A22=M[0]*D;
4147             M[0] = A11; M[1] *= -D;
4148             M[3] *= -D; M[4] = A22;
4149             double b1 = -M[0]*M[2] - M[1]*M[5];
4150             double b2 = -M[3]*M[2] - M[4]*M[5];
4151             M[2] = b1; M[5] = b2;
4152         }
4153     }
4154     matM.convertTo(M0, doubleSupport ? CV_64F : CV_32F);
4155
4156     k.args(ocl::KernelArg::ReadOnly(src), ocl::KernelArg::WriteOnly(dst), ocl::KernelArg::PtrReadOnly(M0),
4157            ocl::KernelArg(0, 0, 0, borderBuf, CV_ELEM_SIZE(sctype)));
4158
4159     size_t globalThreads[2] = { dst.cols, dst.rows };
4160     return k.run(2, globalThreads, NULL, false);
4161 }
4162
4163 #endif
4164
4165 }
4166
4167
4168 void cv::warpAffine( InputArray _src, OutputArray _dst,
4169                      InputArray _M0, Size dsize,
4170                      int flags, int borderType, const Scalar& borderValue )
4171 {
4172     CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
4173                ocl_warpTransform(_src, _dst, _M0, dsize, flags, borderType,
4174                                  borderValue, OCL_OP_AFFINE))
4175
4176     Mat src = _src.getMat(), M0 = _M0.getMat();
4177     _dst.create( dsize.area() == 0 ? src.size() : dsize, src.type() );
4178     Mat dst = _dst.getMat();
4179     CV_Assert( src.cols > 0 && src.rows > 0 );
4180     if( dst.data == src.data )
4181         src = src.clone();
4182
4183     double M[6];
4184     Mat matM(2, 3, CV_64F, M);
4185     int interpolation = flags & INTER_MAX;
4186     if( interpolation == INTER_AREA )
4187         interpolation = INTER_LINEAR;
4188
4189     CV_Assert( (M0.type() == CV_32F || M0.type() == CV_64F) && M0.rows == 2 && M0.cols == 3 );
4190     M0.convertTo(matM, matM.type());
4191
4192 #ifdef HAVE_TEGRA_OPTIMIZATION
4193     if( tegra::warpAffine(src, dst, M, flags, borderType, borderValue) )
4194         return;
4195 #endif
4196
4197     if( !(flags & WARP_INVERSE_MAP) )
4198     {
4199         double D = M[0]*M[4] - M[1]*M[3];
4200         D = D != 0 ? 1./D : 0;
4201         double A11 = M[4]*D, A22=M[0]*D;
4202         M[0] = A11; M[1] *= -D;
4203         M[3] *= -D; M[4] = A22;
4204         double b1 = -M[0]*M[2] - M[1]*M[5];
4205         double b2 = -M[3]*M[2] - M[4]*M[5];
4206         M[2] = b1; M[5] = b2;
4207     }
4208
4209     int x;
4210     AutoBuffer<int> _abdelta(dst.cols*2);
4211     int* adelta = &_abdelta[0], *bdelta = adelta + dst.cols;
4212     const int AB_BITS = MAX(10, (int)INTER_BITS);
4213     const int AB_SCALE = 1 << AB_BITS;
4214 /*
4215 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
4216     int depth = src.depth();
4217     int channels = src.channels();
4218     if( ( depth == CV_8U || depth == CV_16U || depth == CV_32F ) &&
4219         ( channels == 1 || channels == 3 || channels == 4 ) &&
4220         ( borderType == cv::BORDER_TRANSPARENT || ( borderType == cv::BORDER_CONSTANT ) ) )
4221     {
4222         int type = src.type();
4223         ippiWarpAffineBackFunc ippFunc =
4224             type == CV_8UC1 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_8u_C1R :
4225             type == CV_8UC3 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_8u_C3R :
4226             type == CV_8UC4 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_8u_C4R :
4227             type == CV_16UC1 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_16u_C1R :
4228             type == CV_16UC3 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_16u_C3R :
4229             type == CV_16UC4 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_16u_C4R :
4230             type == CV_32FC1 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_32f_C1R :
4231             type == CV_32FC3 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_32f_C3R :
4232             type == CV_32FC4 ? (ippiWarpAffineBackFunc)ippiWarpAffineBack_32f_C4R :
4233             0;
4234         int mode =
4235             flags == INTER_LINEAR ? IPPI_INTER_LINEAR :
4236             flags == INTER_NEAREST ? IPPI_INTER_NN :
4237             flags == INTER_CUBIC ? IPPI_INTER_CUBIC :
4238             0;
4239         if( mode && ippFunc )
4240         {
4241             double coeffs[2][3];
4242             for( int i = 0; i < 2; i++ )
4243             {
4244                 for( int j = 0; j < 3; j++ )
4245                 {
4246                     coeffs[i][j] = matM.at<double>(i, j);
4247                 }
4248             }
4249             bool ok;
4250             Range range(0, dst.rows);
4251             IPPwarpAffineInvoker invoker(src, dst, coeffs, mode, borderType, borderValue, ippFunc, &ok);
4252             parallel_for_(range, invoker, dst.total()/(double)(1<<16));
4253             if( ok )
4254                 return;
4255         }
4256     }
4257 #endif
4258 */
4259     for( x = 0; x < dst.cols; x++ )
4260     {
4261         adelta[x] = saturate_cast<int>(M[0]*x*AB_SCALE);
4262         bdelta[x] = saturate_cast<int>(M[3]*x*AB_SCALE);
4263     }
4264
4265     Range range(0, dst.rows);
4266     warpAffineInvoker invoker(src, dst, interpolation, borderType,
4267                               borderValue, adelta, bdelta, M);
4268     parallel_for_(range, invoker, dst.total()/(double)(1<<16));
4269 }
4270
4271
4272 namespace cv
4273 {
4274
4275 class warpPerspectiveInvoker :
4276     public ParallelLoopBody
4277 {
4278 public:
4279
4280     warpPerspectiveInvoker(const Mat &_src, Mat &_dst, double *_M, int _interpolation,
4281                            int _borderType, const Scalar &_borderValue) :
4282         ParallelLoopBody(), src(_src), dst(_dst), M(_M), interpolation(_interpolation),
4283         borderType(_borderType), borderValue(_borderValue)
4284     {
4285     }
4286
4287     virtual void operator() (const Range& range) const
4288     {
4289         const int BLOCK_SZ = 32;
4290         short XY[BLOCK_SZ*BLOCK_SZ*2], A[BLOCK_SZ*BLOCK_SZ];
4291         int x, y, x1, y1, width = dst.cols, height = dst.rows;
4292
4293         int bh0 = std::min(BLOCK_SZ/2, height);
4294         int bw0 = std::min(BLOCK_SZ*BLOCK_SZ/bh0, width);
4295         bh0 = std::min(BLOCK_SZ*BLOCK_SZ/bw0, height);
4296
4297         for( y = range.start; y < range.end; y += bh0 )
4298         {
4299             for( x = 0; x < width; x += bw0 )
4300             {
4301                 int bw = std::min( bw0, width - x);
4302                 int bh = std::min( bh0, range.end - y); // height
4303
4304                 Mat _XY(bh, bw, CV_16SC2, XY), matA;
4305                 Mat dpart(dst, Rect(x, y, bw, bh));
4306
4307                 for( y1 = 0; y1 < bh; y1++ )
4308                 {
4309                     short* xy = XY + y1*bw*2;
4310                     double X0 = M[0]*x + M[1]*(y + y1) + M[2];
4311                     double Y0 = M[3]*x + M[4]*(y + y1) + M[5];
4312                     double W0 = M[6]*x + M[7]*(y + y1) + M[8];
4313
4314                     if( interpolation == INTER_NEAREST )
4315                         for( x1 = 0; x1 < bw; x1++ )
4316                         {
4317                             double W = W0 + M[6]*x1;
4318                             W = W ? 1./W : 0;
4319                             double fX = std::max((double)INT_MIN, std::min((double)INT_MAX, (X0 + M[0]*x1)*W));
4320                             double fY = std::max((double)INT_MIN, std::min((double)INT_MAX, (Y0 + M[3]*x1)*W));
4321                             int X = saturate_cast<int>(fX);
4322                             int Y = saturate_cast<int>(fY);
4323
4324                             xy[x1*2] = saturate_cast<short>(X);
4325                             xy[x1*2+1] = saturate_cast<short>(Y);
4326                         }
4327                     else
4328                     {
4329                         short* alpha = A + y1*bw;
4330                         for( x1 = 0; x1 < bw; x1++ )
4331                         {
4332                             double W = W0 + M[6]*x1;
4333                             W = W ? INTER_TAB_SIZE/W : 0;
4334                             double fX = std::max((double)INT_MIN, std::min((double)INT_MAX, (X0 + M[0]*x1)*W));
4335                             double fY = std::max((double)INT_MIN, std::min((double)INT_MAX, (Y0 + M[3]*x1)*W));
4336                             int X = saturate_cast<int>(fX);
4337                             int Y = saturate_cast<int>(fY);
4338
4339                             xy[x1*2] = saturate_cast<short>(X >> INTER_BITS);
4340                             xy[x1*2+1] = saturate_cast<short>(Y >> INTER_BITS);
4341                             alpha[x1] = (short)((Y & (INTER_TAB_SIZE-1))*INTER_TAB_SIZE +
4342                                                 (X & (INTER_TAB_SIZE-1)));
4343                         }
4344                     }
4345                 }
4346
4347                 if( interpolation == INTER_NEAREST )
4348                     remap( src, dpart, _XY, Mat(), interpolation, borderType, borderValue );
4349                 else
4350                 {
4351                     Mat _matA(bh, bw, CV_16U, A);
4352                     remap( src, dpart, _XY, _matA, interpolation, borderType, borderValue );
4353                 }
4354             }
4355         }
4356     }
4357
4358 private:
4359     Mat src;
4360     Mat dst;
4361     double* M;
4362     int interpolation, borderType;
4363     Scalar borderValue;
4364 };
4365
4366 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
4367 class IPPwarpPerspectiveInvoker :
4368     public ParallelLoopBody
4369 {
4370 public:
4371     IPPwarpPerspectiveInvoker(Mat &_src, Mat &_dst, double (&_coeffs)[3][3], int &_interpolation, int &_borderType, const Scalar &_borderValue, ippiWarpPerspectiveBackFunc _func, bool *_ok) :
4372       ParallelLoopBody(), src(_src), dst(_dst), mode(_interpolation), coeffs(_coeffs), borderType(_borderType), borderValue(_borderValue), func(_func), ok(_ok)
4373       {
4374           *ok = true;
4375       }
4376
4377       virtual void operator() (const Range& range) const
4378       {
4379           IppiSize srcsize = {src.cols, src.rows};
4380           IppiRect srcroi = {0, 0, src.cols, src.rows};
4381           IppiRect dstroi = {0, range.start, dst.cols, range.end - range.start};
4382           int cnn = src.channels();
4383
4384           if( borderType == BORDER_CONSTANT )
4385           {
4386               IppiSize setSize = {dst.cols, range.end - range.start};
4387               void *dataPointer = dst.data + dst.step[0] * range.start;
4388               if( !IPPSet( borderValue, dataPointer, (int)dst.step[0], setSize, cnn, src.depth() ) )
4389               {
4390                   *ok = false;
4391                   return;
4392               }
4393           }
4394           if( func(src.data, srcsize, (int)src.step[0], srcroi, dst.data, (int)dst.step[0], dstroi, coeffs, mode) < 0)
4395               *ok = false;
4396       }
4397 private:
4398     Mat &src;
4399     Mat &dst;
4400     double (&coeffs)[3][3];
4401     int mode;
4402     int borderType;
4403     const Scalar borderValue;
4404     ippiWarpPerspectiveBackFunc func;
4405     bool *ok;
4406     const IPPwarpPerspectiveInvoker& operator= (const IPPwarpPerspectiveInvoker&);
4407 };
4408 #endif
4409
4410 }
4411
4412 void cv::warpPerspective( InputArray _src, OutputArray _dst, InputArray _M0,
4413                           Size dsize, int flags, int borderType, const Scalar& borderValue )
4414 {
4415     CV_Assert( _src.total() > 0 );
4416
4417     CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
4418                ocl_warpTransform(_src, _dst, _M0, dsize, flags, borderType, borderValue,
4419                               OCL_OP_PERSPECTIVE))
4420
4421     Mat src = _src.getMat(), M0 = _M0.getMat();
4422     _dst.create( dsize.area() == 0 ? src.size() : dsize, src.type() );
4423     Mat dst = _dst.getMat();
4424
4425     if( dst.data == src.data )
4426         src = src.clone();
4427
4428     double M[9];
4429     Mat matM(3, 3, CV_64F, M);
4430     int interpolation = flags & INTER_MAX;
4431     if( interpolation == INTER_AREA )
4432         interpolation = INTER_LINEAR;
4433
4434     CV_Assert( (M0.type() == CV_32F || M0.type() == CV_64F) && M0.rows == 3 && M0.cols == 3 );
4435     M0.convertTo(matM, matM.type());
4436
4437 #ifdef HAVE_TEGRA_OPTIMIZATION
4438     if( tegra::warpPerspective(src, dst, M, flags, borderType, borderValue) )
4439         return;
4440 #endif
4441
4442     if( !(flags & WARP_INVERSE_MAP) )
4443          invert(matM, matM);
4444 /*
4445 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
4446     int depth = src.depth();
4447     int channels = src.channels();
4448     if( ( depth == CV_8U || depth == CV_16U || depth == CV_32F ) &&
4449         ( channels == 1 || channels == 3 || channels == 4 ) &&
4450         ( borderType == cv::BORDER_TRANSPARENT || borderType == cv::BORDER_CONSTANT ) )
4451     {
4452         int type = src.type();
4453         ippiWarpPerspectiveBackFunc ippFunc =
4454             type == CV_8UC1 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_8u_C1R :
4455             type == CV_8UC3 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_8u_C3R :
4456             type == CV_8UC4 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_8u_C4R :
4457             type == CV_16UC1 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_16u_C1R :
4458             type == CV_16UC3 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_16u_C3R :
4459             type == CV_16UC4 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_16u_C4R :
4460             type == CV_32FC1 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_32f_C1R :
4461             type == CV_32FC3 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_32f_C3R :
4462             type == CV_32FC4 ? (ippiWarpPerspectiveBackFunc)ippiWarpPerspectiveBack_32f_C4R :
4463             0;
4464         int mode =
4465             flags == INTER_LINEAR ? IPPI_INTER_LINEAR :
4466             flags == INTER_NEAREST ? IPPI_INTER_NN :
4467             flags == INTER_CUBIC ? IPPI_INTER_CUBIC :
4468             0;
4469         if( mode && ippFunc )
4470         {
4471             double coeffs[3][3];
4472             for( int i = 0; i < 3; i++ )
4473             {
4474                 for( int j = 0; j < 3; j++ )
4475                 {
4476                     coeffs[i][j] = matM.at<double>(i, j);
4477                 }
4478             }
4479             bool ok;
4480             Range range(0, dst.rows);
4481             IPPwarpPerspectiveInvoker invoker(src, dst, coeffs, mode, borderType, borderValue, ippFunc, &ok);
4482             parallel_for_(range, invoker, dst.total()/(double)(1<<16));
4483             if( ok )
4484                 return;
4485         }
4486     }
4487 #endif
4488 */
4489     Range range(0, dst.rows);
4490     warpPerspectiveInvoker invoker(src, dst, M, interpolation, borderType, borderValue);
4491     parallel_for_(range, invoker, dst.total()/(double)(1<<16));
4492 }
4493
4494
4495 cv::Mat cv::getRotationMatrix2D( Point2f center, double angle, double scale )
4496 {
4497     angle *= CV_PI/180;
4498     double alpha = cos(angle)*scale;
4499     double beta = sin(angle)*scale;
4500
4501     Mat M(2, 3, CV_64F);
4502     double* m = (double*)M.data;
4503
4504     m[0] = alpha;
4505     m[1] = beta;
4506     m[2] = (1-alpha)*center.x - beta*center.y;
4507     m[3] = -beta;
4508     m[4] = alpha;
4509     m[5] = beta*center.x + (1-alpha)*center.y;
4510
4511     return M;
4512 }
4513
4514 /* Calculates coefficients of perspective transformation
4515  * which maps (xi,yi) to (ui,vi), (i=1,2,3,4):
4516  *
4517  *      c00*xi + c01*yi + c02
4518  * ui = ---------------------
4519  *      c20*xi + c21*yi + c22
4520  *
4521  *      c10*xi + c11*yi + c12
4522  * vi = ---------------------
4523  *      c20*xi + c21*yi + c22
4524  *
4525  * Coefficients are calculated by solving linear system:
4526  * / x0 y0  1  0  0  0 -x0*u0 -y0*u0 \ /c00\ /u0\
4527  * | x1 y1  1  0  0  0 -x1*u1 -y1*u1 | |c01| |u1|
4528  * | x2 y2  1  0  0  0 -x2*u2 -y2*u2 | |c02| |u2|
4529  * | x3 y3  1  0  0  0 -x3*u3 -y3*u3 |.|c10|=|u3|,
4530  * |  0  0  0 x0 y0  1 -x0*v0 -y0*v0 | |c11| |v0|
4531  * |  0  0  0 x1 y1  1 -x1*v1 -y1*v1 | |c12| |v1|
4532  * |  0  0  0 x2 y2  1 -x2*v2 -y2*v2 | |c20| |v2|
4533  * \  0  0  0 x3 y3  1 -x3*v3 -y3*v3 / \c21/ \v3/
4534  *
4535  * where:
4536  *   cij - matrix coefficients, c22 = 1
4537  */
4538 cv::Mat cv::getPerspectiveTransform( const Point2f src[], const Point2f dst[] )
4539 {
4540     Mat M(3, 3, CV_64F), X(8, 1, CV_64F, M.data);
4541     double a[8][8], b[8];
4542     Mat A(8, 8, CV_64F, a), B(8, 1, CV_64F, b);
4543
4544     for( int i = 0; i < 4; ++i )
4545     {
4546         a[i][0] = a[i+4][3] = src[i].x;
4547         a[i][1] = a[i+4][4] = src[i].y;
4548         a[i][2] = a[i+4][5] = 1;
4549         a[i][3] = a[i][4] = a[i][5] =
4550         a[i+4][0] = a[i+4][1] = a[i+4][2] = 0;
4551         a[i][6] = -src[i].x*dst[i].x;
4552         a[i][7] = -src[i].y*dst[i].x;
4553         a[i+4][6] = -src[i].x*dst[i].y;
4554         a[i+4][7] = -src[i].y*dst[i].y;
4555         b[i] = dst[i].x;
4556         b[i+4] = dst[i].y;
4557     }
4558
4559     solve( A, B, X, DECOMP_SVD );
4560     ((double*)M.data)[8] = 1.;
4561
4562     return M;
4563 }
4564
4565 /* Calculates coefficients of affine transformation
4566  * which maps (xi,yi) to (ui,vi), (i=1,2,3):
4567  *
4568  * ui = c00*xi + c01*yi + c02
4569  *
4570  * vi = c10*xi + c11*yi + c12
4571  *
4572  * Coefficients are calculated by solving linear system:
4573  * / x0 y0  1  0  0  0 \ /c00\ /u0\
4574  * | x1 y1  1  0  0  0 | |c01| |u1|
4575  * | x2 y2  1  0  0  0 | |c02| |u2|
4576  * |  0  0  0 x0 y0  1 | |c10| |v0|
4577  * |  0  0  0 x1 y1  1 | |c11| |v1|
4578  * \  0  0  0 x2 y2  1 / |c12| |v2|
4579  *
4580  * where:
4581  *   cij - matrix coefficients
4582  */
4583
4584 cv::Mat cv::getAffineTransform( const Point2f src[], const Point2f dst[] )
4585 {
4586     Mat M(2, 3, CV_64F), X(6, 1, CV_64F, M.data);
4587     double a[6*6], b[6];
4588     Mat A(6, 6, CV_64F, a), B(6, 1, CV_64F, b);
4589
4590     for( int i = 0; i < 3; i++ )
4591     {
4592         int j = i*12;
4593         int k = i*12+6;
4594         a[j] = a[k+3] = src[i].x;
4595         a[j+1] = a[k+4] = src[i].y;
4596         a[j+2] = a[k+5] = 1;
4597         a[j+3] = a[j+4] = a[j+5] = 0;
4598         a[k] = a[k+1] = a[k+2] = 0;
4599         b[i*2] = dst[i].x;
4600         b[i*2+1] = dst[i].y;
4601     }
4602
4603     solve( A, B, X );
4604     return M;
4605 }
4606
4607 void cv::invertAffineTransform(InputArray _matM, OutputArray __iM)
4608 {
4609     Mat matM = _matM.getMat();
4610     CV_Assert(matM.rows == 2 && matM.cols == 3);
4611     __iM.create(2, 3, matM.type());
4612     Mat _iM = __iM.getMat();
4613
4614     if( matM.type() == CV_32F )
4615     {
4616         const float* M = (const float*)matM.data;
4617         float* iM = (float*)_iM.data;
4618         int step = (int)(matM.step/sizeof(M[0])), istep = (int)(_iM.step/sizeof(iM[0]));
4619
4620         double D = M[0]*M[step+1] - M[1]*M[step];
4621         D = D != 0 ? 1./D : 0;
4622         double A11 = M[step+1]*D, A22 = M[0]*D, A12 = -M[1]*D, A21 = -M[step]*D;
4623         double b1 = -A11*M[2] - A12*M[step+2];
4624         double b2 = -A21*M[2] - A22*M[step+2];
4625
4626         iM[0] = (float)A11; iM[1] = (float)A12; iM[2] = (float)b1;
4627         iM[istep] = (float)A21; iM[istep+1] = (float)A22; iM[istep+2] = (float)b2;
4628     }
4629     else if( matM.type() == CV_64F )
4630     {
4631         const double* M = (const double*)matM.data;
4632         double* iM = (double*)_iM.data;
4633         int step = (int)(matM.step/sizeof(M[0])), istep = (int)(_iM.step/sizeof(iM[0]));
4634
4635         double D = M[0]*M[step+1] - M[1]*M[step];
4636         D = D != 0 ? 1./D : 0;
4637         double A11 = M[step+1]*D, A22 = M[0]*D, A12 = -M[1]*D, A21 = -M[step]*D;
4638         double b1 = -A11*M[2] - A12*M[step+2];
4639         double b2 = -A21*M[2] - A22*M[step+2];
4640
4641         iM[0] = A11; iM[1] = A12; iM[2] = b1;
4642         iM[istep] = A21; iM[istep+1] = A22; iM[istep+2] = b2;
4643     }
4644     else
4645         CV_Error( CV_StsUnsupportedFormat, "" );
4646 }
4647
4648 cv::Mat cv::getPerspectiveTransform(InputArray _src, InputArray _dst)
4649 {
4650     Mat src = _src.getMat(), dst = _dst.getMat();
4651     CV_Assert(src.checkVector(2, CV_32F) == 4 && dst.checkVector(2, CV_32F) == 4);
4652     return getPerspectiveTransform((const Point2f*)src.data, (const Point2f*)dst.data);
4653 }
4654
4655 cv::Mat cv::getAffineTransform(InputArray _src, InputArray _dst)
4656 {
4657     Mat src = _src.getMat(), dst = _dst.getMat();
4658     CV_Assert(src.checkVector(2, CV_32F) == 3 && dst.checkVector(2, CV_32F) == 3);
4659     return getAffineTransform((const Point2f*)src.data, (const Point2f*)dst.data);
4660 }
4661
4662 CV_IMPL void
4663 cvResize( const CvArr* srcarr, CvArr* dstarr, int method )
4664 {
4665     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
4666     CV_Assert( src.type() == dst.type() );
4667     cv::resize( src, dst, dst.size(), (double)dst.cols/src.cols,
4668         (double)dst.rows/src.rows, method );
4669 }
4670
4671
4672 CV_IMPL void
4673 cvWarpAffine( const CvArr* srcarr, CvArr* dstarr, const CvMat* marr,
4674               int flags, CvScalar fillval )
4675 {
4676     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
4677     cv::Mat matrix = cv::cvarrToMat(marr);
4678     CV_Assert( src.type() == dst.type() );
4679     cv::warpAffine( src, dst, matrix, dst.size(), flags,
4680         (flags & CV_WARP_FILL_OUTLIERS) ? cv::BORDER_CONSTANT : cv::BORDER_TRANSPARENT,
4681         fillval );
4682 }
4683
4684 CV_IMPL void
4685 cvWarpPerspective( const CvArr* srcarr, CvArr* dstarr, const CvMat* marr,
4686                    int flags, CvScalar fillval )
4687 {
4688     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
4689     cv::Mat matrix = cv::cvarrToMat(marr);
4690     CV_Assert( src.type() == dst.type() );
4691     cv::warpPerspective( src, dst, matrix, dst.size(), flags,
4692         (flags & CV_WARP_FILL_OUTLIERS) ? cv::BORDER_CONSTANT : cv::BORDER_TRANSPARENT,
4693         fillval );
4694 }
4695
4696 CV_IMPL void
4697 cvRemap( const CvArr* srcarr, CvArr* dstarr,
4698          const CvArr* _mapx, const CvArr* _mapy,
4699          int flags, CvScalar fillval )
4700 {
4701     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr), dst0 = dst;
4702     cv::Mat mapx = cv::cvarrToMat(_mapx), mapy = cv::cvarrToMat(_mapy);
4703     CV_Assert( src.type() == dst.type() && dst.size() == mapx.size() );
4704     cv::remap( src, dst, mapx, mapy, flags & cv::INTER_MAX,
4705         (flags & CV_WARP_FILL_OUTLIERS) ? cv::BORDER_CONSTANT : cv::BORDER_TRANSPARENT,
4706         fillval );
4707     CV_Assert( dst0.data == dst.data );
4708 }
4709
4710
4711 CV_IMPL CvMat*
4712 cv2DRotationMatrix( CvPoint2D32f center, double angle,
4713                     double scale, CvMat* matrix )
4714 {
4715     cv::Mat M0 = cv::cvarrToMat(matrix), M = cv::getRotationMatrix2D(center, angle, scale);
4716     CV_Assert( M.size() == M0.size() );
4717     M.convertTo(M0, M0.type());
4718     return matrix;
4719 }
4720
4721
4722 CV_IMPL CvMat*
4723 cvGetPerspectiveTransform( const CvPoint2D32f* src,
4724                           const CvPoint2D32f* dst,
4725                           CvMat* matrix )
4726 {
4727     cv::Mat M0 = cv::cvarrToMat(matrix),
4728         M = cv::getPerspectiveTransform((const cv::Point2f*)src, (const cv::Point2f*)dst);
4729     CV_Assert( M.size() == M0.size() );
4730     M.convertTo(M0, M0.type());
4731     return matrix;
4732 }
4733
4734
4735 CV_IMPL CvMat*
4736 cvGetAffineTransform( const CvPoint2D32f* src,
4737                           const CvPoint2D32f* dst,
4738                           CvMat* matrix )
4739 {
4740     cv::Mat M0 = cv::cvarrToMat(matrix),
4741         M = cv::getAffineTransform((const cv::Point2f*)src, (const cv::Point2f*)dst);
4742     CV_Assert( M.size() == M0.size() );
4743     M.convertTo(M0, M0.type());
4744     return matrix;
4745 }
4746
4747
4748 CV_IMPL void
4749 cvConvertMaps( const CvArr* arr1, const CvArr* arr2, CvArr* dstarr1, CvArr* dstarr2 )
4750 {
4751     cv::Mat map1 = cv::cvarrToMat(arr1), map2;
4752     cv::Mat dstmap1 = cv::cvarrToMat(dstarr1), dstmap2;
4753
4754     if( arr2 )
4755         map2 = cv::cvarrToMat(arr2);
4756     if( dstarr2 )
4757     {
4758         dstmap2 = cv::cvarrToMat(dstarr2);
4759         if( dstmap2.type() == CV_16SC1 )
4760             dstmap2 = cv::Mat(dstmap2.size(), CV_16UC1, dstmap2.data, dstmap2.step);
4761     }
4762
4763     cv::convertMaps( map1, map2, dstmap1, dstmap2, dstmap1.type(), false );
4764 }
4765
4766 /****************************************************************************************\
4767 *                                   Log-Polar Transform                                  *
4768 \****************************************************************************************/
4769
4770 /* now it is done via Remap; more correct implementation should use
4771    some super-sampling technique outside of the "fovea" circle */
4772 CV_IMPL void
4773 cvLogPolar( const CvArr* srcarr, CvArr* dstarr,
4774             CvPoint2D32f center, double M, int flags )
4775 {
4776     cv::Ptr<CvMat> mapx, mapy;
4777
4778     CvMat srcstub, *src = cvGetMat(srcarr, &srcstub);
4779     CvMat dststub, *dst = cvGetMat(dstarr, &dststub);
4780     CvSize ssize, dsize;
4781
4782     if( !CV_ARE_TYPES_EQ( src, dst ))
4783         CV_Error( CV_StsUnmatchedFormats, "" );
4784
4785     if( M <= 0 )
4786         CV_Error( CV_StsOutOfRange, "M should be >0" );
4787
4788     ssize = cvGetMatSize(src);
4789     dsize = cvGetMatSize(dst);
4790
4791     mapx.reset(cvCreateMat( dsize.height, dsize.width, CV_32F ));
4792     mapy.reset(cvCreateMat( dsize.height, dsize.width, CV_32F ));
4793
4794     if( !(flags & CV_WARP_INVERSE_MAP) )
4795     {
4796         int phi, rho;
4797         cv::AutoBuffer<double> _exp_tab(dsize.width);
4798         double* exp_tab = _exp_tab;
4799
4800         for( rho = 0; rho < dst->width; rho++ )
4801             exp_tab[rho] = std::exp(rho/M);
4802
4803         for( phi = 0; phi < dsize.height; phi++ )
4804         {
4805             double cp = cos(phi*2*CV_PI/dsize.height);
4806             double sp = sin(phi*2*CV_PI/dsize.height);
4807             float* mx = (float*)(mapx->data.ptr + phi*mapx->step);
4808             float* my = (float*)(mapy->data.ptr + phi*mapy->step);
4809
4810             for( rho = 0; rho < dsize.width; rho++ )
4811             {
4812                 double r = exp_tab[rho];
4813                 double x = r*cp + center.x;
4814                 double y = r*sp + center.y;
4815
4816                 mx[rho] = (float)x;
4817                 my[rho] = (float)y;
4818             }
4819         }
4820     }
4821     else
4822     {
4823         int x, y;
4824         CvMat bufx, bufy, bufp, bufa;
4825         double ascale = ssize.height/(2*CV_PI);
4826         cv::AutoBuffer<float> _buf(4*dsize.width);
4827         float* buf = _buf;
4828
4829         bufx = cvMat( 1, dsize.width, CV_32F, buf );
4830         bufy = cvMat( 1, dsize.width, CV_32F, buf + dsize.width );
4831         bufp = cvMat( 1, dsize.width, CV_32F, buf + dsize.width*2 );
4832         bufa = cvMat( 1, dsize.width, CV_32F, buf + dsize.width*3 );
4833
4834         for( x = 0; x < dsize.width; x++ )
4835             bufx.data.fl[x] = (float)x - center.x;
4836
4837         for( y = 0; y < dsize.height; y++ )
4838         {
4839             float* mx = (float*)(mapx->data.ptr + y*mapx->step);
4840             float* my = (float*)(mapy->data.ptr + y*mapy->step);
4841
4842             for( x = 0; x < dsize.width; x++ )
4843                 bufy.data.fl[x] = (float)y - center.y;
4844
4845 #if 1
4846             cvCartToPolar( &bufx, &bufy, &bufp, &bufa );
4847
4848             for( x = 0; x < dsize.width; x++ )
4849                 bufp.data.fl[x] += 1.f;
4850
4851             cvLog( &bufp, &bufp );
4852
4853             for( x = 0; x < dsize.width; x++ )
4854             {
4855                 double rho = bufp.data.fl[x]*M;
4856                 double phi = bufa.data.fl[x]*ascale;
4857
4858                 mx[x] = (float)rho;
4859                 my[x] = (float)phi;
4860             }
4861 #else
4862             for( x = 0; x < dsize.width; x++ )
4863             {
4864                 double xx = bufx.data.fl[x];
4865                 double yy = bufy.data.fl[x];
4866
4867                 double p = log(std::sqrt(xx*xx + yy*yy) + 1.)*M;
4868                 double a = atan2(yy,xx);
4869                 if( a < 0 )
4870                     a = 2*CV_PI + a;
4871                 a *= ascale;
4872
4873                 mx[x] = (float)p;
4874                 my[x] = (float)a;
4875             }
4876 #endif
4877         }
4878     }
4879
4880     cvRemap( src, dst, mapx, mapy, flags, cvScalarAll(0) );
4881 }
4882
4883 void cv::logPolar( InputArray _src, OutputArray _dst,
4884                    Point2f center, double M, int flags )
4885 {
4886     Mat src = _src.getMat();
4887     _dst.create( src.size(), src.type() );
4888     CvMat c_src = src, c_dst = _dst.getMat();
4889     cvLogPolar( &c_src, &c_dst, center, M, flags );
4890 }
4891
4892 /****************************************************************************************
4893                                    Linear-Polar Transform
4894   J.L. Blanco, Apr 2009
4895  ****************************************************************************************/
4896 CV_IMPL
4897 void cvLinearPolar( const CvArr* srcarr, CvArr* dstarr,
4898             CvPoint2D32f center, double maxRadius, int flags )
4899 {
4900     cv::Ptr<CvMat> mapx, mapy;
4901
4902     CvMat srcstub, *src = (CvMat*)srcarr;
4903     CvMat dststub, *dst = (CvMat*)dstarr;
4904     CvSize ssize, dsize;
4905
4906     src = cvGetMat( srcarr, &srcstub,0,0 );
4907     dst = cvGetMat( dstarr, &dststub,0,0 );
4908
4909     if( !CV_ARE_TYPES_EQ( src, dst ))
4910         CV_Error( CV_StsUnmatchedFormats, "" );
4911
4912     ssize.width = src->cols;
4913     ssize.height = src->rows;
4914     dsize.width = dst->cols;
4915     dsize.height = dst->rows;
4916
4917     mapx.reset(cvCreateMat( dsize.height, dsize.width, CV_32F ));
4918     mapy.reset(cvCreateMat( dsize.height, dsize.width, CV_32F ));
4919
4920     if( !(flags & CV_WARP_INVERSE_MAP) )
4921     {
4922         int phi, rho;
4923
4924         for( phi = 0; phi < dsize.height; phi++ )
4925         {
4926             double cp = cos(phi*2*CV_PI/dsize.height);
4927             double sp = sin(phi*2*CV_PI/dsize.height);
4928             float* mx = (float*)(mapx->data.ptr + phi*mapx->step);
4929             float* my = (float*)(mapy->data.ptr + phi*mapy->step);
4930
4931             for( rho = 0; rho < dsize.width; rho++ )
4932             {
4933                 double r = maxRadius*(rho+1)/dsize.width;
4934                 double x = r*cp + center.x;
4935                 double y = r*sp + center.y;
4936
4937                 mx[rho] = (float)x;
4938                 my[rho] = (float)y;
4939             }
4940         }
4941     }
4942     else
4943     {
4944         int x, y;
4945         CvMat bufx, bufy, bufp, bufa;
4946         const double ascale = ssize.height/(2*CV_PI);
4947         const double pscale = ssize.width/maxRadius;
4948
4949         cv::AutoBuffer<float> _buf(4*dsize.width);
4950         float* buf = _buf;
4951
4952         bufx = cvMat( 1, dsize.width, CV_32F, buf );
4953         bufy = cvMat( 1, dsize.width, CV_32F, buf + dsize.width );
4954         bufp = cvMat( 1, dsize.width, CV_32F, buf + dsize.width*2 );
4955         bufa = cvMat( 1, dsize.width, CV_32F, buf + dsize.width*3 );
4956
4957         for( x = 0; x < dsize.width; x++ )
4958             bufx.data.fl[x] = (float)x - center.x;
4959
4960         for( y = 0; y < dsize.height; y++ )
4961         {
4962             float* mx = (float*)(mapx->data.ptr + y*mapx->step);
4963             float* my = (float*)(mapy->data.ptr + y*mapy->step);
4964
4965             for( x = 0; x < dsize.width; x++ )
4966                 bufy.data.fl[x] = (float)y - center.y;
4967
4968             cvCartToPolar( &bufx, &bufy, &bufp, &bufa, 0 );
4969
4970             for( x = 0; x < dsize.width; x++ )
4971                 bufp.data.fl[x] += 1.f;
4972
4973             for( x = 0; x < dsize.width; x++ )
4974             {
4975                 double rho = bufp.data.fl[x]*pscale;
4976                 double phi = bufa.data.fl[x]*ascale;
4977                 mx[x] = (float)rho;
4978                 my[x] = (float)phi;
4979             }
4980         }
4981     }
4982
4983     cvRemap( src, dst, mapx, mapy, flags, cvScalarAll(0) );
4984 }
4985
4986 void cv::linearPolar( InputArray _src, OutputArray _dst,
4987                       Point2f center, double maxRadius, int flags )
4988 {
4989     Mat src = _src.getMat();
4990     _dst.create( src.size(), src.type() );
4991     CvMat c_src = src, c_dst = _dst.getMat();
4992     cvLinearPolar( &c_src, &c_dst, center, maxRadius, flags );
4993 }
4994
4995 /* End of file. */