modules/core/src/matrix.cpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 //
  10 //                           License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
  14 // Copyright (C) 2009, Willow Garage Inc., all rights reserved.
  15 // Third party copyrights are property of their respective owners.
  16 //
  17 // Redistribution and use in source and binary forms, with or without modification,
  18 // are permitted provided that the following conditions are met:
  19 //
  20 //   * Redistribution's of source code must retain the above copyright notice,
  21 //     this list of conditions and the following disclaimer.
  22 //
  23 //   * Redistribution's in binary form must reproduce the above copyright notice,
  24 //     this list of conditions and the following disclaimer in the documentation
  25 //     and/or other materials provided with the distribution.
  26 //
  27 //   * The name of the copyright holders may not be used to endorse or promote products
  28 //     derived from this software without specific prior written permission.
  29 //
  30 // This software is provided by the copyright holders and contributors "as is" and
  31 // any express or implied warranties, including, but not limited to, the implied
  32 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  33 // In no event shall the Intel Corporation or contributors be liable for any direct,
  34 // indirect, incidental, special, exemplary, or consequential damages
  35 // (including, but not limited to, procurement of substitute goods or services;
  36 // loss of use, data, or profits; or business interruption) however caused
  37 // and on any theory of liability, whether in contract, strict liability,
  38 // or tort (including negligence or otherwise) arising in any way out of
  39 // the use of this software, even if advised of the possibility of such damage.
  40 //
  41 //M*/
  42
  43 #include "precomp.hpp"
  44 #include "opencl_kernels_core.hpp"
  45
  46 #include "bufferpool.impl.hpp"
  47
  48 /****************************************************************************************\
  49 *                           [scaled] Identity matrix initialization                      *
  50 \****************************************************************************************/
  51
  52 namespace cv {
  53
  54 void MatAllocator::map(UMatData*, int) const
  55 {
  56 }
  57
  58 void MatAllocator::unmap(UMatData* u) const
  59 {
  60     if(u->urefcount == 0 && u->refcount == 0)
  61     {
  62         deallocate(u);
  63         u = NULL;
  64     }
  65 }
  66
  67 void MatAllocator::download(UMatData* u, void* dstptr,
  68          int dims, const size_t sz[],
  69          const size_t srcofs[], const size_t srcstep[],
  70          const size_t dststep[]) const
  71 {
  72     if(!u)
  73         return;
  74     int isz[CV_MAX_DIM];
  75     uchar* srcptr = u->data;
  76     for( int i = 0; i < dims; i++ )
  77     {
  78         CV_Assert( sz[i] <= (size_t)INT_MAX );
  79         if( sz[i] == 0 )
  80         return;
  81         if( srcofs )
  82         srcptr += srcofs[i]*(i <= dims-2 ? srcstep[i] : 1);
  83         isz[i] = (int)sz[i];
  84     }
  85
  86     Mat src(dims, isz, CV_8U, srcptr, srcstep);
  87     Mat dst(dims, isz, CV_8U, dstptr, dststep);
  88
  89     const Mat* arrays[] = { &src, &dst };
  90     uchar* ptrs[2];
  91     NAryMatIterator it(arrays, ptrs, 2);
  92     size_t j, planesz = it.size;
  93
  94     for( j = 0; j < it.nplanes; j++, ++it )
  95         memcpy(ptrs[1], ptrs[0], planesz);
  96 }
  97
  98
  99 void MatAllocator::upload(UMatData* u, const void* srcptr, int dims, const size_t sz[],
 100                     const size_t dstofs[], const size_t dststep[],
 101                     const size_t srcstep[]) const
 102 {
 103     if(!u)
 104         return;
 105     int isz[CV_MAX_DIM];
 106     uchar* dstptr = u->data;
 107     for( int i = 0; i < dims; i++ )
 108     {
 109         CV_Assert( sz[i] <= (size_t)INT_MAX );
 110         if( sz[i] == 0 )
 111         return;
 112         if( dstofs )
 113         dstptr += dstofs[i]*(i <= dims-2 ? dststep[i] : 1);
 114         isz[i] = (int)sz[i];
 115     }
 116
 117     Mat src(dims, isz, CV_8U, (void*)srcptr, srcstep);
 118     Mat dst(dims, isz, CV_8U, dstptr, dststep);
 119
 120     const Mat* arrays[] = { &src, &dst };
 121     uchar* ptrs[2];
 122     NAryMatIterator it(arrays, ptrs, 2);
 123     size_t j, planesz = it.size;
 124
 125     for( j = 0; j < it.nplanes; j++, ++it )
 126         memcpy(ptrs[1], ptrs[0], planesz);
 127 }
 128
 129 void MatAllocator::copy(UMatData* usrc, UMatData* udst, int dims, const size_t sz[],
 130                   const size_t srcofs[], const size_t srcstep[],
 131                   const size_t dstofs[], const size_t dststep[], bool /*sync*/) const
 132 {
 133     if(!usrc || !udst)
 134         return;
 135     int isz[CV_MAX_DIM];
 136     uchar* srcptr = usrc->data;
 137     uchar* dstptr = udst->data;
 138     for( int i = 0; i < dims; i++ )
 139     {
 140         CV_Assert( sz[i] <= (size_t)INT_MAX );
 141         if( sz[i] == 0 )
 142             return;
 143         if( srcofs )
 144             srcptr += srcofs[i]*(i <= dims-2 ? srcstep[i] : 1);
 145         if( dstofs )
 146             dstptr += dstofs[i]*(i <= dims-2 ? dststep[i] : 1);
 147         isz[i] = (int)sz[i];
 148     }
 149
 150     Mat src(dims, isz, CV_8U, srcptr, srcstep);
 151     Mat dst(dims, isz, CV_8U, dstptr, dststep);
 152
 153     const Mat* arrays[] = { &src, &dst };
 154     uchar* ptrs[2];
 155     NAryMatIterator it(arrays, ptrs, 2);
 156     size_t j, planesz = it.size;
 157
 158     for( j = 0; j < it.nplanes; j++, ++it )
 159         memcpy(ptrs[1], ptrs[0], planesz);
 160 }
 161
 162 BufferPoolController* MatAllocator::getBufferPoolController() const
 163 {
 164     static DummyBufferPoolController dummy;
 165     return &dummy;
 166 }
 167
 168 class StdMatAllocator : public MatAllocator
 169 {
 170 public:
 171     UMatData* allocate(int dims, const int* sizes, int type,
 172                        void* data0, size_t* step, int /*flags*/, UMatUsageFlags /*usageFlags*/) const
 173     {
 174         size_t total = CV_ELEM_SIZE(type);
 175         for( int i = dims-1; i >= 0; i-- )
 176         {
 177             if( step )
 178             {
 179                 if( data0 && step[i] != CV_AUTOSTEP )
 180                 {
 181                     CV_Assert(total <= step[i]);
 182                     total = step[i];
 183                 }
 184                 else
 185                     step[i] = total;
 186             }
 187             total *= sizes[i];
 188         }
 189         uchar* data = data0 ? (uchar*)data0 : (uchar*)fastMalloc(total);
 190         UMatData* u = new UMatData(this);
 191         u->data = u->origdata = data;
 192         u->size = total;
 193         if(data0)
 194             u->flags |= UMatData::USER_ALLOCATED;
 195
 196         return u;
 197     }
 198
 199     bool allocate(UMatData* u, int /*accessFlags*/, UMatUsageFlags /*usageFlags*/) const
 200     {
 201         if(!u) return false;
 202         return true;
 203     }
 204
 205     void deallocate(UMatData* u) const
 206     {
 207         CV_Assert(u->urefcount >= 0);
 208         CV_Assert(u->refcount >= 0);
 209         if(u && u->refcount == 0)
 210         {
 211             if( !(u->flags & UMatData::USER_ALLOCATED) )
 212             {
 213                 fastFree(u->origdata);
 214                 u->origdata = 0;
 215             }
 216             delete u;
 217         }
 218     }
 219 };
 220
 221
 222 MatAllocator* Mat::getStdAllocator()
 223 {
 224     static MatAllocator * allocator = new StdMatAllocator();
 225     return allocator;
 226 }
 227
 228 void swap( Mat& a, Mat& b )
 229 {
 230     std::swap(a.flags, b.flags);
 231     std::swap(a.dims, b.dims);
 232     std::swap(a.rows, b.rows);
 233     std::swap(a.cols, b.cols);
 234     std::swap(a.data, b.data);
 235     std::swap(a.datastart, b.datastart);
 236     std::swap(a.dataend, b.dataend);
 237     std::swap(a.datalimit, b.datalimit);
 238     std::swap(a.allocator, b.allocator);
 239     std::swap(a.u, b.u);
 240
 241     std::swap(a.size.p, b.size.p);
 242     std::swap(a.step.p, b.step.p);
 243     std::swap(a.step.buf[0], b.step.buf[0]);
 244     std::swap(a.step.buf[1], b.step.buf[1]);
 245
 246     if( a.step.p == b.step.buf )
 247     {
 248         a.step.p = a.step.buf;
 249         a.size.p = &a.rows;
 250     }
 251
 252     if( b.step.p == a.step.buf )
 253     {
 254         b.step.p = b.step.buf;
 255         b.size.p = &b.rows;
 256     }
 257 }
 258
 259
 260 static inline void setSize( Mat& m, int _dims, const int* _sz,
 261                             const size_t* _steps, bool autoSteps=false )
 262 {
 263     CV_Assert( 0 <= _dims && _dims <= CV_MAX_DIM );
 264     if( m.dims != _dims )
 265     {
 266         if( m.step.p != m.step.buf )
 267         {
 268             fastFree(m.step.p);
 269             m.step.p = m.step.buf;
 270             m.size.p = &m.rows;
 271         }
 272         if( _dims > 2 )
 273         {
 274             m.step.p = (size_t*)fastMalloc(_dims*sizeof(m.step.p[0]) + (_dims+1)*sizeof(m.size.p[0]));
 275             m.size.p = (int*)(m.step.p + _dims) + 1;
 276             m.size.p[-1] = _dims;
 277             m.rows = m.cols = -1;
 278         }
 279     }
 280
 281     m.dims = _dims;
 282     if( !_sz )
 283         return;
 284
 285     size_t esz = CV_ELEM_SIZE(m.flags), esz1 = CV_ELEM_SIZE1(m.flags), total = esz;
 286     int i;
 287     for( i = _dims-1; i >= 0; i-- )
 288     {
 289         int s = _sz[i];
 290         CV_Assert( s >= 0 );
 291         m.size.p[i] = s;
 292
 293         if( _steps )
 294         {
 295             if (_steps[i] % esz1 != 0)
 296             {
 297                 CV_Error(Error::BadStep, "Step must be a multiple of esz1");
 298             }
 299
 300             m.step.p[i] = i < _dims-1 ? _steps[i] : esz;
 301         }
 302         else if( autoSteps )
 303         {
 304             m.step.p[i] = total;
 305             int64 total1 = (int64)total*s;
 306             if( (uint64)total1 != (size_t)total1 )
 307                 CV_Error( CV_StsOutOfRange, "The total matrix size does not fit to \"size_t\" type" );
 308             total = (size_t)total1;
 309         }
 310     }
 311
 312     if( _dims == 1 )
 313     {
 314         m.dims = 2;
 315         m.cols = 1;
 316         m.step[1] = esz;
 317     }
 318 }
 319
 320 static void updateContinuityFlag(Mat& m)
 321 {
 322     int i, j;
 323     for( i = 0; i < m.dims; i++ )
 324     {
 325         if( m.size[i] > 1 )
 326             break;
 327     }
 328
 329     for( j = m.dims-1; j > i; j-- )
 330     {
 331         if( m.step[j]*m.size[j] < m.step[j-1] )
 332             break;
 333     }
 334
 335     uint64 t = (uint64)m.step[0]*m.size[0];
 336     if( j <= i && t == (size_t)t )
 337         m.flags |= Mat::CONTINUOUS_FLAG;
 338     else
 339         m.flags &= ~Mat::CONTINUOUS_FLAG;
 340 }
 341
 342 static void finalizeHdr(Mat& m)
 343 {
 344     updateContinuityFlag(m);
 345     int d = m.dims;
 346     if( d > 2 )
 347         m.rows = m.cols = -1;
 348     if(m.u)
 349         m.datastart = m.data = m.u->data;
 350     if( m.data )
 351     {
 352         m.datalimit = m.datastart + m.size[0]*m.step[0];
 353         if( m.size[0] > 0 )
 354         {
 355             m.dataend = m.data + m.size[d-1]*m.step[d-1];
 356             for( int i = 0; i < d-1; i++ )
 357                 m.dataend += (m.size[i] - 1)*m.step[i];
 358         }
 359         else
 360             m.dataend = m.datalimit;
 361     }
 362     else
 363         m.dataend = m.datalimit = 0;
 364 }
 365
 366
 367 void Mat::create(int d, const int* _sizes, int _type)
 368 {
 369     int i;
 370     CV_Assert(0 <= d && d <= CV_MAX_DIM && _sizes);
 371     _type = CV_MAT_TYPE(_type);
 372
 373     if( data && (d == dims || (d == 1 && dims <= 2)) && _type == type() )
 374     {
 375         if( d == 2 && rows == _sizes[0] && cols == _sizes[1] )
 376             return;
 377         for( i = 0; i < d; i++ )
 378             if( size[i] != _sizes[i] )
 379                 break;
 380         if( i == d && (d > 1 || size[1] == 1))
 381             return;
 382     }
 383
 384     release();
 385     if( d == 0 )
 386         return;
 387     flags = (_type & CV_MAT_TYPE_MASK) | MAGIC_VAL;
 388     setSize(*this, d, _sizes, 0, true);
 389
 390     if( total() > 0 )
 391     {
 392         MatAllocator *a = allocator, *a0 = getStdAllocator();
 393 #ifdef HAVE_TGPU
 394         if( !a || a == tegra::getAllocator() )
 395             a = tegra::getAllocator(d, _sizes, _type);
 396 #endif
 397         if(!a)
 398             a = a0;
 399         try
 400         {
 401             u = a->allocate(dims, size, _type, 0, step.p, 0, USAGE_DEFAULT);
 402             CV_Assert(u != 0);
 403         }
 404         catch(...)
 405         {
 406             if(a != a0)
 407                 u = a0->allocate(dims, size, _type, 0, step.p, 0, USAGE_DEFAULT);
 408             CV_Assert(u != 0);
 409         }
 410         CV_Assert( step[dims-1] == (size_t)CV_ELEM_SIZE(flags) );
 411     }
 412
 413     addref();
 414     finalizeHdr(*this);
 415 }
 416
 417 void Mat::copySize(const Mat& m)
 418 {
 419     setSize(*this, m.dims, 0, 0);
 420     for( int i = 0; i < dims; i++ )
 421     {
 422         size[i] = m.size[i];
 423         step[i] = m.step[i];
 424     }
 425 }
 426
 427 void Mat::deallocate()
 428 {
 429     if(u)
 430         (u->currAllocator ? u->currAllocator : allocator ? allocator : getStdAllocator())->unmap(u);
 431     u = NULL;
 432 }
 433
 434 Mat::Mat(const Mat& m, const Range& _rowRange, const Range& _colRange)
 435     : flags(MAGIC_VAL), dims(0), rows(0), cols(0), data(0), datastart(0), dataend(0),
 436       datalimit(0), allocator(0), u(0), size(&rows)
 437 {
 438     CV_Assert( m.dims >= 2 );
 439     if( m.dims > 2 )
 440     {
 441         AutoBuffer<Range> rs(m.dims);
 442         rs[0] = _rowRange;
 443         rs[1] = _colRange;
 444         for( int i = 2; i < m.dims; i++ )
 445             rs[i] = Range::all();
 446         *this = m(rs);
 447         return;
 448     }
 449
 450     *this = m;
 451     if( _rowRange != Range::all() && _rowRange != Range(0,rows) )
 452     {
 453         CV_Assert( 0 <= _rowRange.start && _rowRange.start <= _rowRange.end && _rowRange.end <= m.rows );
 454         rows = _rowRange.size();
 455         data += step*_rowRange.start;
 456         flags |= SUBMATRIX_FLAG;
 457     }
 458
 459     if( _colRange != Range::all() && _colRange != Range(0,cols) )
 460     {
 461         CV_Assert( 0 <= _colRange.start && _colRange.start <= _colRange.end && _colRange.end <= m.cols );
 462         cols = _colRange.size();
 463         data += _colRange.start*elemSize();
 464         flags &= cols < m.cols ? ~CONTINUOUS_FLAG : -1;
 465         flags |= SUBMATRIX_FLAG;
 466     }
 467
 468     if( rows == 1 )
 469         flags |= CONTINUOUS_FLAG;
 470
 471     if( rows <= 0 || cols <= 0 )
 472     {
 473         release();
 474         rows = cols = 0;
 475     }
 476 }
 477
 478
 479 Mat::Mat(const Mat& m, const Rect& roi)
 480     : flags(m.flags), dims(2), rows(roi.height), cols(roi.width),
 481     data(m.data + roi.y*m.step[0]),
 482     datastart(m.datastart), dataend(m.dataend), datalimit(m.datalimit),
 483     allocator(m.allocator), u(m.u), size(&rows)
 484 {
 485     CV_Assert( m.dims <= 2 );
 486     flags &= roi.width < m.cols ? ~CONTINUOUS_FLAG : -1;
 487     flags |= roi.height == 1 ? CONTINUOUS_FLAG : 0;
 488
 489     size_t esz = CV_ELEM_SIZE(flags);
 490     data += roi.x*esz;
 491     CV_Assert( 0 <= roi.x && 0 <= roi.width && roi.x + roi.width <= m.cols &&
 492               0 <= roi.y && 0 <= roi.height && roi.y + roi.height <= m.rows );
 493     if( u )
 494         CV_XADD(&u->refcount, 1);
 495     if( roi.width < m.cols || roi.height < m.rows )
 496         flags |= SUBMATRIX_FLAG;
 497
 498     step[0] = m.step[0]; step[1] = esz;
 499
 500     if( rows <= 0 || cols <= 0 )
 501     {
 502         release();
 503         rows = cols = 0;
 504     }
 505 }
 506
 507
 508 Mat::Mat(int _dims, const int* _sizes, int _type, void* _data, const size_t* _steps)
 509     : flags(MAGIC_VAL), dims(0), rows(0), cols(0), data(0), datastart(0), dataend(0),
 510       datalimit(0), allocator(0), u(0), size(&rows)
 511 {
 512     flags |= CV_MAT_TYPE(_type);
 513     datastart = data = (uchar*)_data;
 514     setSize(*this, _dims, _sizes, _steps, true);
 515     finalizeHdr(*this);
 516 }
 517
 518
 519 Mat::Mat(const Mat& m, const Range* ranges)
 520     : flags(MAGIC_VAL), dims(0), rows(0), cols(0), data(0), datastart(0), dataend(0),
 521       datalimit(0), allocator(0), u(0), size(&rows)
 522 {
 523     int i, d = m.dims;
 524
 525     CV_Assert(ranges);
 526     for( i = 0; i < d; i++ )
 527     {
 528         Range r = ranges[i];
 529         CV_Assert( r == Range::all() || (0 <= r.start && r.start < r.end && r.end <= m.size[i]) );
 530     }
 531     *this = m;
 532     for( i = 0; i < d; i++ )
 533     {
 534         Range r = ranges[i];
 535         if( r != Range::all() && r != Range(0, size.p[i]))
 536         {
 537             size.p[i] = r.end - r.start;
 538             data += r.start*step.p[i];
 539             flags |= SUBMATRIX_FLAG;
 540         }
 541     }
 542     updateContinuityFlag(*this);
 543 }
 544
 545
 546 static Mat cvMatNDToMat(const CvMatND* m, bool copyData)
 547 {
 548     Mat thiz;
 549
 550     if( !m )
 551         return thiz;
 552     thiz.datastart = thiz.data = m->data.ptr;
 553     thiz.flags |= CV_MAT_TYPE(m->type);
 554     int _sizes[CV_MAX_DIM];
 555     size_t _steps[CV_MAX_DIM];
 556
 557     int i, d = m->dims;
 558     for( i = 0; i < d; i++ )
 559     {
 560         _sizes[i] = m->dim[i].size;
 561         _steps[i] = m->dim[i].step;
 562     }
 563
 564     setSize(thiz, d, _sizes, _steps);
 565     finalizeHdr(thiz);
 566
 567     if( copyData )
 568     {
 569         Mat temp(thiz);
 570         thiz.release();
 571         temp.copyTo(thiz);
 572     }
 573
 574     return thiz;
 575 }
 576
 577 static Mat cvMatToMat(const CvMat* m, bool copyData)
 578 {
 579     Mat thiz;
 580
 581     if( !m )
 582         return thiz;
 583
 584     if( !copyData )
 585     {
 586         thiz.flags = Mat::MAGIC_VAL + (m->type & (CV_MAT_TYPE_MASK|CV_MAT_CONT_FLAG));
 587         thiz.dims = 2;
 588         thiz.rows = m->rows;
 589         thiz.cols = m->cols;
 590         thiz.datastart = thiz.data = m->data.ptr;
 591         size_t esz = CV_ELEM_SIZE(m->type), minstep = thiz.cols*esz, _step = m->step;
 592         if( _step == 0 )
 593             _step = minstep;
 594         thiz.datalimit = thiz.datastart + _step*thiz.rows;
 595         thiz.dataend = thiz.datalimit - _step + minstep;
 596         thiz.step[0] = _step; thiz.step[1] = esz;
 597     }
 598     else
 599     {
 600         thiz.datastart = thiz.dataend = thiz.data = 0;
 601         Mat(m->rows, m->cols, m->type, m->data.ptr, m->step).copyTo(thiz);
 602     }
 603
 604     return thiz;
 605 }
 606
 607
 608 static Mat iplImageToMat(const IplImage* img, bool copyData)
 609 {
 610     Mat m;
 611
 612     if( !img )
 613         return m;
 614
 615     m.dims = 2;
 616     CV_DbgAssert(CV_IS_IMAGE(img) && img->imageData != 0);
 617
 618     int imgdepth = IPL2CV_DEPTH(img->depth);
 619     size_t esz;
 620     m.step[0] = img->widthStep;
 621
 622     if(!img->roi)
 623     {
 624         CV_Assert(img->dataOrder == IPL_DATA_ORDER_PIXEL);
 625         m.flags = Mat::MAGIC_VAL + CV_MAKETYPE(imgdepth, img->nChannels);
 626         m.rows = img->height;
 627         m.cols = img->width;
 628         m.datastart = m.data = (uchar*)img->imageData;
 629         esz = CV_ELEM_SIZE(m.flags);
 630     }
 631     else
 632     {
 633         CV_Assert(img->dataOrder == IPL_DATA_ORDER_PIXEL || img->roi->coi != 0);
 634         bool selectedPlane = img->roi->coi && img->dataOrder == IPL_DATA_ORDER_PLANE;
 635         m.flags = Mat::MAGIC_VAL + CV_MAKETYPE(imgdepth, selectedPlane ? 1 : img->nChannels);
 636         m.rows = img->roi->height;
 637         m.cols = img->roi->width;
 638         esz = CV_ELEM_SIZE(m.flags);
 639         m.datastart = m.data = (uchar*)img->imageData +
 640             (selectedPlane ? (img->roi->coi - 1)*m.step*img->height : 0) +
 641             img->roi->yOffset*m.step[0] + img->roi->xOffset*esz;
 642     }
 643     m.datalimit = m.datastart + m.step.p[0]*m.rows;
 644     m.dataend = m.datastart + m.step.p[0]*(m.rows-1) + esz*m.cols;
 645     m.flags |= (m.cols*esz == m.step.p[0] || m.rows == 1 ? Mat::CONTINUOUS_FLAG : 0);
 646     m.step[1] = esz;
 647
 648     if( copyData )
 649     {
 650         Mat m2 = m;
 651         m.release();
 652         if( !img->roi || !img->roi->coi ||
 653             img->dataOrder == IPL_DATA_ORDER_PLANE)
 654             m2.copyTo(m);
 655         else
 656         {
 657             int ch[] = {img->roi->coi - 1, 0};
 658             m.create(m2.rows, m2.cols, m2.type());
 659             mixChannels(&m2, 1, &m, 1, ch, 1);
 660         }
 661     }
 662
 663     return m;
 664 }
 665
 666 Mat Mat::diag(int d) const
 667 {
 668     CV_Assert( dims <= 2 );
 669     Mat m = *this;
 670     size_t esz = elemSize();
 671     int len;
 672
 673     if( d >= 0 )
 674     {
 675         len = std::min(cols - d, rows);
 676         m.data += esz*d;
 677     }
 678     else
 679     {
 680         len = std::min(rows + d, cols);
 681         m.data -= step[0]*d;
 682     }
 683     CV_DbgAssert( len > 0 );
 684
 685     m.size[0] = m.rows = len;
 686     m.size[1] = m.cols = 1;
 687     m.step[0] += (len > 1 ? esz : 0);
 688
 689     if( m.rows > 1 )
 690         m.flags &= ~CONTINUOUS_FLAG;
 691     else
 692         m.flags |= CONTINUOUS_FLAG;
 693
 694     if( size() != Size(1,1) )
 695         m.flags |= SUBMATRIX_FLAG;
 696
 697     return m;
 698 }
 699
 700 void Mat::pop_back(size_t nelems)
 701 {
 702     CV_Assert( nelems <= (size_t)size.p[0] );
 703
 704     if( isSubmatrix() )
 705         *this = rowRange(0, size.p[0] - (int)nelems);
 706     else
 707     {
 708         size.p[0] -= (int)nelems;
 709         dataend -= nelems*step.p[0];
 710         /*if( size.p[0] <= 1 )
 711         {
 712             if( dims <= 2 )
 713                 flags |= CONTINUOUS_FLAG;
 714             else
 715                 updateContinuityFlag(*this);
 716         }*/
 717     }
 718 }
 719
 720
 721 void Mat::push_back_(const void* elem)
 722 {
 723     int r = size.p[0];
 724     if( isSubmatrix() || dataend + step.p[0] > datalimit )
 725         reserve( std::max(r + 1, (r*3+1)/2) );
 726
 727     size_t esz = elemSize();
 728     memcpy(data + r*step.p[0], elem, esz);
 729     size.p[0] = r + 1;
 730     dataend += step.p[0];
 731     if( esz < step.p[0] )
 732         flags &= ~CONTINUOUS_FLAG;
 733 }
 734
 735 void Mat::reserve(size_t nelems)
 736 {
 737     const size_t MIN_SIZE = 64;
 738
 739     CV_Assert( (int)nelems >= 0 );
 740     if( !isSubmatrix() && data + step.p[0]*nelems <= datalimit )
 741         return;
 742
 743     int r = size.p[0];
 744
 745     if( (size_t)r >= nelems )
 746         return;
 747
 748     size.p[0] = std::max((int)nelems, 1);
 749     size_t newsize = total()*elemSize();
 750
 751     if( newsize < MIN_SIZE )
 752         size.p[0] = (int)((MIN_SIZE + newsize - 1)*nelems/newsize);
 753
 754     Mat m(dims, size.p, type());
 755     size.p[0] = r;
 756     if( r > 0 )
 757     {
 758         Mat mpart = m.rowRange(0, r);
 759         copyTo(mpart);
 760     }
 761
 762     *this = m;
 763     size.p[0] = r;
 764     dataend = data + step.p[0]*r;
 765 }
 766
 767
 768 void Mat::resize(size_t nelems)
 769 {
 770     int saveRows = size.p[0];
 771     if( saveRows == (int)nelems )
 772         return;
 773     CV_Assert( (int)nelems >= 0 );
 774
 775     if( isSubmatrix() || data + step.p[0]*nelems > datalimit )
 776         reserve(nelems);
 777
 778     size.p[0] = (int)nelems;
 779     dataend += (size.p[0] - saveRows)*step.p[0];
 780
 781     //updateContinuityFlag(*this);
 782 }
 783
 784
 785 void Mat::resize(size_t nelems, const Scalar& s)
 786 {
 787     int saveRows = size.p[0];
 788     resize(nelems);
 789
 790     if( size.p[0] > saveRows )
 791     {
 792         Mat part = rowRange(saveRows, size.p[0]);
 793         part = s;
 794     }
 795 }
 796
 797 void Mat::push_back(const Mat& elems)
 798 {
 799     int r = size.p[0], delta = elems.size.p[0];
 800     if( delta == 0 )
 801         return;
 802     if( this == &elems )
 803     {
 804         Mat tmp = elems;
 805         push_back(tmp);
 806         return;
 807     }
 808     if( !data )
 809     {
 810         *this = elems.clone();
 811         return;
 812     }
 813
 814     size.p[0] = elems.size.p[0];
 815     bool eq = size == elems.size;
 816     size.p[0] = r;
 817     if( !eq )
 818         CV_Error(CV_StsUnmatchedSizes, "");
 819     if( type() != elems.type() )
 820         CV_Error(CV_StsUnmatchedFormats, "");
 821
 822     if( isSubmatrix() || dataend + step.p[0]*delta > datalimit )
 823         reserve( std::max(r + delta, (r*3+1)/2) );
 824
 825     size.p[0] += delta;
 826     dataend += step.p[0]*delta;
 827
 828     //updateContinuityFlag(*this);
 829
 830     if( isContinuous() && elems.isContinuous() )
 831         memcpy(data + r*step.p[0], elems.data, elems.total()*elems.elemSize());
 832     else
 833     {
 834         Mat part = rowRange(r, r + delta);
 835         elems.copyTo(part);
 836     }
 837 }
 838
 839
 840 Mat cvarrToMat(const CvArr* arr, bool copyData,
 841                bool /*allowND*/, int coiMode, AutoBuffer<double>* abuf )
 842 {
 843     if( !arr )
 844         return Mat();
 845     if( CV_IS_MAT_HDR_Z(arr) )
 846         return cvMatToMat((const CvMat*)arr, copyData);
 847     if( CV_IS_MATND(arr) )
 848         return cvMatNDToMat((const CvMatND*)arr, copyData );
 849     if( CV_IS_IMAGE(arr) )
 850     {
 851         const IplImage* iplimg = (const IplImage*)arr;
 852         if( coiMode == 0 && iplimg->roi && iplimg->roi->coi > 0 )
 853             CV_Error(CV_BadCOI, "COI is not supported by the function");
 854         return iplImageToMat(iplimg, copyData);
 855     }
 856     if( CV_IS_SEQ(arr) )
 857     {
 858         CvSeq* seq = (CvSeq*)arr;
 859         int total = seq->total, type = CV_MAT_TYPE(seq->flags), esz = seq->elem_size;
 860         if( total == 0 )
 861             return Mat();
 862         CV_Assert(total > 0 && CV_ELEM_SIZE(seq->flags) == esz);
 863         if(!copyData && seq->first->next == seq->first)
 864             return Mat(total, 1, type, seq->first->data);
 865         if( abuf )
 866         {
 867             abuf->allocate(((size_t)total*esz + sizeof(double)-1)/sizeof(double));
 868             double* bufdata = *abuf;
 869             cvCvtSeqToArray(seq, bufdata, CV_WHOLE_SEQ);
 870             return Mat(total, 1, type, bufdata);
 871         }
 872
 873         Mat buf(total, 1, type);
 874         cvCvtSeqToArray(seq, buf.data, CV_WHOLE_SEQ);
 875         return buf;
 876     }
 877     CV_Error(CV_StsBadArg, "Unknown array type");
 878     return Mat();
 879 }
 880
 881 void Mat::locateROI( Size& wholeSize, Point& ofs ) const
 882 {
 883     CV_Assert( dims <= 2 && step[0] > 0 );
 884     size_t esz = elemSize(), minstep;
 885     ptrdiff_t delta1 = data - datastart, delta2 = dataend - datastart;
 886
 887     if( delta1 == 0 )
 888         ofs.x = ofs.y = 0;
 889     else
 890     {
 891         ofs.y = (int)(delta1/step[0]);
 892         ofs.x = (int)((delta1 - step[0]*ofs.y)/esz);
 893         CV_DbgAssert( data == datastart + ofs.y*step[0] + ofs.x*esz );
 894     }
 895     minstep = (ofs.x + cols)*esz;
 896     wholeSize.height = (int)((delta2 - minstep)/step[0] + 1);
 897     wholeSize.height = std::max(wholeSize.height, ofs.y + rows);
 898     wholeSize.width = (int)((delta2 - step*(wholeSize.height-1))/esz);
 899     wholeSize.width = std::max(wholeSize.width, ofs.x + cols);
 900 }
 901
 902 Mat& Mat::adjustROI( int dtop, int dbottom, int dleft, int dright )
 903 {
 904     CV_Assert( dims <= 2 && step[0] > 0 );
 905     Size wholeSize; Point ofs;
 906     size_t esz = elemSize();
 907     locateROI( wholeSize, ofs );
 908     int row1 = std::max(ofs.y - dtop, 0), row2 = std::min(ofs.y + rows + dbottom, wholeSize.height);
 909     int col1 = std::max(ofs.x - dleft, 0), col2 = std::min(ofs.x + cols + dright, wholeSize.width);
 910     data += (row1 - ofs.y)*step + (col1 - ofs.x)*esz;
 911     rows = row2 - row1; cols = col2 - col1;
 912     size.p[0] = rows; size.p[1] = cols;
 913     if( esz*cols == step[0] || rows == 1 )
 914         flags |= CONTINUOUS_FLAG;
 915     else
 916         flags &= ~CONTINUOUS_FLAG;
 917     return *this;
 918 }
 919
 920 }
 921
 922 void cv::extractImageCOI(const CvArr* arr, OutputArray _ch, int coi)
 923 {
 924     Mat mat = cvarrToMat(arr, false, true, 1);
 925     _ch.create(mat.dims, mat.size, mat.depth());
 926     Mat ch = _ch.getMat();
 927     if(coi < 0)
 928     {
 929         CV_Assert( CV_IS_IMAGE(arr) );
 930         coi = cvGetImageCOI((const IplImage*)arr)-1;
 931     }
 932     CV_Assert(0 <= coi && coi < mat.channels());
 933     int _pairs[] = { coi, 0 };
 934     mixChannels( &mat, 1, &ch, 1, _pairs, 1 );
 935 }
 936
 937 void cv::insertImageCOI(InputArray _ch, CvArr* arr, int coi)
 938 {
 939     Mat ch = _ch.getMat(), mat = cvarrToMat(arr, false, true, 1);
 940     if(coi < 0)
 941     {
 942         CV_Assert( CV_IS_IMAGE(arr) );
 943         coi = cvGetImageCOI((const IplImage*)arr)-1;
 944     }
 945     CV_Assert(ch.size == mat.size && ch.depth() == mat.depth() && 0 <= coi && coi < mat.channels());
 946     int _pairs[] = { 0, coi };
 947     mixChannels( &ch, 1, &mat, 1, _pairs, 1 );
 948 }
 949
 950 namespace cv
 951 {
 952
 953 Mat Mat::reshape(int new_cn, int new_rows) const
 954 {
 955     int cn = channels();
 956     Mat hdr = *this;
 957
 958     if( dims > 2 && new_rows == 0 && new_cn != 0 && size[dims-1]*cn % new_cn == 0 )
 959     {
 960         hdr.flags = (hdr.flags & ~CV_MAT_CN_MASK) | ((new_cn-1) << CV_CN_SHIFT);
 961         hdr.step[dims-1] = CV_ELEM_SIZE(hdr.flags);
 962         hdr.size[dims-1] = hdr.size[dims-1]*cn / new_cn;
 963         return hdr;
 964     }
 965
 966     CV_Assert( dims <= 2 );
 967
 968     if( new_cn == 0 )
 969         new_cn = cn;
 970
 971     int total_width = cols * cn;
 972
 973     if( (new_cn > total_width || total_width % new_cn != 0) && new_rows == 0 )
 974         new_rows = rows * total_width / new_cn;
 975
 976     if( new_rows != 0 && new_rows != rows )
 977     {
 978         int total_size = total_width * rows;
 979         if( !isContinuous() )
 980             CV_Error( CV_BadStep,
 981             "The matrix is not continuous, thus its number of rows can not be changed" );
 982
 983         if( (unsigned)new_rows > (unsigned)total_size )
 984             CV_Error( CV_StsOutOfRange, "Bad new number of rows" );
 985
 986         total_width = total_size / new_rows;
 987
 988         if( total_width * new_rows != total_size )
 989             CV_Error( CV_StsBadArg, "The total number of matrix elements "
 990                                     "is not divisible by the new number of rows" );
 991
 992         hdr.rows = new_rows;
 993         hdr.step[0] = total_width * elemSize1();
 994     }
 995
 996     int new_width = total_width / new_cn;
 997
 998     if( new_width * new_cn != total_width )
 999         CV_Error( CV_BadNumChannels,
1000         "The total width is not divisible by the new number of channels" );
1001
1002     hdr.cols = new_width;
1003     hdr.flags = (hdr.flags & ~CV_MAT_CN_MASK) | ((new_cn-1) << CV_CN_SHIFT);
1004     hdr.step[1] = CV_ELEM_SIZE(hdr.flags);
1005     return hdr;
1006 }
1007
1008 Mat Mat::diag(const Mat& d)
1009 {
1010     CV_Assert( d.cols == 1 || d.rows == 1 );
1011     int len = d.rows + d.cols - 1;
1012     Mat m(len, len, d.type(), Scalar(0));
1013     Mat md = m.diag();
1014     if( d.cols == 1 )
1015         d.copyTo(md);
1016     else
1017         transpose(d, md);
1018     return m;
1019 }
1020
1021 int Mat::checkVector(int _elemChannels, int _depth, bool _requireContinuous) const
1022 {
1023     return (depth() == _depth || _depth <= 0) &&
1024         (isContinuous() || !_requireContinuous) &&
1025         ((dims == 2 && (((rows == 1 || cols == 1) && channels() == _elemChannels) ||
1026                         (cols == _elemChannels && channels() == 1))) ||
1027         (dims == 3 && channels() == 1 && size.p[2] == _elemChannels && (size.p[0] == 1 || size.p[1] == 1) &&
1028          (isContinuous() || step.p[1] == step.p[2]*size.p[2])))
1029     ? (int)(total()*channels()/_elemChannels) : -1;
1030 }
1031
1032
1033 void scalarToRawData(const Scalar& s, void* _buf, int type, int unroll_to)
1034 {
1035     int i, depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
1036     CV_Assert(cn <= 4);
1037     switch(depth)
1038     {
1039     case CV_8U:
1040         {
1041         uchar* buf = (uchar*)_buf;
1042         for(i = 0; i < cn; i++)
1043             buf[i] = saturate_cast<uchar>(s.val[i]);
1044         for(; i < unroll_to; i++)
1045             buf[i] = buf[i-cn];
1046         }
1047         break;
1048     case CV_8S:
1049         {
1050         schar* buf = (schar*)_buf;
1051         for(i = 0; i < cn; i++)
1052             buf[i] = saturate_cast<schar>(s.val[i]);
1053         for(; i < unroll_to; i++)
1054             buf[i] = buf[i-cn];
1055         }
1056         break;
1057     case CV_16U:
1058         {
1059         ushort* buf = (ushort*)_buf;
1060         for(i = 0; i < cn; i++)
1061             buf[i] = saturate_cast<ushort>(s.val[i]);
1062         for(; i < unroll_to; i++)
1063             buf[i] = buf[i-cn];
1064         }
1065         break;
1066     case CV_16S:
1067         {
1068         short* buf = (short*)_buf;
1069         for(i = 0; i < cn; i++)
1070             buf[i] = saturate_cast<short>(s.val[i]);
1071         for(; i < unroll_to; i++)
1072             buf[i] = buf[i-cn];
1073         }
1074         break;
1075     case CV_32S:
1076         {
1077         int* buf = (int*)_buf;
1078         for(i = 0; i < cn; i++)
1079             buf[i] = saturate_cast<int>(s.val[i]);
1080         for(; i < unroll_to; i++)
1081             buf[i] = buf[i-cn];
1082         }
1083         break;
1084     case CV_32F:
1085         {
1086         float* buf = (float*)_buf;
1087         for(i = 0; i < cn; i++)
1088             buf[i] = saturate_cast<float>(s.val[i]);
1089         for(; i < unroll_to; i++)
1090             buf[i] = buf[i-cn];
1091         }
1092         break;
1093     case CV_64F:
1094         {
1095         double* buf = (double*)_buf;
1096         for(i = 0; i < cn; i++)
1097             buf[i] = saturate_cast<double>(s.val[i]);
1098         for(; i < unroll_to; i++)
1099             buf[i] = buf[i-cn];
1100         break;
1101         }
1102     default:
1103         CV_Error(CV_StsUnsupportedFormat,"");
1104     }
1105 }
1106
1107
1108 /*************************************************************************************************\
1109                                         Input/Output Array
1110 \*************************************************************************************************/
1111
1112 Mat _InputArray::getMat(int i) const
1113 {
1114     int k = kind();
1115     int accessFlags = flags & ACCESS_MASK;
1116
1117     if( k == MAT )
1118     {
1119         const Mat* m = (const Mat*)obj;
1120         if( i < 0 )
1121             return *m;
1122         return m->row(i);
1123     }
1124
1125     if( k == UMAT )
1126     {
1127         const UMat* m = (const UMat*)obj;
1128         if( i < 0 )
1129             return m->getMat(accessFlags);
1130         return m->getMat(accessFlags).row(i);
1131     }
1132
1133     if( k == EXPR )
1134     {
1135         CV_Assert( i < 0 );
1136         return (Mat)*((const MatExpr*)obj);
1137     }
1138
1139     if( k == MATX )
1140     {
1141         CV_Assert( i < 0 );
1142         return Mat(sz, flags, obj);
1143     }
1144
1145     if( k == STD_VECTOR )
1146     {
1147         CV_Assert( i < 0 );
1148         int t = CV_MAT_TYPE(flags);
1149         const std::vector<uchar>& v = *(const std::vector<uchar>*)obj;
1150
1151         return !v.empty() ? Mat(size(), t, (void*)&v[0]) : Mat();
1152     }
1153
1154     if( k == NONE )
1155         return Mat();
1156
1157     if( k == STD_VECTOR_VECTOR )
1158     {
1159         int t = type(i);
1160         const std::vector<std::vector<uchar> >& vv = *(const std::vector<std::vector<uchar> >*)obj;
1161         CV_Assert( 0 <= i && i < (int)vv.size() );
1162         const std::vector<uchar>& v = vv[i];
1163
1164         return !v.empty() ? Mat(size(i), t, (void*)&v[0]) : Mat();
1165     }
1166
1167     if( k == STD_VECTOR_MAT )
1168     {
1169         const std::vector<Mat>& v = *(const std::vector<Mat>*)obj;
1170         CV_Assert( 0 <= i && i < (int)v.size() );
1171
1172         return v[i];
1173     }
1174
1175     if( k == STD_VECTOR_UMAT )
1176     {
1177         const std::vector<UMat>& v = *(const std::vector<UMat>*)obj;
1178         CV_Assert( 0 <= i && i < (int)v.size() );
1179
1180         return v[i].getMat(accessFlags);
1181     }
1182
1183     if( k == OPENGL_BUFFER )
1184     {
1185         CV_Assert( i < 0 );
1186         CV_Error(cv::Error::StsNotImplemented, "You should explicitly call mapHost/unmapHost methods for ogl::Buffer object");
1187         return Mat();
1188     }
1189
1190     if( k == GPU_MAT )
1191     {
1192         CV_Assert( i < 0 );
1193         CV_Error(cv::Error::StsNotImplemented, "You should explicitly call download method for cuda::GpuMat object");
1194         return Mat();
1195     }
1196
1197     if( k == CUDA_MEM )
1198     {
1199         CV_Assert( i < 0 );
1200
1201         const cuda::CudaMem* cuda_mem = (const cuda::CudaMem*)obj;
1202
1203         return cuda_mem->createMatHeader();
1204     }
1205
1206     CV_Error(Error::StsNotImplemented, "Unknown/unsupported array type");
1207     return Mat();
1208 }
1209
1210 UMat _InputArray::getUMat(int i) const
1211 {
1212     int k = kind();
1213     int accessFlags = flags & ACCESS_MASK;
1214
1215     if( k == UMAT )
1216     {
1217         const UMat* m = (const UMat*)obj;
1218         if( i < 0 )
1219             return *m;
1220         return m->row(i);
1221     }
1222
1223     if( k == STD_VECTOR_UMAT )
1224     {
1225         const std::vector<UMat>& v = *(const std::vector<UMat>*)obj;
1226         CV_Assert( 0 <= i && i < (int)v.size() );
1227
1228         return v[i];
1229     }
1230
1231     if( k == MAT )
1232     {
1233         const Mat* m = (const Mat*)obj;
1234         if( i < 0 )
1235             return m->getUMat(accessFlags);
1236         return m->row(i).getUMat(accessFlags);
1237     }
1238
1239     return getMat(i).getUMat(accessFlags);
1240 }
1241
1242 void _InputArray::getMatVector(std::vector<Mat>& mv) const
1243 {
1244     int k = kind();
1245     int accessFlags = flags & ACCESS_MASK;
1246
1247     if( k == MAT )
1248     {
1249         const Mat& m = *(const Mat*)obj;
1250         int i, n = (int)m.size[0];
1251         mv.resize(n);
1252
1253         for( i = 0; i < n; i++ )
1254             mv[i] = m.dims == 2 ? Mat(1, m.cols, m.type(), (void*)m.ptr(i)) :
1255                 Mat(m.dims-1, &m.size[1], m.type(), (void*)m.ptr(i), &m.step[1]);
1256         return;
1257     }
1258
1259     if( k == EXPR )
1260     {
1261         Mat m = *(const MatExpr*)obj;
1262         int i, n = m.size[0];
1263         mv.resize(n);
1264
1265         for( i = 0; i < n; i++ )
1266             mv[i] = m.row(i);
1267         return;
1268     }
1269
1270     if( k == MATX )
1271     {
1272         size_t i, n = sz.height, esz = CV_ELEM_SIZE(flags);
1273         mv.resize(n);
1274
1275         for( i = 0; i < n; i++ )
1276             mv[i] = Mat(1, sz.width, CV_MAT_TYPE(flags), (uchar*)obj + esz*sz.width*i);
1277         return;
1278     }
1279
1280     if( k == STD_VECTOR )
1281     {
1282         const std::vector<uchar>& v = *(const std::vector<uchar>*)obj;
1283
1284         size_t i, n = v.size(), esz = CV_ELEM_SIZE(flags);
1285         int t = CV_MAT_DEPTH(flags), cn = CV_MAT_CN(flags);
1286         mv.resize(n);
1287
1288         for( i = 0; i < n; i++ )
1289             mv[i] = Mat(1, cn, t, (void*)(&v[0] + esz*i));
1290         return;
1291     }
1292
1293     if( k == NONE )
1294     {
1295         mv.clear();
1296         return;
1297     }
1298
1299     if( k == STD_VECTOR_VECTOR )
1300     {
1301         const std::vector<std::vector<uchar> >& vv = *(const std::vector<std::vector<uchar> >*)obj;
1302         int i, n = (int)vv.size();
1303         int t = CV_MAT_TYPE(flags);
1304         mv.resize(n);
1305
1306         for( i = 0; i < n; i++ )
1307         {
1308             const std::vector<uchar>& v = vv[i];
1309             mv[i] = Mat(size(i), t, (void*)&v[0]);
1310         }
1311         return;
1312     }
1313
1314     if( k == STD_VECTOR_MAT )
1315     {
1316         const std::vector<Mat>& v = *(const std::vector<Mat>*)obj;
1317         size_t i, n = v.size();
1318         mv.resize(n);
1319
1320         for( i = 0; i < n; i++ )
1321             mv[i] = v[i];
1322         return;
1323     }
1324
1325     if( k == STD_VECTOR_UMAT )
1326     {
1327         const std::vector<UMat>& v = *(const std::vector<UMat>*)obj;
1328         size_t i, n = v.size();
1329         mv.resize(n);
1330
1331         for( i = 0; i < n; i++ )
1332             mv[i] = v[i].getMat(accessFlags);
1333         return;
1334     }
1335
1336     CV_Error(Error::StsNotImplemented, "Unknown/unsupported array type");
1337 }
1338
1339 void _InputArray::getUMatVector(std::vector<UMat>& umv) const
1340 {
1341     int k = kind();
1342     int accessFlags = flags & ACCESS_MASK;
1343
1344     if( k == NONE )
1345     {
1346         umv.clear();
1347         return;
1348     }
1349
1350     if( k == STD_VECTOR_MAT )
1351     {
1352         const std::vector<Mat>& v = *(const std::vector<Mat>*)obj;
1353         size_t i, n = v.size();
1354         umv.resize(n);
1355
1356         for( i = 0; i < n; i++ )
1357             umv[i] = v[i].getUMat(accessFlags);
1358         return;
1359     }
1360
1361     if( k == STD_VECTOR_UMAT )
1362     {
1363         const std::vector<UMat>& v = *(const std::vector<UMat>*)obj;
1364         size_t i, n = v.size();
1365         umv.resize(n);
1366
1367         for( i = 0; i < n; i++ )
1368             umv[i] = v[i];
1369         return;
1370     }
1371
1372     if( k == UMAT )
1373     {
1374         UMat& v = *(UMat*)obj;
1375         umv.resize(1);
1376         umv[0] = v;
1377         return;
1378     }
1379     if( k == MAT )
1380     {
1381         Mat& v = *(Mat*)obj;
1382         umv.resize(1);
1383         umv[0] = v.getUMat(accessFlags);
1384         return;
1385     }
1386
1387     CV_Error(Error::StsNotImplemented, "Unknown/unsupported array type");
1388 }
1389
1390 cuda::GpuMat _InputArray::getGpuMat() const
1391 {
1392     int k = kind();
1393
1394     if (k == GPU_MAT)
1395     {
1396         const cuda::GpuMat* d_mat = (const cuda::GpuMat*)obj;
1397         return *d_mat;
1398     }
1399
1400     if (k == CUDA_MEM)
1401     {
1402         const cuda::CudaMem* cuda_mem = (const cuda::CudaMem*)obj;
1403         return cuda_mem->createGpuMatHeader();
1404     }
1405
1406     if (k == OPENGL_BUFFER)
1407     {
1408         CV_Error(cv::Error::StsNotImplemented, "You should explicitly call mapDevice/unmapDevice methods for ogl::Buffer object");
1409         return cuda::GpuMat();
1410     }
1411
1412     if (k == NONE)
1413         return cuda::GpuMat();
1414
1415     CV_Error(cv::Error::StsNotImplemented, "getGpuMat is available only for cuda::GpuMat and cuda::CudaMem");
1416     return cuda::GpuMat();
1417 }
1418
1419 ogl::Buffer _InputArray::getOGlBuffer() const
1420 {
1421     int k = kind();
1422
1423     CV_Assert(k == OPENGL_BUFFER);
1424
1425     const ogl::Buffer* gl_buf = (const ogl::Buffer*)obj;
1426     return *gl_buf;
1427 }
1428
1429 int _InputArray::kind() const
1430 {
1431     return flags & KIND_MASK;
1432 }
1433
1434 int _InputArray::rows(int i) const
1435 {
1436     return size(i).height;
1437 }
1438
1439 int _InputArray::cols(int i) const
1440 {
1441     return size(i).width;
1442 }
1443
1444 Size _InputArray::size(int i) const
1445 {
1446     int k = kind();
1447
1448     if( k == MAT )
1449     {
1450         CV_Assert( i < 0 );
1451         return ((const Mat*)obj)->size();
1452     }
1453
1454     if( k == EXPR )
1455     {
1456         CV_Assert( i < 0 );
1457         return ((const MatExpr*)obj)->size();
1458     }
1459
1460     if( k == UMAT )
1461     {
1462         CV_Assert( i < 0 );
1463         return ((const UMat*)obj)->size();
1464     }
1465
1466     if( k == MATX )
1467     {
1468         CV_Assert( i < 0 );
1469         return sz;
1470     }
1471
1472     if( k == STD_VECTOR )
1473     {
1474         CV_Assert( i < 0 );
1475         const std::vector<uchar>& v = *(const std::vector<uchar>*)obj;
1476         const std::vector<int>& iv = *(const std::vector<int>*)obj;
1477         size_t szb = v.size(), szi = iv.size();
1478         return szb == szi ? Size((int)szb, 1) : Size((int)(szb/CV_ELEM_SIZE(flags)), 1);
1479     }
1480
1481     if( k == NONE )
1482         return Size();
1483
1484     if( k == STD_VECTOR_VECTOR )
1485     {
1486         const std::vector<std::vector<uchar> >& vv = *(const std::vector<std::vector<uchar> >*)obj;
1487         if( i < 0 )
1488             return vv.empty() ? Size() : Size((int)vv.size(), 1);
1489         CV_Assert( i < (int)vv.size() );
1490         const std::vector<std::vector<int> >& ivv = *(const std::vector<std::vector<int> >*)obj;
1491
1492         size_t szb = vv[i].size(), szi = ivv[i].size();
1493         return szb == szi ? Size((int)szb, 1) : Size((int)(szb/CV_ELEM_SIZE(flags)), 1);
1494     }
1495
1496     if( k == STD_VECTOR_MAT )
1497     {
1498         const std::vector<Mat>& vv = *(const std::vector<Mat>*)obj;
1499         if( i < 0 )
1500             return vv.empty() ? Size() : Size((int)vv.size(), 1);
1501         CV_Assert( i < (int)vv.size() );
1502
1503         return vv[i].size();
1504     }
1505
1506     if( k == STD_VECTOR_UMAT )
1507     {
1508         const std::vector<UMat>& vv = *(const std::vector<UMat>*)obj;
1509         if( i < 0 )
1510             return vv.empty() ? Size() : Size((int)vv.size(), 1);
1511         CV_Assert( i < (int)vv.size() );
1512
1513         return vv[i].size();
1514     }
1515
1516     if( k == OPENGL_BUFFER )
1517     {
1518         CV_Assert( i < 0 );
1519         const ogl::Buffer* buf = (const ogl::Buffer*)obj;
1520         return buf->size();
1521     }
1522
1523     if( k == GPU_MAT )
1524     {
1525         CV_Assert( i < 0 );
1526         const cuda::GpuMat* d_mat = (const cuda::GpuMat*)obj;
1527         return d_mat->size();
1528     }
1529
1530     CV_Assert( k == CUDA_MEM );
1531     //if( k == CUDA_MEM )
1532     {
1533         CV_Assert( i < 0 );
1534         const cuda::CudaMem* cuda_mem = (const cuda::CudaMem*)obj;
1535         return cuda_mem->size();
1536     }
1537 }
1538
1539 int _InputArray::sizend(int* arrsz, int i) const
1540 {
1541     int j, d=0, k = kind();
1542
1543     if( k == NONE )
1544         ;
1545     else if( k == MAT )
1546     {
1547         CV_Assert( i < 0 );
1548         const Mat& m = *(const Mat*)obj;
1549         d = m.dims;
1550         if(arrsz)
1551             for(j = 0; j < d; j++)
1552                 arrsz[j] = m.size.p[j];
1553     }
1554     else if( k == UMAT )
1555     {
1556         CV_Assert( i < 0 );
1557         const UMat& m = *(const UMat*)obj;
1558         d = m.dims;
1559         if(arrsz)
1560             for(j = 0; j < d; j++)
1561                 arrsz[j] = m.size.p[j];
1562     }
1563     else if( k == STD_VECTOR_MAT && i >= 0 )
1564     {
1565         const std::vector<Mat>& vv = *(const std::vector<Mat>*)obj;
1566         CV_Assert( i < (int)vv.size() );
1567         const Mat& m = vv[i];
1568         d = m.dims;
1569         if(arrsz)
1570             for(j = 0; j < d; j++)
1571                 arrsz[j] = m.size.p[j];
1572     }
1573     else if( k == STD_VECTOR_UMAT && i >= 0 )
1574     {
1575         const std::vector<UMat>& vv = *(const std::vector<UMat>*)obj;
1576         CV_Assert( i < (int)vv.size() );
1577         const UMat& m = vv[i];
1578         d = m.dims;
1579         if(arrsz)
1580             for(j = 0; j < d; j++)
1581                 arrsz[j] = m.size.p[j];
1582     }
1583     else
1584     {
1585         Size sz2d = size(i);
1586         d = 2;
1587         if(arrsz)
1588         {
1589             arrsz[0] = sz2d.height;
1590             arrsz[1] = sz2d.width;
1591         }
1592     }
1593
1594     return d;
1595 }
1596
1597 bool _InputArray::sameSize(const _InputArray& arr) const
1598 {
1599     int k1 = kind(), k2 = arr.kind();
1600     Size sz1;
1601
1602     if( k1 == MAT )
1603     {
1604         const Mat* m = ((const Mat*)obj);
1605         if( k2 == MAT )
1606             return m->size == ((const Mat*)arr.obj)->size;
1607         if( k2 == UMAT )
1608             return m->size == ((const UMat*)arr.obj)->size;
1609         if( m->dims > 2 )
1610             return false;
1611         sz1 = m->size();
1612     }
1613     else if( k1 == UMAT )
1614     {
1615         const UMat* m = ((const UMat*)obj);
1616         if( k2 == MAT )
1617             return m->size == ((const Mat*)arr.obj)->size;
1618         if( k2 == UMAT )
1619             return m->size == ((const UMat*)arr.obj)->size;
1620         if( m->dims > 2 )
1621             return false;
1622         sz1 = m->size();
1623     }
1624     else
1625         sz1 = size();
1626     if( arr.dims() > 2 )
1627         return false;
1628     return sz1 == arr.size();
1629 }
1630
1631 int _InputArray::dims(int i) const
1632 {
1633     int k = kind();
1634
1635     if( k == MAT )
1636     {
1637         CV_Assert( i < 0 );
1638         return ((const Mat*)obj)->dims;
1639     }
1640
1641     if( k == EXPR )
1642     {
1643         CV_Assert( i < 0 );
1644         return ((const MatExpr*)obj)->a.dims;
1645     }
1646
1647     if( k == UMAT )
1648     {
1649         CV_Assert( i < 0 );
1650         return ((const UMat*)obj)->dims;
1651     }
1652
1653     if( k == MATX )
1654     {
1655         CV_Assert( i < 0 );
1656         return 2;
1657     }
1658
1659     if( k == STD_VECTOR )
1660     {
1661         CV_Assert( i < 0 );
1662         return 2;
1663     }
1664
1665     if( k == NONE )
1666         return 0;
1667
1668     if( k == STD_VECTOR_VECTOR )
1669     {
1670         const std::vector<std::vector<uchar> >& vv = *(const std::vector<std::vector<uchar> >*)obj;
1671         if( i < 0 )
1672             return 1;
1673         CV_Assert( i < (int)vv.size() );
1674         return 2;
1675     }
1676
1677     if( k == STD_VECTOR_MAT )
1678     {
1679         const std::vector<Mat>& vv = *(const std::vector<Mat>*)obj;
1680         if( i < 0 )
1681             return 1;
1682         CV_Assert( i < (int)vv.size() );
1683
1684         return vv[i].dims;
1685     }
1686
1687     if( k == STD_VECTOR_UMAT )
1688     {
1689         const std::vector<UMat>& vv = *(const std::vector<UMat>*)obj;
1690         if( i < 0 )
1691             return 1;
1692         CV_Assert( i < (int)vv.size() );
1693
1694         return vv[i].dims;
1695     }
1696
1697     if( k == OPENGL_BUFFER )
1698     {
1699         CV_Assert( i < 0 );
1700         return 2;
1701     }
1702
1703     if( k == GPU_MAT )
1704     {
1705         CV_Assert( i < 0 );
1706         return 2;
1707     }
1708
1709     CV_Assert( k == CUDA_MEM );
1710     //if( k == CUDA_MEM )
1711     {
1712         CV_Assert( i < 0 );
1713         return 2;
1714     }
1715 }
1716
1717 size_t _InputArray::total(int i) const
1718 {
1719     int k = kind();
1720
1721     if( k == MAT )
1722     {
1723         CV_Assert( i < 0 );
1724         return ((const Mat*)obj)->total();
1725     }
1726
1727     if( k == UMAT )
1728     {
1729         CV_Assert( i < 0 );
1730         return ((const UMat*)obj)->total();
1731     }
1732
1733     if( k == STD_VECTOR_MAT )
1734     {
1735         const std::vector<Mat>& vv = *(const std::vector<Mat>*)obj;
1736         if( i < 0 )
1737             return vv.size();
1738
1739         CV_Assert( i < (int)vv.size() );
1740         return vv[i].total();
1741     }
1742
1743     if( k == STD_VECTOR_UMAT )
1744     {
1745         const std::vector<UMat>& vv = *(const std::vector<UMat>*)obj;
1746         if( i < 0 )
1747             return vv.size();
1748
1749         CV_Assert( i < (int)vv.size() );
1750         return vv[i].total();
1751     }
1752
1753     return size(i).area();
1754 }
1755
1756 int _InputArray::type(int i) const
1757 {
1758     int k = kind();
1759
1760     if( k == MAT )
1761         return ((const Mat*)obj)->type();
1762
1763     if( k == UMAT )
1764         return ((const UMat*)obj)->type();
1765
1766     if( k == EXPR )
1767         return ((const MatExpr*)obj)->type();
1768
1769     if( k == MATX || k == STD_VECTOR || k == STD_VECTOR_VECTOR )
1770         return CV_MAT_TYPE(flags);
1771
1772     if( k == NONE )
1773         return -1;
1774
1775     if( k == STD_VECTOR_UMAT )
1776     {
1777         const std::vector<UMat>& vv = *(const std::vector<UMat>*)obj;
1778         if( vv.empty() )
1779         {
1780             CV_Assert((flags & FIXED_TYPE) != 0);
1781             return CV_MAT_TYPE(flags);
1782         }
1783         CV_Assert( i < (int)vv.size() );
1784         return vv[i >= 0 ? i : 0].type();
1785     }
1786
1787     if( k == STD_VECTOR_MAT )
1788     {
1789         const std::vector<Mat>& vv = *(const std::vector<Mat>*)obj;
1790         if( vv.empty() )
1791         {
1792             CV_Assert((flags & FIXED_TYPE) != 0);
1793             return CV_MAT_TYPE(flags);
1794         }
1795         CV_Assert( i < (int)vv.size() );
1796         return vv[i >= 0 ? i : 0].type();
1797     }
1798
1799     if( k == OPENGL_BUFFER )
1800         return ((const ogl::Buffer*)obj)->type();
1801
1802     if( k == GPU_MAT )
1803         return ((const cuda::GpuMat*)obj)->type();
1804
1805     CV_Assert( k == CUDA_MEM );
1806     //if( k == CUDA_MEM )
1807         return ((const cuda::CudaMem*)obj)->type();
1808 }
1809
1810 int _InputArray::depth(int i) const
1811 {
1812     return CV_MAT_DEPTH(type(i));
1813 }
1814
1815 int _InputArray::channels(int i) const
1816 {
1817     return CV_MAT_CN(type(i));
1818 }
1819
1820 bool _InputArray::empty() const
1821 {
1822     int k = kind();
1823
1824     if( k == MAT )
1825         return ((const Mat*)obj)->empty();
1826
1827     if( k == UMAT )
1828         return ((const UMat*)obj)->empty();
1829
1830     if( k == EXPR )
1831         return false;
1832
1833     if( k == MATX )
1834         return false;
1835
1836     if( k == STD_VECTOR )
1837     {
1838         const std::vector<uchar>& v = *(const std::vector<uchar>*)obj;
1839         return v.empty();
1840     }
1841
1842     if( k == NONE )
1843         return true;
1844
1845     if( k == STD_VECTOR_VECTOR )
1846     {
1847         const std::vector<std::vector<uchar> >& vv = *(const std::vector<std::vector<uchar> >*)obj;
1848         return vv.empty();
1849     }
1850
1851     if( k == STD_VECTOR_MAT )
1852     {
1853         const std::vector<Mat>& vv = *(const std::vector<Mat>*)obj;
1854         return vv.empty();
1855     }
1856
1857     if( k == STD_VECTOR_UMAT )
1858     {
1859         const std::vector<UMat>& vv = *(const std::vector<UMat>*)obj;
1860         return vv.empty();
1861     }
1862
1863     if( k == OPENGL_BUFFER )
1864         return ((const ogl::Buffer*)obj)->empty();
1865
1866     if( k == GPU_MAT )
1867         return ((const cuda::GpuMat*)obj)->empty();
1868
1869     CV_Assert( k == CUDA_MEM );
1870     //if( k == CUDA_MEM )
1871         return ((const cuda::CudaMem*)obj)->empty();
1872 }
1873
1874 bool _InputArray::isContinuous(int i) const
1875 {
1876     int k = kind();
1877
1878     if( k == MAT )
1879         return i < 0 ? ((const Mat*)obj)->isContinuous() : true;
1880
1881     if( k == UMAT )
1882         return i < 0 ? ((const UMat*)obj)->isContinuous() : true;
1883
1884     if( k == EXPR || k == MATX || k == STD_VECTOR || k == NONE || k == STD_VECTOR_VECTOR)
1885         return true;
1886
1887     if( k == STD_VECTOR_MAT )
1888     {
1889         const std::vector<Mat>& vv = *(const std::vector<Mat>*)obj;
1890         CV_Assert((size_t)i < vv.size());
1891         return vv[i].isContinuous();
1892     }
1893
1894     if( k == STD_VECTOR_UMAT )
1895     {
1896         const std::vector<UMat>& vv = *(const std::vector<UMat>*)obj;
1897         CV_Assert((size_t)i < vv.size());
1898         return vv[i].isContinuous();
1899     }
1900
1901     CV_Error(CV_StsNotImplemented, "Unknown/unsupported array type");
1902     return false;
1903 }
1904
1905 bool _InputArray::isSubmatrix(int i) const
1906 {
1907     int k = kind();
1908
1909     if( k == MAT )
1910         return i < 0 ? ((const Mat*)obj)->isSubmatrix() : false;
1911
1912     if( k == UMAT )
1913         return i < 0 ? ((const UMat*)obj)->isSubmatrix() : false;
1914
1915     if( k == EXPR || k == MATX || k == STD_VECTOR || k == NONE || k == STD_VECTOR_VECTOR)
1916         return false;
1917
1918     if( k == STD_VECTOR_MAT )
1919     {
1920         const std::vector<Mat>& vv = *(const std::vector<Mat>*)obj;
1921         CV_Assert((size_t)i < vv.size());
1922         return vv[i].isSubmatrix();
1923     }
1924
1925     if( k == STD_VECTOR_UMAT )
1926     {
1927         const std::vector<UMat>& vv = *(const std::vector<UMat>*)obj;
1928         CV_Assert((size_t)i < vv.size());
1929         return vv[i].isSubmatrix();
1930     }
1931
1932     CV_Error(CV_StsNotImplemented, "");
1933     return false;
1934 }
1935
1936 size_t _InputArray::offset(int i) const
1937 {
1938     int k = kind();
1939
1940     if( k == MAT )
1941     {
1942         CV_Assert( i < 0 );
1943         const Mat * const m = ((const Mat*)obj);
1944         return (size_t)(m->data - m->datastart);
1945     }
1946
1947     if( k == UMAT )
1948     {
1949         CV_Assert( i < 0 );
1950         return ((const UMat*)obj)->offset;
1951     }
1952
1953     if( k == EXPR || k == MATX || k == STD_VECTOR || k == NONE || k == STD_VECTOR_VECTOR)
1954         return 0;
1955
1956     if( k == STD_VECTOR_MAT )
1957     {
1958         const std::vector<Mat>& vv = *(const std::vector<Mat>*)obj;
1959         if( i < 0 )
1960             return 1;
1961         CV_Assert( i < (int)vv.size() );
1962
1963         return (size_t)(vv[i].data - vv[i].datastart);
1964     }
1965
1966     if( k == STD_VECTOR_UMAT )
1967     {
1968         const std::vector<UMat>& vv = *(const std::vector<UMat>*)obj;
1969         CV_Assert((size_t)i < vv.size());
1970         return vv[i].offset;
1971     }
1972
1973     if( k == GPU_MAT )
1974     {
1975         CV_Assert( i < 0 );
1976         const cuda::GpuMat * const m = ((const cuda::GpuMat*)obj);
1977         return (size_t)(m->data - m->datastart);
1978     }
1979
1980     CV_Error(Error::StsNotImplemented, "");
1981     return 0;
1982 }
1983
1984 size_t _InputArray::step(int i) const
1985 {
1986     int k = kind();
1987
1988     if( k == MAT )
1989     {
1990         CV_Assert( i < 0 );
1991         return ((const Mat*)obj)->step;
1992     }
1993
1994     if( k == UMAT )
1995     {
1996         CV_Assert( i < 0 );
1997         return ((const UMat*)obj)->step;
1998     }
1999
2000     if( k == EXPR || k == MATX || k == STD_VECTOR || k == NONE || k == STD_VECTOR_VECTOR)
2001         return 0;
2002
2003     if( k == STD_VECTOR_MAT )
2004     {
2005         const std::vector<Mat>& vv = *(const std::vector<Mat>*)obj;
2006         if( i < 0 )
2007             return 1;
2008         CV_Assert( i < (int)vv.size() );
2009         return vv[i].step;
2010     }
2011
2012     if( k == STD_VECTOR_UMAT )
2013     {
2014         const std::vector<UMat>& vv = *(const std::vector<UMat>*)obj;
2015         CV_Assert((size_t)i < vv.size());
2016         return vv[i].step;
2017     }
2018
2019     if( k == GPU_MAT )
2020     {
2021         CV_Assert( i < 0 );
2022         return ((const cuda::GpuMat*)obj)->step;
2023     }
2024
2025     CV_Error(Error::StsNotImplemented, "");
2026     return 0;
2027 }
2028
2029 void _InputArray::copyTo(const _OutputArray& arr) const
2030 {
2031     int k = kind();
2032
2033     if( k == NONE )
2034         arr.release();
2035     else if( k == MAT || k == MATX || k == STD_VECTOR )
2036     {
2037         Mat m = getMat();
2038         m.copyTo(arr);
2039     }
2040     else if( k == EXPR )
2041     {
2042         const MatExpr& e = *((MatExpr*)obj);
2043         if( arr.kind() == MAT )
2044             arr.getMatRef() = e;
2045         else
2046             Mat(e).copyTo(arr);
2047     }
2048     else if( k == UMAT )
2049         ((UMat*)obj)->copyTo(arr);
2050     else
2051         CV_Error(Error::StsNotImplemented, "");
2052 }
2053
2054 void _InputArray::copyTo(const _OutputArray& arr, const _InputArray & mask) const
2055 {
2056     int k = kind();
2057
2058     if( k == NONE )
2059         arr.release();
2060     else if( k == MAT || k == MATX || k == STD_VECTOR )
2061     {
2062         Mat m = getMat();
2063         m.copyTo(arr, mask);
2064     }
2065     else if( k == UMAT )
2066         ((UMat*)obj)->copyTo(arr, mask);
2067     else
2068         CV_Error(Error::StsNotImplemented, "");
2069 }
2070
2071 bool _OutputArray::fixedSize() const
2072 {
2073     return (flags & FIXED_SIZE) == FIXED_SIZE;
2074 }
2075
2076 bool _OutputArray::fixedType() const
2077 {
2078     return (flags & FIXED_TYPE) == FIXED_TYPE;
2079 }
2080
2081 void _OutputArray::create(Size _sz, int mtype, int i, bool allowTransposed, int fixedDepthMask) const
2082 {
2083     int k = kind();
2084     if( k == MAT && i < 0 && !allowTransposed && fixedDepthMask == 0 )
2085     {
2086         CV_Assert(!fixedSize() || ((Mat*)obj)->size.operator()() == _sz);
2087         CV_Assert(!fixedType() || ((Mat*)obj)->type() == mtype);
2088         ((Mat*)obj)->create(_sz, mtype);
2089         return;
2090     }
2091     if( k == UMAT && i < 0 && !allowTransposed && fixedDepthMask == 0 )
2092     {
2093         CV_Assert(!fixedSize() || ((UMat*)obj)->size.operator()() == _sz);
2094         CV_Assert(!fixedType() || ((UMat*)obj)->type() == mtype);
2095         ((UMat*)obj)->create(_sz, mtype);
2096         return;
2097     }
2098     if( k == GPU_MAT && i < 0 && !allowTransposed && fixedDepthMask == 0 )
2099     {
2100         CV_Assert(!fixedSize() || ((cuda::GpuMat*)obj)->size() == _sz);
2101         CV_Assert(!fixedType() || ((cuda::GpuMat*)obj)->type() == mtype);
2102         ((cuda::GpuMat*)obj)->create(_sz, mtype);
2103         return;
2104     }
2105     if( k == OPENGL_BUFFER && i < 0 && !allowTransposed && fixedDepthMask == 0 )
2106     {
2107         CV_Assert(!fixedSize() || ((ogl::Buffer*)obj)->size() == _sz);
2108         CV_Assert(!fixedType() || ((ogl::Buffer*)obj)->type() == mtype);
2109         ((ogl::Buffer*)obj)->create(_sz, mtype);
2110         return;
2111     }
2112     if( k == CUDA_MEM && i < 0 && !allowTransposed && fixedDepthMask == 0 )
2113     {
2114         CV_Assert(!fixedSize() || ((cuda::CudaMem*)obj)->size() == _sz);
2115         CV_Assert(!fixedType() || ((cuda::CudaMem*)obj)->type() == mtype);
2116         ((cuda::CudaMem*)obj)->create(_sz, mtype);
2117         return;
2118     }
2119     int sizes[] = {_sz.height, _sz.width};
2120     create(2, sizes, mtype, i, allowTransposed, fixedDepthMask);
2121 }
2122
2123 void _OutputArray::create(int _rows, int _cols, int mtype, int i, bool allowTransposed, int fixedDepthMask) const
2124 {
2125     int k = kind();
2126     if( k == MAT && i < 0 && !allowTransposed && fixedDepthMask == 0 )
2127     {
2128         CV_Assert(!fixedSize() || ((Mat*)obj)->size.operator()() == Size(_cols, _rows));
2129         CV_Assert(!fixedType() || ((Mat*)obj)->type() == mtype);
2130         ((Mat*)obj)->create(_rows, _cols, mtype);
2131         return;
2132     }
2133     if( k == UMAT && i < 0 && !allowTransposed && fixedDepthMask == 0 )
2134     {
2135         CV_Assert(!fixedSize() || ((UMat*)obj)->size.operator()() == Size(_cols, _rows));
2136         CV_Assert(!fixedType() || ((UMat*)obj)->type() == mtype);
2137         ((UMat*)obj)->create(_rows, _cols, mtype);
2138         return;
2139     }
2140     if( k == GPU_MAT && i < 0 && !allowTransposed && fixedDepthMask == 0 )
2141     {
2142         CV_Assert(!fixedSize() || ((cuda::GpuMat*)obj)->size() == Size(_cols, _rows));
2143         CV_Assert(!fixedType() || ((cuda::GpuMat*)obj)->type() == mtype);
2144         ((cuda::GpuMat*)obj)->create(_rows, _cols, mtype);
2145         return;
2146     }
2147     if( k == OPENGL_BUFFER && i < 0 && !allowTransposed && fixedDepthMask == 0 )
2148     {
2149         CV_Assert(!fixedSize() || ((ogl::Buffer*)obj)->size() == Size(_cols, _rows));
2150         CV_Assert(!fixedType() || ((ogl::Buffer*)obj)->type() == mtype);
2151         ((ogl::Buffer*)obj)->create(_rows, _cols, mtype);
2152         return;
2153     }
2154     if( k == CUDA_MEM && i < 0 && !allowTransposed && fixedDepthMask == 0 )
2155     {
2156         CV_Assert(!fixedSize() || ((cuda::CudaMem*)obj)->size() == Size(_cols, _rows));
2157         CV_Assert(!fixedType() || ((cuda::CudaMem*)obj)->type() == mtype);
2158         ((cuda::CudaMem*)obj)->create(_rows, _cols, mtype);
2159         return;
2160     }
2161     int sizes[] = {_rows, _cols};
2162     create(2, sizes, mtype, i, allowTransposed, fixedDepthMask);
2163 }
2164
2165 void _OutputArray::create(int d, const int* sizes, int mtype, int i,
2166                           bool allowTransposed, int fixedDepthMask) const
2167 {
2168     int k = kind();
2169     mtype = CV_MAT_TYPE(mtype);
2170
2171     if( k == MAT )
2172     {
2173         CV_Assert( i < 0 );
2174         Mat& m = *(Mat*)obj;
2175         if( allowTransposed )
2176         {
2177             if( !m.isContinuous() )
2178             {
2179                 CV_Assert(!fixedType() && !fixedSize());
2180                 m.release();
2181             }
2182
2183             if( d == 2 && m.dims == 2 && m.data &&
2184                 m.type() == mtype && m.rows == sizes[1] && m.cols == sizes[0] )
2185                 return;
2186         }
2187
2188         if(fixedType())
2189         {
2190             if(CV_MAT_CN(mtype) == m.channels() && ((1 << CV_MAT_TYPE(flags)) & fixedDepthMask) != 0 )
2191                 mtype = m.type();
2192             else
2193                 CV_Assert(CV_MAT_TYPE(mtype) == m.type());
2194         }
2195         if(fixedSize())
2196         {
2197             CV_Assert(m.dims == d);
2198             for(int j = 0; j < d; ++j)
2199                 CV_Assert(m.size[j] == sizes[j]);
2200         }
2201         m.create(d, sizes, mtype);
2202         return;
2203     }
2204
2205     if( k == UMAT )
2206     {
2207         CV_Assert( i < 0 );
2208         UMat& m = *(UMat*)obj;
2209         if( allowTransposed )
2210         {
2211             if( !m.isContinuous() )
2212             {
2213                 CV_Assert(!fixedType() && !fixedSize());
2214                 m.release();
2215             }
2216
2217             if( d == 2 && m.dims == 2 && !m.empty() &&
2218                 m.type() == mtype && m.rows == sizes[1] && m.cols == sizes[0] )
2219                 return;
2220         }
2221
2222         if(fixedType())
2223         {
2224             if(CV_MAT_CN(mtype) == m.channels() && ((1 << CV_MAT_TYPE(flags)) & fixedDepthMask) != 0 )
2225                 mtype = m.type();
2226             else
2227                 CV_Assert(CV_MAT_TYPE(mtype) == m.type());
2228         }
2229         if(fixedSize())
2230         {
2231             CV_Assert(m.dims == d);
2232             for(int j = 0; j < d; ++j)
2233                 CV_Assert(m.size[j] == sizes[j]);
2234         }
2235         m.create(d, sizes, mtype);
2236         return;
2237     }
2238
2239     if( k == MATX )
2240     {
2241         CV_Assert( i < 0 );
2242         int type0 = CV_MAT_TYPE(flags);
2243         CV_Assert( mtype == type0 || (CV_MAT_CN(mtype) == 1 && ((1 << type0) & fixedDepthMask) != 0) );
2244         CV_Assert( d == 2 && ((sizes[0] == sz.height && sizes[1] == sz.width) ||
2245                                  (allowTransposed && sizes[0] == sz.width && sizes[1] == sz.height)));
2246         return;
2247     }
2248
2249     if( k == STD_VECTOR || k == STD_VECTOR_VECTOR )
2250     {
2251         CV_Assert( d == 2 && (sizes[0] == 1 || sizes[1] == 1 || sizes[0]*sizes[1] == 0) );
2252         size_t len = sizes[0]*sizes[1] > 0 ? sizes[0] + sizes[1] - 1 : 0;
2253         std::vector<uchar>* v = (std::vector<uchar>*)obj;
2254
2255         if( k == STD_VECTOR_VECTOR )
2256         {
2257             std::vector<std::vector<uchar> >& vv = *(std::vector<std::vector<uchar> >*)obj;
2258             if( i < 0 )
2259             {
2260                 CV_Assert(!fixedSize() || len == vv.size());
2261                 vv.resize(len);
2262                 return;
2263             }
2264             CV_Assert( i < (int)vv.size() );
2265             v = &vv[i];
2266         }
2267         else
2268             CV_Assert( i < 0 );
2269
2270         int type0 = CV_MAT_TYPE(flags);
2271         CV_Assert( mtype == type0 || (CV_MAT_CN(mtype) == CV_MAT_CN(type0) && ((1 << type0) & fixedDepthMask) != 0) );
2272
2273         int esz = CV_ELEM_SIZE(type0);
2274         CV_Assert(!fixedSize() || len == ((std::vector<uchar>*)v)->size() / esz);
2275         switch( esz )
2276         {
2277         case 1:
2278             ((std::vector<uchar>*)v)->resize(len);
2279             break;
2280         case 2:
2281             ((std::vector<Vec2b>*)v)->resize(len);
2282             break;
2283         case 3:
2284             ((std::vector<Vec3b>*)v)->resize(len);
2285             break;
2286         case 4:
2287             ((std::vector<int>*)v)->resize(len);
2288             break;
2289         case 6:
2290             ((std::vector<Vec3s>*)v)->resize(len);
2291             break;
2292         case 8:
2293             ((std::vector<Vec2i>*)v)->resize(len);
2294             break;
2295         case 12:
2296             ((std::vector<Vec3i>*)v)->resize(len);
2297             break;
2298         case 16:
2299             ((std::vector<Vec4i>*)v)->resize(len);
2300             break;
2301         case 24:
2302             ((std::vector<Vec6i>*)v)->resize(len);
2303             break;
2304         case 32:
2305             ((std::vector<Vec8i>*)v)->resize(len);
2306             break;
2307         case 36:
2308             ((std::vector<Vec<int, 9> >*)v)->resize(len);
2309             break;
2310         case 48:
2311             ((std::vector<Vec<int, 12> >*)v)->resize(len);
2312             break;
2313         case 64:
2314             ((std::vector<Vec<int, 16> >*)v)->resize(len);
2315             break;
2316         case 128:
2317             ((std::vector<Vec<int, 32> >*)v)->resize(len);
2318             break;
2319         case 256:
2320             ((std::vector<Vec<int, 64> >*)v)->resize(len);
2321             break;
2322         case 512:
2323             ((std::vector<Vec<int, 128> >*)v)->resize(len);
2324             break;
2325         default:
2326             CV_Error_(CV_StsBadArg, ("Vectors with element size %d are not supported. Please, modify OutputArray::create()\n", esz));
2327         }
2328         return;
2329     }
2330
2331     if( k == NONE )
2332     {
2333         CV_Error(CV_StsNullPtr, "create() called for the missing output array" );
2334         return;
2335     }
2336
2337     if( k == STD_VECTOR_MAT )
2338     {
2339         std::vector<Mat>& v = *(std::vector<Mat>*)obj;
2340
2341         if( i < 0 )
2342         {
2343             CV_Assert( d == 2 && (sizes[0] == 1 || sizes[1] == 1 || sizes[0]*sizes[1] == 0) );
2344             size_t len = sizes[0]*sizes[1] > 0 ? sizes[0] + sizes[1] - 1 : 0, len0 = v.size();
2345
2346             CV_Assert(!fixedSize() || len == len0);
2347             v.resize(len);
2348             if( fixedType() )
2349             {
2350                 int _type = CV_MAT_TYPE(flags);
2351                 for( size_t j = len0; j < len; j++ )
2352                 {
2353                     if( v[j].type() == _type )
2354                         continue;
2355                     CV_Assert( v[j].empty() );
2356                     v[j].flags = (v[j].flags & ~CV_MAT_TYPE_MASK) | _type;
2357                 }
2358             }
2359             return;
2360         }
2361
2362         CV_Assert( i < (int)v.size() );
2363         Mat& m = v[i];
2364
2365         if( allowTransposed )
2366         {
2367             if( !m.isContinuous() )
2368             {
2369                 CV_Assert(!fixedType() && !fixedSize());
2370                 m.release();
2371             }
2372
2373             if( d == 2 && m.dims == 2 && m.data &&
2374                 m.type() == mtype && m.rows == sizes[1] && m.cols == sizes[0] )
2375                 return;
2376         }
2377
2378         if(fixedType())
2379         {
2380             if(CV_MAT_CN(mtype) == m.channels() && ((1 << CV_MAT_TYPE(flags)) & fixedDepthMask) != 0 )
2381                 mtype = m.type();
2382             else
2383                 CV_Assert(CV_MAT_TYPE(mtype) == m.type());
2384         }
2385         if(fixedSize())
2386         {
2387             CV_Assert(m.dims == d);
2388             for(int j = 0; j < d; ++j)
2389                 CV_Assert(m.size[j] == sizes[j]);
2390         }
2391
2392         m.create(d, sizes, mtype);
2393         return;
2394     }
2395
2396     if( k == STD_VECTOR_UMAT )
2397     {
2398         std::vector<UMat>& v = *(std::vector<UMat>*)obj;
2399
2400         if( i < 0 )
2401         {
2402             CV_Assert( d == 2 && (sizes[0] == 1 || sizes[1] == 1 || sizes[0]*sizes[1] == 0) );
2403             size_t len = sizes[0]*sizes[1] > 0 ? sizes[0] + sizes[1] - 1 : 0, len0 = v.size();
2404
2405             CV_Assert(!fixedSize() || len == len0);
2406             v.resize(len);
2407             if( fixedType() )
2408             {
2409                 int _type = CV_MAT_TYPE(flags);
2410                 for( size_t j = len0; j < len; j++ )
2411                 {
2412                     if( v[j].type() == _type )
2413                         continue;
2414                     CV_Assert( v[j].empty() );
2415                     v[j].flags = (v[j].flags & ~CV_MAT_TYPE_MASK) | _type;
2416                 }
2417             }
2418             return;
2419         }
2420
2421         CV_Assert( i < (int)v.size() );
2422         UMat& m = v[i];
2423
2424         if( allowTransposed )
2425         {
2426             if( !m.isContinuous() )
2427             {
2428                 CV_Assert(!fixedType() && !fixedSize());
2429                 m.release();
2430             }
2431
2432             if( d == 2 && m.dims == 2 && m.u &&
2433                 m.type() == mtype && m.rows == sizes[1] && m.cols == sizes[0] )
2434                 return;
2435         }
2436
2437         if(fixedType())
2438         {
2439             if(CV_MAT_CN(mtype) == m.channels() && ((1 << CV_MAT_TYPE(flags)) & fixedDepthMask) != 0 )
2440                 mtype = m.type();
2441             else
2442                 CV_Assert(CV_MAT_TYPE(mtype) == m.type());
2443         }
2444         if(fixedSize())
2445         {
2446             CV_Assert(m.dims == d);
2447             for(int j = 0; j < d; ++j)
2448                 CV_Assert(m.size[j] == sizes[j]);
2449         }
2450
2451         m.create(d, sizes, mtype);
2452         return;
2453     }
2454
2455     CV_Error(Error::StsNotImplemented, "Unknown/unsupported array type");
2456 }
2457
2458 void _OutputArray::createSameSize(const _InputArray& arr, int mtype) const
2459 {
2460     int arrsz[CV_MAX_DIM], d = arr.sizend(arrsz);
2461     create(d, arrsz, mtype);
2462 }
2463
2464 void _OutputArray::release() const
2465 {
2466     CV_Assert(!fixedSize());
2467
2468     int k = kind();
2469
2470     if( k == MAT )
2471     {
2472         ((Mat*)obj)->release();
2473         return;
2474     }
2475
2476     if( k == UMAT )
2477     {
2478         ((UMat*)obj)->release();
2479         return;
2480     }
2481
2482     if( k == GPU_MAT )
2483     {
2484         ((cuda::GpuMat*)obj)->release();
2485         return;
2486     }
2487
2488     if( k == CUDA_MEM )
2489     {
2490         ((cuda::CudaMem*)obj)->release();
2491         return;
2492     }
2493
2494     if( k == OPENGL_BUFFER )
2495     {
2496         ((ogl::Buffer*)obj)->release();
2497         return;
2498     }
2499
2500     if( k == NONE )
2501         return;
2502
2503     if( k == STD_VECTOR )
2504     {
2505         create(Size(), CV_MAT_TYPE(flags));
2506         return;
2507     }
2508
2509     if( k == STD_VECTOR_VECTOR )
2510     {
2511         ((std::vector<std::vector<uchar> >*)obj)->clear();
2512         return;
2513     }
2514
2515     if( k == STD_VECTOR_MAT )
2516     {
2517         ((std::vector<Mat>*)obj)->clear();
2518         return;
2519     }
2520
2521     if( k == STD_VECTOR_UMAT )
2522     {
2523         ((std::vector<UMat>*)obj)->clear();
2524         return;
2525     }
2526
2527     CV_Error(Error::StsNotImplemented, "Unknown/unsupported array type");
2528 }
2529
2530 void _OutputArray::clear() const
2531 {
2532     int k = kind();
2533
2534     if( k == MAT )
2535     {
2536         CV_Assert(!fixedSize());
2537         ((Mat*)obj)->resize(0);
2538         return;
2539     }
2540
2541     release();
2542 }
2543
2544 bool _OutputArray::needed() const
2545 {
2546     return kind() != NONE;
2547 }
2548
2549 Mat& _OutputArray::getMatRef(int i) const
2550 {
2551     int k = kind();
2552     if( i < 0 )
2553     {
2554         CV_Assert( k == MAT );
2555         return *(Mat*)obj;
2556     }
2557     else
2558     {
2559         CV_Assert( k == STD_VECTOR_MAT );
2560         std::vector<Mat>& v = *(std::vector<Mat>*)obj;
2561         CV_Assert( i < (int)v.size() );
2562         return v[i];
2563     }
2564 }
2565
2566 UMat& _OutputArray::getUMatRef(int i) const
2567 {
2568     int k = kind();
2569     if( i < 0 )
2570     {
2571         CV_Assert( k == UMAT );
2572         return *(UMat*)obj;
2573     }
2574     else
2575     {
2576         CV_Assert( k == STD_VECTOR_UMAT );
2577         std::vector<UMat>& v = *(std::vector<UMat>*)obj;
2578         CV_Assert( i < (int)v.size() );
2579         return v[i];
2580     }
2581 }
2582
2583 cuda::GpuMat& _OutputArray::getGpuMatRef() const
2584 {
2585     int k = kind();
2586     CV_Assert( k == GPU_MAT );
2587     return *(cuda::GpuMat*)obj;
2588 }
2589
2590 ogl::Buffer& _OutputArray::getOGlBufferRef() const
2591 {
2592     int k = kind();
2593     CV_Assert( k == OPENGL_BUFFER );
2594     return *(ogl::Buffer*)obj;
2595 }
2596
2597 cuda::CudaMem& _OutputArray::getCudaMemRef() const
2598 {
2599     int k = kind();
2600     CV_Assert( k == CUDA_MEM );
2601     return *(cuda::CudaMem*)obj;
2602 }
2603
2604 void _OutputArray::setTo(const _InputArray& arr, const _InputArray & mask) const
2605 {
2606     int k = kind();
2607
2608     if( k == NONE )
2609         ;
2610     else if( k == MAT || k == MATX || k == STD_VECTOR )
2611     {
2612         Mat m = getMat();
2613         m.setTo(arr, mask);
2614     }
2615     else if( k == UMAT )
2616         ((UMat*)obj)->setTo(arr, mask);
2617     else if( k == GPU_MAT )
2618     {
2619         Mat value = arr.getMat();
2620         CV_Assert( checkScalar(value, type(), arr.kind(), _InputArray::GPU_MAT) );
2621         ((cuda::GpuMat*)obj)->setTo(Scalar(Vec<double, 4>((double *)value.data)), mask);
2622     }
2623     else
2624         CV_Error(Error::StsNotImplemented, "");
2625 }
2626
2627
2628 void _OutputArray::assign(const UMat& u) const
2629 {
2630     int k = kind();
2631     if (k == UMAT)
2632     {
2633         *(UMat*)obj = u;
2634     }
2635     else if (k == MAT)
2636     {
2637         u.copyTo(*(Mat*)obj); // TODO check u.getMat()
2638     }
2639     else
2640     {
2641         CV_Error(Error::StsNotImplemented, "");
2642     }
2643 }
2644
2645
2646 void _OutputArray::assign(const Mat& m) const
2647 {
2648     int k = kind();
2649     if (k == UMAT)
2650     {
2651         m.copyTo(*(UMat*)obj); // TODO check m.getUMat()
2652     }
2653     else if (k == MAT)
2654     {
2655         *(Mat*)obj = m;
2656     }
2657     else
2658     {
2659         CV_Error(Error::StsNotImplemented, "");
2660     }
2661 }
2662
2663
2664 static _InputOutputArray _none;
2665 InputOutputArray noArray() { return _none; }
2666
2667 }
2668
2669 /*************************************************************************************************\
2670                                         Matrix Operations
2671 \*************************************************************************************************/
2672
2673 void cv::hconcat(const Mat* src, size_t nsrc, OutputArray _dst)
2674 {
2675     if( nsrc == 0 || !src )
2676     {
2677         _dst.release();
2678         return;
2679     }
2680
2681     int totalCols = 0, cols = 0;
2682     size_t i;
2683     for( i = 0; i < nsrc; i++ )
2684     {
2685         CV_Assert( src[i].dims <= 2 &&
2686                    src[i].rows == src[0].rows &&
2687                    src[i].type() == src[0].type());
2688         totalCols += src[i].cols;
2689     }
2690     _dst.create( src[0].rows, totalCols, src[0].type());
2691     Mat dst = _dst.getMat();
2692     for( i = 0; i < nsrc; i++ )
2693     {
2694         Mat dpart = dst(Rect(cols, 0, src[i].cols, src[i].rows));
2695         src[i].copyTo(dpart);
2696         cols += src[i].cols;
2697     }
2698 }
2699
2700 void cv::hconcat(InputArray src1, InputArray src2, OutputArray dst)
2701 {
2702     Mat src[] = {src1.getMat(), src2.getMat()};
2703     hconcat(src, 2, dst);
2704 }
2705
2706 void cv::hconcat(InputArray _src, OutputArray dst)
2707 {
2708     std::vector<Mat> src;
2709     _src.getMatVector(src);
2710     hconcat(!src.empty() ? &src[0] : 0, src.size(), dst);
2711 }
2712
2713 void cv::vconcat(const Mat* src, size_t nsrc, OutputArray _dst)
2714 {
2715     if( nsrc == 0 || !src )
2716     {
2717         _dst.release();
2718         return;
2719     }
2720
2721     int totalRows = 0, rows = 0;
2722     size_t i;
2723     for( i = 0; i < nsrc; i++ )
2724     {
2725         CV_Assert(src[i].dims <= 2 &&
2726                   src[i].cols == src[0].cols &&
2727                   src[i].type() == src[0].type());
2728         totalRows += src[i].rows;
2729     }
2730     _dst.create( totalRows, src[0].cols, src[0].type());
2731     Mat dst = _dst.getMat();
2732     for( i = 0; i < nsrc; i++ )
2733     {
2734         Mat dpart(dst, Rect(0, rows, src[i].cols, src[i].rows));
2735         src[i].copyTo(dpart);
2736         rows += src[i].rows;
2737     }
2738 }
2739
2740 void cv::vconcat(InputArray src1, InputArray src2, OutputArray dst)
2741 {
2742     Mat src[] = {src1.getMat(), src2.getMat()};
2743     vconcat(src, 2, dst);
2744 }
2745
2746 void cv::vconcat(InputArray _src, OutputArray dst)
2747 {
2748     std::vector<Mat> src;
2749     _src.getMatVector(src);
2750     vconcat(!src.empty() ? &src[0] : 0, src.size(), dst);
2751 }
2752
2753 //////////////////////////////////////// set identity ////////////////////////////////////////////
2754
2755 #ifdef HAVE_OPENCL
2756
2757 namespace cv {
2758
2759 static bool ocl_setIdentity( InputOutputArray _m, const Scalar& s )
2760 {
2761     int type = _m.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type), kercn = cn;
2762     if (cn == 1)
2763     {
2764         kercn = std::min(ocl::predictOptimalVectorWidth(_m), 4);
2765         if (kercn != 4)
2766             kercn = 1;
2767     }
2768     int sctype = CV_MAKE_TYPE(depth, cn == 3 ? 4 : cn),
2769             rowsPerWI = ocl::Device::getDefault().isIntel() ? 4 : 1;
2770
2771     ocl::Kernel k("setIdentity", ocl::core::set_identity_oclsrc,
2772                   format("-D T=%s -D T1=%s -D cn=%d -D ST=%s -D kercn=%d -D rowsPerWI=%d",
2773                          ocl::memopTypeToStr(CV_MAKE_TYPE(depth, kercn)),
2774                          ocl::memopTypeToStr(depth), cn,
2775                          ocl::memopTypeToStr(sctype),
2776                          kercn, rowsPerWI));
2777     if (k.empty())
2778         return false;
2779
2780     UMat m = _m.getUMat();
2781     k.args(ocl::KernelArg::WriteOnly(m, cn, kercn),
2782            ocl::KernelArg::Constant(Mat(1, 1, sctype, s)));
2783
2784     size_t globalsize[2] = { m.cols * cn / kercn, (m.rows + rowsPerWI - 1) / rowsPerWI };
2785     return k.run(2, globalsize, NULL, false);
2786 }
2787
2788 }
2789
2790 #endif
2791
2792 void cv::setIdentity( InputOutputArray _m, const Scalar& s )
2793 {
2794     CV_Assert( _m.dims() <= 2 );
2795
2796     CV_OCL_RUN(_m.isUMat(),
2797                ocl_setIdentity(_m, s))
2798
2799     Mat m = _m.getMat();
2800     int i, j, rows = m.rows, cols = m.cols, type = m.type();
2801
2802     if( type == CV_32FC1 )
2803     {
2804         float* data = (float*)m.data;
2805         float val = (float)s[0];
2806         size_t step = m.step/sizeof(data[0]);
2807
2808         for( i = 0; i < rows; i++, data += step )
2809         {
2810             for( j = 0; j < cols; j++ )
2811                 data[j] = 0;
2812             if( i < cols )
2813                 data[i] = val;
2814         }
2815     }
2816     else if( type == CV_64FC1 )
2817     {
2818         double* data = (double*)m.data;
2819         double val = s[0];
2820         size_t step = m.step/sizeof(data[0]);
2821
2822         for( i = 0; i < rows; i++, data += step )
2823         {
2824             for( j = 0; j < cols; j++ )
2825                 data[j] = j == i ? val : 0;
2826         }
2827     }
2828     else
2829     {
2830         m = Scalar(0);
2831         m.diag() = s;
2832     }
2833 }
2834
2835 //////////////////////////////////////////// trace ///////////////////////////////////////////
2836
2837 cv::Scalar cv::trace( InputArray _m )
2838 {
2839     Mat m = _m.getMat();
2840     CV_Assert( m.dims <= 2 );
2841     int i, type = m.type();
2842     int nm = std::min(m.rows, m.cols);
2843
2844     if( type == CV_32FC1 )
2845     {
2846         const float* ptr = (const float*)m.data;
2847         size_t step = m.step/sizeof(ptr[0]) + 1;
2848         double _s = 0;
2849         for( i = 0; i < nm; i++ )
2850             _s += ptr[i*step];
2851         return _s;
2852     }
2853
2854     if( type == CV_64FC1 )
2855     {
2856         const double* ptr = (const double*)m.data;
2857         size_t step = m.step/sizeof(ptr[0]) + 1;
2858         double _s = 0;
2859         for( i = 0; i < nm; i++ )
2860             _s += ptr[i*step];
2861         return _s;
2862     }
2863
2864     return cv::sum(m.diag());
2865 }
2866
2867 ////////////////////////////////////// transpose /////////////////////////////////////////
2868
2869 namespace cv
2870 {
2871
2872 template<typename T> static void
2873 transpose_( const uchar* src, size_t sstep, uchar* dst, size_t dstep, Size sz )
2874 {
2875     int i=0, j, m = sz.width, n = sz.height;
2876
2877     #if CV_ENABLE_UNROLLED
2878     for(; i <= m - 4; i += 4 )
2879     {
2880         T* d0 = (T*)(dst + dstep*i);
2881         T* d1 = (T*)(dst + dstep*(i+1));
2882         T* d2 = (T*)(dst + dstep*(i+2));
2883         T* d3 = (T*)(dst + dstep*(i+3));
2884
2885         for( j = 0; j <= n - 4; j += 4 )
2886         {
2887             const T* s0 = (const T*)(src + i*sizeof(T) + sstep*j);
2888             const T* s1 = (const T*)(src + i*sizeof(T) + sstep*(j+1));
2889             const T* s2 = (const T*)(src + i*sizeof(T) + sstep*(j+2));
2890             const T* s3 = (const T*)(src + i*sizeof(T) + sstep*(j+3));
2891
2892             d0[j] = s0[0]; d0[j+1] = s1[0]; d0[j+2] = s2[0]; d0[j+3] = s3[0];
2893             d1[j] = s0[1]; d1[j+1] = s1[1]; d1[j+2] = s2[1]; d1[j+3] = s3[1];
2894             d2[j] = s0[2]; d2[j+1] = s1[2]; d2[j+2] = s2[2]; d2[j+3] = s3[2];
2895             d3[j] = s0[3]; d3[j+1] = s1[3]; d3[j+2] = s2[3]; d3[j+3] = s3[3];
2896         }
2897
2898         for( ; j < n; j++ )
2899         {
2900             const T* s0 = (const T*)(src + i*sizeof(T) + j*sstep);
2901             d0[j] = s0[0]; d1[j] = s0[1]; d2[j] = s0[2]; d3[j] = s0[3];
2902         }
2903     }
2904     #endif
2905     for( ; i < m; i++ )
2906     {
2907         T* d0 = (T*)(dst + dstep*i);
2908         j = 0;
2909         #if CV_ENABLE_UNROLLED
2910         for(; j <= n - 4; j += 4 )
2911         {
2912             const T* s0 = (const T*)(src + i*sizeof(T) + sstep*j);
2913             const T* s1 = (const T*)(src + i*sizeof(T) + sstep*(j+1));
2914             const T* s2 = (const T*)(src + i*sizeof(T) + sstep*(j+2));
2915             const T* s3 = (const T*)(src + i*sizeof(T) + sstep*(j+3));
2916
2917             d0[j] = s0[0]; d0[j+1] = s1[0]; d0[j+2] = s2[0]; d0[j+3] = s3[0];
2918         }
2919         #endif
2920         for( ; j < n; j++ )
2921         {
2922             const T* s0 = (const T*)(src + i*sizeof(T) + j*sstep);
2923             d0[j] = s0[0];
2924         }
2925     }
2926 }
2927
2928 template<typename T> static void
2929 transposeI_( uchar* data, size_t step, int n )
2930 {
2931     int i, j;
2932     for( i = 0; i < n; i++ )
2933     {
2934         T* row = (T*)(data + step*i);
2935         uchar* data1 = data + i*sizeof(T);
2936         for( j = i+1; j < n; j++ )
2937             std::swap( row[j], *(T*)(data1 + step*j) );
2938     }
2939 }
2940
2941 typedef void (*TransposeFunc)( const uchar* src, size_t sstep, uchar* dst, size_t dstep, Size sz );
2942 typedef void (*TransposeInplaceFunc)( uchar* data, size_t step, int n );
2943
2944 #define DEF_TRANSPOSE_FUNC(suffix, type) \
2945 static void transpose_##suffix( const uchar* src, size_t sstep, uchar* dst, size_t dstep, Size sz ) \
2946 { transpose_<type>(src, sstep, dst, dstep, sz); } \
2947 \
2948 static void transposeI_##suffix( uchar* data, size_t step, int n ) \
2949 { transposeI_<type>(data, step, n); }
2950
2951 DEF_TRANSPOSE_FUNC(8u, uchar)
2952 DEF_TRANSPOSE_FUNC(16u, ushort)
2953 DEF_TRANSPOSE_FUNC(8uC3, Vec3b)
2954 DEF_TRANSPOSE_FUNC(32s, int)
2955 DEF_TRANSPOSE_FUNC(16uC3, Vec3s)
2956 DEF_TRANSPOSE_FUNC(32sC2, Vec2i)
2957 DEF_TRANSPOSE_FUNC(32sC3, Vec3i)
2958 DEF_TRANSPOSE_FUNC(32sC4, Vec4i)
2959 DEF_TRANSPOSE_FUNC(32sC6, Vec6i)
2960 DEF_TRANSPOSE_FUNC(32sC8, Vec8i)
2961
2962 static TransposeFunc transposeTab[] =
2963 {
2964     0, transpose_8u, transpose_16u, transpose_8uC3, transpose_32s, 0, transpose_16uC3, 0,
2965     transpose_32sC2, 0, 0, 0, transpose_32sC3, 0, 0, 0, transpose_32sC4,
2966     0, 0, 0, 0, 0, 0, 0, transpose_32sC6, 0, 0, 0, 0, 0, 0, 0, transpose_32sC8
2967 };
2968
2969 static TransposeInplaceFunc transposeInplaceTab[] =
2970 {
2971     0, transposeI_8u, transposeI_16u, transposeI_8uC3, transposeI_32s, 0, transposeI_16uC3, 0,
2972     transposeI_32sC2, 0, 0, 0, transposeI_32sC3, 0, 0, 0, transposeI_32sC4,
2973     0, 0, 0, 0, 0, 0, 0, transposeI_32sC6, 0, 0, 0, 0, 0, 0, 0, transposeI_32sC8
2974 };
2975
2976 #ifdef HAVE_OPENCL
2977
2978 static inline int divUp(int a, int b)
2979 {
2980     return (a + b - 1) / b;
2981 }
2982
2983 static bool ocl_transpose( InputArray _src, OutputArray _dst )
2984 {
2985     const ocl::Device & dev = ocl::Device::getDefault();
2986     const int TILE_DIM = 32, BLOCK_ROWS = 8;
2987     int type = _src.type(), cn = CV_MAT_CN(type), depth = CV_MAT_DEPTH(type),
2988         rowsPerWI = dev.isIntel() ? 4 : 1;
2989
2990     UMat src = _src.getUMat();
2991     _dst.create(src.cols, src.rows, type);
2992     UMat dst = _dst.getUMat();
2993
2994     String kernelName("transpose");
2995     bool inplace = dst.u == src.u;
2996
2997     if (inplace)
2998     {
2999         CV_Assert(dst.cols == dst.rows);
3000         kernelName += "_inplace";
3001     }
3002
3003     ocl::Kernel k(kernelName.c_str(), ocl::core::transpose_oclsrc,
3004                   format("-D T=%s -D T1=%s -D cn=%d -D TILE_DIM=%d -D BLOCK_ROWS=%d -D rowsPerWI=%d",
3005                          ocl::memopTypeToStr(type), ocl::memopTypeToStr(depth),
3006                          cn, TILE_DIM, BLOCK_ROWS, rowsPerWI));
3007     if (k.empty())
3008         return false;
3009
3010     if (inplace)
3011         k.args(ocl::KernelArg::ReadWriteNoSize(dst), dst.rows);
3012     else
3013         k.args(ocl::KernelArg::ReadOnly(src),
3014                ocl::KernelArg::WriteOnlyNoSize(dst));
3015
3016     size_t localsize[2]  = { TILE_DIM, BLOCK_ROWS };
3017     size_t globalsize[2] = { src.cols, inplace ? (src.rows + rowsPerWI - 1) / rowsPerWI : (divUp(src.rows, TILE_DIM) * BLOCK_ROWS) };
3018
3019     if (inplace && dev.isIntel())
3020     {
3021         localsize[0] = 16;
3022         localsize[1] = dev.maxWorkGroupSize() / localsize[0];
3023     }
3024
3025     return k.run(2, globalsize, localsize, false);
3026 }
3027
3028 #endif
3029
3030 }
3031
3032 void cv::transpose( InputArray _src, OutputArray _dst )
3033 {
3034     int type = _src.type(), esz = CV_ELEM_SIZE(type);
3035     CV_Assert( _src.dims() <= 2 && esz <= 32 );
3036
3037     CV_OCL_RUN(_dst.isUMat(),
3038                ocl_transpose(_src, _dst))
3039
3040     Mat src = _src.getMat();
3041     if( src.empty() )
3042     {
3043         _dst.release();
3044         return;
3045     }
3046
3047     _dst.create(src.cols, src.rows, src.type());
3048     Mat dst = _dst.getMat();
3049
3050     // handle the case of single-column/single-row matrices, stored in STL vectors.
3051     if( src.rows != dst.cols || src.cols != dst.rows )
3052     {
3053         CV_Assert( src.size() == dst.size() && (src.cols == 1 || src.rows == 1) );
3054         src.copyTo(dst);
3055         return;
3056     }
3057
3058 #if defined HAVE_IPP
3059     typedef IppStatus (CV_STDCALL * ippiTranspose)(const void * pSrc, int srcStep, void * pDst, int dstStep, IppiSize roiSize);
3060     typedef IppStatus (CV_STDCALL * ippiTransposeI)(const void * pSrcDst, int srcDstStep, IppiSize roiSize);
3061     ippiTranspose ippFunc = 0;
3062     ippiTransposeI ippFuncI = 0;
3063
3064     if (dst.data == src.data && dst.cols == dst.rows)
3065     {
3066         CV_SUPPRESS_DEPRECATED_START
3067         ippFuncI =
3068             type == CV_8UC1 ? (ippiTransposeI)ippiTranspose_8u_C1IR :
3069             type == CV_8UC3 ? (ippiTransposeI)ippiTranspose_8u_C3IR :
3070             type == CV_8UC4 ? (ippiTransposeI)ippiTranspose_8u_C4IR :
3071             type == CV_16UC1 ? (ippiTransposeI)ippiTranspose_16u_C1IR :
3072             type == CV_16UC3 ? (ippiTransposeI)ippiTranspose_16u_C3IR :
3073             type == CV_16UC4 ? (ippiTransposeI)ippiTranspose_16u_C4IR :
3074             type == CV_16SC1 ? (ippiTransposeI)ippiTranspose_16s_C1IR :
3075             type == CV_16SC3 ? (ippiTransposeI)ippiTranspose_16s_C3IR :
3076             type == CV_16SC4 ? (ippiTransposeI)ippiTranspose_16s_C4IR :
3077             type == CV_32SC1 ? (ippiTransposeI)ippiTranspose_32s_C1IR :
3078             type == CV_32SC3 ? (ippiTransposeI)ippiTranspose_32s_C3IR :
3079             type == CV_32SC4 ? (ippiTransposeI)ippiTranspose_32s_C4IR :
3080             type == CV_32FC1 ? (ippiTransposeI)ippiTranspose_32f_C1IR :
3081             type == CV_32FC3 ? (ippiTransposeI)ippiTranspose_32f_C3IR :
3082             type == CV_32FC4 ? (ippiTransposeI)ippiTranspose_32f_C4IR : 0;
3083         CV_SUPPRESS_DEPRECATED_END
3084     }
3085     else
3086     {
3087         ippFunc =
3088             type == CV_8UC1 ? (ippiTranspose)ippiTranspose_8u_C1R :
3089             type == CV_8UC3 ? (ippiTranspose)ippiTranspose_8u_C3R :
3090             type == CV_8UC4 ? (ippiTranspose)ippiTranspose_8u_C4R :
3091             type == CV_16UC1 ? (ippiTranspose)ippiTranspose_16u_C1R :
3092             type == CV_16UC3 ? (ippiTranspose)ippiTranspose_16u_C3R :
3093             type == CV_16UC4 ? (ippiTranspose)ippiTranspose_16u_C4R :
3094             type == CV_16SC1 ? (ippiTranspose)ippiTranspose_16s_C1R :
3095             type == CV_16SC3 ? (ippiTranspose)ippiTranspose_16s_C3R :
3096             type == CV_16SC4 ? (ippiTranspose)ippiTranspose_16s_C4R :
3097             type == CV_32SC1 ? (ippiTranspose)ippiTranspose_32s_C1R :
3098             type == CV_32SC3 ? (ippiTranspose)ippiTranspose_32s_C3R :
3099             type == CV_32SC4 ? (ippiTranspose)ippiTranspose_32s_C4R :
3100             type == CV_32FC1 ? (ippiTranspose)ippiTranspose_32f_C1R :
3101             type == CV_32FC3 ? (ippiTranspose)ippiTranspose_32f_C3R :
3102             type == CV_32FC4 ? (ippiTranspose)ippiTranspose_32f_C4R : 0;
3103     }
3104
3105     IppiSize roiSize = { src.cols, src.rows };
3106     if (ippFunc != 0)
3107     {
3108         if (ippFunc(src.data, (int)src.step, dst.data, (int)dst.step, roiSize) >= 0)
3109             return;
3110         setIppErrorStatus();
3111     }
3112     else if (ippFuncI != 0)
3113     {
3114         if (ippFuncI(dst.data, (int)dst.step, roiSize) >= 0)
3115             return;
3116         setIppErrorStatus();
3117     }
3118 #endif
3119
3120     if( dst.data == src.data )
3121     {
3122         TransposeInplaceFunc func = transposeInplaceTab[esz];
3123         CV_Assert( func != 0 );
3124         CV_Assert( dst.cols == dst.rows );
3125         func( dst.data, dst.step, dst.rows );
3126     }
3127     else
3128     {
3129         TransposeFunc func = transposeTab[esz];
3130         CV_Assert( func != 0 );
3131         func( src.data, src.step, dst.data, dst.step, src.size() );
3132     }
3133 }
3134
3135
3136 ////////////////////////////////////// completeSymm /////////////////////////////////////////
3137
3138 void cv::completeSymm( InputOutputArray _m, bool LtoR )
3139 {
3140     Mat m = _m.getMat();
3141     size_t step = m.step, esz = m.elemSize();
3142     CV_Assert( m.dims <= 2 && m.rows == m.cols );
3143
3144     int rows = m.rows;
3145     int j0 = 0, j1 = rows;
3146
3147     uchar* data = m.data;
3148     for( int i = 0; i < rows; i++ )
3149     {
3150         if( !LtoR ) j1 = i; else j0 = i+1;
3151         for( int j = j0; j < j1; j++ )
3152             memcpy(data + (i*step + j*esz), data + (j*step + i*esz), esz);
3153     }
3154 }
3155
3156
3157 cv::Mat cv::Mat::cross(InputArray _m) const
3158 {
3159     Mat m = _m.getMat();
3160     int tp = type(), d = CV_MAT_DEPTH(tp);
3161     CV_Assert( dims <= 2 && m.dims <= 2 && size() == m.size() && tp == m.type() &&
3162         ((rows == 3 && cols == 1) || (cols*channels() == 3 && rows == 1)));
3163     Mat result(rows, cols, tp);
3164
3165     if( d == CV_32F )
3166     {
3167         const float *a = (const float*)data, *b = (const float*)m.data;
3168         float* c = (float*)result.data;
3169         size_t lda = rows > 1 ? step/sizeof(a[0]) : 1;
3170         size_t ldb = rows > 1 ? m.step/sizeof(b[0]) : 1;
3171
3172         c[0] = a[lda] * b[ldb*2] - a[lda*2] * b[ldb];
3173         c[1] = a[lda*2] * b[0] - a[0] * b[ldb*2];
3174         c[2] = a[0] * b[ldb] - a[lda] * b[0];
3175     }
3176     else if( d == CV_64F )
3177     {
3178         const double *a = (const double*)data, *b = (const double*)m.data;
3179         double* c = (double*)result.data;
3180         size_t lda = rows > 1 ? step/sizeof(a[0]) : 1;
3181         size_t ldb = rows > 1 ? m.step/sizeof(b[0]) : 1;
3182
3183         c[0] = a[lda] * b[ldb*2] - a[lda*2] * b[ldb];
3184         c[1] = a[lda*2] * b[0] - a[0] * b[ldb*2];
3185         c[2] = a[0] * b[ldb] - a[lda] * b[0];
3186     }
3187
3188     return result;
3189 }
3190
3191
3192 ////////////////////////////////////////// reduce ////////////////////////////////////////////
3193
3194 namespace cv
3195 {
3196
3197 template<typename T, typename ST, class Op> static void
3198 reduceR_( const Mat& srcmat, Mat& dstmat )
3199 {
3200     typedef typename Op::rtype WT;
3201     Size size = srcmat.size();
3202     size.width *= srcmat.channels();
3203     AutoBuffer<WT> buffer(size.width);
3204     WT* buf = buffer;
3205     ST* dst = (ST*)dstmat.data;
3206     const T* src = (const T*)srcmat.data;
3207     size_t srcstep = srcmat.step/sizeof(src[0]);
3208     int i;
3209     Op op;
3210
3211     for( i = 0; i < size.width; i++ )
3212         buf[i] = src[i];
3213
3214     for( ; --size.height; )
3215     {
3216         src += srcstep;
3217         i = 0;
3218         #if CV_ENABLE_UNROLLED
3219         for(; i <= size.width - 4; i += 4 )
3220         {
3221             WT s0, s1;
3222             s0 = op(buf[i], (WT)src[i]);
3223             s1 = op(buf[i+1], (WT)src[i+1]);
3224             buf[i] = s0; buf[i+1] = s1;
3225
3226             s0 = op(buf[i+2], (WT)src[i+2]);
3227             s1 = op(buf[i+3], (WT)src[i+3]);
3228             buf[i+2] = s0; buf[i+3] = s1;
3229         }
3230         #endif
3231         for( ; i < size.width; i++ )
3232             buf[i] = op(buf[i], (WT)src[i]);
3233     }
3234
3235     for( i = 0; i < size.width; i++ )
3236         dst[i] = (ST)buf[i];
3237 }
3238
3239
3240 template<typename T, typename ST, class Op> static void
3241 reduceC_( const Mat& srcmat, Mat& dstmat )
3242 {
3243     typedef typename Op::rtype WT;
3244     Size size = srcmat.size();
3245     int i, k, cn = srcmat.channels();
3246     size.width *= cn;
3247     Op op;
3248
3249     for( int y = 0; y < size.height; y++ )
3250     {
3251         const T* src = (const T*)(srcmat.data + srcmat.step*y);
3252         ST* dst = (ST*)(dstmat.data + dstmat.step*y);
3253         if( size.width == cn )
3254             for( k = 0; k < cn; k++ )
3255                 dst[k] = src[k];
3256         else
3257         {
3258             for( k = 0; k < cn; k++ )
3259             {
3260                 WT a0 = src[k], a1 = src[k+cn];
3261                 for( i = 2*cn; i <= size.width - 4*cn; i += 4*cn )
3262                 {
3263                     a0 = op(a0, (WT)src[i+k]);
3264                     a1 = op(a1, (WT)src[i+k+cn]);
3265                     a0 = op(a0, (WT)src[i+k+cn*2]);
3266                     a1 = op(a1, (WT)src[i+k+cn*3]);
3267                 }
3268
3269                 for( ; i < size.width; i += cn )
3270                 {
3271                     a0 = op(a0, (WT)src[i+k]);
3272                 }
3273                 a0 = op(a0, a1);
3274               dst[k] = (ST)a0;
3275             }
3276         }
3277     }
3278 }
3279
3280 typedef void (*ReduceFunc)( const Mat& src, Mat& dst );
3281
3282 }
3283
3284 #define reduceSumR8u32s  reduceR_<uchar, int,   OpAdd<int> >
3285 #define reduceSumR8u32f  reduceR_<uchar, float, OpAdd<int> >
3286 #define reduceSumR8u64f  reduceR_<uchar, double,OpAdd<int> >
3287 #define reduceSumR16u32f reduceR_<ushort,float, OpAdd<float> >
3288 #define reduceSumR16u64f reduceR_<ushort,double,OpAdd<double> >
3289 #define reduceSumR16s32f reduceR_<short, float, OpAdd<float> >
3290 #define reduceSumR16s64f reduceR_<short, double,OpAdd<double> >
3291 #define reduceSumR32f32f reduceR_<float, float, OpAdd<float> >
3292 #define reduceSumR32f64f reduceR_<float, double,OpAdd<double> >
3293 #define reduceSumR64f64f reduceR_<double,double,OpAdd<double> >
3294
3295 #define reduceMaxR8u  reduceR_<uchar, uchar, OpMax<uchar> >
3296 #define reduceMaxR16u reduceR_<ushort,ushort,OpMax<ushort> >
3297 #define reduceMaxR16s reduceR_<short, short, OpMax<short> >
3298 #define reduceMaxR32f reduceR_<float, float, OpMax<float> >
3299 #define reduceMaxR64f reduceR_<double,double,OpMax<double> >
3300
3301 #define reduceMinR8u  reduceR_<uchar, uchar, OpMin<uchar> >
3302 #define reduceMinR16u reduceR_<ushort,ushort,OpMin<ushort> >
3303 #define reduceMinR16s reduceR_<short, short, OpMin<short> >
3304 #define reduceMinR32f reduceR_<float, float, OpMin<float> >
3305 #define reduceMinR64f reduceR_<double,double,OpMin<double> >
3306
3307 #if IPP_VERSION_X100 > 0
3308
3309 static inline void reduceSumC_8u16u16s32f_64f(const cv::Mat& srcmat, cv::Mat& dstmat)
3310 {
3311     cv::Size size = srcmat.size();
3312     IppiSize roisize = { size.width, 1 };
3313     int sstep = (int)srcmat.step, stype = srcmat.type(),
3314             sdepth = CV_MAT_DEPTH(stype), ddepth = dstmat.depth();
3315
3316     typedef IppStatus (CV_STDCALL * ippiSum)(const void * pSrc, int srcStep, IppiSize roiSize, Ipp64f* pSum);
3317     typedef IppStatus (CV_STDCALL * ippiSumHint)(const void * pSrc, int srcStep, IppiSize roiSize, Ipp64f* pSum, IppHintAlgorithm hint);
3318     ippiSum ippFunc = 0;
3319     ippiSumHint ippFuncHint = 0;
3320     cv::ReduceFunc func = 0;
3321
3322     if (ddepth == CV_64F)
3323     {
3324         ippFunc =
3325             stype == CV_8UC1 ? (ippiSum)ippiSum_8u_C1R :
3326             stype == CV_8UC3 ? (ippiSum)ippiSum_8u_C3R :
3327             stype == CV_8UC4 ? (ippiSum)ippiSum_8u_C4R :
3328             stype == CV_16UC1 ? (ippiSum)ippiSum_16u_C1R :
3329             stype == CV_16UC3 ? (ippiSum)ippiSum_16u_C3R :
3330             stype == CV_16UC4 ? (ippiSum)ippiSum_16u_C4R :
3331             stype == CV_16SC1 ? (ippiSum)ippiSum_16s_C1R :
3332             stype == CV_16SC3 ? (ippiSum)ippiSum_16s_C3R :
3333             stype == CV_16SC4 ? (ippiSum)ippiSum_16s_C4R : 0;
3334         ippFuncHint =
3335             stype == CV_32FC1 ? (ippiSumHint)ippiSum_32f_C1R :
3336             stype == CV_32FC3 ? (ippiSumHint)ippiSum_32f_C3R :
3337             stype == CV_32FC4 ? (ippiSumHint)ippiSum_32f_C4R : 0;
3338         func =
3339             sdepth == CV_8U ? (cv::ReduceFunc)cv::reduceC_<uchar, double,   cv::OpAdd<double> > :
3340             sdepth == CV_16U ? (cv::ReduceFunc)cv::reduceC_<ushort, double,   cv::OpAdd<double> > :
3341             sdepth == CV_16S ? (cv::ReduceFunc)cv::reduceC_<short, double,   cv::OpAdd<double> > :
3342             sdepth == CV_32F ? (cv::ReduceFunc)cv::reduceC_<float, double,   cv::OpAdd<double> > : 0;
3343     }
3344     CV_Assert(!(ippFunc && ippFuncHint) && func);
3345
3346     if (ippFunc)
3347     {
3348         for (int y = 0; y < size.height; ++y)
3349             if (ippFunc(srcmat.data + sstep * y, sstep, roisize, dstmat.ptr<Ipp64f>(y)) < 0)
3350             {
3351                 setIppErrorStatus();
3352                 cv::Mat dstroi = dstmat.rowRange(y, y + 1);
3353                 func(srcmat.rowRange(y, y + 1), dstroi);
3354             }
3355         return;
3356     }
3357     else if (ippFuncHint)
3358     {
3359         for (int y = 0; y < size.height; ++y)
3360             if (ippFuncHint(srcmat.data + sstep * y, sstep, roisize, dstmat.ptr<Ipp64f>(y), ippAlgHintAccurate) < 0)
3361             {
3362                 setIppErrorStatus();
3363                 cv::Mat dstroi = dstmat.rowRange(y, y + 1);
3364                 func(srcmat.rowRange(y, y + 1), dstroi);
3365             }
3366         return;
3367     }
3368
3369     func(srcmat, dstmat);
3370 }
3371
3372 #endif
3373
3374 #define reduceSumC8u32s  reduceC_<uchar, int,   OpAdd<int> >
3375 #define reduceSumC8u32f  reduceC_<uchar, float, OpAdd<int> >
3376 #define reduceSumC16u32f reduceC_<ushort,float, OpAdd<float> >
3377 #define reduceSumC16s32f reduceC_<short, float, OpAdd<float> >
3378 #define reduceSumC32f32f reduceC_<float, float, OpAdd<float> >
3379 #define reduceSumC64f64f reduceC_<double,double,OpAdd<double> >
3380
3381 #if IPP_VERSION_X100 > 0
3382 #define reduceSumC8u64f  reduceSumC_8u16u16s32f_64f
3383 #define reduceSumC16u64f reduceSumC_8u16u16s32f_64f
3384 #define reduceSumC16s64f reduceSumC_8u16u16s32f_64f
3385 #define reduceSumC32f64f reduceSumC_8u16u16s32f_64f
3386 #else
3387 #define reduceSumC8u64f  reduceC_<uchar, double,OpAdd<int> >
3388 #define reduceSumC16u64f reduceC_<ushort,double,OpAdd<double> >
3389 #define reduceSumC16s64f reduceC_<short, double,OpAdd<double> >
3390 #define reduceSumC32f64f reduceC_<float, double,OpAdd<double> >
3391 #endif
3392
3393 #if IPP_VERSION_X100 > 0
3394 #define REDUCE_OP(favor, optype, type1, type2) \
3395 static inline void reduce##optype##C##favor(const cv::Mat& srcmat, cv::Mat& dstmat) \
3396 { \
3397     typedef Ipp##favor IppType; \
3398     cv::Size size = srcmat.size(); \
3399     IppiSize roisize = ippiSize(size.width, 1);\
3400     int sstep = (int)srcmat.step; \
3401      \
3402     if (srcmat.channels() == 1) \
3403     { \
3404         for (int y = 0; y < size.height; ++y) \
3405             if (ippi##optype##_##favor##_C1R(srcmat.ptr<IppType>(y), sstep, roisize, dstmat.ptr<IppType>(y)) < 0) \
3406             { \
3407                 setIppErrorStatus(); \
3408                 cv::Mat dstroi = dstmat.rowRange(y, y + 1); \
3409                 cv::reduceC_ < type1, type2, cv::Op##optype < type2 > >(srcmat.rowRange(y, y + 1), dstroi); \
3410             } \
3411         return; \
3412     } \
3413     cv::reduceC_ < type1, type2, cv::Op##optype < type2 > >(srcmat, dstmat); \
3414 }
3415 #endif
3416
3417 #if IPP_VERSION_X100 > 0
3418 REDUCE_OP(8u, Max, uchar, uchar)
3419 REDUCE_OP(16u, Max, ushort, ushort)
3420 REDUCE_OP(16s, Max, short, short)
3421 REDUCE_OP(32f, Max, float, float)
3422 #else
3423 #define reduceMaxC8u  reduceC_<uchar, uchar, OpMax<uchar> >
3424 #define reduceMaxC16u reduceC_<ushort,ushort,OpMax<ushort> >
3425 #define reduceMaxC16s reduceC_<short, short, OpMax<short> >
3426 #define reduceMaxC32f reduceC_<float, float, OpMax<float> >
3427 #endif
3428 #define reduceMaxC64f reduceC_<double,double,OpMax<double> >
3429
3430 #if IPP_VERSION_X100 > 0
3431 REDUCE_OP(8u, Min, uchar, uchar)
3432 REDUCE_OP(16u, Min, ushort, ushort)
3433 REDUCE_OP(16s, Min, short, short)
3434 REDUCE_OP(32f, Min, float, float)
3435 #else
3436 #define reduceMinC8u  reduceC_<uchar, uchar, OpMin<uchar> >
3437 #define reduceMinC16u reduceC_<ushort,ushort,OpMin<ushort> >
3438 #define reduceMinC16s reduceC_<short, short, OpMin<short> >
3439 #define reduceMinC32f reduceC_<float, float, OpMin<float> >
3440 #endif
3441 #define reduceMinC64f reduceC_<double,double,OpMin<double> >
3442
3443 #ifdef HAVE_OPENCL
3444
3445 namespace cv {
3446
3447 static bool ocl_reduce(InputArray _src, OutputArray _dst,
3448                        int dim, int op, int op0, int stype, int dtype)
3449 {
3450     const int min_opt_cols = 128, buf_cols = 32;
3451     int sdepth = CV_MAT_DEPTH(stype), cn = CV_MAT_CN(stype),
3452             ddepth = CV_MAT_DEPTH(dtype), ddepth0 = ddepth;
3453     const ocl::Device &defDev = ocl::Device::getDefault();
3454     bool doubleSupport = defDev.doubleFPConfig() > 0;
3455
3456     size_t wgs = defDev.maxWorkGroupSize();
3457     bool useOptimized = 1 == dim && _src.cols() > min_opt_cols && (wgs >= buf_cols);
3458
3459     if (!doubleSupport && (sdepth == CV_64F || ddepth == CV_64F))
3460         return false;
3461
3462     if ((op == CV_REDUCE_SUM && sdepth == CV_32F) || op == CV_REDUCE_MIN || op == CV_REDUCE_MAX)
3463         return false;
3464
3465     if (op == CV_REDUCE_AVG)
3466     {
3467         if (sdepth < CV_32S && ddepth < CV_32S)
3468             ddepth = CV_32S;
3469     }
3470
3471     const char * const ops[4] = { "OCL_CV_REDUCE_SUM", "OCL_CV_REDUCE_AVG",
3472                                   "OCL_CV_REDUCE_MAX", "OCL_CV_REDUCE_MIN" };
3473     int wdepth = std::max(ddepth, CV_32F);
3474     if (useOptimized)
3475     {
3476         size_t tileHeight = (size_t)(wgs / buf_cols);
3477         if (defDev.isIntel())
3478         {
3479             static const size_t maxItemInGroupCount = 16;
3480             tileHeight = min(tileHeight, defDev.localMemSize() / buf_cols / CV_ELEM_SIZE(CV_MAKETYPE(wdepth, cn)) / maxItemInGroupCount);
3481         }
3482         char cvt[3][40];
3483         cv::String build_opt = format("-D OP_REDUCE_PRE -D BUF_COLS=%d -D TILE_HEIGHT=%d -D %s -D dim=1"
3484                                             " -D cn=%d -D ddepth=%d"
3485                                             " -D srcT=%s -D bufT=%s -D dstT=%s"
3486                                             " -D convertToWT=%s -D convertToBufT=%s -D convertToDT=%s%s",
3487                                             buf_cols, tileHeight, ops[op], cn, ddepth,
3488                                             ocl::typeToStr(sdepth),
3489                                             ocl::typeToStr(ddepth),
3490                                             ocl::typeToStr(ddepth0),
3491                                             ocl::convertTypeStr(ddepth, wdepth, 1, cvt[0]),
3492                                             ocl::convertTypeStr(sdepth, ddepth, 1, cvt[1]),
3493                                             ocl::convertTypeStr(wdepth, ddepth0, 1, cvt[2]),
3494                                             doubleSupport ? " -D DOUBLE_SUPPORT" : "");
3495         ocl::Kernel k("reduce_horz_opt", ocl::core::reduce2_oclsrc, build_opt);
3496         if (k.empty())
3497             return false;
3498         UMat src = _src.getUMat();
3499         Size dsize(1, src.rows);
3500         _dst.create(dsize, dtype);
3501         UMat dst = _dst.getUMat();
3502
3503         if (op0 == CV_REDUCE_AVG)
3504             k.args(ocl::KernelArg::ReadOnly(src),
3505                       ocl::KernelArg::WriteOnlyNoSize(dst), 1.0f / src.cols);
3506         else
3507             k.args(ocl::KernelArg::ReadOnly(src),
3508                       ocl::KernelArg::WriteOnlyNoSize(dst));
3509
3510         size_t localSize[2] = { buf_cols, tileHeight};
3511         size_t globalSize[2] = { buf_cols, src.rows };
3512         return k.run(2, globalSize, localSize, false);
3513     }
3514     else
3515     {
3516         char cvt[2][40];
3517         cv::String build_opt = format("-D %s -D dim=%d -D cn=%d -D ddepth=%d"
3518                                       " -D srcT=%s -D dstT=%s -D dstT0=%s -D convertToWT=%s"
3519                                       " -D convertToDT=%s -D convertToDT0=%s%s",
3520                                       ops[op], dim, cn, ddepth, ocl::typeToStr(useOptimized ? ddepth : sdepth),
3521                                       ocl::typeToStr(ddepth), ocl::typeToStr(ddepth0),
3522                                       ocl::convertTypeStr(ddepth, wdepth, 1, cvt[0]),
3523                                       ocl::convertTypeStr(sdepth, ddepth, 1, cvt[0]),
3524                                       ocl::convertTypeStr(wdepth, ddepth0, 1, cvt[1]),
3525                                       doubleSupport ? " -D DOUBLE_SUPPORT" : "");
3526
3527         ocl::Kernel k("reduce", ocl::core::reduce2_oclsrc, build_opt);
3528         if (k.empty())
3529             return false;
3530
3531         UMat src = _src.getUMat();
3532         Size dsize(dim == 0 ? src.cols : 1, dim == 0 ? 1 : src.rows);
3533         _dst.create(dsize, dtype);
3534         UMat dst = _dst.getUMat();
3535
3536         ocl::KernelArg srcarg = ocl::KernelArg::ReadOnly(src),
3537                 temparg = ocl::KernelArg::WriteOnlyNoSize(dst);
3538
3539         if (op0 == CV_REDUCE_AVG)
3540             k.args(srcarg, temparg, 1.0f / (dim == 0 ? src.rows : src.cols));
3541         else
3542             k.args(srcarg, temparg);
3543
3544         size_t globalsize = std::max(dsize.width, dsize.height);
3545         return k.run(1, &globalsize, NULL, false);
3546     }
3547 }
3548
3549 }
3550
3551 #endif
3552
3553 void cv::reduce(InputArray _src, OutputArray _dst, int dim, int op, int dtype)
3554 {
3555     CV_Assert( _src.dims() <= 2 );
3556     int op0 = op;
3557     int stype = _src.type(), sdepth = CV_MAT_DEPTH(stype), cn = CV_MAT_CN(stype);
3558     if( dtype < 0 )
3559         dtype = _dst.fixedType() ? _dst.type() : stype;
3560     dtype = CV_MAKETYPE(dtype >= 0 ? dtype : stype, cn);
3561     int ddepth = CV_MAT_DEPTH(dtype);
3562
3563     CV_Assert( cn == CV_MAT_CN(dtype) );
3564     CV_Assert( op == CV_REDUCE_SUM || op == CV_REDUCE_MAX ||
3565                op == CV_REDUCE_MIN || op == CV_REDUCE_AVG );
3566
3567     CV_OCL_RUN(_dst.isUMat(),
3568                ocl_reduce(_src, _dst, dim, op, op0, stype, dtype))
3569
3570     Mat src = _src.getMat();
3571     _dst.create(dim == 0 ? 1 : src.rows, dim == 0 ? src.cols : 1, dtype);
3572     Mat dst = _dst.getMat(), temp = dst;
3573
3574     if( op == CV_REDUCE_AVG )
3575     {
3576         op = CV_REDUCE_SUM;
3577         if( sdepth < CV_32S && ddepth < CV_32S )
3578         {
3579             temp.create(dst.rows, dst.cols, CV_32SC(cn));
3580             ddepth = CV_32S;
3581         }
3582     }
3583
3584     ReduceFunc func = 0;
3585     if( dim == 0 )
3586     {
3587         if( op == CV_REDUCE_SUM )
3588         {
3589             if(sdepth == CV_8U && ddepth == CV_32S)
3590                 func = GET_OPTIMIZED(reduceSumR8u32s);
3591             else if(sdepth == CV_8U && ddepth == CV_32F)
3592                 func = reduceSumR8u32f;
3593             else if(sdepth == CV_8U && ddepth == CV_64F)
3594                 func = reduceSumR8u64f;
3595             else if(sdepth == CV_16U && ddepth == CV_32F)
3596                 func = reduceSumR16u32f;
3597             else if(sdepth == CV_16U && ddepth == CV_64F)
3598                 func = reduceSumR16u64f;
3599             else if(sdepth == CV_16S && ddepth == CV_32F)
3600                 func = reduceSumR16s32f;
3601             else if(sdepth == CV_16S && ddepth == CV_64F)
3602                 func = reduceSumR16s64f;
3603             else if(sdepth == CV_32F && ddepth == CV_32F)
3604                 func = GET_OPTIMIZED(reduceSumR32f32f);
3605             else if(sdepth == CV_32F && ddepth == CV_64F)
3606                 func = reduceSumR32f64f;
3607             else if(sdepth == CV_64F && ddepth == CV_64F)
3608                 func = reduceSumR64f64f;
3609         }
3610         else if(op == CV_REDUCE_MAX)
3611         {
3612             if(sdepth == CV_8U && ddepth == CV_8U)
3613                 func = GET_OPTIMIZED(reduceMaxR8u);
3614             else if(sdepth == CV_16U && ddepth == CV_16U)
3615                 func = reduceMaxR16u;
3616             else if(sdepth == CV_16S && ddepth == CV_16S)
3617                 func = reduceMaxR16s;
3618             else if(sdepth == CV_32F && ddepth == CV_32F)
3619                 func = GET_OPTIMIZED(reduceMaxR32f);
3620             else if(sdepth == CV_64F && ddepth == CV_64F)
3621                 func = reduceMaxR64f;
3622         }
3623         else if(op == CV_REDUCE_MIN)
3624         {
3625             if(sdepth == CV_8U && ddepth == CV_8U)
3626                 func = GET_OPTIMIZED(reduceMinR8u);
3627             else if(sdepth == CV_16U && ddepth == CV_16U)
3628                 func = reduceMinR16u;
3629             else if(sdepth == CV_16S && ddepth == CV_16S)
3630                 func = reduceMinR16s;
3631             else if(sdepth == CV_32F && ddepth == CV_32F)
3632                 func = GET_OPTIMIZED(reduceMinR32f);
3633             else if(sdepth == CV_64F && ddepth == CV_64F)
3634                 func = reduceMinR64f;
3635         }
3636     }
3637     else
3638     {
3639         if(op == CV_REDUCE_SUM)
3640         {
3641             if(sdepth == CV_8U && ddepth == CV_32S)
3642                 func = GET_OPTIMIZED(reduceSumC8u32s);
3643             else if(sdepth == CV_8U && ddepth == CV_32F)
3644                 func = reduceSumC8u32f;
3645             else if(sdepth == CV_8U && ddepth == CV_64F)
3646                 func = reduceSumC8u64f;
3647             else if(sdepth == CV_16U && ddepth == CV_32F)
3648                 func = reduceSumC16u32f;
3649             else if(sdepth == CV_16U && ddepth == CV_64F)
3650                 func = reduceSumC16u64f;
3651             else if(sdepth == CV_16S && ddepth == CV_32F)
3652                 func = reduceSumC16s32f;
3653             else if(sdepth == CV_16S && ddepth == CV_64F)
3654                 func = reduceSumC16s64f;
3655             else if(sdepth == CV_32F && ddepth == CV_32F)
3656                 func = GET_OPTIMIZED(reduceSumC32f32f);
3657             else if(sdepth == CV_32F && ddepth == CV_64F)
3658                 func = reduceSumC32f64f;
3659             else if(sdepth == CV_64F && ddepth == CV_64F)
3660                 func = reduceSumC64f64f;
3661         }
3662         else if(op == CV_REDUCE_MAX)
3663         {
3664             if(sdepth == CV_8U && ddepth == CV_8U)
3665                 func = GET_OPTIMIZED(reduceMaxC8u);
3666             else if(sdepth == CV_16U && ddepth == CV_16U)
3667                 func = reduceMaxC16u;
3668             else if(sdepth == CV_16S && ddepth == CV_16S)
3669                 func = reduceMaxC16s;
3670             else if(sdepth == CV_32F && ddepth == CV_32F)
3671                 func = GET_OPTIMIZED(reduceMaxC32f);
3672             else if(sdepth == CV_64F && ddepth == CV_64F)
3673                 func = reduceMaxC64f;
3674         }
3675         else if(op == CV_REDUCE_MIN)
3676         {
3677             if(sdepth == CV_8U && ddepth == CV_8U)
3678                 func = GET_OPTIMIZED(reduceMinC8u);
3679             else if(sdepth == CV_16U && ddepth == CV_16U)
3680                 func = reduceMinC16u;
3681             else if(sdepth == CV_16S && ddepth == CV_16S)
3682                 func = reduceMinC16s;
3683             else if(sdepth == CV_32F && ddepth == CV_32F)
3684                 func = GET_OPTIMIZED(reduceMinC32f);
3685             else if(sdepth == CV_64F && ddepth == CV_64F)
3686                 func = reduceMinC64f;
3687         }
3688     }
3689
3690     if( !func )
3691         CV_Error( CV_StsUnsupportedFormat,
3692                   "Unsupported combination of input and output array formats" );
3693
3694     func( src, temp );
3695
3696     if( op0 == CV_REDUCE_AVG )
3697         temp.convertTo(dst, dst.type(), 1./(dim == 0 ? src.rows : src.cols));
3698 }
3699
3700
3701 //////////////////////////////////////// sort ///////////////////////////////////////////
3702
3703 namespace cv
3704 {
3705
3706 #if IPP_VERSION_X100 > 0
3707 #define USE_IPP_SORT
3708
3709 typedef IppStatus (CV_STDCALL * IppSortFunc)(void *, int);
3710 typedef IppSortFunc IppFlipFunc;
3711
3712 static IppSortFunc getSortFunc(int depth, bool sortDescending)
3713 {
3714     if (!sortDescending)
3715         return depth == CV_8U ? (IppSortFunc)ippsSortAscend_8u_I :
3716             /*depth == CV_16U ? (IppSortFunc)ippsSortAscend_16u_I :
3717             depth == CV_16S ? (IppSortFunc)ippsSortAscend_16s_I :
3718             depth == CV_32S ? (IppSortFunc)ippsSortAscend_32s_I :
3719             depth == CV_32F ? (IppSortFunc)ippsSortAscend_32f_I :
3720             depth == CV_64F ? (IppSortFunc)ippsSortAscend_64f_I :*/ 0;
3721     else
3722         return depth == CV_8U ? (IppSortFunc)ippsSortDescend_8u_I :
3723             /*depth == CV_16U ? (IppSortFunc)ippsSortDescend_16u_I :
3724             depth == CV_16S ? (IppSortFunc)ippsSortDescend_16s_I :
3725             depth == CV_32S ? (IppSortFunc)ippsSortDescend_32s_I :
3726             depth == CV_32F ? (IppSortFunc)ippsSortDescend_32f_I :
3727             depth == CV_64F ? (IppSortFunc)ippsSortDescend_64f_I :*/ 0;
3728 }
3729
3730 static IppFlipFunc getFlipFunc(int depth)
3731 {
3732     CV_SUPPRESS_DEPRECATED_START
3733     return
3734             depth == CV_8U || depth == CV_8S ? (IppFlipFunc)ippsFlip_8u_I :
3735             depth == CV_16U || depth == CV_16S ? (IppFlipFunc)ippsFlip_16u_I :
3736             depth == CV_32S || depth == CV_32F ? (IppFlipFunc)ippsFlip_32f_I :
3737             depth == CV_64F ? (IppFlipFunc)ippsFlip_64f_I : 0;
3738     CV_SUPPRESS_DEPRECATED_END
3739 }
3740
3741
3742 #endif
3743
3744 template<typename T> static void sort_( const Mat& src, Mat& dst, int flags )
3745 {
3746     AutoBuffer<T> buf;
3747     T* bptr;
3748     int i, j, n, len;
3749     bool sortRows = (flags & 1) == CV_SORT_EVERY_ROW;
3750     bool inplace = src.data == dst.data;
3751     bool sortDescending = (flags & CV_SORT_DESCENDING) != 0;
3752
3753     if( sortRows )
3754         n = src.rows, len = src.cols;
3755     else
3756     {
3757         n = src.cols, len = src.rows;
3758         buf.allocate(len);
3759     }
3760     bptr = (T*)buf;
3761
3762 #ifdef USE_IPP_SORT
3763     int depth = src.depth();
3764     IppSortFunc ippSortFunc = getSortFunc(depth, sortDescending);
3765     IppFlipFunc ippFlipFunc = getFlipFunc(depth);
3766 #endif
3767
3768     for( i = 0; i < n; i++ )
3769     {
3770         T* ptr = bptr;
3771         if( sortRows )
3772         {
3773             T* dptr = (T*)(dst.data + dst.step*i);
3774             if( !inplace )
3775             {
3776                 const T* sptr = (const T*)(src.data + src.step*i);
3777                 memcpy(dptr, sptr, sizeof(T) * len);
3778             }
3779             ptr = dptr;
3780         }
3781         else
3782         {
3783             for( j = 0; j < len; j++ )
3784                 ptr[j] = ((const T*)(src.data + src.step*j))[i];
3785         }
3786
3787 #ifdef USE_IPP_SORT
3788         if (!ippSortFunc || ippSortFunc(ptr, len) < 0)
3789 #endif
3790         {
3791 #ifdef USE_IPP_SORT
3792             if (depth == CV_8U)
3793                 setIppErrorStatus();
3794 #endif
3795             std::sort( ptr, ptr + len );
3796             if( sortDescending )
3797             {
3798 #ifdef USE_IPP_SORT
3799                 if (!ippFlipFunc || ippFlipFunc(ptr, len) < 0)
3800 #endif
3801                 {
3802 #ifdef USE_IPP_SORT
3803                     setIppErrorStatus();
3804 #endif
3805                     for( j = 0; j < len/2; j++ )
3806                         std::swap(ptr[j], ptr[len-1-j]);
3807                 }
3808             }
3809         }
3810
3811         if( !sortRows )
3812             for( j = 0; j < len; j++ )
3813                 ((T*)(dst.data + dst.step*j))[i] = ptr[j];
3814     }
3815 }
3816
3817 template<typename _Tp> class LessThanIdx
3818 {
3819 public:
3820     LessThanIdx( const _Tp* _arr ) : arr(_arr) {}
3821     bool operator()(int a, int b) const { return arr[a] < arr[b]; }
3822     const _Tp* arr;
3823 };
3824
3825 #if defined USE_IPP_SORT && 0
3826
3827 typedef IppStatus (CV_STDCALL *IppSortIndexFunc)(void *, int *, int);
3828
3829 static IppSortIndexFunc getSortIndexFunc(int depth, bool sortDescending)
3830 {
3831     if (!sortDescending)
3832         return depth == CV_8U ? (IppSortIndexFunc)ippsSortIndexAscend_8u_I :
3833             depth == CV_16U ? (IppSortIndexFunc)ippsSortIndexAscend_16u_I :
3834             depth == CV_16S ? (IppSortIndexFunc)ippsSortIndexAscend_16s_I :
3835             depth == CV_32S ? (IppSortIndexFunc)ippsSortIndexAscend_32s_I :
3836             depth == CV_32F ? (IppSortIndexFunc)ippsSortIndexAscend_32f_I :
3837             depth == CV_64F ? (IppSortIndexFunc)ippsSortIndexAscend_64f_I : 0;
3838     else
3839         return depth == CV_8U ? (IppSortIndexFunc)ippsSortIndexDescend_8u_I :
3840             depth == CV_16U ? (IppSortIndexFunc)ippsSortIndexDescend_16u_I :
3841             depth == CV_16S ? (IppSortIndexFunc)ippsSortIndexDescend_16s_I :
3842             depth == CV_32S ? (IppSortIndexFunc)ippsSortIndexDescend_32s_I :
3843             depth == CV_32F ? (IppSortIndexFunc)ippsSortIndexDescend_32f_I :
3844             depth == CV_64F ? (IppSortIndexFunc)ippsSortIndexDescend_64f_I : 0;
3845 }
3846
3847 #endif
3848
3849 template<typename T> static void sortIdx_( const Mat& src, Mat& dst, int flags )
3850 {
3851     AutoBuffer<T> buf;
3852     AutoBuffer<int> ibuf;
3853     T* bptr;
3854     int* _iptr;
3855     int i, j, n, len;
3856     bool sortRows = (flags & 1) == CV_SORT_EVERY_ROW;
3857     bool sortDescending = (flags & CV_SORT_DESCENDING) != 0;
3858
3859     CV_Assert( src.data != dst.data );
3860
3861     if( sortRows )
3862         n = src.rows, len = src.cols;
3863     else
3864     {
3865         n = src.cols, len = src.rows;
3866         buf.allocate(len);
3867         ibuf.allocate(len);
3868     }
3869     bptr = (T*)buf;
3870     _iptr = (int*)ibuf;
3871
3872 #if defined USE_IPP_SORT && 0
3873     int depth = src.depth();
3874     IppSortIndexFunc ippFunc = getSortIndexFunc(depth, sortDescending);
3875     IppFlipFunc ippFlipFunc = getFlipFunc(depth);
3876 #endif
3877
3878     for( i = 0; i < n; i++ )
3879     {
3880         T* ptr = bptr;
3881         int* iptr = _iptr;
3882
3883         if( sortRows )
3884         {
3885             ptr = (T*)(src.data + src.step*i);
3886             iptr = (int*)(dst.data + dst.step*i);
3887         }
3888         else
3889         {
3890             for( j = 0; j < len; j++ )
3891                 ptr[j] = ((const T*)(src.data + src.step*j))[i];
3892         }
3893         for( j = 0; j < len; j++ )
3894             iptr[j] = j;
3895
3896 #if defined USE_IPP_SORT && 0
3897         if (sortRows || !ippFunc || ippFunc(ptr, iptr, len) < 0)
3898 #endif
3899         {
3900 #if defined USE_IPP_SORT && 0
3901             setIppErrorStatus();
3902 #endif
3903             std::sort( iptr, iptr + len, LessThanIdx<T>(ptr) );
3904             if( sortDescending )
3905             {
3906 #if defined USE_IPP_SORT && 0
3907                 if (!ippFlipFunc || ippFlipFunc(iptr, len) < 0)
3908 #endif
3909                 {
3910 #if defined USE_IPP_SORT && 0
3911                     setIppErrorStatus();
3912 #endif
3913                     for( j = 0; j < len/2; j++ )
3914                         std::swap(iptr[j], iptr[len-1-j]);
3915                 }
3916             }
3917         }
3918
3919         if( !sortRows )
3920             for( j = 0; j < len; j++ )
3921                 ((int*)(dst.data + dst.step*j))[i] = iptr[j];
3922     }
3923 }
3924
3925 typedef void (*SortFunc)(const Mat& src, Mat& dst, int flags);
3926
3927 }
3928
3929 void cv::sort( InputArray _src, OutputArray _dst, int flags )
3930 {
3931     static SortFunc tab[] =
3932     {
3933         sort_<uchar>, sort_<schar>, sort_<ushort>, sort_<short>,
3934         sort_<int>, sort_<float>, sort_<double>, 0
3935     };
3936     Mat src = _src.getMat();
3937     SortFunc func = tab[src.depth()];
3938     CV_Assert( src.dims <= 2 && src.channels() == 1 && func != 0 );
3939     _dst.create( src.size(), src.type() );
3940     Mat dst = _dst.getMat();
3941     func( src, dst, flags );
3942 }
3943
3944 void cv::sortIdx( InputArray _src, OutputArray _dst, int flags )
3945 {
3946     static SortFunc tab[] =
3947     {
3948         sortIdx_<uchar>, sortIdx_<schar>, sortIdx_<ushort>, sortIdx_<short>,
3949         sortIdx_<int>, sortIdx_<float>, sortIdx_<double>, 0
3950     };
3951     Mat src = _src.getMat();
3952     SortFunc func = tab[src.depth()];
3953     CV_Assert( src.dims <= 2 && src.channels() == 1 && func != 0 );
3954
3955     Mat dst = _dst.getMat();
3956     if( dst.data == src.data )
3957         _dst.release();
3958     _dst.create( src.size(), CV_32S );
3959     dst = _dst.getMat();
3960     func( src, dst, flags );
3961 }
3962
3963
3964 ////////////////////////////////////////// kmeans ////////////////////////////////////////////
3965
3966 namespace cv
3967 {
3968
3969 static void generateRandomCenter(const std::vector<Vec2f>& box, float* center, RNG& rng)
3970 {
3971     size_t j, dims = box.size();
3972     float margin = 1.f/dims;
3973     for( j = 0; j < dims; j++ )
3974         center[j] = ((float)rng*(1.f+margin*2.f)-margin)*(box[j][1] - box[j][0]) + box[j][0];
3975 }
3976
3977 class KMeansPPDistanceComputer : public ParallelLoopBody
3978 {
3979 public:
3980     KMeansPPDistanceComputer( float *_tdist2,
3981                               const float *_data,
3982                               const float *_dist,
3983                               int _dims,
3984                               size_t _step,
3985                               size_t _stepci )
3986         : tdist2(_tdist2),
3987           data(_data),
3988           dist(_dist),
3989           dims(_dims),
3990           step(_step),
3991           stepci(_stepci) { }
3992
3993     void operator()( const cv::Range& range ) const
3994     {
3995         const int begin = range.start;
3996         const int end = range.end;
3997
3998         for ( int i = begin; i<end; i++ )
3999         {
4000             tdist2[i] = std::min(normL2Sqr_(data + step*i, data + stepci, dims), dist[i]);
4001         }
4002     }
4003
4004 private:
4005     KMeansPPDistanceComputer& operator=(const KMeansPPDistanceComputer&); // to quiet MSVC
4006
4007     float *tdist2;
4008     const float *data;
4009     const float *dist;
4010     const int dims;
4011     const size_t step;
4012     const size_t stepci;
4013 };
4014
4015 /*
4016 k-means center initialization using the following algorithm:
4017 Arthur & Vassilvitskii (2007) k-means++: The Advantages of Careful Seeding
4018 */
4019 static void generateCentersPP(const Mat& _data, Mat& _out_centers,
4020                               int K, RNG& rng, int trials)
4021 {
4022     int i, j, k, dims = _data.cols, N = _data.rows;
4023     const float* data = _data.ptr<float>(0);
4024     size_t step = _data.step/sizeof(data[0]);
4025     std::vector<int> _centers(K);
4026     int* centers = &_centers[0];
4027     std::vector<float> _dist(N*3);
4028     float* dist = &_dist[0], *tdist = dist + N, *tdist2 = tdist + N;
4029     double sum0 = 0;
4030
4031     centers[0] = (unsigned)rng % N;
4032
4033     for( i = 0; i < N; i++ )
4034     {
4035         dist[i] = normL2Sqr_(data + step*i, data + step*centers[0], dims);
4036         sum0 += dist[i];
4037     }
4038
4039     for( k = 1; k < K; k++ )
4040     {
4041         double bestSum = DBL_MAX;
4042         int bestCenter = -1;
4043
4044         for( j = 0; j < trials; j++ )
4045         {
4046             double p = (double)rng*sum0, s = 0;
4047             for( i = 0; i < N-1; i++ )
4048                 if( (p -= dist[i]) <= 0 )
4049                     break;
4050             int ci = i;
4051
4052             parallel_for_(Range(0, N),
4053                          KMeansPPDistanceComputer(tdist2, data, dist, dims, step, step*ci));
4054             for( i = 0; i < N; i++ )
4055             {
4056                 s += tdist2[i];
4057             }
4058
4059             if( s < bestSum )
4060             {
4061                 bestSum = s;
4062                 bestCenter = ci;
4063                 std::swap(tdist, tdist2);
4064             }
4065         }
4066         centers[k] = bestCenter;
4067         sum0 = bestSum;
4068         std::swap(dist, tdist);
4069     }
4070
4071     for( k = 0; k < K; k++ )
4072     {
4073         const float* src = data + step*centers[k];
4074         float* dst = _out_centers.ptr<float>(k);
4075         for( j = 0; j < dims; j++ )
4076             dst[j] = src[j];
4077     }
4078 }
4079
4080 class KMeansDistanceComputer : public ParallelLoopBody
4081 {
4082 public:
4083     KMeansDistanceComputer( double *_distances,
4084                             int *_labels,
4085                             const Mat& _data,
4086                             const Mat& _centers )
4087         : distances(_distances),
4088           labels(_labels),
4089           data(_data),
4090           centers(_centers)
4091     {
4092     }
4093
4094     void operator()( const Range& range ) const
4095     {
4096         const int begin = range.start;
4097         const int end = range.end;
4098         const int K = centers.rows;
4099         const int dims = centers.cols;
4100
4101         const float *sample;
4102         for( int i = begin; i<end; ++i)
4103         {
4104             sample = data.ptr<float>(i);
4105             int k_best = 0;
4106             double min_dist = DBL_MAX;
4107
4108             for( int k = 0; k < K; k++ )
4109             {
4110                 const float* center = centers.ptr<float>(k);
4111                 const double dist = normL2Sqr_(sample, center, dims);
4112
4113                 if( min_dist > dist )
4114                 {
4115                     min_dist = dist;
4116                     k_best = k;
4117                 }
4118             }
4119
4120             distances[i] = min_dist;
4121             labels[i] = k_best;
4122         }
4123     }
4124
4125 private:
4126     KMeansDistanceComputer& operator=(const KMeansDistanceComputer&); // to quiet MSVC
4127
4128     double *distances;
4129     int *labels;
4130     const Mat& data;
4131     const Mat& centers;
4132 };
4133
4134 }
4135
4136 double cv::kmeans( InputArray _data, int K,
4137                    InputOutputArray _bestLabels,
4138                    TermCriteria criteria, int attempts,
4139                    int flags, OutputArray _centers )
4140 {
4141     const int SPP_TRIALS = 3;
4142     Mat data0 = _data.getMat();
4143     bool isrow = data0.rows == 1 && data0.channels() > 1;
4144     int N = !isrow ? data0.rows : data0.cols;
4145     int dims = (!isrow ? data0.cols : 1)*data0.channels();
4146     int type = data0.depth();
4147
4148     attempts = std::max(attempts, 1);
4149     CV_Assert( data0.dims <= 2 && type == CV_32F && K > 0 );
4150     CV_Assert( N >= K );
4151
4152     Mat data(N, dims, CV_32F, data0.data, isrow ? dims * sizeof(float) : static_cast<size_t>(data0.step));
4153
4154     _bestLabels.create(N, 1, CV_32S, -1, true);
4155
4156     Mat _labels, best_labels = _bestLabels.getMat();
4157     if( flags & CV_KMEANS_USE_INITIAL_LABELS )
4158     {
4159         CV_Assert( (best_labels.cols == 1 || best_labels.rows == 1) &&
4160                   best_labels.cols*best_labels.rows == N &&
4161                   best_labels.type() == CV_32S &&
4162                   best_labels.isContinuous());
4163         best_labels.copyTo(_labels);
4164     }
4165     else
4166     {
4167         if( !((best_labels.cols == 1 || best_labels.rows == 1) &&
4168              best_labels.cols*best_labels.rows == N &&
4169             best_labels.type() == CV_32S &&
4170             best_labels.isContinuous()))
4171             best_labels.create(N, 1, CV_32S);
4172         _labels.create(best_labels.size(), best_labels.type());
4173     }
4174     int* labels = _labels.ptr<int>();
4175
4176     Mat centers(K, dims, type), old_centers(K, dims, type), temp(1, dims, type);
4177     std::vector<int> counters(K);
4178     std::vector<Vec2f> _box(dims);
4179     Vec2f* box = &_box[0];
4180     double best_compactness = DBL_MAX, compactness = 0;
4181     RNG& rng = theRNG();
4182     int a, iter, i, j, k;
4183
4184     if( criteria.type & TermCriteria::EPS )
4185         criteria.epsilon = std::max(criteria.epsilon, 0.);
4186     else
4187         criteria.epsilon = FLT_EPSILON;
4188     criteria.epsilon *= criteria.epsilon;
4189
4190     if( criteria.type & TermCriteria::COUNT )
4191         criteria.maxCount = std::min(std::max(criteria.maxCount, 2), 100);
4192     else
4193         criteria.maxCount = 100;
4194
4195     if( K == 1 )
4196     {
4197         attempts = 1;
4198         criteria.maxCount = 2;
4199     }
4200
4201     const float* sample = data.ptr<float>(0);
4202     for( j = 0; j < dims; j++ )
4203         box[j] = Vec2f(sample[j], sample[j]);
4204
4205     for( i = 1; i < N; i++ )
4206     {
4207         sample = data.ptr<float>(i);
4208         for( j = 0; j < dims; j++ )
4209         {
4210             float v = sample[j];
4211             box[j][0] = std::min(box[j][0], v);
4212             box[j][1] = std::max(box[j][1], v);
4213         }
4214     }
4215
4216     for( a = 0; a < attempts; a++ )
4217     {
4218         double max_center_shift = DBL_MAX;
4219         for( iter = 0;; )
4220         {
4221             swap(centers, old_centers);
4222
4223             if( iter == 0 && (a > 0 || !(flags & KMEANS_USE_INITIAL_LABELS)) )
4224             {
4225                 if( flags & KMEANS_PP_CENTERS )
4226                     generateCentersPP(data, centers, K, rng, SPP_TRIALS);
4227                 else
4228                 {
4229                     for( k = 0; k < K; k++ )
4230                         generateRandomCenter(_box, centers.ptr<float>(k), rng);
4231                 }
4232             }
4233             else
4234             {
4235                 if( iter == 0 && a == 0 && (flags & KMEANS_USE_INITIAL_LABELS) )
4236                 {
4237                     for( i = 0; i < N; i++ )
4238                         CV_Assert( (unsigned)labels[i] < (unsigned)K );
4239                 }
4240
4241                 // compute centers
4242                 centers = Scalar(0);
4243                 for( k = 0; k < K; k++ )
4244                     counters[k] = 0;
4245
4246                 for( i = 0; i < N; i++ )
4247                 {
4248                     sample = data.ptr<float>(i);
4249                     k = labels[i];
4250                     float* center = centers.ptr<float>(k);
4251                     j=0;
4252                     #if CV_ENABLE_UNROLLED
4253                     for(; j <= dims - 4; j += 4 )
4254                     {
4255                         float t0 = center[j] + sample[j];
4256                         float t1 = center[j+1] + sample[j+1];
4257
4258                         center[j] = t0;
4259                         center[j+1] = t1;
4260
4261                         t0 = center[j+2] + sample[j+2];
4262                         t1 = center[j+3] + sample[j+3];
4263
4264                         center[j+2] = t0;
4265                         center[j+3] = t1;
4266                     }
4267                     #endif
4268                     for( ; j < dims; j++ )
4269                         center[j] += sample[j];
4270                     counters[k]++;
4271                 }
4272
4273                 if( iter > 0 )
4274                     max_center_shift = 0;
4275
4276                 for( k = 0; k < K; k++ )
4277                 {
4278                     if( counters[k] != 0 )
4279                         continue;
4280
4281                     // if some cluster appeared to be empty then:
4282                     //   1. find the biggest cluster
4283                     //   2. find the farthest from the center point in the biggest cluster
4284                     //   3. exclude the farthest point from the biggest cluster and form a new 1-point cluster.
4285                     int max_k = 0;
4286                     for( int k1 = 1; k1 < K; k1++ )
4287                     {
4288                         if( counters[max_k] < counters[k1] )
4289                             max_k = k1;
4290                     }
4291
4292                     double max_dist = 0;
4293                     int farthest_i = -1;
4294                     float* new_center = centers.ptr<float>(k);
4295                     float* old_center = centers.ptr<float>(max_k);
4296                     float* _old_center = temp.ptr<float>(); // normalized
4297                     float scale = 1.f/counters[max_k];
4298                     for( j = 0; j < dims; j++ )
4299                         _old_center[j] = old_center[j]*scale;
4300
4301                     for( i = 0; i < N; i++ )
4302                     {
4303                         if( labels[i] != max_k )
4304                             continue;
4305                         sample = data.ptr<float>(i);
4306                         double dist = normL2Sqr_(sample, _old_center, dims);
4307
4308                         if( max_dist <= dist )
4309                         {
4310                             max_dist = dist;
4311                             farthest_i = i;
4312                         }
4313                     }
4314
4315                     counters[max_k]--;
4316                     counters[k]++;
4317                     labels[farthest_i] = k;
4318                     sample = data.ptr<float>(farthest_i);
4319
4320                     for( j = 0; j < dims; j++ )
4321                     {
4322                         old_center[j] -= sample[j];
4323                         new_center[j] += sample[j];
4324                     }
4325                 }
4326
4327                 for( k = 0; k < K; k++ )
4328                 {
4329                     float* center = centers.ptr<float>(k);
4330                     CV_Assert( counters[k] != 0 );
4331
4332                     float scale = 1.f/counters[k];
4333                     for( j = 0; j < dims; j++ )
4334                         center[j] *= scale;
4335
4336                     if( iter > 0 )
4337                     {
4338                         double dist = 0;
4339                         const float* old_center = old_centers.ptr<float>(k);
4340                         for( j = 0; j < dims; j++ )
4341                         {
4342                             double t = center[j] - old_center[j];
4343                             dist += t*t;
4344                         }
4345                         max_center_shift = std::max(max_center_shift, dist);
4346                     }
4347                 }
4348             }
4349
4350             if( ++iter == MAX(criteria.maxCount, 2) || max_center_shift <= criteria.epsilon )
4351                 break;
4352
4353             // assign labels
4354             Mat dists(1, N, CV_64F);
4355             double* dist = dists.ptr<double>(0);
4356             parallel_for_(Range(0, N),
4357                          KMeansDistanceComputer(dist, labels, data, centers));
4358             compactness = 0;
4359             for( i = 0; i < N; i++ )
4360             {
4361                 compactness += dist[i];
4362             }
4363         }
4364
4365         if( compactness < best_compactness )
4366         {
4367             best_compactness = compactness;
4368             if( _centers.needed() )
4369                 centers.copyTo(_centers);
4370             _labels.copyTo(best_labels);
4371         }
4372     }
4373
4374     return best_compactness;
4375 }
4376
4377
4378 CV_IMPL void cvSetIdentity( CvArr* arr, CvScalar value )
4379 {
4380     cv::Mat m = cv::cvarrToMat(arr);
4381     cv::setIdentity(m, value);
4382 }
4383
4384
4385 CV_IMPL CvScalar cvTrace( const CvArr* arr )
4386 {
4387     return cv::trace(cv::cvarrToMat(arr));
4388 }
4389
4390
4391 CV_IMPL void cvTranspose( const CvArr* srcarr, CvArr* dstarr )
4392 {
4393     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
4394
4395     CV_Assert( src.rows == dst.cols && src.cols == dst.rows && src.type() == dst.type() );
4396     transpose( src, dst );
4397 }
4398
4399
4400 CV_IMPL void cvCompleteSymm( CvMat* matrix, int LtoR )
4401 {
4402     cv::Mat m = cv::cvarrToMat(matrix);
4403     cv::completeSymm( m, LtoR != 0 );
4404 }
4405
4406
4407 CV_IMPL void cvCrossProduct( const CvArr* srcAarr, const CvArr* srcBarr, CvArr* dstarr )
4408 {
4409     cv::Mat srcA = cv::cvarrToMat(srcAarr), dst = cv::cvarrToMat(dstarr);
4410
4411     CV_Assert( srcA.size() == dst.size() && srcA.type() == dst.type() );
4412     srcA.cross(cv::cvarrToMat(srcBarr)).copyTo(dst);
4413 }
4414
4415
4416 CV_IMPL void
4417 cvReduce( const CvArr* srcarr, CvArr* dstarr, int dim, int op )
4418 {
4419     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
4420
4421     if( dim < 0 )
4422         dim = src.rows > dst.rows ? 0 : src.cols > dst.cols ? 1 : dst.cols == 1;
4423
4424     if( dim > 1 )
4425         CV_Error( CV_StsOutOfRange, "The reduced dimensionality index is out of range" );
4426
4427     if( (dim == 0 && (dst.cols != src.cols || dst.rows != 1)) ||
4428         (dim == 1 && (dst.rows != src.rows || dst.cols != 1)) )
4429         CV_Error( CV_StsBadSize, "The output array size is incorrect" );
4430
4431     if( src.channels() != dst.channels() )
4432         CV_Error( CV_StsUnmatchedFormats, "Input and output arrays must have the same number of channels" );
4433
4434     cv::reduce(src, dst, dim, op, dst.type());
4435 }
4436
4437
4438 CV_IMPL CvArr*
4439 cvRange( CvArr* arr, double start, double end )
4440 {
4441     int ok = 0;
4442
4443     CvMat stub, *mat = (CvMat*)arr;
4444     double delta;
4445     int type, step;
4446     double val = start;
4447     int i, j;
4448     int rows, cols;
4449
4450     if( !CV_IS_MAT(mat) )
4451         mat = cvGetMat( mat, &stub);
4452
4453     rows = mat->rows;
4454     cols = mat->cols;
4455     type = CV_MAT_TYPE(mat->type);
4456     delta = (end-start)/(rows*cols);
4457
4458     if( CV_IS_MAT_CONT(mat->type) )
4459     {
4460         cols *= rows;
4461         rows = 1;
4462         step = 1;
4463     }
4464     else
4465         step = mat->step / CV_ELEM_SIZE(type);
4466
4467     if( type == CV_32SC1 )
4468     {
4469         int* idata = mat->data.i;
4470         int ival = cvRound(val), idelta = cvRound(delta);
4471
4472         if( fabs(val - ival) < DBL_EPSILON &&
4473             fabs(delta - idelta) < DBL_EPSILON )
4474         {
4475             for( i = 0; i < rows; i++, idata += step )
4476                 for( j = 0; j < cols; j++, ival += idelta )
4477                     idata[j] = ival;
4478         }
4479         else
4480         {
4481             for( i = 0; i < rows; i++, idata += step )
4482                 for( j = 0; j < cols; j++, val += delta )
4483                     idata[j] = cvRound(val);
4484         }
4485     }
4486     else if( type == CV_32FC1 )
4487     {
4488         float* fdata = mat->data.fl;
4489         for( i = 0; i < rows; i++, fdata += step )
4490             for( j = 0; j < cols; j++, val += delta )
4491                 fdata[j] = (float)val;
4492     }
4493     else
4494         CV_Error( CV_StsUnsupportedFormat, "The function only supports 32sC1 and 32fC1 datatypes" );
4495
4496     ok = 1;
4497     return ok ? arr : 0;
4498 }
4499
4500
4501 CV_IMPL void
4502 cvSort( const CvArr* _src, CvArr* _dst, CvArr* _idx, int flags )
4503 {
4504     cv::Mat src = cv::cvarrToMat(_src);
4505
4506     if( _idx )
4507     {
4508         cv::Mat idx0 = cv::cvarrToMat(_idx), idx = idx0;
4509         CV_Assert( src.size() == idx.size() && idx.type() == CV_32S && src.data != idx.data );
4510         cv::sortIdx( src, idx, flags );
4511         CV_Assert( idx0.data == idx.data );
4512     }
4513
4514     if( _dst )
4515     {
4516         cv::Mat dst0 = cv::cvarrToMat(_dst), dst = dst0;
4517         CV_Assert( src.size() == dst.size() && src.type() == dst.type() );
4518         cv::sort( src, dst, flags );
4519         CV_Assert( dst0.data == dst.data );
4520     }
4521 }
4522
4523
4524 CV_IMPL int
4525 cvKMeans2( const CvArr* _samples, int cluster_count, CvArr* _labels,
4526            CvTermCriteria termcrit, int attempts, CvRNG*,
4527            int flags, CvArr* _centers, double* _compactness )
4528 {
4529     cv::Mat data = cv::cvarrToMat(_samples), labels = cv::cvarrToMat(_labels), centers;
4530     if( _centers )
4531     {
4532         centers = cv::cvarrToMat(_centers);
4533
4534         centers = centers.reshape(1);
4535         data = data.reshape(1);
4536
4537         CV_Assert( !centers.empty() );
4538         CV_Assert( centers.rows == cluster_count );
4539         CV_Assert( centers.cols == data.cols );
4540         CV_Assert( centers.depth() == data.depth() );
4541     }
4542     CV_Assert( labels.isContinuous() && labels.type() == CV_32S &&
4543         (labels.cols == 1 || labels.rows == 1) &&
4544         labels.cols + labels.rows - 1 == data.rows );
4545
4546     double compactness = cv::kmeans(data, cluster_count, labels, termcrit, attempts,
4547                                     flags, _centers ? cv::_OutputArray(centers) : cv::_OutputArray() );
4548     if( _compactness )
4549         *_compactness = compactness;
4550     return 1;
4551 }
4552
4553 ///////////////////////////// n-dimensional matrices ////////////////////////////
4554
4555 namespace cv
4556 {
4557
4558 Mat Mat::reshape(int _cn, int _newndims, const int* _newsz) const
4559 {
4560     if(_newndims == dims)
4561     {
4562         if(_newsz == 0)
4563             return reshape(_cn);
4564         if(_newndims == 2)
4565             return reshape(_cn, _newsz[0]);
4566     }
4567
4568     CV_Error(CV_StsNotImplemented, "");
4569     // TBD
4570     return Mat();
4571 }
4572
4573 NAryMatIterator::NAryMatIterator()
4574     : arrays(0), planes(0), ptrs(0), narrays(0), nplanes(0), size(0), iterdepth(0), idx(0)
4575 {
4576 }
4577
4578 NAryMatIterator::NAryMatIterator(const Mat** _arrays, Mat* _planes, int _narrays)
4579 : arrays(0), planes(0), ptrs(0), narrays(0), nplanes(0), size(0), iterdepth(0), idx(0)
4580 {
4581     init(_arrays, _planes, 0, _narrays);
4582 }
4583
4584 NAryMatIterator::NAryMatIterator(const Mat** _arrays, uchar** _ptrs, int _narrays)
4585     : arrays(0), planes(0), ptrs(0), narrays(0), nplanes(0), size(0), iterdepth(0), idx(0)
4586 {
4587     init(_arrays, 0, _ptrs, _narrays);
4588 }
4589
4590 void NAryMatIterator::init(const Mat** _arrays, Mat* _planes, uchar** _ptrs, int _narrays)
4591 {
4592     CV_Assert( _arrays && (_ptrs || _planes) );
4593     int i, j, d1=0, i0 = -1, d = -1;
4594
4595     arrays = _arrays;
4596     ptrs = _ptrs;
4597     planes = _planes;
4598     narrays = _narrays;
4599     nplanes = 0;
4600     size = 0;
4601
4602     if( narrays < 0 )
4603     {
4604         for( i = 0; _arrays[i] != 0; i++ )
4605             ;
4606         narrays = i;
4607         CV_Assert(narrays <= 1000);
4608     }
4609
4610     iterdepth = 0;
4611
4612     for( i = 0; i < narrays; i++ )
4613     {
4614         CV_Assert(arrays[i] != 0);
4615         const Mat& A = *arrays[i];
4616         if( ptrs )
4617             ptrs[i] = A.data;
4618
4619         if( !A.data )
4620             continue;
4621
4622         if( i0 < 0 )
4623         {
4624             i0 = i;
4625             d = A.dims;
4626
4627             // find the first dimensionality which is different from 1;
4628             // in any of the arrays the first "d1" step do not affect the continuity
4629             for( d1 = 0; d1 < d; d1++ )
4630                 if( A.size[d1] > 1 )
4631                     break;
4632         }
4633         else
4634             CV_Assert( A.size == arrays[i0]->size );
4635
4636         if( !A.isContinuous() )
4637         {
4638             CV_Assert( A.step[d-1] == A.elemSize() );
4639             for( j = d-1; j > d1; j-- )
4640                 if( A.step[j]*A.size[j] < A.step[j-1] )
4641                     break;
4642             iterdepth = std::max(iterdepth, j);
4643         }
4644     }
4645
4646     if( i0 >= 0 )
4647     {
4648         size = arrays[i0]->size[d-1];
4649         for( j = d-1; j > iterdepth; j-- )
4650         {
4651             int64 total1 = (int64)size*arrays[i0]->size[j-1];
4652             if( total1 != (int)total1 )
4653                 break;
4654             size = (int)total1;
4655         }
4656
4657         iterdepth = j;
4658         if( iterdepth == d1 )
4659             iterdepth = 0;
4660
4661         nplanes = 1;
4662         for( j = iterdepth-1; j >= 0; j-- )
4663             nplanes *= arrays[i0]->size[j];
4664     }
4665     else
4666         iterdepth = 0;
4667
4668     idx = 0;
4669
4670     if( !planes )
4671         return;
4672
4673     for( i = 0; i < narrays; i++ )
4674     {
4675         CV_Assert(arrays[i] != 0);
4676         const Mat& A = *arrays[i];
4677
4678         if( !A.data )
4679         {
4680             planes[i] = Mat();
4681             continue;
4682         }
4683
4684         planes[i] = Mat(1, (int)size, A.type(), A.data);
4685     }
4686 }
4687
4688
4689 NAryMatIterator& NAryMatIterator::operator ++()
4690 {
4691     if( idx >= nplanes-1 )
4692         return *this;
4693     ++idx;
4694
4695     if( iterdepth == 1 )
4696     {
4697         if( ptrs )
4698         {
4699             for( int i = 0; i < narrays; i++ )
4700             {
4701                 if( !ptrs[i] )
4702                     continue;
4703                 ptrs[i] = arrays[i]->data + arrays[i]->step[0]*idx;
4704             }
4705         }
4706         if( planes )
4707         {
4708             for( int i = 0; i < narrays; i++ )
4709             {
4710                 if( !planes[i].data )
4711                     continue;
4712                 planes[i].data = arrays[i]->data + arrays[i]->step[0]*idx;
4713             }
4714         }
4715     }
4716     else
4717     {
4718         for( int i = 0; i < narrays; i++ )
4719         {
4720             const Mat& A = *arrays[i];
4721             if( !A.data )
4722                 continue;
4723             int _idx = (int)idx;
4724             uchar* data = A.data;
4725             for( int j = iterdepth-1; j >= 0 && _idx > 0; j-- )
4726             {
4727                 int szi = A.size[j], t = _idx/szi;
4728                 data += (_idx - t * szi)*A.step[j];
4729                 _idx = t;
4730             }
4731             if( ptrs )
4732                 ptrs[i] = data;
4733             if( planes )
4734                 planes[i].data = data;
4735         }
4736     }
4737
4738     return *this;
4739 }
4740
4741 NAryMatIterator NAryMatIterator::operator ++(int)
4742 {
4743     NAryMatIterator it = *this;
4744     ++*this;
4745     return it;
4746 }
4747
4748 ///////////////////////////////////////////////////////////////////////////
4749 //                              MatConstIterator                         //
4750 ///////////////////////////////////////////////////////////////////////////
4751
4752 Point MatConstIterator::pos() const
4753 {
4754     if( !m )
4755         return Point();
4756     CV_DbgAssert(m->dims <= 2);
4757
4758     ptrdiff_t ofs = ptr - m->data;
4759     int y = (int)(ofs/m->step[0]);
4760     return Point((int)((ofs - y*m->step[0])/elemSize), y);
4761 }
4762
4763 void MatConstIterator::pos(int* _idx) const
4764 {
4765     CV_Assert(m != 0 && _idx);
4766     ptrdiff_t ofs = ptr - m->data;
4767     for( int i = 0; i < m->dims; i++ )
4768     {
4769         size_t s = m->step[i], v = ofs/s;
4770         ofs -= v*s;
4771         _idx[i] = (int)v;
4772     }
4773 }
4774
4775 ptrdiff_t MatConstIterator::lpos() const
4776 {
4777     if(!m)
4778         return 0;
4779     if( m->isContinuous() )
4780         return (ptr - sliceStart)/elemSize;
4781     ptrdiff_t ofs = ptr - m->data;
4782     int i, d = m->dims;
4783     if( d == 2 )
4784     {
4785         ptrdiff_t y = ofs/m->step[0];
4786         return y*m->cols + (ofs - y*m->step[0])/elemSize;
4787     }
4788     ptrdiff_t result = 0;
4789     for( i = 0; i < d; i++ )
4790     {
4791         size_t s = m->step[i], v = ofs/s;
4792         ofs -= v*s;
4793         result = result*m->size[i] + v;
4794     }
4795     return result;
4796 }
4797
4798 void MatConstIterator::seek(ptrdiff_t ofs, bool relative)
4799 {
4800     if( m->isContinuous() )
4801     {
4802         ptr = (relative ? ptr : sliceStart) + ofs*elemSize;
4803         if( ptr < sliceStart )
4804             ptr = sliceStart;
4805         else if( ptr > sliceEnd )
4806             ptr = sliceEnd;
4807         return;
4808     }
4809
4810     int d = m->dims;
4811     if( d == 2 )
4812     {
4813         ptrdiff_t ofs0, y;
4814         if( relative )
4815         {
4816             ofs0 = ptr - m->data;
4817             y = ofs0/m->step[0];
4818             ofs += y*m->cols + (ofs0 - y*m->step[0])/elemSize;
4819         }
4820         y = ofs/m->cols;
4821         int y1 = std::min(std::max((int)y, 0), m->rows-1);
4822         sliceStart = m->data + y1*m->step[0];
4823         sliceEnd = sliceStart + m->cols*elemSize;
4824         ptr = y < 0 ? sliceStart : y >= m->rows ? sliceEnd :
4825             sliceStart + (ofs - y*m->cols)*elemSize;
4826         return;
4827     }
4828
4829     if( relative )
4830         ofs += lpos();
4831
4832     if( ofs < 0 )
4833         ofs = 0;
4834
4835     int szi = m->size[d-1];
4836     ptrdiff_t t = ofs/szi;
4837     int v = (int)(ofs - t*szi);
4838     ofs = t;
4839     ptr = m->data + v*elemSize;
4840     sliceStart = m->data;
4841
4842     for( int i = d-2; i >= 0; i-- )
4843     {
4844         szi = m->size[i];
4845         t = ofs/szi;
4846         v = (int)(ofs - t*szi);
4847         ofs = t;
4848         sliceStart += v*m->step[i];
4849     }
4850
4851     sliceEnd = sliceStart + m->size[d-1]*elemSize;
4852     if( ofs > 0 )
4853         ptr = sliceEnd;
4854     else
4855         ptr = sliceStart + (ptr - m->data);
4856 }
4857
4858 void MatConstIterator::seek(const int* _idx, bool relative)
4859 {
4860     int i, d = m->dims;
4861     ptrdiff_t ofs = 0;
4862     if( !_idx )
4863         ;
4864     else if( d == 2 )
4865         ofs = _idx[0]*m->size[1] + _idx[1];
4866     else
4867     {
4868         for( i = 0; i < d; i++ )
4869             ofs = ofs*m->size[i] + _idx[i];
4870     }
4871     seek(ofs, relative);
4872 }
4873
4874 //////////////////////////////// SparseMat ////////////////////////////////
4875
4876 template<typename T1, typename T2> void
4877 convertData_(const void* _from, void* _to, int cn)
4878 {
4879     const T1* from = (const T1*)_from;
4880     T2* to = (T2*)_to;
4881     if( cn == 1 )
4882         *to = saturate_cast<T2>(*from);
4883     else
4884         for( int i = 0; i < cn; i++ )
4885             to[i] = saturate_cast<T2>(from[i]);
4886 }
4887
4888 template<typename T1, typename T2> void
4889 convertScaleData_(const void* _from, void* _to, int cn, double alpha, double beta)
4890 {
4891     const T1* from = (const T1*)_from;
4892     T2* to = (T2*)_to;
4893     if( cn == 1 )
4894         *to = saturate_cast<T2>(*from*alpha + beta);
4895     else
4896         for( int i = 0; i < cn; i++ )
4897             to[i] = saturate_cast<T2>(from[i]*alpha + beta);
4898 }
4899
4900 typedef void (*ConvertData)(const void* from, void* to, int cn);
4901 typedef void (*ConvertScaleData)(const void* from, void* to, int cn, double alpha, double beta);
4902
4903 static ConvertData getConvertElem(int fromType, int toType)
4904 {
4905     static ConvertData tab[][8] =
4906     {{ convertData_<uchar, uchar>, convertData_<uchar, schar>,
4907       convertData_<uchar, ushort>, convertData_<uchar, short>,
4908       convertData_<uchar, int>, convertData_<uchar, float>,
4909       convertData_<uchar, double>, 0 },
4910
4911     { convertData_<schar, uchar>, convertData_<schar, schar>,
4912       convertData_<schar, ushort>, convertData_<schar, short>,
4913       convertData_<schar, int>, convertData_<schar, float>,
4914       convertData_<schar, double>, 0 },
4915
4916     { convertData_<ushort, uchar>, convertData_<ushort, schar>,
4917       convertData_<ushort, ushort>, convertData_<ushort, short>,
4918       convertData_<ushort, int>, convertData_<ushort, float>,
4919       convertData_<ushort, double>, 0 },
4920
4921     { convertData_<short, uchar>, convertData_<short, schar>,
4922       convertData_<short, ushort>, convertData_<short, short>,
4923       convertData_<short, int>, convertData_<short, float>,
4924       convertData_<short, double>, 0 },
4925
4926     { convertData_<int, uchar>, convertData_<int, schar>,
4927       convertData_<int, ushort>, convertData_<int, short>,
4928       convertData_<int, int>, convertData_<int, float>,
4929       convertData_<int, double>, 0 },
4930
4931     { convertData_<float, uchar>, convertData_<float, schar>,
4932       convertData_<float, ushort>, convertData_<float, short>,
4933       convertData_<float, int>, convertData_<float, float>,
4934       convertData_<float, double>, 0 },
4935
4936     { convertData_<double, uchar>, convertData_<double, schar>,
4937       convertData_<double, ushort>, convertData_<double, short>,
4938       convertData_<double, int>, convertData_<double, float>,
4939       convertData_<double, double>, 0 },
4940
4941     { 0, 0, 0, 0, 0, 0, 0, 0 }};
4942
4943     ConvertData func = tab[CV_MAT_DEPTH(fromType)][CV_MAT_DEPTH(toType)];
4944     CV_Assert( func != 0 );
4945     return func;
4946 }
4947
4948 static ConvertScaleData getConvertScaleElem(int fromType, int toType)
4949 {
4950     static ConvertScaleData tab[][8] =
4951     {{ convertScaleData_<uchar, uchar>, convertScaleData_<uchar, schar>,
4952       convertScaleData_<uchar, ushort>, convertScaleData_<uchar, short>,
4953       convertScaleData_<uchar, int>, convertScaleData_<uchar, float>,
4954       convertScaleData_<uchar, double>, 0 },
4955
4956     { convertScaleData_<schar, uchar>, convertScaleData_<schar, schar>,
4957       convertScaleData_<schar, ushort>, convertScaleData_<schar, short>,
4958       convertScaleData_<schar, int>, convertScaleData_<schar, float>,
4959       convertScaleData_<schar, double>, 0 },
4960
4961     { convertScaleData_<ushort, uchar>, convertScaleData_<ushort, schar>,
4962       convertScaleData_<ushort, ushort>, convertScaleData_<ushort, short>,
4963       convertScaleData_<ushort, int>, convertScaleData_<ushort, float>,
4964       convertScaleData_<ushort, double>, 0 },
4965
4966     { convertScaleData_<short, uchar>, convertScaleData_<short, schar>,
4967       convertScaleData_<short, ushort>, convertScaleData_<short, short>,
4968       convertScaleData_<short, int>, convertScaleData_<short, float>,
4969       convertScaleData_<short, double>, 0 },
4970
4971     { convertScaleData_<int, uchar>, convertScaleData_<int, schar>,
4972       convertScaleData_<int, ushort>, convertScaleData_<int, short>,
4973       convertScaleData_<int, int>, convertScaleData_<int, float>,
4974       convertScaleData_<int, double>, 0 },
4975
4976     { convertScaleData_<float, uchar>, convertScaleData_<float, schar>,
4977       convertScaleData_<float, ushort>, convertScaleData_<float, short>,
4978       convertScaleData_<float, int>, convertScaleData_<float, float>,
4979       convertScaleData_<float, double>, 0 },
4980
4981     { convertScaleData_<double, uchar>, convertScaleData_<double, schar>,
4982       convertScaleData_<double, ushort>, convertScaleData_<double, short>,
4983       convertScaleData_<double, int>, convertScaleData_<double, float>,
4984       convertScaleData_<double, double>, 0 },
4985
4986     { 0, 0, 0, 0, 0, 0, 0, 0 }};
4987
4988     ConvertScaleData func = tab[CV_MAT_DEPTH(fromType)][CV_MAT_DEPTH(toType)];
4989     CV_Assert( func != 0 );
4990     return func;
4991 }
4992
4993 enum { HASH_SIZE0 = 8 };
4994
4995 static inline void copyElem(const uchar* from, uchar* to, size_t elemSize)
4996 {
4997     size_t i;
4998     for( i = 0; i + sizeof(int) <= elemSize; i += sizeof(int) )
4999         *(int*)(to + i) = *(const int*)(from + i);
5000     for( ; i < elemSize; i++ )
5001         to[i] = from[i];
5002 }
5003
5004 static inline bool isZeroElem(const uchar* data, size_t elemSize)
5005 {
5006     size_t i;
5007     for( i = 0; i + sizeof(int) <= elemSize; i += sizeof(int) )
5008         if( *(int*)(data + i) != 0 )
5009             return false;
5010     for( ; i < elemSize; i++ )
5011         if( data[i] != 0 )
5012             return false;
5013     return true;
5014 }
5015
5016 SparseMat::Hdr::Hdr( int _dims, const int* _sizes, int _type )
5017 {
5018     refcount = 1;
5019
5020     dims = _dims;
5021     valueOffset = (int)alignSize(sizeof(SparseMat::Node) +
5022         sizeof(int)*std::max(dims - CV_MAX_DIM, 0), CV_ELEM_SIZE1(_type));
5023     nodeSize = alignSize(valueOffset +
5024         CV_ELEM_SIZE(_type), (int)sizeof(size_t));
5025
5026     int i;
5027     for( i = 0; i < dims; i++ )
5028         size[i] = _sizes[i];
5029     for( ; i < CV_MAX_DIM; i++ )
5030         size[i] = 0;
5031     clear();
5032 }
5033
5034 void SparseMat::Hdr::clear()
5035 {
5036     hashtab.clear();
5037     hashtab.resize(HASH_SIZE0);
5038     pool.clear();
5039     pool.resize(nodeSize);
5040     nodeCount = freeList = 0;
5041 }
5042
5043
5044 SparseMat::SparseMat(const Mat& m)
5045 : flags(MAGIC_VAL), hdr(0)
5046 {
5047     create( m.dims, m.size, m.type() );
5048
5049     int i, idx[CV_MAX_DIM] = {0}, d = m.dims, lastSize = m.size[d - 1];
5050     size_t esz = m.elemSize();
5051     uchar* dptr = m.data;
5052
5053     for(;;)
5054     {
5055         for( i = 0; i < lastSize; i++, dptr += esz )
5056         {
5057             if( isZeroElem(dptr, esz) )
5058                 continue;
5059             idx[d-1] = i;
5060             uchar* to = newNode(idx, hash(idx));
5061             copyElem( dptr, to, esz );
5062         }
5063
5064         for( i = d - 2; i >= 0; i-- )
5065         {
5066             dptr += m.step[i] - m.size[i+1]*m.step[i+1];
5067             if( ++idx[i] < m.size[i] )
5068                 break;
5069             idx[i] = 0;
5070         }
5071         if( i < 0 )
5072             break;
5073     }
5074 }
5075
5076 void SparseMat::create(int d, const int* _sizes, int _type)
5077 {
5078     int i;
5079     CV_Assert( _sizes && 0 < d && d <= CV_MAX_DIM );
5080     for( i = 0; i < d; i++ )
5081         CV_Assert( _sizes[i] > 0 );
5082     _type = CV_MAT_TYPE(_type);
5083     if( hdr && _type == type() && hdr->dims == d && hdr->refcount == 1 )
5084     {
5085         for( i = 0; i < d; i++ )
5086             if( _sizes[i] != hdr->size[i] )
5087                 break;
5088         if( i == d )
5089         {
5090             clear();
5091             return;
5092         }
5093     }
5094     release();
5095     flags = MAGIC_VAL | _type;
5096     hdr = new Hdr(d, _sizes, _type);
5097 }
5098
5099 void SparseMat::copyTo( SparseMat& m ) const
5100 {
5101     if( hdr == m.hdr )
5102         return;
5103     if( !hdr )
5104     {
5105         m.release();
5106         return;
5107     }
5108     m.create( hdr->dims, hdr->size, type() );
5109     SparseMatConstIterator from = begin();
5110     size_t i, N = nzcount(), esz = elemSize();
5111
5112     for( i = 0; i < N; i++, ++from )
5113     {
5114         const Node* n = from.node();
5115         uchar* to = m.newNode(n->idx, n->hashval);
5116         copyElem( from.ptr, to, esz );
5117     }
5118 }
5119
5120 void SparseMat::copyTo( Mat& m ) const
5121 {
5122     CV_Assert( hdr );
5123     m.create( dims(), hdr->size, type() );
5124     m = Scalar(0);
5125
5126     SparseMatConstIterator from = begin();
5127     size_t i, N = nzcount(), esz = elemSize();
5128
5129     for( i = 0; i < N; i++, ++from )
5130     {
5131         const Node* n = from.node();
5132         copyElem( from.ptr, m.ptr(n->idx), esz);
5133     }
5134 }
5135
5136
5137 void SparseMat::convertTo( SparseMat& m, int rtype, double alpha ) const
5138 {
5139     int cn = channels();
5140     if( rtype < 0 )
5141         rtype = type();
5142     rtype = CV_MAKETYPE(rtype, cn);
5143     if( hdr == m.hdr && rtype != type()  )
5144     {
5145         SparseMat temp;
5146         convertTo(temp, rtype, alpha);
5147         m = temp;
5148         return;
5149     }
5150
5151     CV_Assert(hdr != 0);
5152     if( hdr != m.hdr )
5153         m.create( hdr->dims, hdr->size, rtype );
5154
5155     SparseMatConstIterator from = begin();
5156     size_t i, N = nzcount();
5157
5158     if( alpha == 1 )
5159     {
5160         ConvertData cvtfunc = getConvertElem(type(), rtype);
5161         for( i = 0; i < N; i++, ++from )
5162         {
5163             const Node* n = from.node();
5164             uchar* to = hdr == m.hdr ? from.ptr : m.newNode(n->idx, n->hashval);
5165             cvtfunc( from.ptr, to, cn );
5166         }
5167     }
5168     else
5169     {
5170         ConvertScaleData cvtfunc = getConvertScaleElem(type(), rtype);
5171         for( i = 0; i < N; i++, ++from )
5172         {
5173             const Node* n = from.node();
5174             uchar* to = hdr == m.hdr ? from.ptr : m.newNode(n->idx, n->hashval);
5175             cvtfunc( from.ptr, to, cn, alpha, 0 );
5176         }
5177     }
5178 }
5179
5180
5181 void SparseMat::convertTo( Mat& m, int rtype, double alpha, double beta ) const
5182 {
5183     int cn = channels();
5184     if( rtype < 0 )
5185         rtype = type();
5186     rtype = CV_MAKETYPE(rtype, cn);
5187
5188     CV_Assert( hdr );
5189     m.create( dims(), hdr->size, rtype );
5190     m = Scalar(beta);
5191
5192     SparseMatConstIterator from = begin();
5193     size_t i, N = nzcount();
5194
5195     if( alpha == 1 && beta == 0 )
5196     {
5197         ConvertData cvtfunc = getConvertElem(type(), rtype);
5198         for( i = 0; i < N; i++, ++from )
5199         {
5200             const Node* n = from.node();
5201             uchar* to = m.ptr(n->idx);
5202             cvtfunc( from.ptr, to, cn );
5203         }
5204     }
5205     else
5206     {
5207         ConvertScaleData cvtfunc = getConvertScaleElem(type(), rtype);
5208         for( i = 0; i < N; i++, ++from )
5209         {
5210             const Node* n = from.node();
5211             uchar* to = m.ptr(n->idx);
5212             cvtfunc( from.ptr, to, cn, alpha, beta );
5213         }
5214     }
5215 }
5216
5217 void SparseMat::clear()
5218 {
5219     if( hdr )
5220         hdr->clear();
5221 }
5222
5223 uchar* SparseMat::ptr(int i0, bool createMissing, size_t* hashval)
5224 {
5225     CV_Assert( hdr && hdr->dims == 1 );
5226     size_t h = hashval ? *hashval : hash(i0);
5227     size_t hidx = h & (hdr->hashtab.size() - 1), nidx = hdr->hashtab[hidx];
5228     uchar* pool = &hdr->pool[0];
5229     while( nidx != 0 )
5230     {
5231         Node* elem = (Node*)(pool + nidx);
5232         if( elem->hashval == h && elem->idx[0] == i0 )
5233             return &value<uchar>(elem);
5234         nidx = elem->next;
5235     }
5236
5237     if( createMissing )
5238     {
5239         int idx[] = { i0 };
5240         return newNode( idx, h );
5241     }
5242     return 0;
5243 }
5244
5245 uchar* SparseMat::ptr(int i0, int i1, bool createMissing, size_t* hashval)
5246 {
5247     CV_Assert( hdr && hdr->dims == 2 );
5248     size_t h = hashval ? *hashval : hash(i0, i1);
5249     size_t hidx = h & (hdr->hashtab.size() - 1), nidx = hdr->hashtab[hidx];
5250     uchar* pool = &hdr->pool[0];
5251     while( nidx != 0 )
5252     {
5253         Node* elem = (Node*)(pool + nidx);
5254         if( elem->hashval == h && elem->idx[0] == i0 && elem->idx[1] == i1 )
5255             return &value<uchar>(elem);
5256         nidx = elem->next;
5257     }
5258
5259     if( createMissing )
5260     {
5261         int idx[] = { i0, i1 };
5262         return newNode( idx, h );
5263     }
5264     return 0;
5265 }
5266
5267 uchar* SparseMat::ptr(int i0, int i1, int i2, bool createMissing, size_t* hashval)
5268 {
5269     CV_Assert( hdr && hdr->dims == 3 );
5270     size_t h = hashval ? *hashval : hash(i0, i1, i2);
5271     size_t hidx = h & (hdr->hashtab.size() - 1), nidx = hdr->hashtab[hidx];
5272     uchar* pool = &hdr->pool[0];
5273     while( nidx != 0 )
5274     {
5275         Node* elem = (Node*)(pool + nidx);
5276         if( elem->hashval == h && elem->idx[0] == i0 &&
5277             elem->idx[1] == i1 && elem->idx[2] == i2 )
5278             return &value<uchar>(elem);
5279         nidx = elem->next;
5280     }
5281
5282     if( createMissing )
5283     {
5284         int idx[] = { i0, i1, i2 };
5285         return newNode( idx, h );
5286     }
5287     return 0;
5288 }
5289
5290 uchar* SparseMat::ptr(const int* idx, bool createMissing, size_t* hashval)
5291 {
5292     CV_Assert( hdr );
5293     int i, d = hdr->dims;
5294     size_t h = hashval ? *hashval : hash(idx);
5295     size_t hidx = h & (hdr->hashtab.size() - 1), nidx = hdr->hashtab[hidx];
5296     uchar* pool = &hdr->pool[0];
5297     while( nidx != 0 )
5298     {
5299         Node* elem = (Node*)(pool + nidx);
5300         if( elem->hashval == h )
5301         {
5302             for( i = 0; i < d; i++ )
5303                 if( elem->idx[i] != idx[i] )
5304                     break;
5305             if( i == d )
5306                 return &value<uchar>(elem);
5307         }
5308         nidx = elem->next;
5309     }
5310
5311     return createMissing ? newNode(idx, h) : 0;
5312 }
5313
5314 void SparseMat::erase(int i0, int i1, size_t* hashval)
5315 {
5316     CV_Assert( hdr && hdr->dims == 2 );
5317     size_t h = hashval ? *hashval : hash(i0, i1);
5318     size_t hidx = h & (hdr->hashtab.size() - 1), nidx = hdr->hashtab[hidx], previdx=0;
5319     uchar* pool = &hdr->pool[0];
5320     while( nidx != 0 )
5321     {
5322         Node* elem = (Node*)(pool + nidx);
5323         if( elem->hashval == h && elem->idx[0] == i0 && elem->idx[1] == i1 )
5324             break;
5325         previdx = nidx;
5326         nidx = elem->next;
5327     }
5328
5329     if( nidx )
5330         removeNode(hidx, nidx, previdx);
5331 }
5332
5333 void SparseMat::erase(int i0, int i1, int i2, size_t* hashval)
5334 {
5335     CV_Assert( hdr && hdr->dims == 3 );
5336     size_t h = hashval ? *hashval : hash(i0, i1, i2);
5337     size_t hidx = h & (hdr->hashtab.size() - 1), nidx = hdr->hashtab[hidx], previdx=0;
5338     uchar* pool = &hdr->pool[0];
5339     while( nidx != 0 )
5340     {
5341         Node* elem = (Node*)(pool + nidx);
5342         if( elem->hashval == h && elem->idx[0] == i0 &&
5343             elem->idx[1] == i1 && elem->idx[2] == i2 )
5344             break;
5345         previdx = nidx;
5346         nidx = elem->next;
5347     }
5348
5349     if( nidx )
5350         removeNode(hidx, nidx, previdx);
5351 }
5352
5353 void SparseMat::erase(const int* idx, size_t* hashval)
5354 {
5355     CV_Assert( hdr );
5356     int i, d = hdr->dims;
5357     size_t h = hashval ? *hashval : hash(idx);
5358     size_t hidx = h & (hdr->hashtab.size() - 1), nidx = hdr->hashtab[hidx], previdx=0;
5359     uchar* pool = &hdr->pool[0];
5360     while( nidx != 0 )
5361     {
5362         Node* elem = (Node*)(pool + nidx);
5363         if( elem->hashval == h )
5364         {
5365             for( i = 0; i < d; i++ )
5366                 if( elem->idx[i] != idx[i] )
5367                     break;
5368             if( i == d )
5369                 break;
5370         }
5371         previdx = nidx;
5372         nidx = elem->next;
5373     }
5374
5375     if( nidx )
5376         removeNode(hidx, nidx, previdx);
5377 }
5378
5379 void SparseMat::resizeHashTab(size_t newsize)
5380 {
5381     newsize = std::max(newsize, (size_t)8);
5382     if((newsize & (newsize-1)) != 0)
5383         newsize = (size_t)1 << cvCeil(std::log((double)newsize)/CV_LOG2);
5384
5385     size_t i, hsize = hdr->hashtab.size();
5386     std::vector<size_t> _newh(newsize);
5387     size_t* newh = &_newh[0];
5388     for( i = 0; i < newsize; i++ )
5389         newh[i] = 0;
5390     uchar* pool = &hdr->pool[0];
5391     for( i = 0; i < hsize; i++ )
5392     {
5393         size_t nidx = hdr->hashtab[i];
5394         while( nidx )
5395         {
5396             Node* elem = (Node*)(pool + nidx);
5397             size_t next = elem->next;
5398             size_t newhidx = elem->hashval & (newsize - 1);
5399             elem->next = newh[newhidx];
5400             newh[newhidx] = nidx;
5401             nidx = next;
5402         }
5403     }
5404     hdr->hashtab = _newh;
5405 }
5406
5407 uchar* SparseMat::newNode(const int* idx, size_t hashval)
5408 {
5409     const int HASH_MAX_FILL_FACTOR=3;
5410     assert(hdr);
5411     size_t hsize = hdr->hashtab.size();
5412     if( ++hdr->nodeCount > hsize*HASH_MAX_FILL_FACTOR )
5413     {
5414         resizeHashTab(std::max(hsize*2, (size_t)8));
5415         hsize = hdr->hashtab.size();
5416     }
5417
5418     if( !hdr->freeList )
5419     {
5420         size_t i, nsz = hdr->nodeSize, psize = hdr->pool.size(),
5421             newpsize = std::max(psize*2, 8*nsz);
5422         hdr->pool.resize(newpsize);
5423         uchar* pool = &hdr->pool[0];
5424         hdr->freeList = std::max(psize, nsz);
5425         for( i = hdr->freeList; i < newpsize - nsz; i += nsz )
5426             ((Node*)(pool + i))->next = i + nsz;
5427         ((Node*)(pool + i))->next = 0;
5428     }
5429     size_t nidx = hdr->freeList;
5430     Node* elem = (Node*)&hdr->pool[nidx];
5431     hdr->freeList = elem->next;
5432     elem->hashval = hashval;
5433     size_t hidx = hashval & (hsize - 1);
5434     elem->next = hdr->hashtab[hidx];
5435     hdr->hashtab[hidx] = nidx;
5436
5437     int i, d = hdr->dims;
5438     for( i = 0; i < d; i++ )
5439         elem->idx[i] = idx[i];
5440     size_t esz = elemSize();
5441     uchar* p = &value<uchar>(elem);
5442     if( esz == sizeof(float) )
5443         *((float*)p) = 0.f;
5444     else if( esz == sizeof(double) )
5445         *((double*)p) = 0.;
5446     else
5447         memset(p, 0, esz);
5448
5449     return p;
5450 }
5451
5452
5453 void SparseMat::removeNode(size_t hidx, size_t nidx, size_t previdx)
5454 {
5455     Node* n = node(nidx);
5456     if( previdx )
5457     {
5458         Node* prev = node(previdx);
5459         prev->next = n->next;
5460     }
5461     else
5462         hdr->hashtab[hidx] = n->next;
5463     n->next = hdr->freeList;
5464     hdr->freeList = nidx;
5465     --hdr->nodeCount;
5466 }
5467
5468
5469 SparseMatConstIterator::SparseMatConstIterator(const SparseMat* _m)
5470 : m((SparseMat*)_m), hashidx(0), ptr(0)
5471 {
5472     if(!_m || !_m->hdr)
5473         return;
5474     SparseMat::Hdr& hdr = *m->hdr;
5475     const std::vector<size_t>& htab = hdr.hashtab;
5476     size_t i, hsize = htab.size();
5477     for( i = 0; i < hsize; i++ )
5478     {
5479         size_t nidx = htab[i];
5480         if( nidx )
5481         {
5482             hashidx = i;
5483             ptr = &hdr.pool[nidx] + hdr.valueOffset;
5484             return;
5485         }
5486     }
5487 }
5488
5489 SparseMatConstIterator& SparseMatConstIterator::operator ++()
5490 {
5491     if( !ptr || !m || !m->hdr )
5492         return *this;
5493     SparseMat::Hdr& hdr = *m->hdr;
5494     size_t next = ((const SparseMat::Node*)(ptr - hdr.valueOffset))->next;
5495     if( next )
5496     {
5497         ptr = &hdr.pool[next] + hdr.valueOffset;
5498         return *this;
5499     }
5500     size_t i = hashidx + 1, sz = hdr.hashtab.size();
5501     for( ; i < sz; i++ )
5502     {
5503         size_t nidx = hdr.hashtab[i];
5504         if( nidx )
5505         {
5506             hashidx = i;
5507             ptr = &hdr.pool[nidx] + hdr.valueOffset;
5508             return *this;
5509         }
5510     }
5511     hashidx = sz;
5512     ptr = 0;
5513     return *this;
5514 }
5515
5516
5517 double norm( const SparseMat& src, int normType )
5518 {
5519     SparseMatConstIterator it = src.begin();
5520
5521     size_t i, N = src.nzcount();
5522     normType &= NORM_TYPE_MASK;
5523     int type = src.type();
5524     double result = 0;
5525
5526     CV_Assert( normType == NORM_INF || normType == NORM_L1 || normType == NORM_L2 );
5527
5528     if( type == CV_32F )
5529     {
5530         if( normType == NORM_INF )
5531             for( i = 0; i < N; i++, ++it )
5532                 result = std::max(result, std::abs((double)it.value<float>()));
5533         else if( normType == NORM_L1 )
5534             for( i = 0; i < N; i++, ++it )
5535                 result += std::abs(it.value<float>());
5536         else
5537             for( i = 0; i < N; i++, ++it )
5538             {
5539                 double v = it.value<float>();
5540                 result += v*v;
5541             }
5542     }
5543     else if( type == CV_64F )
5544     {
5545         if( normType == NORM_INF )
5546             for( i = 0; i < N; i++, ++it )
5547                 result = std::max(result, std::abs(it.value<double>()));
5548         else if( normType == NORM_L1 )
5549             for( i = 0; i < N; i++, ++it )
5550                 result += std::abs(it.value<double>());
5551         else
5552             for( i = 0; i < N; i++, ++it )
5553             {
5554                 double v = it.value<double>();
5555                 result += v*v;
5556             }
5557     }
5558     else
5559         CV_Error( CV_StsUnsupportedFormat, "Only 32f and 64f are supported" );
5560
5561     if( normType == NORM_L2 )
5562         result = std::sqrt(result);
5563     return result;
5564 }
5565
5566 void minMaxLoc( const SparseMat& src, double* _minval, double* _maxval, int* _minidx, int* _maxidx )
5567 {
5568     SparseMatConstIterator it = src.begin();
5569     size_t i, N = src.nzcount(), d = src.hdr ? src.hdr->dims : 0;
5570     int type = src.type();
5571     const int *minidx = 0, *maxidx = 0;
5572
5573     if( type == CV_32F )
5574     {
5575         float minval = FLT_MAX, maxval = -FLT_MAX;
5576         for( i = 0; i < N; i++, ++it )
5577         {
5578             float v = it.value<float>();
5579             if( v < minval )
5580             {
5581                 minval = v;
5582                 minidx = it.node()->idx;
5583             }
5584             if( v > maxval )
5585             {
5586                 maxval = v;
5587                 maxidx = it.node()->idx;
5588             }
5589         }
5590         if( _minval )
5591             *_minval = minval;
5592         if( _maxval )
5593             *_maxval = maxval;
5594     }
5595     else if( type == CV_64F )
5596     {
5597         double minval = DBL_MAX, maxval = -DBL_MAX;
5598         for( i = 0; i < N; i++, ++it )
5599         {
5600             double v = it.value<float>();
5601             if( v < minval )
5602             {
5603                 minval = v;
5604                 minidx = it.node()->idx;
5605             }
5606             if( v > maxval )
5607             {
5608                 maxval = v;
5609                 maxidx = it.node()->idx;
5610             }
5611         }
5612         if( _minval )
5613             *_minval = minval;
5614         if( _maxval )
5615             *_maxval = maxval;
5616     }
5617     else
5618         CV_Error( CV_StsUnsupportedFormat, "Only 32f and 64f are supported" );
5619
5620     if( _minidx )
5621         for( i = 0; i < d; i++ )
5622             _minidx[i] = minidx[i];
5623     if( _maxidx )
5624         for( i = 0; i < d; i++ )
5625             _maxidx[i] = maxidx[i];
5626 }
5627
5628
5629 void normalize( const SparseMat& src, SparseMat& dst, double a, int norm_type )
5630 {
5631     double scale = 1;
5632     if( norm_type == CV_L2 || norm_type == CV_L1 || norm_type == CV_C )
5633     {
5634         scale = norm( src, norm_type );
5635         scale = scale > DBL_EPSILON ? a/scale : 0.;
5636     }
5637     else
5638         CV_Error( CV_StsBadArg, "Unknown/unsupported norm type" );
5639
5640     src.convertTo( dst, -1, scale );
5641 }
5642
5643 ////////////////////// RotatedRect //////////////////////
5644
5645 RotatedRect::RotatedRect(const Point2f& _point1, const Point2f& _point2, const Point2f& _point3)
5646 {
5647     Point2f _center = 0.5f * (_point1 + _point3);
5648     Vec2f vecs[2];
5649     vecs[0] = Vec2f(_point1 - _point2);
5650     vecs[1] = Vec2f(_point2 - _point3);
5651     // check that given sides are perpendicular
5652     CV_Assert( abs(vecs[0].dot(vecs[1])) / (norm(vecs[0]) * norm(vecs[1])) <= FLT_EPSILON );
5653
5654     // wd_i stores which vector (0,1) or (1,2) will make the width
5655     // One of them will definitely have slope within -1 to 1
5656     int wd_i = 0;
5657     if( abs(vecs[1][1]) < abs(vecs[1][0]) ) wd_i = 1;
5658     int ht_i = (wd_i + 1) % 2;
5659
5660     float _angle = atan(vecs[wd_i][1] / vecs[wd_i][0]) * 180.0f / (float) CV_PI;
5661     float _width = (float) norm(vecs[wd_i]);
5662     float _height = (float) norm(vecs[ht_i]);
5663
5664     center = _center;
5665     size = Size2f(_width, _height);
5666     angle = _angle;
5667 }
5668
5669 void RotatedRect::points(Point2f pt[]) const
5670 {
5671     double _angle = angle*CV_PI/180.;
5672     float b = (float)cos(_angle)*0.5f;
5673     float a = (float)sin(_angle)*0.5f;
5674
5675     pt[0].x = center.x - a*size.height - b*size.width;
5676     pt[0].y = center.y + b*size.height - a*size.width;
5677     pt[1].x = center.x + a*size.height - b*size.width;
5678     pt[1].y = center.y - b*size.height - a*size.width;
5679     pt[2].x = 2*center.x - pt[0].x;
5680     pt[2].y = 2*center.y - pt[0].y;
5681     pt[3].x = 2*center.x - pt[1].x;
5682     pt[3].y = 2*center.y - pt[1].y;
5683 }
5684
5685 Rect RotatedRect::boundingRect() const
5686 {
5687     Point2f pt[4];
5688     points(pt);
5689     Rect r(cvFloor(std::min(std::min(std::min(pt[0].x, pt[1].x), pt[2].x), pt[3].x)),
5690            cvFloor(std::min(std::min(std::min(pt[0].y, pt[1].y), pt[2].y), pt[3].y)),
5691            cvCeil(std::max(std::max(std::max(pt[0].x, pt[1].x), pt[2].x), pt[3].x)),
5692            cvCeil(std::max(std::max(std::max(pt[0].y, pt[1].y), pt[2].y), pt[3].y)));
5693     r.width -= r.x - 1;
5694     r.height -= r.y - 1;
5695     return r;
5696 }
5697
5698 }
5699
5700 // glue
5701
5702 CvMatND::CvMatND(const cv::Mat& m)
5703 {
5704     cvInitMatNDHeader(this, m.dims, m.size, m.type(), m.data );
5705     int i, d = m.dims;
5706     for( i = 0; i < d; i++ )
5707         dim[i].step = (int)m.step[i];
5708     type |= m.flags & cv::Mat::CONTINUOUS_FLAG;
5709 }
5710
5711 _IplImage::_IplImage(const cv::Mat& m)
5712 {
5713     CV_Assert( m.dims <= 2 );
5714     cvInitImageHeader(this, m.size(), cvIplDepth(m.flags), m.channels());
5715     cvSetData(this, m.data, (int)m.step[0]);
5716 }
5717
5718 CvSparseMat* cvCreateSparseMat(const cv::SparseMat& sm)
5719 {
5720     if( !sm.hdr )
5721         return 0;
5722
5723     CvSparseMat* m = cvCreateSparseMat(sm.hdr->dims, sm.hdr->size, sm.type());
5724
5725     cv::SparseMatConstIterator from = sm.begin();
5726     size_t i, N = sm.nzcount(), esz = sm.elemSize();
5727
5728     for( i = 0; i < N; i++, ++from )
5729     {
5730         const cv::SparseMat::Node* n = from.node();
5731         uchar* to = cvPtrND(m, n->idx, 0, -2, 0);
5732         cv::copyElem(from.ptr, to, esz);
5733     }
5734     return m;
5735 }
5736
5737 void CvSparseMat::copyToSparseMat(cv::SparseMat& m) const
5738 {
5739     m.create( dims, &size[0], type );
5740
5741     CvSparseMatIterator it;
5742     CvSparseNode* n = cvInitSparseMatIterator(this, &it);
5743     size_t esz = m.elemSize();
5744
5745     for( ; n != 0; n = cvGetNextSparseNode(&it) )
5746     {
5747         const int* idx = CV_NODE_IDX(this, n);
5748         uchar* to = m.newNode(idx, m.hash(idx));
5749         cv::copyElem((const uchar*)CV_NODE_VAL(this, n), to, esz);
5750     }
5751 }
5752
5753
5754 /* End of file. */