modules/imgproc/src/histogram.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 //            Intel License Agreement
  11 //        For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
  14 // Third party copyrights are property of their respective owners.
  15 //
  16 // Redistribution and use in source and binary forms, with or without modification,
  17 // are permitted provided that the following conditions are met:
  18 //
  19 //   * Redistribution's of source code must retain the above copyright notice,
  20 //     this list of conditions and the following disclaimer.
  21 //
  22 //   * Redistribution's in binary form must reproduce the above copyright notice,
  23 //     this list of conditions and the following disclaimer in the documentation
  24 //     and/or other materials provided with the distribution.
  25 //
  26 //   * The name of Intel Corporation may not be used to endorse or promote products
  27 //     derived from this software without specific prior written permission.
  28 //
  29 // This software is provided by the copyright holders and contributors "as is" and
  30 // any express or implied warranties, including, but not limited to, the implied
  31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  32 // In no event shall the Intel Corporation or contributors be liable for any direct,
  33 // indirect, incidental, special, exemplary, or consequential damages
  34 // (including, but not limited to, procurement of substitute goods or services;
  35 // loss of use, data, or profits; or business interruption) however caused
  36 // and on any theory of liability, whether in contract, strict liability,
  37 // or tort (including negligence or otherwise) arising in any way out of
  38 // the use of this software, even if advised of the possibility of such damage.
  39 //
  40 //M*/
  41
  42 #include "precomp.hpp"
  43 #include "opencl_kernels_imgproc.hpp"
  44
  45 namespace cv
  46 {
  47
  48 ////////////////// Helper functions //////////////////////
  49
  50 static const size_t OUT_OF_RANGE = (size_t)1 << (sizeof(size_t)*8 - 2);
  51
  52 static void
  53 calcHistLookupTables_8u( const Mat& hist, const SparseMat& shist,
  54                          int dims, const float** ranges, const double* uniranges,
  55                          bool uniform, bool issparse, std::vector<size_t>& _tab )
  56 {
  57     const int low = 0, high = 256;
  58     int i, j;
  59     _tab.resize((high-low)*dims);
  60     size_t* tab = &_tab[0];
  61
  62     if( uniform )
  63     {
  64         for( i = 0; i < dims; i++ )
  65         {
  66             double a = uniranges[i*2];
  67             double b = uniranges[i*2+1];
  68             int sz = !issparse ? hist.size[i] : shist.size(i);
  69             size_t step = !issparse ? hist.step[i] : 1;
  70
  71             for( j = low; j < high; j++ )
  72             {
  73                 int idx = cvFloor(j*a + b);
  74                 size_t written_idx;
  75                 if( (unsigned)idx < (unsigned)sz )
  76                     written_idx = idx*step;
  77                 else
  78                     written_idx = OUT_OF_RANGE;
  79
  80                 tab[i*(high - low) + j - low] = written_idx;
  81             }
  82         }
  83     }
  84     else
  85     {
  86         for( i = 0; i < dims; i++ )
  87         {
  88             int limit = std::min(cvCeil(ranges[i][0]), high);
  89             int idx = -1, sz = !issparse ? hist.size[i] : shist.size(i);
  90             size_t written_idx = OUT_OF_RANGE;
  91             size_t step = !issparse ? hist.step[i] : 1;
  92
  93             for(j = low;;)
  94             {
  95                 for( ; j < limit; j++ )
  96                     tab[i*(high - low) + j - low] = written_idx;
  97
  98                 if( (unsigned)(++idx) < (unsigned)sz )
  99                 {
 100                     limit = std::min(cvCeil(ranges[i][idx+1]), high);
 101                     written_idx = idx*step;
 102                 }
 103                 else
 104                 {
 105                     for( ; j < high; j++ )
 106                         tab[i*(high - low) + j - low] = OUT_OF_RANGE;
 107                     break;
 108                 }
 109             }
 110         }
 111     }
 112 }
 113
 114
 115 static void histPrepareImages( const Mat* images, int nimages, const int* channels,
 116                                const Mat& mask, int dims, const int* histSize,
 117                                const float** ranges, bool uniform,
 118                                std::vector<uchar*>& ptrs, std::vector<int>& deltas,
 119                                Size& imsize, std::vector<double>& uniranges )
 120 {
 121     int i, j, c;
 122     CV_Assert( channels != 0 || nimages == dims );
 123
 124     imsize = images[0].size();
 125     int depth = images[0].depth(), esz1 = (int)images[0].elemSize1();
 126     bool isContinuous = true;
 127
 128     ptrs.resize(dims + 1);
 129     deltas.resize((dims + 1)*2);
 130
 131     for( i = 0; i < dims; i++ )
 132     {
 133         if(!channels)
 134         {
 135             j = i;
 136             c = 0;
 137             CV_Assert( images[j].channels() == 1 );
 138         }
 139         else
 140         {
 141             c = channels[i];
 142             CV_Assert( c >= 0 );
 143             for( j = 0; j < nimages; c -= images[j].channels(), j++ )
 144                 if( c < images[j].channels() )
 145                     break;
 146             CV_Assert( j < nimages );
 147         }
 148
 149         CV_Assert( images[j].size() == imsize && images[j].depth() == depth );
 150         if( !images[j].isContinuous() )
 151             isContinuous = false;
 152         ptrs[i] = images[j].data + c*esz1;
 153         deltas[i*2] = images[j].channels();
 154         deltas[i*2+1] = (int)(images[j].step/esz1 - imsize.width*deltas[i*2]);
 155     }
 156
 157     if( !mask.empty() )
 158     {
 159         CV_Assert( mask.size() == imsize && mask.channels() == 1 );
 160         isContinuous = isContinuous && mask.isContinuous();
 161         ptrs[dims] = mask.data;
 162         deltas[dims*2] = 1;
 163         deltas[dims*2 + 1] = (int)(mask.step/mask.elemSize1());
 164     }
 165
 166 #ifndef HAVE_TBB
 167     if( isContinuous )
 168     {
 169         imsize.width *= imsize.height;
 170         imsize.height = 1;
 171     }
 172 #endif
 173
 174     if( !ranges )
 175     {
 176         CV_Assert( depth == CV_8U );
 177
 178         uniranges.resize( dims*2 );
 179         for( i = 0; i < dims; i++ )
 180         {
 181             uniranges[i*2] = histSize[i]/256.;
 182             uniranges[i*2+1] = 0;
 183         }
 184     }
 185     else if( uniform )
 186     {
 187         uniranges.resize( dims*2 );
 188         for( i = 0; i < dims; i++ )
 189         {
 190             CV_Assert( ranges[i] && ranges[i][0] < ranges[i][1] );
 191             double low = ranges[i][0], high = ranges[i][1];
 192             double t = histSize[i]/(high - low);
 193             uniranges[i*2] = t;
 194             uniranges[i*2+1] = -t*low;
 195         }
 196     }
 197     else
 198     {
 199         for( i = 0; i < dims; i++ )
 200         {
 201             size_t n = histSize[i];
 202             for(size_t k = 0; k < n; k++ )
 203                 CV_Assert( ranges[i][k] < ranges[i][k+1] );
 204         }
 205     }
 206 }
 207
 208
 209 ////////////////////////////////// C A L C U L A T E    H I S T O G R A M ////////////////////////////////////
 210 #ifdef HAVE_TBB
 211 enum {one = 1, two, three}; // array elements number
 212
 213 template<typename T>
 214 class calcHist1D_Invoker
 215 {
 216 public:
 217     calcHist1D_Invoker( const std::vector<uchar*>& _ptrs, const std::vector<int>& _deltas,
 218                         Mat& hist, const double* _uniranges, int sz, int dims,
 219                         Size& imageSize )
 220         : mask_(_ptrs[dims]),
 221           mstep_(_deltas[dims*2 + 1]),
 222           imageWidth_(imageSize.width),
 223           histogramSize_(hist.size()), histogramType_(hist.type()),
 224           globalHistogram_((tbb::atomic<int>*)hist.data)
 225     {
 226         p_[0] = ((T**)&_ptrs[0])[0];
 227         step_[0] = (&_deltas[0])[1];
 228         d_[0] = (&_deltas[0])[0];
 229         a_[0] = (&_uniranges[0])[0];
 230         b_[0] = (&_uniranges[0])[1];
 231         size_[0] = sz;
 232     }
 233
 234     void operator()( const BlockedRange& range ) const
 235     {
 236         T* p0 = p_[0] + range.begin() * (step_[0] + imageWidth_*d_[0]);
 237         uchar* mask = mask_ + range.begin()*mstep_;
 238
 239         for( int row = range.begin(); row < range.end(); row++, p0 += step_[0] )
 240         {
 241             if( !mask_ )
 242             {
 243                 for( int x = 0; x < imageWidth_; x++, p0 += d_[0] )
 244                 {
 245                     int idx = cvFloor(*p0*a_[0] + b_[0]);
 246                     if( (unsigned)idx < (unsigned)size_[0] )
 247                     {
 248                         globalHistogram_[idx].fetch_and_add(1);
 249                     }
 250                 }
 251             }
 252             else
 253             {
 254                 for( int x = 0; x < imageWidth_; x++, p0 += d_[0] )
 255                 {
 256                     if( mask[x] )
 257                     {
 258                         int idx = cvFloor(*p0*a_[0] + b_[0]);
 259                         if( (unsigned)idx < (unsigned)size_[0] )
 260                         {
 261                             globalHistogram_[idx].fetch_and_add(1);
 262                         }
 263                     }
 264                 }
 265                 mask += mstep_;
 266             }
 267         }
 268     }
 269
 270 private:
 271     calcHist1D_Invoker operator=(const calcHist1D_Invoker&);
 272
 273     T* p_[one];
 274     uchar* mask_;
 275     int step_[one];
 276     int d_[one];
 277     int mstep_;
 278     double a_[one];
 279     double b_[one];
 280     int size_[one];
 281     int imageWidth_;
 282     Size histogramSize_;
 283     int histogramType_;
 284     tbb::atomic<int>* globalHistogram_;
 285 };
 286
 287 template<typename T>
 288 class calcHist2D_Invoker
 289 {
 290 public:
 291     calcHist2D_Invoker( const std::vector<uchar*>& _ptrs, const std::vector<int>& _deltas,
 292                         Mat& hist, const double* _uniranges, const int* size,
 293                         int dims, Size& imageSize, size_t* hstep )
 294         : mask_(_ptrs[dims]),
 295           mstep_(_deltas[dims*2 + 1]),
 296           imageWidth_(imageSize.width),
 297           histogramSize_(hist.size()), histogramType_(hist.type()),
 298           globalHistogram_(hist.data)
 299     {
 300         p_[0] = ((T**)&_ptrs[0])[0]; p_[1] = ((T**)&_ptrs[0])[1];
 301         step_[0] = (&_deltas[0])[1]; step_[1] = (&_deltas[0])[3];
 302         d_[0] = (&_deltas[0])[0];    d_[1] = (&_deltas[0])[2];
 303         a_[0] = (&_uniranges[0])[0]; a_[1] = (&_uniranges[0])[2];
 304         b_[0] = (&_uniranges[0])[1]; b_[1] = (&_uniranges[0])[3];
 305         size_[0] = size[0];          size_[1] = size[1];
 306         hstep_[0] = hstep[0];
 307     }
 308
 309     void operator()(const BlockedRange& range) const
 310     {
 311         T* p0 = p_[0] + range.begin()*(step_[0] + imageWidth_*d_[0]);
 312         T* p1 = p_[1] + range.begin()*(step_[1] + imageWidth_*d_[1]);
 313         uchar* mask = mask_ + range.begin()*mstep_;
 314
 315         for( int row = range.begin(); row < range.end(); row++, p0 += step_[0], p1 += step_[1] )
 316         {
 317             if( !mask_ )
 318             {
 319                 for( int x = 0; x < imageWidth_; x++, p0 += d_[0], p1 += d_[1] )
 320                 {
 321                     int idx0 = cvFloor(*p0*a_[0] + b_[0]);
 322                     int idx1 = cvFloor(*p1*a_[1] + b_[1]);
 323                     if( (unsigned)idx0 < (unsigned)size_[0] && (unsigned)idx1 < (unsigned)size_[1] )
 324                         ( (tbb::atomic<int>*)(globalHistogram_ + hstep_[0]*idx0) )[idx1].fetch_and_add(1);
 325                 }
 326             }
 327             else
 328             {
 329                 for( int x = 0; x < imageWidth_; x++, p0 += d_[0], p1 += d_[1] )
 330                 {
 331                     if( mask[x] )
 332                     {
 333                         int idx0 = cvFloor(*p0*a_[0] + b_[0]);
 334                         int idx1 = cvFloor(*p1*a_[1] + b_[1]);
 335                         if( (unsigned)idx0 < (unsigned)size_[0] && (unsigned)idx1 < (unsigned)size_[1] )
 336                             ((tbb::atomic<int>*)(globalHistogram_ + hstep_[0]*idx0))[idx1].fetch_and_add(1);
 337                     }
 338                 }
 339                 mask += mstep_;
 340             }
 341         }
 342     }
 343
 344 private:
 345     calcHist2D_Invoker operator=(const calcHist2D_Invoker&);
 346
 347     T* p_[two];
 348     uchar* mask_;
 349     int step_[two];
 350     int d_[two];
 351     int mstep_;
 352     double a_[two];
 353     double b_[two];
 354     int size_[two];
 355     const int imageWidth_;
 356     size_t hstep_[one];
 357     Size histogramSize_;
 358     int histogramType_;
 359     uchar* globalHistogram_;
 360 };
 361
 362
 363 template<typename T>
 364 class calcHist3D_Invoker
 365 {
 366 public:
 367     calcHist3D_Invoker( const std::vector<uchar*>& _ptrs, const std::vector<int>& _deltas,
 368                         Size imsize, Mat& hist, const double* uniranges, int _dims,
 369                         size_t* hstep, int* size )
 370         : mask_(_ptrs[_dims]),
 371           mstep_(_deltas[_dims*2 + 1]),
 372           imageWidth_(imsize.width),
 373           globalHistogram_(hist.data)
 374     {
 375         p_[0] = ((T**)&_ptrs[0])[0]; p_[1] = ((T**)&_ptrs[0])[1]; p_[2] = ((T**)&_ptrs[0])[2];
 376         step_[0] = (&_deltas[0])[1]; step_[1] = (&_deltas[0])[3]; step_[2] = (&_deltas[0])[5];
 377         d_[0] = (&_deltas[0])[0];    d_[1] = (&_deltas[0])[2];    d_[2] = (&_deltas[0])[4];
 378         a_[0] = uniranges[0];        a_[1] = uniranges[2];        a_[2] = uniranges[4];
 379         b_[0] = uniranges[1];        b_[1] = uniranges[3];        b_[2] = uniranges[5];
 380         size_[0] = size[0];          size_[1] = size[1];          size_[2] = size[2];
 381         hstep_[0] = hstep[0];        hstep_[1] = hstep[1];
 382     }
 383
 384     void operator()( const BlockedRange& range ) const
 385     {
 386         T* p0 = p_[0] + range.begin()*(imageWidth_*d_[0] + step_[0]);
 387         T* p1 = p_[1] + range.begin()*(imageWidth_*d_[1] + step_[1]);
 388         T* p2 = p_[2] + range.begin()*(imageWidth_*d_[2] + step_[2]);
 389         uchar* mask = mask_ + range.begin()*mstep_;
 390
 391         for( int i = range.begin(); i < range.end(); i++, p0 += step_[0], p1 += step_[1], p2 += step_[2] )
 392         {
 393             if( !mask_ )
 394             {
 395                 for( int x = 0; x < imageWidth_; x++, p0 += d_[0], p1 += d_[1], p2 += d_[2] )
 396                 {
 397                     int idx0 = cvFloor(*p0*a_[0] + b_[0]);
 398                     int idx1 = cvFloor(*p1*a_[1] + b_[1]);
 399                     int idx2 = cvFloor(*p2*a_[2] + b_[2]);
 400                     if( (unsigned)idx0 < (unsigned)size_[0] &&
 401                             (unsigned)idx1 < (unsigned)size_[1] &&
 402                             (unsigned)idx2 < (unsigned)size_[2] )
 403                     {
 404                         ( (tbb::atomic<int>*)(globalHistogram_ + hstep_[0]*idx0 + hstep_[1]*idx1) )[idx2].fetch_and_add(1);
 405                     }
 406                 }
 407             }
 408             else
 409             {
 410                 for( int x = 0; x < imageWidth_; x++, p0 += d_[0], p1 += d_[1], p2 += d_[2] )
 411                 {
 412                     if( mask[x] )
 413                     {
 414                         int idx0 = cvFloor(*p0*a_[0] + b_[0]);
 415                         int idx1 = cvFloor(*p1*a_[1] + b_[1]);
 416                         int idx2 = cvFloor(*p2*a_[2] + b_[2]);
 417                         if( (unsigned)idx0 < (unsigned)size_[0] &&
 418                                 (unsigned)idx1 < (unsigned)size_[1] &&
 419                                 (unsigned)idx2 < (unsigned)size_[2] )
 420                         {
 421                             ( (tbb::atomic<int>*)(globalHistogram_ + hstep_[0]*idx0 + hstep_[1]*idx1) )[idx2].fetch_and_add(1);
 422                         }
 423                     }
 424                 }
 425                 mask += mstep_;
 426             }
 427         }
 428     }
 429
 430     static bool isFit( const Mat& histogram, const Size imageSize )
 431     {
 432         return ( imageSize.width * imageSize.height >= 320*240
 433                  && histogram.total() >= 8*8*8 );
 434     }
 435
 436 private:
 437     calcHist3D_Invoker operator=(const calcHist3D_Invoker&);
 438
 439     T* p_[three];
 440     uchar* mask_;
 441     int step_[three];
 442     int d_[three];
 443     const int mstep_;
 444     double a_[three];
 445     double b_[three];
 446     int size_[three];
 447     int imageWidth_;
 448     size_t hstep_[two];
 449     uchar* globalHistogram_;
 450 };
 451
 452 class CalcHist1D_8uInvoker
 453 {
 454 public:
 455     CalcHist1D_8uInvoker( const std::vector<uchar*>& ptrs, const std::vector<int>& deltas,
 456                           Size imsize, Mat& hist, int dims, const std::vector<size_t>& tab,
 457                           tbb::mutex* lock )
 458         : mask_(ptrs[dims]),
 459           mstep_(deltas[dims*2 + 1]),
 460           imageWidth_(imsize.width),
 461           imageSize_(imsize),
 462           histSize_(hist.size()), histType_(hist.type()),
 463           tab_((size_t*)&tab[0]),
 464           histogramWriteLock_(lock),
 465           globalHistogram_(hist.data)
 466     {
 467         p_[0] = (&ptrs[0])[0];
 468         step_[0] = (&deltas[0])[1];
 469         d_[0] = (&deltas[0])[0];
 470     }
 471
 472     void operator()( const BlockedRange& range ) const
 473     {
 474         int localHistogram[256] = { 0, };
 475         uchar* mask = mask_;
 476         uchar* p0 = p_[0];
 477         int x;
 478         tbb::mutex::scoped_lock lock;
 479
 480         if( !mask_ )
 481         {
 482             int n = (imageWidth_ - 4) / 4 + 1;
 483             int tail = imageWidth_ - n*4;
 484
 485             int xN = 4*n;
 486             p0 += (xN*d_[0] + tail*d_[0] + step_[0]) * range.begin();
 487         }
 488         else
 489         {
 490             p0 += (imageWidth_*d_[0] + step_[0]) * range.begin();
 491             mask += mstep_*range.begin();
 492         }
 493
 494         for( int i = range.begin(); i < range.end(); i++, p0 += step_[0] )
 495         {
 496             if( !mask_ )
 497             {
 498                 if( d_[0] == 1 )
 499                 {
 500                     for( x = 0; x <= imageWidth_ - 4; x += 4 )
 501                     {
 502                         int t0 = p0[x], t1 = p0[x+1];
 503                         localHistogram[t0]++; localHistogram[t1]++;
 504                         t0 = p0[x+2]; t1 = p0[x+3];
 505                         localHistogram[t0]++; localHistogram[t1]++;
 506                     }
 507                     p0 += x;
 508                 }
 509                 else
 510                 {
 511                     for( x = 0; x <= imageWidth_ - 4; x += 4 )
 512                     {
 513                         int t0 = p0[0], t1 = p0[d_[0]];
 514                         localHistogram[t0]++; localHistogram[t1]++;
 515                         p0 += d_[0]*2;
 516                         t0 = p0[0]; t1 = p0[d_[0]];
 517                         localHistogram[t0]++; localHistogram[t1]++;
 518                         p0 += d_[0]*2;
 519                     }
 520                 }
 521
 522                 for( ; x < imageWidth_; x++, p0 += d_[0] )
 523                 {
 524                     localHistogram[*p0]++;
 525                 }
 526             }
 527             else
 528             {
 529                 for( x = 0; x < imageWidth_; x++, p0 += d_[0] )
 530                 {
 531                     if( mask[x] )
 532                     {
 533                         localHistogram[*p0]++;
 534                     }
 535                 }
 536                 mask += mstep_;
 537             }
 538         }
 539
 540         lock.acquire(*histogramWriteLock_);
 541         for(int i = 0; i < 256; i++ )
 542         {
 543             size_t hidx = tab_[i];
 544             if( hidx < OUT_OF_RANGE )
 545             {
 546                 *(int*)((globalHistogram_ + hidx)) += localHistogram[i];
 547             }
 548         }
 549         lock.release();
 550     }
 551
 552     static bool isFit( const Mat& histogram, const Size imageSize )
 553     {
 554         return ( histogram.total() >= 8
 555                 && imageSize.width * imageSize.height >= 160*120 );
 556     }
 557
 558 private:
 559     uchar* p_[one];
 560     uchar* mask_;
 561     int mstep_;
 562     int step_[one];
 563     int d_[one];
 564     int imageWidth_;
 565     Size imageSize_;
 566     Size histSize_;
 567     int histType_;
 568     size_t* tab_;
 569     tbb::mutex* histogramWriteLock_;
 570     uchar* globalHistogram_;
 571 };
 572
 573 class CalcHist2D_8uInvoker
 574 {
 575 public:
 576     CalcHist2D_8uInvoker( const std::vector<uchar*>& _ptrs, const std::vector<int>& _deltas,
 577                           Size imsize, Mat& hist, int dims, const std::vector<size_t>& _tab,
 578                           tbb::mutex* lock )
 579         : mask_(_ptrs[dims]),
 580           mstep_(_deltas[dims*2 + 1]),
 581           imageWidth_(imsize.width),
 582           histSize_(hist.size()), histType_(hist.type()),
 583           tab_((size_t*)&_tab[0]),
 584           histogramWriteLock_(lock),
 585           globalHistogram_(hist.data)
 586     {
 587         p_[0] = (uchar*)(&_ptrs[0])[0]; p_[1] = (uchar*)(&_ptrs[0])[1];
 588         step_[0] = (&_deltas[0])[1];    step_[1] = (&_deltas[0])[3];
 589         d_[0] = (&_deltas[0])[0];       d_[1] = (&_deltas[0])[2];
 590     }
 591
 592     void operator()( const BlockedRange& range ) const
 593     {
 594         uchar* p0 = p_[0] + range.begin()*(step_[0] + imageWidth_*d_[0]);
 595         uchar* p1 = p_[1] + range.begin()*(step_[1] + imageWidth_*d_[1]);
 596         uchar* mask = mask_ + range.begin()*mstep_;
 597
 598         Mat localHist = Mat::zeros(histSize_, histType_);
 599         uchar* localHistData = localHist.data;
 600         tbb::mutex::scoped_lock lock;
 601
 602         for(int i = range.begin(); i < range.end(); i++, p0 += step_[0], p1 += step_[1])
 603         {
 604             if( !mask_ )
 605             {
 606                 for( int x = 0; x < imageWidth_; x++, p0 += d_[0], p1 += d_[1] )
 607                 {
 608                     size_t idx = tab_[*p0] + tab_[*p1 + 256];
 609                     if( idx < OUT_OF_RANGE )
 610                     {
 611                         ++*(int*)(localHistData + idx);
 612                     }
 613                 }
 614             }
 615             else
 616             {
 617                 for( int x = 0; x < imageWidth_; x++, p0 += d_[0], p1 += d_[1] )
 618                 {
 619                     size_t idx;
 620                     if( mask[x] && (idx = tab_[*p0] + tab_[*p1 + 256]) < OUT_OF_RANGE )
 621                     {
 622                         ++*(int*)(localHistData + idx);
 623                     }
 624                 }
 625                 mask += mstep_;
 626             }
 627         }
 628
 629         lock.acquire(*histogramWriteLock_);
 630         for(int i = 0; i < histSize_.width*histSize_.height; i++)
 631         {
 632             ((int*)globalHistogram_)[i] += ((int*)localHistData)[i];
 633         }
 634         lock.release();
 635     }
 636
 637     static bool isFit( const Mat& histogram, const Size imageSize )
 638     {
 639         return ( (histogram.total() > 4*4 &&  histogram.total() <= 116*116
 640                   && imageSize.width * imageSize.height >= 320*240)
 641                  || (histogram.total() > 116*116 && imageSize.width * imageSize.height >= 1280*720) );
 642     }
 643
 644 private:
 645     uchar* p_[two];
 646     uchar* mask_;
 647     int step_[two];
 648     int d_[two];
 649     int mstep_;
 650     int imageWidth_;
 651     Size histSize_;
 652     int histType_;
 653     size_t* tab_;
 654     tbb::mutex* histogramWriteLock_;
 655     uchar* globalHistogram_;
 656 };
 657
 658 class CalcHist3D_8uInvoker
 659 {
 660 public:
 661     CalcHist3D_8uInvoker( const std::vector<uchar*>& _ptrs, const std::vector<int>& _deltas,
 662                           Size imsize, Mat& hist, int dims, const std::vector<size_t>& tab )
 663         : mask_(_ptrs[dims]),
 664           mstep_(_deltas[dims*2 + 1]),
 665           histogramSize_(hist.size.p), histogramType_(hist.type()),
 666           imageWidth_(imsize.width),
 667           tab_((size_t*)&tab[0]),
 668           globalHistogram_(hist.data)
 669     {
 670         p_[0] = (uchar*)(&_ptrs[0])[0]; p_[1] = (uchar*)(&_ptrs[0])[1]; p_[2] = (uchar*)(&_ptrs[0])[2];
 671         step_[0] = (&_deltas[0])[1];    step_[1] = (&_deltas[0])[3];    step_[2] = (&_deltas[0])[5];
 672         d_[0] = (&_deltas[0])[0];       d_[1] = (&_deltas[0])[2];       d_[2] = (&_deltas[0])[4];
 673     }
 674
 675     void operator()( const BlockedRange& range ) const
 676     {
 677         uchar* p0 = p_[0] + range.begin()*(step_[0] + imageWidth_*d_[0]);
 678         uchar* p1 = p_[1] + range.begin()*(step_[1] + imageWidth_*d_[1]);
 679         uchar* p2 = p_[2] + range.begin()*(step_[2] + imageWidth_*d_[2]);
 680         uchar* mask = mask_ + range.begin()*mstep_;
 681
 682         for(int i = range.begin(); i < range.end(); i++, p0 += step_[0], p1 += step_[1], p2 += step_[2] )
 683         {
 684             if( !mask_ )
 685             {
 686                 for( int x = 0; x < imageWidth_; x++, p0 += d_[0], p1 += d_[1], p2 += d_[2] )
 687                 {
 688                     size_t idx = tab_[*p0] + tab_[*p1 + 256] + tab_[*p2 + 512];
 689                     if( idx < OUT_OF_RANGE )
 690                     {
 691                         ( *(tbb::atomic<int>*)(globalHistogram_ + idx) ).fetch_and_add(1);
 692                     }
 693                 }
 694             }
 695             else
 696             {
 697                 for( int x = 0; x < imageWidth_; x++, p0 += d_[0], p1 += d_[1], p2 += d_[2] )
 698                 {
 699                     size_t idx;
 700                     if( mask[x] && (idx = tab_[*p0] + tab_[*p1 + 256] + tab_[*p2 + 512]) < OUT_OF_RANGE )
 701                     {
 702                         (*(tbb::atomic<int>*)(globalHistogram_ + idx)).fetch_and_add(1);
 703                     }
 704                 }
 705                 mask += mstep_;
 706             }
 707         }
 708     }
 709
 710     static bool isFit( const Mat& histogram, const Size imageSize )
 711     {
 712         return ( histogram.total() >= 128*128*128
 713                  && imageSize.width * imageSize.width >= 320*240 );
 714     }
 715
 716 private:
 717     uchar* p_[three];
 718     uchar* mask_;
 719     int mstep_;
 720     int step_[three];
 721     int d_[three];
 722     int* histogramSize_;
 723     int histogramType_;
 724     int imageWidth_;
 725     size_t* tab_;
 726     uchar* globalHistogram_;
 727 };
 728
 729 static void
 730 callCalcHist2D_8u( std::vector<uchar*>& _ptrs, const std::vector<int>& _deltas,
 731                    Size imsize, Mat& hist, int dims,  std::vector<size_t>& _tab )
 732 {
 733     int grainSize = imsize.height / tbb::task_scheduler_init::default_num_threads();
 734     tbb::mutex histogramWriteLock;
 735
 736     CalcHist2D_8uInvoker body(_ptrs, _deltas, imsize, hist, dims, _tab, &histogramWriteLock);
 737     parallel_for(BlockedRange(0, imsize.height, grainSize), body);
 738 }
 739
 740 static void
 741 callCalcHist3D_8u( std::vector<uchar*>& _ptrs, const std::vector<int>& _deltas,
 742                    Size imsize, Mat& hist, int dims,  std::vector<size_t>& _tab )
 743 {
 744     CalcHist3D_8uInvoker body(_ptrs, _deltas, imsize, hist, dims, _tab);
 745     parallel_for(BlockedRange(0, imsize.height), body);
 746 }
 747 #endif
 748
 749 template<typename T> static void
 750 calcHist_( std::vector<uchar*>& _ptrs, const std::vector<int>& _deltas,
 751            Size imsize, Mat& hist, int dims, const float** _ranges,
 752            const double* _uniranges, bool uniform )
 753 {
 754     T** ptrs = (T**)&_ptrs[0];
 755     const int* deltas = &_deltas[0];
 756     uchar* H = hist.ptr();
 757     int i, x;
 758     const uchar* mask = _ptrs[dims];
 759     int mstep = _deltas[dims*2 + 1];
 760     int size[CV_MAX_DIM];
 761     size_t hstep[CV_MAX_DIM];
 762
 763     for( i = 0; i < dims; i++ )
 764     {
 765         size[i] = hist.size[i];
 766         hstep[i] = hist.step[i];
 767     }
 768
 769     if( uniform )
 770     {
 771         const double* uniranges = &_uniranges[0];
 772
 773         if( dims == 1 )
 774         {
 775 #ifdef HAVE_TBB
 776             calcHist1D_Invoker<T> body(_ptrs, _deltas, hist, _uniranges, size[0], dims, imsize);
 777             parallel_for(BlockedRange(0, imsize.height), body);
 778 #else
 779             double a = uniranges[0], b = uniranges[1];
 780             int sz = size[0], d0 = deltas[0], step0 = deltas[1];
 781             const T* p0 = (const T*)ptrs[0];
 782
 783             for( ; imsize.height--; p0 += step0, mask += mstep )
 784             {
 785                 if( !mask )
 786                     for( x = 0; x < imsize.width; x++, p0 += d0 )
 787                     {
 788                         int idx = cvFloor(*p0*a + b);
 789                         if( (unsigned)idx < (unsigned)sz )
 790                             ((int*)H)[idx]++;
 791                     }
 792                 else
 793                     for( x = 0; x < imsize.width; x++, p0 += d0 )
 794                         if( mask[x] )
 795                         {
 796                             int idx = cvFloor(*p0*a + b);
 797                             if( (unsigned)idx < (unsigned)sz )
 798                                 ((int*)H)[idx]++;
 799                         }
 800             }
 801 #endif //HAVE_TBB
 802             return;
 803         }
 804         else if( dims == 2 )
 805         {
 806 #ifdef HAVE_TBB
 807             calcHist2D_Invoker<T> body(_ptrs, _deltas, hist, _uniranges, size, dims, imsize, hstep);
 808             parallel_for(BlockedRange(0, imsize.height), body);
 809 #else
 810             double a0 = uniranges[0], b0 = uniranges[1], a1 = uniranges[2], b1 = uniranges[3];
 811             int sz0 = size[0], sz1 = size[1];
 812             int d0 = deltas[0], step0 = deltas[1],
 813                 d1 = deltas[2], step1 = deltas[3];
 814             size_t hstep0 = hstep[0];
 815             const T* p0 = (const T*)ptrs[0];
 816             const T* p1 = (const T*)ptrs[1];
 817
 818             for( ; imsize.height--; p0 += step0, p1 += step1, mask += mstep )
 819             {
 820                 if( !mask )
 821                     for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
 822                     {
 823                         int idx0 = cvFloor(*p0*a0 + b0);
 824                         int idx1 = cvFloor(*p1*a1 + b1);
 825                         if( (unsigned)idx0 < (unsigned)sz0 && (unsigned)idx1 < (unsigned)sz1 )
 826                             ((int*)(H + hstep0*idx0))[idx1]++;
 827                     }
 828                 else
 829                     for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
 830                         if( mask[x] )
 831                         {
 832                             int idx0 = cvFloor(*p0*a0 + b0);
 833                             int idx1 = cvFloor(*p1*a1 + b1);
 834                             if( (unsigned)idx0 < (unsigned)sz0 && (unsigned)idx1 < (unsigned)sz1 )
 835                                 ((int*)(H + hstep0*idx0))[idx1]++;
 836                         }
 837             }
 838 #endif //HAVE_TBB
 839             return;
 840         }
 841         else if( dims == 3 )
 842         {
 843 #ifdef HAVE_TBB
 844             if( calcHist3D_Invoker<T>::isFit(hist, imsize) )
 845             {
 846                 calcHist3D_Invoker<T> body(_ptrs, _deltas, imsize, hist, uniranges, dims, hstep, size);
 847                 parallel_for(BlockedRange(0, imsize.height), body);
 848                 return;
 849             }
 850 #endif
 851             double a0 = uniranges[0], b0 = uniranges[1],
 852                    a1 = uniranges[2], b1 = uniranges[3],
 853                    a2 = uniranges[4], b2 = uniranges[5];
 854             int sz0 = size[0], sz1 = size[1], sz2 = size[2];
 855             int d0 = deltas[0], step0 = deltas[1],
 856                 d1 = deltas[2], step1 = deltas[3],
 857                 d2 = deltas[4], step2 = deltas[5];
 858             size_t hstep0 = hstep[0], hstep1 = hstep[1];
 859             const T* p0 = (const T*)ptrs[0];
 860             const T* p1 = (const T*)ptrs[1];
 861             const T* p2 = (const T*)ptrs[2];
 862
 863             for( ; imsize.height--; p0 += step0, p1 += step1, p2 += step2, mask += mstep )
 864             {
 865                 if( !mask )
 866                     for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1, p2 += d2 )
 867                     {
 868                         int idx0 = cvFloor(*p0*a0 + b0);
 869                         int idx1 = cvFloor(*p1*a1 + b1);
 870                         int idx2 = cvFloor(*p2*a2 + b2);
 871                         if( (unsigned)idx0 < (unsigned)sz0 &&
 872                             (unsigned)idx1 < (unsigned)sz1 &&
 873                             (unsigned)idx2 < (unsigned)sz2 )
 874                             ((int*)(H + hstep0*idx0 + hstep1*idx1))[idx2]++;
 875                     }
 876                 else
 877                     for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1, p2 += d2 )
 878                         if( mask[x] )
 879                         {
 880                             int idx0 = cvFloor(*p0*a0 + b0);
 881                             int idx1 = cvFloor(*p1*a1 + b1);
 882                             int idx2 = cvFloor(*p2*a2 + b2);
 883                             if( (unsigned)idx0 < (unsigned)sz0 &&
 884                                (unsigned)idx1 < (unsigned)sz1 &&
 885                                (unsigned)idx2 < (unsigned)sz2 )
 886                                 ((int*)(H + hstep0*idx0 + hstep1*idx1))[idx2]++;
 887                         }
 888             }
 889         }
 890         else
 891         {
 892             for( ; imsize.height--; mask += mstep )
 893             {
 894                 if( !mask )
 895                     for( x = 0; x < imsize.width; x++ )
 896                     {
 897                         uchar* Hptr = H;
 898                         for( i = 0; i < dims; i++ )
 899                         {
 900                             int idx = cvFloor(*ptrs[i]*uniranges[i*2] + uniranges[i*2+1]);
 901                             if( (unsigned)idx >= (unsigned)size[i] )
 902                                 break;
 903                             ptrs[i] += deltas[i*2];
 904                             Hptr += idx*hstep[i];
 905                         }
 906
 907                         if( i == dims )
 908                             ++*((int*)Hptr);
 909                         else
 910                             for( ; i < dims; i++ )
 911                                 ptrs[i] += deltas[i*2];
 912                     }
 913                 else
 914                     for( x = 0; x < imsize.width; x++ )
 915                     {
 916                         uchar* Hptr = H;
 917                         i = 0;
 918                         if( mask[x] )
 919                             for( ; i < dims; i++ )
 920                             {
 921                                 int idx = cvFloor(*ptrs[i]*uniranges[i*2] + uniranges[i*2+1]);
 922                                 if( (unsigned)idx >= (unsigned)size[i] )
 923                                     break;
 924                                 ptrs[i] += deltas[i*2];
 925                                 Hptr += idx*hstep[i];
 926                             }
 927
 928                         if( i == dims )
 929                             ++*((int*)Hptr);
 930                         else
 931                             for( ; i < dims; i++ )
 932                                 ptrs[i] += deltas[i*2];
 933                     }
 934                 for( i = 0; i < dims; i++ )
 935                     ptrs[i] += deltas[i*2 + 1];
 936             }
 937         }
 938     }
 939     else
 940     {
 941         // non-uniform histogram
 942         const float* ranges[CV_MAX_DIM];
 943         for( i = 0; i < dims; i++ )
 944             ranges[i] = &_ranges[i][0];
 945
 946         for( ; imsize.height--; mask += mstep )
 947         {
 948             for( x = 0; x < imsize.width; x++ )
 949             {
 950                 uchar* Hptr = H;
 951                 i = 0;
 952
 953                 if( !mask || mask[x] )
 954                     for( ; i < dims; i++ )
 955                     {
 956                         float v = (float)*ptrs[i];
 957                         const float* R = ranges[i];
 958                         int idx = -1, sz = size[i];
 959
 960                         while( v >= R[idx+1] && ++idx < sz )
 961                             ; // nop
 962
 963                         if( (unsigned)idx >= (unsigned)sz )
 964                             break;
 965
 966                         ptrs[i] += deltas[i*2];
 967                         Hptr += idx*hstep[i];
 968                     }
 969
 970                 if( i == dims )
 971                     ++*((int*)Hptr);
 972                 else
 973                     for( ; i < dims; i++ )
 974                         ptrs[i] += deltas[i*2];
 975             }
 976
 977             for( i = 0; i < dims; i++ )
 978                 ptrs[i] += deltas[i*2 + 1];
 979         }
 980     }
 981 }
 982
 983
 984 static void
 985 calcHist_8u( std::vector<uchar*>& _ptrs, const std::vector<int>& _deltas,
 986              Size imsize, Mat& hist, int dims, const float** _ranges,
 987              const double* _uniranges, bool uniform )
 988 {
 989     uchar** ptrs = &_ptrs[0];
 990     const int* deltas = &_deltas[0];
 991     uchar* H = hist.ptr();
 992     int x;
 993     const uchar* mask = _ptrs[dims];
 994     int mstep = _deltas[dims*2 + 1];
 995     std::vector<size_t> _tab;
 996
 997     calcHistLookupTables_8u( hist, SparseMat(), dims, _ranges, _uniranges, uniform, false, _tab );
 998     const size_t* tab = &_tab[0];
 999
1000     if( dims == 1 )
1001     {
1002 #ifdef HAVE_TBB
1003         if( CalcHist1D_8uInvoker::isFit(hist, imsize) )
1004         {
1005             int treadsNumber = tbb::task_scheduler_init::default_num_threads();
1006             int grainSize = imsize.height/treadsNumber;
1007             tbb::mutex histogramWriteLock;
1008
1009             CalcHist1D_8uInvoker body(_ptrs, _deltas, imsize, hist, dims, _tab, &histogramWriteLock);
1010             parallel_for(BlockedRange(0, imsize.height, grainSize), body);
1011             return;
1012         }
1013 #endif
1014         int d0 = deltas[0], step0 = deltas[1];
1015         int matH[256] = { 0, };
1016         const uchar* p0 = (const uchar*)ptrs[0];
1017
1018         for( ; imsize.height--; p0 += step0, mask += mstep )
1019         {
1020             if( !mask )
1021             {
1022                 if( d0 == 1 )
1023                 {
1024                     for( x = 0; x <= imsize.width - 4; x += 4 )
1025                     {
1026                         int t0 = p0[x], t1 = p0[x+1];
1027                         matH[t0]++; matH[t1]++;
1028                         t0 = p0[x+2]; t1 = p0[x+3];
1029                         matH[t0]++; matH[t1]++;
1030                     }
1031                     p0 += x;
1032                 }
1033                 else
1034                     for( x = 0; x <= imsize.width - 4; x += 4 )
1035                     {
1036                         int t0 = p0[0], t1 = p0[d0];
1037                         matH[t0]++; matH[t1]++;
1038                         p0 += d0*2;
1039                         t0 = p0[0]; t1 = p0[d0];
1040                         matH[t0]++; matH[t1]++;
1041                         p0 += d0*2;
1042                     }
1043
1044                 for( ; x < imsize.width; x++, p0 += d0 )
1045                     matH[*p0]++;
1046             }
1047             else
1048                 for( x = 0; x < imsize.width; x++, p0 += d0 )
1049                     if( mask[x] )
1050                         matH[*p0]++;
1051         }
1052
1053         for(int i = 0; i < 256; i++ )
1054         {
1055             size_t hidx = tab[i];
1056             if( hidx < OUT_OF_RANGE )
1057                 *(int*)(H + hidx) += matH[i];
1058         }
1059     }
1060     else if( dims == 2 )
1061     {
1062 #ifdef HAVE_TBB
1063         if( CalcHist2D_8uInvoker::isFit(hist, imsize) )
1064         {
1065             callCalcHist2D_8u(_ptrs, _deltas, imsize, hist, dims, _tab);
1066             return;
1067         }
1068 #endif
1069         int d0 = deltas[0], step0 = deltas[1],
1070             d1 = deltas[2], step1 = deltas[3];
1071         const uchar* p0 = (const uchar*)ptrs[0];
1072         const uchar* p1 = (const uchar*)ptrs[1];
1073
1074         for( ; imsize.height--; p0 += step0, p1 += step1, mask += mstep )
1075         {
1076             if( !mask )
1077                 for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
1078                 {
1079                     size_t idx = tab[*p0] + tab[*p1 + 256];
1080                     if( idx < OUT_OF_RANGE )
1081                         ++*(int*)(H + idx);
1082                 }
1083             else
1084                 for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
1085                 {
1086                     size_t idx;
1087                     if( mask[x] && (idx = tab[*p0] + tab[*p1 + 256]) < OUT_OF_RANGE )
1088                         ++*(int*)(H + idx);
1089                 }
1090         }
1091     }
1092     else if( dims == 3 )
1093     {
1094 #ifdef HAVE_TBB
1095         if( CalcHist3D_8uInvoker::isFit(hist, imsize) )
1096         {
1097             callCalcHist3D_8u(_ptrs, _deltas, imsize, hist, dims, _tab);
1098             return;
1099         }
1100 #endif
1101         int d0 = deltas[0], step0 = deltas[1],
1102             d1 = deltas[2], step1 = deltas[3],
1103             d2 = deltas[4], step2 = deltas[5];
1104
1105         const uchar* p0 = (const uchar*)ptrs[0];
1106         const uchar* p1 = (const uchar*)ptrs[1];
1107         const uchar* p2 = (const uchar*)ptrs[2];
1108
1109         for( ; imsize.height--; p0 += step0, p1 += step1, p2 += step2, mask += mstep )
1110         {
1111             if( !mask )
1112                 for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1, p2 += d2 )
1113                 {
1114                     size_t idx = tab[*p0] + tab[*p1 + 256] + tab[*p2 + 512];
1115                     if( idx < OUT_OF_RANGE )
1116                         ++*(int*)(H + idx);
1117                 }
1118             else
1119                 for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1, p2 += d2 )
1120                 {
1121                     size_t idx;
1122                     if( mask[x] && (idx = tab[*p0] + tab[*p1 + 256] + tab[*p2 + 512]) < OUT_OF_RANGE )
1123                         ++*(int*)(H + idx);
1124                 }
1125         }
1126     }
1127     else
1128     {
1129         for( ; imsize.height--; mask += mstep )
1130         {
1131             if( !mask )
1132                 for( x = 0; x < imsize.width; x++ )
1133                 {
1134                     uchar* Hptr = H;
1135                     int i = 0;
1136                     for( ; i < dims; i++ )
1137                     {
1138                         size_t idx = tab[*ptrs[i] + i*256];
1139                         if( idx >= OUT_OF_RANGE )
1140                             break;
1141                         Hptr += idx;
1142                         ptrs[i] += deltas[i*2];
1143                     }
1144
1145                     if( i == dims )
1146                         ++*((int*)Hptr);
1147                     else
1148                         for( ; i < dims; i++ )
1149                             ptrs[i] += deltas[i*2];
1150                 }
1151             else
1152                 for( x = 0; x < imsize.width; x++ )
1153                 {
1154                     uchar* Hptr = H;
1155                     int i = 0;
1156                     if( mask[x] )
1157                         for( ; i < dims; i++ )
1158                         {
1159                             size_t idx = tab[*ptrs[i] + i*256];
1160                             if( idx >= OUT_OF_RANGE )
1161                                 break;
1162                             Hptr += idx;
1163                             ptrs[i] += deltas[i*2];
1164                         }
1165
1166                     if( i == dims )
1167                         ++*((int*)Hptr);
1168                     else
1169                         for( ; i < dims; i++ )
1170                             ptrs[i] += deltas[i*2];
1171                 }
1172             for(int i = 0; i < dims; i++ )
1173                 ptrs[i] += deltas[i*2 + 1];
1174         }
1175     }
1176 }
1177
1178 #ifdef HAVE_IPP
1179
1180 class IPPCalcHistInvoker :
1181     public ParallelLoopBody
1182 {
1183 public:
1184     IPPCalcHistInvoker(const Mat & _src, Mat & _hist, AutoBuffer<Ipp32s> & _levels, Ipp32s _histSize, Ipp32s _low, Ipp32s _high, bool * _ok) :
1185         ParallelLoopBody(), src(&_src), hist(&_hist), levels(&_levels), histSize(_histSize), low(_low), high(_high), ok(_ok)
1186     {
1187         *ok = true;
1188     }
1189
1190     virtual void operator() (const Range & range) const
1191     {
1192         Mat phist(hist->size(), hist->type(), Scalar::all(0));
1193
1194         IppStatus status = ippiHistogramEven_8u_C1R(
1195             src->ptr(range.start), (int)src->step, ippiSize(src->cols, range.end - range.start),
1196             phist.ptr<Ipp32s>(), (Ipp32s *)*levels, histSize, low, high);
1197
1198         if (status < 0)
1199         {
1200             *ok = false;
1201             return;
1202         }
1203         CV_IMPL_ADD(CV_IMPL_IPP|CV_IMPL_MT);
1204
1205         for (int i = 0; i < histSize; ++i)
1206             CV_XADD((int *)(hist->data + i * hist->step), *(int *)(phist.data + i * phist.step));
1207     }
1208
1209 private:
1210     const Mat * src;
1211     Mat * hist;
1212     AutoBuffer<Ipp32s> * levels;
1213     Ipp32s histSize, low, high;
1214     bool * ok;
1215
1216     const IPPCalcHistInvoker & operator = (const IPPCalcHistInvoker & );
1217 };
1218
1219 #endif
1220
1221 }
1222
1223 void cv::calcHist( const Mat* images, int nimages, const int* channels,
1224                    InputArray _mask, OutputArray _hist, int dims, const int* histSize,
1225                    const float** ranges, bool uniform, bool accumulate )
1226 {
1227     Mat mask = _mask.getMat();
1228
1229     CV_Assert(dims > 0 && histSize);
1230
1231     const uchar* const histdata = _hist.getMat().ptr();
1232     _hist.create(dims, histSize, CV_32F);
1233     Mat hist = _hist.getMat(), ihist = hist;
1234     ihist.flags = (ihist.flags & ~CV_MAT_TYPE_MASK)|CV_32S;
1235
1236 #ifdef HAVE_IPP
1237     CV_IPP_CHECK()
1238     {
1239         if (nimages == 1 && images[0].type() == CV_8UC1 && dims == 1 && channels &&
1240                 channels[0] == 0 && mask.empty() && images[0].dims <= 2 &&
1241                 !accumulate && uniform)
1242         {
1243             ihist.setTo(Scalar::all(0));
1244             AutoBuffer<Ipp32s> levels(histSize[0] + 1);
1245
1246             bool ok = true;
1247             const Mat & src = images[0];
1248             int nstripes = std::min<int>(8, static_cast<int>(src.total() / (1 << 16)));
1249 #ifdef HAVE_CONCURRENCY
1250             nstripes = 1;
1251 #endif
1252             IPPCalcHistInvoker invoker(src, ihist, levels, histSize[0] + 1, (Ipp32s)ranges[0][0], (Ipp32s)ranges[0][1], &ok);
1253             Range range(0, src.rows);
1254             parallel_for_(range, invoker, nstripes);
1255
1256             if (ok)
1257             {
1258                 ihist.convertTo(hist, CV_32F);
1259                 CV_IMPL_ADD(CV_IMPL_IPP|CV_IMPL_MT);
1260                 return;
1261             }
1262             setIppErrorStatus();
1263         }
1264     }
1265 #endif
1266
1267     if( !accumulate || histdata != hist.data )
1268         hist = Scalar(0.);
1269     else
1270         hist.convertTo(ihist, CV_32S);
1271
1272     std::vector<uchar*> ptrs;
1273     std::vector<int> deltas;
1274     std::vector<double> uniranges;
1275     Size imsize;
1276
1277     CV_Assert( mask.empty() || mask.type() == CV_8UC1 );
1278     histPrepareImages( images, nimages, channels, mask, dims, hist.size, ranges,
1279                        uniform, ptrs, deltas, imsize, uniranges );
1280     const double* _uniranges = uniform ? &uniranges[0] : 0;
1281
1282     int depth = images[0].depth();
1283
1284     if( depth == CV_8U )
1285         calcHist_8u(ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform );
1286     else if( depth == CV_16U )
1287         calcHist_<ushort>(ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform );
1288     else if( depth == CV_32F )
1289         calcHist_<float>(ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform );
1290     else
1291         CV_Error(CV_StsUnsupportedFormat, "");
1292
1293     ihist.convertTo(hist, CV_32F);
1294 }
1295
1296
1297 namespace cv
1298 {
1299
1300 template<typename T> static void
1301 calcSparseHist_( std::vector<uchar*>& _ptrs, const std::vector<int>& _deltas,
1302                  Size imsize, SparseMat& hist, int dims, const float** _ranges,
1303                  const double* _uniranges, bool uniform )
1304 {
1305     T** ptrs = (T**)&_ptrs[0];
1306     const int* deltas = &_deltas[0];
1307     int i, x;
1308     const uchar* mask = _ptrs[dims];
1309     int mstep = _deltas[dims*2 + 1];
1310     const int* size = hist.hdr->size;
1311     int idx[CV_MAX_DIM];
1312
1313     if( uniform )
1314     {
1315         const double* uniranges = &_uniranges[0];
1316
1317         for( ; imsize.height--; mask += mstep )
1318         {
1319             for( x = 0; x < imsize.width; x++ )
1320             {
1321                 i = 0;
1322                 if( !mask || mask[x] )
1323                     for( ; i < dims; i++ )
1324                     {
1325                         idx[i] = cvFloor(*ptrs[i]*uniranges[i*2] + uniranges[i*2+1]);
1326                         if( (unsigned)idx[i] >= (unsigned)size[i] )
1327                             break;
1328                         ptrs[i] += deltas[i*2];
1329                     }
1330
1331                 if( i == dims )
1332                     ++*(int*)hist.ptr(idx, true);
1333                 else
1334                     for( ; i < dims; i++ )
1335                         ptrs[i] += deltas[i*2];
1336             }
1337             for( i = 0; i < dims; i++ )
1338                 ptrs[i] += deltas[i*2 + 1];
1339         }
1340     }
1341     else
1342     {
1343         // non-uniform histogram
1344         const float* ranges[CV_MAX_DIM];
1345         for( i = 0; i < dims; i++ )
1346             ranges[i] = &_ranges[i][0];
1347
1348         for( ; imsize.height--; mask += mstep )
1349         {
1350             for( x = 0; x < imsize.width; x++ )
1351             {
1352                 i = 0;
1353
1354                 if( !mask || mask[x] )
1355                     for( ; i < dims; i++ )
1356                     {
1357                         float v = (float)*ptrs[i];
1358                         const float* R = ranges[i];
1359                         int j = -1, sz = size[i];
1360
1361                         while( v >= R[j+1] && ++j < sz )
1362                             ; // nop
1363
1364                         if( (unsigned)j >= (unsigned)sz )
1365                             break;
1366                         ptrs[i] += deltas[i*2];
1367                         idx[i] = j;
1368                     }
1369
1370                 if( i == dims )
1371                     ++*(int*)hist.ptr(idx, true);
1372                 else
1373                     for( ; i < dims; i++ )
1374                         ptrs[i] += deltas[i*2];
1375             }
1376
1377             for( i = 0; i < dims; i++ )
1378                 ptrs[i] += deltas[i*2 + 1];
1379         }
1380     }
1381 }
1382
1383
1384 static void
1385 calcSparseHist_8u( std::vector<uchar*>& _ptrs, const std::vector<int>& _deltas,
1386                    Size imsize, SparseMat& hist, int dims, const float** _ranges,
1387                    const double* _uniranges, bool uniform )
1388 {
1389     uchar** ptrs = (uchar**)&_ptrs[0];
1390     const int* deltas = &_deltas[0];
1391     int x;
1392     const uchar* mask = _ptrs[dims];
1393     int mstep = _deltas[dims*2 + 1];
1394     int idx[CV_MAX_DIM];
1395     std::vector<size_t> _tab;
1396
1397     calcHistLookupTables_8u( Mat(), hist, dims, _ranges, _uniranges, uniform, true, _tab );
1398     const size_t* tab = &_tab[0];
1399
1400     for( ; imsize.height--; mask += mstep )
1401     {
1402         for( x = 0; x < imsize.width; x++ )
1403         {
1404             int i = 0;
1405             if( !mask || mask[x] )
1406                 for( ; i < dims; i++ )
1407                 {
1408                     size_t hidx = tab[*ptrs[i] + i*256];
1409                     if( hidx >= OUT_OF_RANGE )
1410                         break;
1411                     ptrs[i] += deltas[i*2];
1412                     idx[i] = (int)hidx;
1413                 }
1414
1415             if( i == dims )
1416                 ++*(int*)hist.ptr(idx,true);
1417             else
1418                 for( ; i < dims; i++ )
1419                     ptrs[i] += deltas[i*2];
1420         }
1421         for(int i = 0; i < dims; i++ )
1422             ptrs[i] += deltas[i*2 + 1];
1423     }
1424 }
1425
1426
1427 static void calcHist( const Mat* images, int nimages, const int* channels,
1428                       const Mat& mask, SparseMat& hist, int dims, const int* histSize,
1429                       const float** ranges, bool uniform, bool accumulate, bool keepInt )
1430 {
1431     size_t i, N;
1432
1433     if( !accumulate )
1434         hist.create(dims, histSize, CV_32F);
1435     else
1436     {
1437         SparseMatIterator it = hist.begin();
1438         for( i = 0, N = hist.nzcount(); i < N; i++, ++it )
1439         {
1440             Cv32suf* val = (Cv32suf*)it.ptr;
1441             val->i = cvRound(val->f);
1442         }
1443     }
1444
1445     std::vector<uchar*> ptrs;
1446     std::vector<int> deltas;
1447     std::vector<double> uniranges;
1448     Size imsize;
1449
1450     CV_Assert( mask.empty() || mask.type() == CV_8UC1 );
1451     histPrepareImages( images, nimages, channels, mask, dims, hist.hdr->size, ranges,
1452                        uniform, ptrs, deltas, imsize, uniranges );
1453     const double* _uniranges = uniform ? &uniranges[0] : 0;
1454
1455     int depth = images[0].depth();
1456     if( depth == CV_8U )
1457         calcSparseHist_8u(ptrs, deltas, imsize, hist, dims, ranges, _uniranges, uniform );
1458     else if( depth == CV_16U )
1459         calcSparseHist_<ushort>(ptrs, deltas, imsize, hist, dims, ranges, _uniranges, uniform );
1460     else if( depth == CV_32F )
1461         calcSparseHist_<float>(ptrs, deltas, imsize, hist, dims, ranges, _uniranges, uniform );
1462     else
1463         CV_Error(CV_StsUnsupportedFormat, "");
1464
1465     if( !keepInt )
1466     {
1467         SparseMatIterator it = hist.begin();
1468         for( i = 0, N = hist.nzcount(); i < N; i++, ++it )
1469         {
1470             Cv32suf* val = (Cv32suf*)it.ptr;
1471             val->f = (float)val->i;
1472         }
1473     }
1474 }
1475
1476 #ifdef HAVE_OPENCL
1477
1478 enum
1479 {
1480     BINS = 256
1481 };
1482
1483 static bool ocl_calcHist1(InputArray _src, OutputArray _hist, int ddepth = CV_32S)
1484 {
1485     const ocl::Device & dev = ocl::Device::getDefault();
1486     int compunits = dev.maxComputeUnits();
1487     size_t wgs = dev.maxWorkGroupSize();
1488     Size size = _src.size();
1489     bool use16 = size.width % 16 == 0 && _src.offset() % 16 == 0 && _src.step() % 16 == 0;
1490     int kercn = dev.isAMD() && use16 ? 16 : std::min(4, ocl::predictOptimalVectorWidth(_src));
1491
1492     ocl::Kernel k1("calculate_histogram", ocl::imgproc::histogram_oclsrc,
1493                    format("-D BINS=%d -D HISTS_COUNT=%d -D WGS=%d -D kercn=%d -D T=%s%s",
1494                           BINS, compunits, wgs, kercn,
1495                           kercn == 4 ? "int" : ocl::typeToStr(CV_8UC(kercn)),
1496                           _src.isContinuous() ? " -D HAVE_SRC_CONT" : ""));
1497     if (k1.empty())
1498         return false;
1499
1500     _hist.create(BINS, 1, ddepth);
1501     UMat src = _src.getUMat(), ghist(1, BINS * compunits, CV_32SC1),
1502             hist = _hist.getUMat();
1503
1504     k1.args(ocl::KernelArg::ReadOnly(src),
1505             ocl::KernelArg::PtrWriteOnly(ghist), (int)src.total());
1506
1507     size_t globalsize = compunits * wgs;
1508     if (!k1.run(1, &globalsize, &wgs, false))
1509         return false;
1510
1511     char cvt[40];
1512     ocl::Kernel k2("merge_histogram", ocl::imgproc::histogram_oclsrc,
1513                    format("-D BINS=%d -D HISTS_COUNT=%d -D WGS=%d -D convertToHT=%s -D HT=%s",
1514                           BINS, compunits, (int)wgs, ocl::convertTypeStr(CV_32S, ddepth, 1, cvt),
1515                           ocl::typeToStr(ddepth)));
1516     if (k2.empty())
1517         return false;
1518
1519     k2.args(ocl::KernelArg::PtrReadOnly(ghist),
1520             ocl::KernelArg::WriteOnlyNoSize(hist));
1521
1522     return k2.run(1, &wgs, &wgs, false);
1523 }
1524
1525 static bool ocl_calcHist(InputArrayOfArrays images, OutputArray hist)
1526 {
1527     std::vector<UMat> v;
1528     images.getUMatVector(v);
1529
1530     return ocl_calcHist1(v[0], hist, CV_32F);
1531 }
1532
1533 #endif
1534
1535 }
1536
1537 void cv::calcHist( const Mat* images, int nimages, const int* channels,
1538                InputArray _mask, SparseMat& hist, int dims, const int* histSize,
1539                const float** ranges, bool uniform, bool accumulate )
1540 {
1541     Mat mask = _mask.getMat();
1542     calcHist( images, nimages, channels, mask, hist, dims, histSize,
1543               ranges, uniform, accumulate, false );
1544 }
1545
1546
1547 void cv::calcHist( InputArrayOfArrays images, const std::vector<int>& channels,
1548                    InputArray mask, OutputArray hist,
1549                    const std::vector<int>& histSize,
1550                    const std::vector<float>& ranges,
1551                    bool accumulate )
1552 {
1553     CV_OCL_RUN(images.total() == 1 && channels.size() == 1 && images.channels(0) == 1 &&
1554                channels[0] == 0 && images.isUMatVector() && mask.empty() && !accumulate &&
1555                histSize.size() == 1 && histSize[0] == BINS && ranges.size() == 2 &&
1556                ranges[0] == 0 && ranges[1] == BINS,
1557                ocl_calcHist(images, hist))
1558
1559     int i, dims = (int)histSize.size(), rsz = (int)ranges.size(), csz = (int)channels.size();
1560     int nimages = (int)images.total();
1561
1562     CV_Assert(nimages > 0 && dims > 0);
1563     CV_Assert(rsz == dims*2 || (rsz == 0 && images.depth(0) == CV_8U));
1564     CV_Assert(csz == 0 || csz == dims);
1565     float* _ranges[CV_MAX_DIM];
1566     if( rsz > 0 )
1567     {
1568         for( i = 0; i < rsz/2; i++ )
1569             _ranges[i] = (float*)&ranges[i*2];
1570     }
1571
1572     AutoBuffer<Mat> buf(nimages);
1573     for( i = 0; i < nimages; i++ )
1574         buf[i] = images.getMat(i);
1575
1576     calcHist(&buf[0], nimages, csz ? &channels[0] : 0,
1577             mask, hist, dims, &histSize[0], rsz ? (const float**)_ranges : 0,
1578             true, accumulate);
1579 }
1580
1581
1582 /////////////////////////////////////// B A C K   P R O J E C T ////////////////////////////////////
1583
1584 namespace cv
1585 {
1586
1587 template<typename T, typename BT> static void
1588 calcBackProj_( std::vector<uchar*>& _ptrs, const std::vector<int>& _deltas,
1589                Size imsize, const Mat& hist, int dims, const float** _ranges,
1590                const double* _uniranges, float scale, bool uniform )
1591 {
1592     T** ptrs = (T**)&_ptrs[0];
1593     const int* deltas = &_deltas[0];
1594     const uchar* H = hist.ptr();
1595     int i, x;
1596     BT* bproj = (BT*)_ptrs[dims];
1597     int bpstep = _deltas[dims*2 + 1];
1598     int size[CV_MAX_DIM];
1599     size_t hstep[CV_MAX_DIM];
1600
1601     for( i = 0; i < dims; i++ )
1602     {
1603         size[i] = hist.size[i];
1604         hstep[i] = hist.step[i];
1605     }
1606
1607     if( uniform )
1608     {
1609         const double* uniranges = &_uniranges[0];
1610
1611         if( dims == 1 )
1612         {
1613             double a = uniranges[0], b = uniranges[1];
1614             int sz = size[0], d0 = deltas[0], step0 = deltas[1];
1615             const T* p0 = (const T*)ptrs[0];
1616
1617             for( ; imsize.height--; p0 += step0, bproj += bpstep )
1618             {
1619                 for( x = 0; x < imsize.width; x++, p0 += d0 )
1620                 {
1621                     int idx = cvFloor(*p0*a + b);
1622                     bproj[x] = (unsigned)idx < (unsigned)sz ? saturate_cast<BT>(((const float*)H)[idx]*scale) : 0;
1623                 }
1624             }
1625         }
1626         else if( dims == 2 )
1627         {
1628             double a0 = uniranges[0], b0 = uniranges[1],
1629                    a1 = uniranges[2], b1 = uniranges[3];
1630             int sz0 = size[0], sz1 = size[1];
1631             int d0 = deltas[0], step0 = deltas[1],
1632                 d1 = deltas[2], step1 = deltas[3];
1633             size_t hstep0 = hstep[0];
1634             const T* p0 = (const T*)ptrs[0];
1635             const T* p1 = (const T*)ptrs[1];
1636
1637             for( ; imsize.height--; p0 += step0, p1 += step1, bproj += bpstep )
1638             {
1639                 for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
1640                 {
1641                     int idx0 = cvFloor(*p0*a0 + b0);
1642                     int idx1 = cvFloor(*p1*a1 + b1);
1643                     bproj[x] = (unsigned)idx0 < (unsigned)sz0 &&
1644                                (unsigned)idx1 < (unsigned)sz1 ?
1645                         saturate_cast<BT>(((const float*)(H + hstep0*idx0))[idx1]*scale) : 0;
1646                 }
1647             }
1648         }
1649         else if( dims == 3 )
1650         {
1651             double a0 = uniranges[0], b0 = uniranges[1],
1652                    a1 = uniranges[2], b1 = uniranges[3],
1653                    a2 = uniranges[4], b2 = uniranges[5];
1654             int sz0 = size[0], sz1 = size[1], sz2 = size[2];
1655             int d0 = deltas[0], step0 = deltas[1],
1656                 d1 = deltas[2], step1 = deltas[3],
1657                 d2 = deltas[4], step2 = deltas[5];
1658             size_t hstep0 = hstep[0], hstep1 = hstep[1];
1659             const T* p0 = (const T*)ptrs[0];
1660             const T* p1 = (const T*)ptrs[1];
1661             const T* p2 = (const T*)ptrs[2];
1662
1663             for( ; imsize.height--; p0 += step0, p1 += step1, p2 += step2, bproj += bpstep )
1664             {
1665                 for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1, p2 += d2 )
1666                 {
1667                     int idx0 = cvFloor(*p0*a0 + b0);
1668                     int idx1 = cvFloor(*p1*a1 + b1);
1669                     int idx2 = cvFloor(*p2*a2 + b2);
1670                     bproj[x] = (unsigned)idx0 < (unsigned)sz0 &&
1671                                (unsigned)idx1 < (unsigned)sz1 &&
1672                                (unsigned)idx2 < (unsigned)sz2 ?
1673                         saturate_cast<BT>(((const float*)(H + hstep0*idx0 + hstep1*idx1))[idx2]*scale) : 0;
1674                 }
1675             }
1676         }
1677         else
1678         {
1679             for( ; imsize.height--; bproj += bpstep )
1680             {
1681                 for( x = 0; x < imsize.width; x++ )
1682                 {
1683                     const uchar* Hptr = H;
1684                     for( i = 0; i < dims; i++ )
1685                     {
1686                         int idx = cvFloor(*ptrs[i]*uniranges[i*2] + uniranges[i*2+1]);
1687                         if( (unsigned)idx >= (unsigned)size[i] || (_ranges && *ptrs[i] >= _ranges[i][1]))
1688                             break;
1689                         ptrs[i] += deltas[i*2];
1690                         Hptr += idx*hstep[i];
1691                     }
1692
1693                     if( i == dims )
1694                         bproj[x] = saturate_cast<BT>(*(const float*)Hptr*scale);
1695                     else
1696                     {
1697                         bproj[x] = 0;
1698                         for( ; i < dims; i++ )
1699                             ptrs[i] += deltas[i*2];
1700                     }
1701                 }
1702                 for( i = 0; i < dims; i++ )
1703                     ptrs[i] += deltas[i*2 + 1];
1704             }
1705         }
1706     }
1707     else
1708     {
1709         // non-uniform histogram
1710         const float* ranges[CV_MAX_DIM];
1711         for( i = 0; i < dims; i++ )
1712             ranges[i] = &_ranges[i][0];
1713
1714         for( ; imsize.height--; bproj += bpstep )
1715         {
1716             for( x = 0; x < imsize.width; x++ )
1717             {
1718                 const uchar* Hptr = H;
1719                 for( i = 0; i < dims; i++ )
1720                 {
1721                     float v = (float)*ptrs[i];
1722                     const float* R = ranges[i];
1723                     int idx = -1, sz = size[i];
1724
1725                     while( v >= R[idx+1] && ++idx < sz )
1726                         ; // nop
1727
1728                     if( (unsigned)idx >= (unsigned)sz )
1729                         break;
1730
1731                     ptrs[i] += deltas[i*2];
1732                     Hptr += idx*hstep[i];
1733                 }
1734
1735                 if( i == dims )
1736                     bproj[x] = saturate_cast<BT>(*(const float*)Hptr*scale);
1737                 else
1738                 {
1739                     bproj[x] = 0;
1740                     for( ; i < dims; i++ )
1741                         ptrs[i] += deltas[i*2];
1742                 }
1743             }
1744
1745             for( i = 0; i < dims; i++ )
1746                 ptrs[i] += deltas[i*2 + 1];
1747         }
1748     }
1749 }
1750
1751
1752 static void
1753 calcBackProj_8u( std::vector<uchar*>& _ptrs, const std::vector<int>& _deltas,
1754                  Size imsize, const Mat& hist, int dims, const float** _ranges,
1755                  const double* _uniranges, float scale, bool uniform )
1756 {
1757     uchar** ptrs = &_ptrs[0];
1758     const int* deltas = &_deltas[0];
1759     const uchar* H = hist.ptr();
1760     int i, x;
1761     uchar* bproj = _ptrs[dims];
1762     int bpstep = _deltas[dims*2 + 1];
1763     std::vector<size_t> _tab;
1764
1765     calcHistLookupTables_8u( hist, SparseMat(), dims, _ranges, _uniranges, uniform, false, _tab );
1766     const size_t* tab = &_tab[0];
1767
1768     if( dims == 1 )
1769     {
1770         int d0 = deltas[0], step0 = deltas[1];
1771         uchar matH[256] = {0};
1772         const uchar* p0 = (const uchar*)ptrs[0];
1773
1774         for( i = 0; i < 256; i++ )
1775         {
1776             size_t hidx = tab[i];
1777             if( hidx < OUT_OF_RANGE )
1778                 matH[i] = saturate_cast<uchar>(*(float*)(H + hidx)*scale);
1779         }
1780
1781         for( ; imsize.height--; p0 += step0, bproj += bpstep )
1782         {
1783             if( d0 == 1 )
1784             {
1785                 for( x = 0; x <= imsize.width - 4; x += 4 )
1786                 {
1787                     uchar t0 = matH[p0[x]], t1 = matH[p0[x+1]];
1788                     bproj[x] = t0; bproj[x+1] = t1;
1789                     t0 = matH[p0[x+2]]; t1 = matH[p0[x+3]];
1790                     bproj[x+2] = t0; bproj[x+3] = t1;
1791                 }
1792                 p0 += x;
1793             }
1794             else
1795                 for( x = 0; x <= imsize.width - 4; x += 4 )
1796                 {
1797                     uchar t0 = matH[p0[0]], t1 = matH[p0[d0]];
1798                     bproj[x] = t0; bproj[x+1] = t1;
1799                     p0 += d0*2;
1800                     t0 = matH[p0[0]]; t1 = matH[p0[d0]];
1801                     bproj[x+2] = t0; bproj[x+3] = t1;
1802                     p0 += d0*2;
1803                 }
1804
1805             for( ; x < imsize.width; x++, p0 += d0 )
1806                 bproj[x] = matH[*p0];
1807         }
1808     }
1809     else if( dims == 2 )
1810     {
1811         int d0 = deltas[0], step0 = deltas[1],
1812             d1 = deltas[2], step1 = deltas[3];
1813         const uchar* p0 = (const uchar*)ptrs[0];
1814         const uchar* p1 = (const uchar*)ptrs[1];
1815
1816         for( ; imsize.height--; p0 += step0, p1 += step1, bproj += bpstep )
1817         {
1818             for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
1819             {
1820                 size_t idx = tab[*p0] + tab[*p1 + 256];
1821                 bproj[x] = idx < OUT_OF_RANGE ? saturate_cast<uchar>(*(const float*)(H + idx)*scale) : 0;
1822             }
1823         }
1824     }
1825     else if( dims == 3 )
1826     {
1827         int d0 = deltas[0], step0 = deltas[1],
1828         d1 = deltas[2], step1 = deltas[3],
1829         d2 = deltas[4], step2 = deltas[5];
1830         const uchar* p0 = (const uchar*)ptrs[0];
1831         const uchar* p1 = (const uchar*)ptrs[1];
1832         const uchar* p2 = (const uchar*)ptrs[2];
1833
1834         for( ; imsize.height--; p0 += step0, p1 += step1, p2 += step2, bproj += bpstep )
1835         {
1836             for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1, p2 += d2 )
1837             {
1838                 size_t idx = tab[*p0] + tab[*p1 + 256] + tab[*p2 + 512];
1839                 bproj[x] = idx < OUT_OF_RANGE ? saturate_cast<uchar>(*(const float*)(H + idx)*scale) : 0;
1840             }
1841         }
1842     }
1843     else
1844     {
1845         for( ; imsize.height--; bproj += bpstep )
1846         {
1847             for( x = 0; x < imsize.width; x++ )
1848             {
1849                 const uchar* Hptr = H;
1850                 for( i = 0; i < dims; i++ )
1851                 {
1852                     size_t idx = tab[*ptrs[i] + i*256];
1853                     if( idx >= OUT_OF_RANGE )
1854                         break;
1855                     ptrs[i] += deltas[i*2];
1856                     Hptr += idx;
1857                 }
1858
1859                 if( i == dims )
1860                     bproj[x] = saturate_cast<uchar>(*(const float*)Hptr*scale);
1861                 else
1862                 {
1863                     bproj[x] = 0;
1864                     for( ; i < dims; i++ )
1865                         ptrs[i] += deltas[i*2];
1866                 }
1867             }
1868             for( i = 0; i < dims; i++ )
1869                 ptrs[i] += deltas[i*2 + 1];
1870         }
1871     }
1872 }
1873
1874 }
1875
1876 void cv::calcBackProject( const Mat* images, int nimages, const int* channels,
1877                           InputArray _hist, OutputArray _backProject,
1878                           const float** ranges, double scale, bool uniform )
1879 {
1880     Mat hist = _hist.getMat();
1881     std::vector<uchar*> ptrs;
1882     std::vector<int> deltas;
1883     std::vector<double> uniranges;
1884     Size imsize;
1885     int dims = hist.dims == 2 && hist.size[1] == 1 ? 1 : hist.dims;
1886
1887     CV_Assert( dims > 0 && !hist.empty() );
1888     _backProject.create( images[0].size(), images[0].depth() );
1889     Mat backProject = _backProject.getMat();
1890     histPrepareImages( images, nimages, channels, backProject, dims, hist.size, ranges,
1891                        uniform, ptrs, deltas, imsize, uniranges );
1892     const double* _uniranges = uniform ? &uniranges[0] : 0;
1893
1894     int depth = images[0].depth();
1895     if( depth == CV_8U )
1896         calcBackProj_8u(ptrs, deltas, imsize, hist, dims, ranges, _uniranges, (float)scale, uniform);
1897     else if( depth == CV_16U )
1898         calcBackProj_<ushort, ushort>(ptrs, deltas, imsize, hist, dims, ranges, _uniranges, (float)scale, uniform );
1899     else if( depth == CV_32F )
1900         calcBackProj_<float, float>(ptrs, deltas, imsize, hist, dims, ranges, _uniranges, (float)scale, uniform );
1901     else
1902         CV_Error(CV_StsUnsupportedFormat, "");
1903 }
1904
1905
1906 namespace cv
1907 {
1908
1909 template<typename T, typename BT> static void
1910 calcSparseBackProj_( std::vector<uchar*>& _ptrs, const std::vector<int>& _deltas,
1911                      Size imsize, const SparseMat& hist, int dims, const float** _ranges,
1912                      const double* _uniranges, float scale, bool uniform )
1913 {
1914     T** ptrs = (T**)&_ptrs[0];
1915     const int* deltas = &_deltas[0];
1916     int i, x;
1917     BT* bproj = (BT*)_ptrs[dims];
1918     int bpstep = _deltas[dims*2 + 1];
1919     const int* size = hist.hdr->size;
1920     int idx[CV_MAX_DIM];
1921     const SparseMat_<float>& hist_ = (const SparseMat_<float>&)hist;
1922
1923     if( uniform )
1924     {
1925         const double* uniranges = &_uniranges[0];
1926         for( ; imsize.height--; bproj += bpstep )
1927         {
1928             for( x = 0; x < imsize.width; x++ )
1929             {
1930                 for( i = 0; i < dims; i++ )
1931                 {
1932                     idx[i] = cvFloor(*ptrs[i]*uniranges[i*2] + uniranges[i*2+1]);
1933                     if( (unsigned)idx[i] >= (unsigned)size[i] )
1934                         break;
1935                     ptrs[i] += deltas[i*2];
1936                 }
1937
1938                 if( i == dims )
1939                     bproj[x] = saturate_cast<BT>(hist_(idx)*scale);
1940                 else
1941                 {
1942                     bproj[x] = 0;
1943                     for( ; i < dims; i++ )
1944                         ptrs[i] += deltas[i*2];
1945                 }
1946             }
1947             for( i = 0; i < dims; i++ )
1948                 ptrs[i] += deltas[i*2 + 1];
1949         }
1950     }
1951     else
1952     {
1953         // non-uniform histogram
1954         const float* ranges[CV_MAX_DIM];
1955         for( i = 0; i < dims; i++ )
1956             ranges[i] = &_ranges[i][0];
1957
1958         for( ; imsize.height--; bproj += bpstep )
1959         {
1960             for( x = 0; x < imsize.width; x++ )
1961             {
1962                 for( i = 0; i < dims; i++ )
1963                 {
1964                     float v = (float)*ptrs[i];
1965                     const float* R = ranges[i];
1966                     int j = -1, sz = size[i];
1967
1968                     while( v >= R[j+1] && ++j < sz )
1969                         ; // nop
1970
1971                     if( (unsigned)j >= (unsigned)sz )
1972                         break;
1973                     idx[i] = j;
1974                     ptrs[i] += deltas[i*2];
1975                 }
1976
1977                 if( i == dims )
1978                     bproj[x] = saturate_cast<BT>(hist_(idx)*scale);
1979                 else
1980                 {
1981                     bproj[x] = 0;
1982                     for( ; i < dims; i++ )
1983                         ptrs[i] += deltas[i*2];
1984                 }
1985             }
1986
1987             for( i = 0; i < dims; i++ )
1988                 ptrs[i] += deltas[i*2 + 1];
1989         }
1990     }
1991 }
1992
1993
1994 static void
1995 calcSparseBackProj_8u( std::vector<uchar*>& _ptrs, const std::vector<int>& _deltas,
1996                        Size imsize, const SparseMat& hist, int dims, const float** _ranges,
1997                        const double* _uniranges, float scale, bool uniform )
1998 {
1999     uchar** ptrs = &_ptrs[0];
2000     const int* deltas = &_deltas[0];
2001     int i, x;
2002     uchar* bproj = _ptrs[dims];
2003     int bpstep = _deltas[dims*2 + 1];
2004     std::vector<size_t> _tab;
2005     int idx[CV_MAX_DIM];
2006
2007     calcHistLookupTables_8u( Mat(), hist, dims, _ranges, _uniranges, uniform, true, _tab );
2008     const size_t* tab = &_tab[0];
2009
2010     for( ; imsize.height--; bproj += bpstep )
2011     {
2012         for( x = 0; x < imsize.width; x++ )
2013         {
2014             for( i = 0; i < dims; i++ )
2015             {
2016                 size_t hidx = tab[*ptrs[i] + i*256];
2017                 if( hidx >= OUT_OF_RANGE )
2018                     break;
2019                 idx[i] = (int)hidx;
2020                 ptrs[i] += deltas[i*2];
2021             }
2022
2023             if( i == dims )
2024                 bproj[x] = saturate_cast<uchar>(hist.value<float>(idx)*scale);
2025             else
2026             {
2027                 bproj[x] = 0;
2028                 for( ; i < dims; i++ )
2029                     ptrs[i] += deltas[i*2];
2030             }
2031         }
2032         for( i = 0; i < dims; i++ )
2033             ptrs[i] += deltas[i*2 + 1];
2034     }
2035 }
2036
2037 }
2038
2039 void cv::calcBackProject( const Mat* images, int nimages, const int* channels,
2040                           const SparseMat& hist, OutputArray _backProject,
2041                           const float** ranges, double scale, bool uniform )
2042 {
2043     std::vector<uchar*> ptrs;
2044     std::vector<int> deltas;
2045     std::vector<double> uniranges;
2046     Size imsize;
2047     int dims = hist.dims();
2048
2049     CV_Assert( dims > 0 );
2050     _backProject.create( images[0].size(), images[0].depth() );
2051     Mat backProject = _backProject.getMat();
2052     histPrepareImages( images, nimages, channels, backProject,
2053                        dims, hist.hdr->size, ranges,
2054                        uniform, ptrs, deltas, imsize, uniranges );
2055     const double* _uniranges = uniform ? &uniranges[0] : 0;
2056
2057     int depth = images[0].depth();
2058     if( depth == CV_8U )
2059         calcSparseBackProj_8u(ptrs, deltas, imsize, hist, dims, ranges,
2060                               _uniranges, (float)scale, uniform);
2061     else if( depth == CV_16U )
2062         calcSparseBackProj_<ushort, ushort>(ptrs, deltas, imsize, hist, dims, ranges,
2063                                           _uniranges, (float)scale, uniform );
2064     else if( depth == CV_32F )
2065         calcSparseBackProj_<float, float>(ptrs, deltas, imsize, hist, dims, ranges,
2066                                           _uniranges, (float)scale, uniform );
2067     else
2068         CV_Error(CV_StsUnsupportedFormat, "");
2069 }
2070
2071 #ifdef HAVE_OPENCL
2072
2073 namespace cv {
2074
2075 static void getUMatIndex(const std::vector<UMat> & um, int cn, int & idx, int & cnidx)
2076 {
2077     int totalChannels = 0;
2078     for (size_t i = 0, size = um.size(); i < size; ++i)
2079     {
2080         int ccn = um[i].channels();
2081         totalChannels += ccn;
2082
2083         if (totalChannels == cn)
2084         {
2085             idx = (int)(i + 1);
2086             cnidx = 0;
2087             return;
2088         }
2089         else if (totalChannels > cn)
2090         {
2091             idx = (int)i;
2092             cnidx = i == 0 ? cn : (cn - totalChannels + ccn);
2093             return;
2094         }
2095     }
2096
2097     idx = cnidx = -1;
2098 }
2099
2100 static bool ocl_calcBackProject( InputArrayOfArrays _images, std::vector<int> channels,
2101                                  InputArray _hist, OutputArray _dst,
2102                                  const std::vector<float>& ranges,
2103                                  float scale, size_t histdims )
2104 {
2105     std::vector<UMat> images;
2106     _images.getUMatVector(images);
2107
2108     size_t nimages = images.size(), totalcn = images[0].channels();
2109
2110     CV_Assert(nimages > 0);
2111     Size size = images[0].size();
2112     int depth = images[0].depth();
2113
2114     //kernels are valid for this type only
2115     if (depth != CV_8U)
2116         return false;
2117
2118     for (size_t i = 1; i < nimages; ++i)
2119     {
2120         const UMat & m = images[i];
2121         totalcn += m.channels();
2122         CV_Assert(size == m.size() && depth == m.depth());
2123     }
2124
2125     std::sort(channels.begin(), channels.end());
2126     for (size_t i = 0; i < histdims; ++i)
2127         CV_Assert(channels[i] < (int)totalcn);
2128
2129     if (histdims == 1)
2130     {
2131         int idx, cnidx;
2132         getUMatIndex(images, channels[0], idx, cnidx);
2133         CV_Assert(idx >= 0);
2134         UMat im = images[idx];
2135
2136         String opts = format("-D histdims=1 -D scn=%d", im.channels());
2137         ocl::Kernel lutk("calcLUT", ocl::imgproc::calc_back_project_oclsrc, opts);
2138         if (lutk.empty())
2139             return false;
2140
2141         size_t lsize = 256;
2142         UMat lut(1, (int)lsize, CV_32SC1), hist = _hist.getUMat(), uranges(ranges, true);
2143
2144         lutk.args(ocl::KernelArg::ReadOnlyNoSize(hist), hist.rows,
2145                   ocl::KernelArg::PtrWriteOnly(lut), scale, ocl::KernelArg::PtrReadOnly(uranges));
2146         if (!lutk.run(1, &lsize, NULL, false))
2147             return false;
2148
2149         ocl::Kernel mapk("LUT", ocl::imgproc::calc_back_project_oclsrc, opts);
2150         if (mapk.empty())
2151             return false;
2152
2153         _dst.create(size, depth);
2154         UMat dst = _dst.getUMat();
2155
2156         im.offset += cnidx;
2157         mapk.args(ocl::KernelArg::ReadOnlyNoSize(im), ocl::KernelArg::PtrReadOnly(lut),
2158                   ocl::KernelArg::WriteOnly(dst));
2159
2160         size_t globalsize[2] = { size.width, size.height };
2161         return mapk.run(2, globalsize, NULL, false);
2162     }
2163     else if (histdims == 2)
2164     {
2165         int idx0, idx1, cnidx0, cnidx1;
2166         getUMatIndex(images, channels[0], idx0, cnidx0);
2167         getUMatIndex(images, channels[1], idx1, cnidx1);
2168         CV_Assert(idx0 >= 0 && idx1 >= 0);
2169         UMat im0 = images[idx0], im1 = images[idx1];
2170
2171         // Lut for the first dimension
2172         String opts = format("-D histdims=2 -D scn1=%d -D scn2=%d", im0.channels(), im1.channels());
2173         ocl::Kernel lutk1("calcLUT", ocl::imgproc::calc_back_project_oclsrc, opts);
2174         if (lutk1.empty())
2175             return false;
2176
2177         size_t lsize = 256;
2178         UMat lut(1, (int)lsize<<1, CV_32SC1), uranges(ranges, true), hist = _hist.getUMat();
2179
2180         lutk1.args(hist.rows, ocl::KernelArg::PtrWriteOnly(lut), (int)0, ocl::KernelArg::PtrReadOnly(uranges), (int)0);
2181         if (!lutk1.run(1, &lsize, NULL, false))
2182             return false;
2183
2184         // lut for the second dimension
2185         ocl::Kernel lutk2("calcLUT", ocl::imgproc::calc_back_project_oclsrc, opts);
2186         if (lutk2.empty())
2187             return false;
2188
2189         lut.offset += lsize * sizeof(int);
2190         lutk2.args(hist.cols, ocl::KernelArg::PtrWriteOnly(lut), (int)256, ocl::KernelArg::PtrReadOnly(uranges), (int)2);
2191         if (!lutk2.run(1, &lsize, NULL, false))
2192             return false;
2193
2194         // perform lut
2195         ocl::Kernel mapk("LUT", ocl::imgproc::calc_back_project_oclsrc, opts);
2196         if (mapk.empty())
2197             return false;
2198
2199         _dst.create(size, depth);
2200         UMat dst = _dst.getUMat();
2201
2202         im0.offset += cnidx0;
2203         im1.offset += cnidx1;
2204         mapk.args(ocl::KernelArg::ReadOnlyNoSize(im0), ocl::KernelArg::ReadOnlyNoSize(im1),
2205                ocl::KernelArg::ReadOnlyNoSize(hist), ocl::KernelArg::PtrReadOnly(lut), scale, ocl::KernelArg::WriteOnly(dst));
2206
2207         size_t globalsize[2] = { size.width, size.height };
2208         return mapk.run(2, globalsize, NULL, false);
2209     }
2210     return false;
2211 }
2212
2213 }
2214
2215 #endif
2216
2217 void cv::calcBackProject( InputArrayOfArrays images, const std::vector<int>& channels,
2218                           InputArray hist, OutputArray dst,
2219                           const std::vector<float>& ranges,
2220                           double scale )
2221 {
2222 #ifdef HAVE_OPENCL
2223     Size histSize = hist.size();
2224     bool _1D = histSize.height == 1 || histSize.width == 1;
2225     size_t histdims = _1D ? 1 : hist.dims();
2226 #endif
2227
2228     CV_OCL_RUN(dst.isUMat() && hist.type() == CV_32FC1 &&
2229                histdims <= 2 && ranges.size() == histdims * 2 && histdims == channels.size(),
2230                ocl_calcBackProject(images, channels, hist, dst, ranges, (float)scale, histdims))
2231
2232     Mat H0 = hist.getMat(), H;
2233     int hcn = H0.channels();
2234
2235     if( hcn > 1 )
2236     {
2237         CV_Assert( H0.isContinuous() );
2238         int hsz[CV_CN_MAX+1];
2239         memcpy(hsz, &H0.size[0], H0.dims*sizeof(hsz[0]));
2240         hsz[H0.dims] = hcn;
2241         H = Mat(H0.dims+1, hsz, H0.depth(), H0.ptr());
2242     }
2243     else
2244         H = H0;
2245
2246     bool _1d = H.rows == 1 || H.cols == 1;
2247     int i, dims = H.dims, rsz = (int)ranges.size(), csz = (int)channels.size();
2248     int nimages = (int)images.total();
2249
2250     CV_Assert(nimages > 0);
2251     CV_Assert(rsz == dims*2 || (rsz == 2 && _1d) || (rsz == 0 && images.depth(0) == CV_8U));
2252     CV_Assert(csz == 0 || csz == dims || (csz == 1 && _1d));
2253
2254     float* _ranges[CV_MAX_DIM];
2255     if( rsz > 0 )
2256     {
2257         for( i = 0; i < rsz/2; i++ )
2258             _ranges[i] = (float*)&ranges[i*2];
2259     }
2260
2261     AutoBuffer<Mat> buf(nimages);
2262     for( i = 0; i < nimages; i++ )
2263         buf[i] = images.getMat(i);
2264
2265     calcBackProject(&buf[0], nimages, csz ? &channels[0] : 0,
2266         hist, dst, rsz ? (const float**)_ranges : 0, scale, true);
2267 }
2268
2269
2270 ////////////////// C O M P A R E   H I S T O G R A M S ////////////////////////
2271
2272 double cv::compareHist( InputArray _H1, InputArray _H2, int method )
2273 {
2274     Mat H1 = _H1.getMat(), H2 = _H2.getMat();
2275     const Mat* arrays[] = {&H1, &H2, 0};
2276     Mat planes[2];
2277     NAryMatIterator it(arrays, planes);
2278     double result = 0;
2279     int j, len = (int)it.size;
2280
2281     CV_Assert( H1.type() == H2.type() && H1.depth() == CV_32F );
2282
2283     double s1 = 0, s2 = 0, s11 = 0, s12 = 0, s22 = 0;
2284
2285     CV_Assert( it.planes[0].isContinuous() && it.planes[1].isContinuous() );
2286
2287     for( size_t i = 0; i < it.nplanes; i++, ++it )
2288     {
2289         const float* h1 = it.planes[0].ptr<float>();
2290         const float* h2 = it.planes[1].ptr<float>();
2291         len = it.planes[0].rows*it.planes[0].cols*H1.channels();
2292
2293         if( (method == CV_COMP_CHISQR) || (method == CV_COMP_CHISQR_ALT))
2294         {
2295             for( j = 0; j < len; j++ )
2296             {
2297                 double a = h1[j] - h2[j];
2298                 double b = (method == CV_COMP_CHISQR) ? h1[j] : h1[j] + h2[j];
2299                 if( fabs(b) > DBL_EPSILON )
2300                     result += a*a/b;
2301             }
2302         }
2303         else if( method == CV_COMP_CORREL )
2304         {
2305             for( j = 0; j < len; j++ )
2306             {
2307                 double a = h1[j];
2308                 double b = h2[j];
2309
2310                 s12 += a*b;
2311                 s1 += a;
2312                 s11 += a*a;
2313                 s2 += b;
2314                 s22 += b*b;
2315             }
2316         }
2317         else if( method == CV_COMP_INTERSECT )
2318         {
2319             j = 0;
2320             #if CV_NEON
2321             float32x4_t v_result = vdupq_n_f32(0.0f);
2322             for( ; j <= len - 4; j += 4 )
2323                 v_result = vaddq_f32(v_result, vminq_f32(vld1q_f32(h1 + j), vld1q_f32(h2 + j)));
2324             float CV_DECL_ALIGNED(16) ar[4];
2325             vst1q_f32(ar, v_result);
2326             result += ar[0] + ar[1] + ar[2] + ar[3];
2327             #endif
2328             for( ; j < len; j++ )
2329                 result += std::min(h1[j], h2[j]);
2330         }
2331         else if( method == CV_COMP_BHATTACHARYYA )
2332         {
2333             for( j = 0; j < len; j++ )
2334             {
2335                 double a = h1[j];
2336                 double b = h2[j];
2337                 result += std::sqrt(a*b);
2338                 s1 += a;
2339                 s2 += b;
2340             }
2341         }
2342         else if( method == CV_COMP_KL_DIV )
2343         {
2344             for( j = 0; j < len; j++ )
2345             {
2346                 double p = h1[j];
2347                 double q = h2[j];
2348                 if( fabs(p) <= DBL_EPSILON ) {
2349                     continue;
2350                 }
2351                 if(  fabs(q) <= DBL_EPSILON ) {
2352                     q = 1e-10;
2353                 }
2354                 result += p * std::log( p / q );
2355             }
2356         }
2357         else
2358             CV_Error( CV_StsBadArg, "Unknown comparison method" );
2359     }
2360
2361     if( method == CV_COMP_CHISQR_ALT )
2362         result *= 2;
2363     else if( method == CV_COMP_CORREL )
2364     {
2365         size_t total = H1.total();
2366         double scale = 1./total;
2367         double num = s12 - s1*s2*scale;
2368         double denom2 = (s11 - s1*s1*scale)*(s22 - s2*s2*scale);
2369         result = std::abs(denom2) > DBL_EPSILON ? num/std::sqrt(denom2) : 1.;
2370     }
2371     else if( method == CV_COMP_BHATTACHARYYA )
2372     {
2373         s1 *= s2;
2374         s1 = fabs(s1) > FLT_EPSILON ? 1./std::sqrt(s1) : 1.;
2375         result = std::sqrt(std::max(1. - result*s1, 0.));
2376     }
2377
2378     return result;
2379 }
2380
2381
2382 double cv::compareHist( const SparseMat& H1, const SparseMat& H2, int method )
2383 {
2384     double result = 0;
2385     int i, dims = H1.dims();
2386
2387     CV_Assert( dims > 0 && dims == H2.dims() && H1.type() == H2.type() && H1.type() == CV_32F );
2388     for( i = 0; i < dims; i++ )
2389         CV_Assert( H1.size(i) == H2.size(i) );
2390
2391     const SparseMat *PH1 = &H1, *PH2 = &H2;
2392     if( PH1->nzcount() > PH2->nzcount() && method != CV_COMP_CHISQR && method != CV_COMP_CHISQR_ALT && method != CV_COMP_KL_DIV )
2393         std::swap(PH1, PH2);
2394
2395     SparseMatConstIterator it = PH1->begin();
2396     int N1 = (int)PH1->nzcount(), N2 = (int)PH2->nzcount();
2397
2398     if( (method == CV_COMP_CHISQR) || (method == CV_COMP_CHISQR_ALT) )
2399     {
2400         for( i = 0; i < N1; i++, ++it )
2401         {
2402             float v1 = it.value<float>();
2403             const SparseMat::Node* node = it.node();
2404             float v2 = PH2->value<float>(node->idx, (size_t*)&node->hashval);
2405             double a = v1 - v2;
2406             double b = (method == CV_COMP_CHISQR) ? v1 : v1 + v2;
2407             if( fabs(b) > DBL_EPSILON )
2408                 result += a*a/b;
2409         }
2410     }
2411     else if( method == CV_COMP_CORREL )
2412     {
2413         double s1 = 0, s2 = 0, s11 = 0, s12 = 0, s22 = 0;
2414
2415         for( i = 0; i < N1; i++, ++it )
2416         {
2417             double v1 = it.value<float>();
2418             const SparseMat::Node* node = it.node();
2419             s12 += v1*PH2->value<float>(node->idx, (size_t*)&node->hashval);
2420             s1 += v1;
2421             s11 += v1*v1;
2422         }
2423
2424         it = PH2->begin();
2425         for( i = 0; i < N2; i++, ++it )
2426         {
2427             double v2 = it.value<float>();
2428             s2 += v2;
2429             s22 += v2*v2;
2430         }
2431
2432         size_t total = 1;
2433         for( i = 0; i < H1.dims(); i++ )
2434             total *= H1.size(i);
2435         double scale = 1./total;
2436         double num = s12 - s1*s2*scale;
2437         double denom2 = (s11 - s1*s1*scale)*(s22 - s2*s2*scale);
2438         result = std::abs(denom2) > DBL_EPSILON ? num/std::sqrt(denom2) : 1.;
2439     }
2440     else if( method == CV_COMP_INTERSECT )
2441     {
2442         for( i = 0; i < N1; i++, ++it )
2443         {
2444             float v1 = it.value<float>();
2445             const SparseMat::Node* node = it.node();
2446             float v2 = PH2->value<float>(node->idx, (size_t*)&node->hashval);
2447             if( v2 )
2448                 result += std::min(v1, v2);
2449         }
2450     }
2451     else if( method == CV_COMP_BHATTACHARYYA )
2452     {
2453         double s1 = 0, s2 = 0;
2454
2455         for( i = 0; i < N1; i++, ++it )
2456         {
2457             double v1 = it.value<float>();
2458             const SparseMat::Node* node = it.node();
2459             double v2 = PH2->value<float>(node->idx, (size_t*)&node->hashval);
2460             result += std::sqrt(v1*v2);
2461             s1 += v1;
2462         }
2463
2464         it = PH2->begin();
2465         for( i = 0; i < N2; i++, ++it )
2466             s2 += it.value<float>();
2467
2468         s1 *= s2;
2469         s1 = fabs(s1) > FLT_EPSILON ? 1./std::sqrt(s1) : 1.;
2470         result = std::sqrt(std::max(1. - result*s1, 0.));
2471     }
2472     else if( method == CV_COMP_KL_DIV )
2473     {
2474         for( i = 0; i < N1; i++, ++it )
2475         {
2476             double v1 = it.value<float>();
2477             const SparseMat::Node* node = it.node();
2478             double v2 = PH2->value<float>(node->idx, (size_t*)&node->hashval);
2479             if( !v2 )
2480                 v2 = 1e-10;
2481             result += v1 * std::log( v1 / v2 );
2482         }
2483     }
2484     else
2485         CV_Error( CV_StsBadArg, "Unknown comparison method" );
2486
2487     if( method == CV_COMP_CHISQR_ALT )
2488         result *= 2;
2489
2490     return result;
2491 }
2492
2493
2494 const int CV_HIST_DEFAULT_TYPE = CV_32F;
2495
2496 /* Creates new histogram */
2497 CvHistogram *
2498 cvCreateHist( int dims, int *sizes, CvHistType type, float** ranges, int uniform )
2499 {
2500     CvHistogram *hist = 0;
2501
2502     if( (unsigned)dims > CV_MAX_DIM )
2503         CV_Error( CV_BadOrder, "Number of dimensions is out of range" );
2504
2505     if( !sizes )
2506         CV_Error( CV_HeaderIsNull, "Null <sizes> pointer" );
2507
2508     hist = (CvHistogram *)cvAlloc( sizeof( CvHistogram ));
2509     hist->type = CV_HIST_MAGIC_VAL + ((int)type & 1);
2510     if (uniform) hist->type|= CV_HIST_UNIFORM_FLAG;
2511     hist->thresh2 = 0;
2512     hist->bins = 0;
2513     if( type == CV_HIST_ARRAY )
2514     {
2515         hist->bins = cvInitMatNDHeader( &hist->mat, dims, sizes,
2516                                         CV_HIST_DEFAULT_TYPE );
2517         cvCreateData( hist->bins );
2518     }
2519     else if( type == CV_HIST_SPARSE )
2520         hist->bins = cvCreateSparseMat( dims, sizes, CV_HIST_DEFAULT_TYPE );
2521     else
2522         CV_Error( CV_StsBadArg, "Invalid histogram type" );
2523
2524     if( ranges )
2525         cvSetHistBinRanges( hist, ranges, uniform );
2526
2527     return hist;
2528 }
2529
2530
2531 /* Creates histogram wrapping header for given array */
2532 CV_IMPL CvHistogram*
2533 cvMakeHistHeaderForArray( int dims, int *sizes, CvHistogram *hist,
2534                           float *data, float **ranges, int uniform )
2535 {
2536     if( !hist )
2537         CV_Error( CV_StsNullPtr, "Null histogram header pointer" );
2538
2539     if( !data )
2540         CV_Error( CV_StsNullPtr, "Null data pointer" );
2541
2542     hist->thresh2 = 0;
2543     hist->type = CV_HIST_MAGIC_VAL;
2544     hist->bins = cvInitMatNDHeader( &hist->mat, dims, sizes, CV_HIST_DEFAULT_TYPE, data );
2545
2546     if( ranges )
2547     {
2548         if( !uniform )
2549             CV_Error( CV_StsBadArg, "Only uniform bin ranges can be used here "
2550                                     "(to avoid memory allocation)" );
2551         cvSetHistBinRanges( hist, ranges, uniform );
2552     }
2553
2554     return hist;
2555 }
2556
2557
2558 CV_IMPL void
2559 cvReleaseHist( CvHistogram **hist )
2560 {
2561     if( !hist )
2562         CV_Error( CV_StsNullPtr, "" );
2563
2564     if( *hist )
2565     {
2566         CvHistogram* temp = *hist;
2567
2568         if( !CV_IS_HIST(temp))
2569             CV_Error( CV_StsBadArg, "Invalid histogram header" );
2570         *hist = 0;
2571
2572         if( CV_IS_SPARSE_HIST( temp ))
2573             cvReleaseSparseMat( (CvSparseMat**)&temp->bins );
2574         else
2575         {
2576             cvReleaseData( temp->bins );
2577             temp->bins = 0;
2578         }
2579
2580         if( temp->thresh2 )
2581             cvFree( &temp->thresh2 );
2582         cvFree( &temp );
2583     }
2584 }
2585
2586 CV_IMPL void
2587 cvClearHist( CvHistogram *hist )
2588 {
2589     if( !CV_IS_HIST(hist) )
2590         CV_Error( CV_StsBadArg, "Invalid histogram header" );
2591     cvZero( hist->bins );
2592 }
2593
2594
2595 // Clears histogram bins that are below than threshold
2596 CV_IMPL void
2597 cvThreshHist( CvHistogram* hist, double thresh )
2598 {
2599     if( !CV_IS_HIST(hist) )
2600         CV_Error( CV_StsBadArg, "Invalid histogram header" );
2601
2602     if( !CV_IS_SPARSE_MAT(hist->bins) )
2603     {
2604         CvMat mat;
2605         cvGetMat( hist->bins, &mat, 0, 1 );
2606         cvThreshold( &mat, &mat, thresh, 0, CV_THRESH_TOZERO );
2607     }
2608     else
2609     {
2610         CvSparseMat* mat = (CvSparseMat*)hist->bins;
2611         CvSparseMatIterator iterator;
2612         CvSparseNode *node;
2613
2614         for( node = cvInitSparseMatIterator( mat, &iterator );
2615              node != 0; node = cvGetNextSparseNode( &iterator ))
2616         {
2617             float* val = (float*)CV_NODE_VAL( mat, node );
2618             if( *val <= thresh )
2619                 *val = 0;
2620         }
2621     }
2622 }
2623
2624
2625 // Normalizes histogram (make sum of the histogram bins == factor)
2626 CV_IMPL void
2627 cvNormalizeHist( CvHistogram* hist, double factor )
2628 {
2629     double sum = 0;
2630
2631     if( !CV_IS_HIST(hist) )
2632         CV_Error( CV_StsBadArg, "Invalid histogram header" );
2633
2634     if( !CV_IS_SPARSE_HIST(hist) )
2635     {
2636         CvMat mat;
2637         cvGetMat( hist->bins, &mat, 0, 1 );
2638         sum = cvSum( &mat ).val[0];
2639         if( fabs(sum) < DBL_EPSILON )
2640             sum = 1;
2641         cvScale( &mat, &mat, factor/sum, 0 );
2642     }
2643     else
2644     {
2645         CvSparseMat* mat = (CvSparseMat*)hist->bins;
2646         CvSparseMatIterator iterator;
2647         CvSparseNode *node;
2648         float scale;
2649
2650         for( node = cvInitSparseMatIterator( mat, &iterator );
2651              node != 0; node = cvGetNextSparseNode( &iterator ))
2652         {
2653             sum += *(float*)CV_NODE_VAL(mat,node);
2654         }
2655
2656         if( fabs(sum) < DBL_EPSILON )
2657             sum = 1;
2658         scale = (float)(factor/sum);
2659
2660         for( node = cvInitSparseMatIterator( mat, &iterator );
2661              node != 0; node = cvGetNextSparseNode( &iterator ))
2662         {
2663             *(float*)CV_NODE_VAL(mat,node) *= scale;
2664         }
2665     }
2666 }
2667
2668
2669 // Retrieves histogram global min, max and their positions
2670 CV_IMPL void
2671 cvGetMinMaxHistValue( const CvHistogram* hist,
2672                       float *value_min, float* value_max,
2673                       int* idx_min, int* idx_max )
2674 {
2675     double minVal, maxVal;
2676     int dims, size[CV_MAX_DIM];
2677
2678     if( !CV_IS_HIST(hist) )
2679         CV_Error( CV_StsBadArg, "Invalid histogram header" );
2680
2681     dims = cvGetDims( hist->bins, size );
2682
2683     if( !CV_IS_SPARSE_HIST(hist) )
2684     {
2685         CvMat mat;
2686         CvPoint minPt, maxPt;
2687
2688         cvGetMat( hist->bins, &mat, 0, 1 );
2689         cvMinMaxLoc( &mat, &minVal, &maxVal, &minPt, &maxPt );
2690
2691         if( dims == 1 )
2692         {
2693             if( idx_min )
2694                 *idx_min = minPt.y + minPt.x;
2695             if( idx_max )
2696                 *idx_max = maxPt.y + maxPt.x;
2697         }
2698         else if( dims == 2 )
2699         {
2700             if( idx_min )
2701                 idx_min[0] = minPt.y, idx_min[1] = minPt.x;
2702             if( idx_max )
2703                 idx_max[0] = maxPt.y, idx_max[1] = maxPt.x;
2704         }
2705         else if( idx_min || idx_max )
2706         {
2707             int imin = minPt.y*mat.cols + minPt.x;
2708             int imax = maxPt.y*mat.cols + maxPt.x;
2709
2710             for(int i = dims - 1; i >= 0; i-- )
2711             {
2712                 if( idx_min )
2713                 {
2714                     int t = imin / size[i];
2715                     idx_min[i] = imin - t*size[i];
2716                     imin = t;
2717                 }
2718
2719                 if( idx_max )
2720                 {
2721                     int t = imax / size[i];
2722                     idx_max[i] = imax - t*size[i];
2723                     imax = t;
2724                 }
2725             }
2726         }
2727     }
2728     else
2729     {
2730         CvSparseMat* mat = (CvSparseMat*)hist->bins;
2731         CvSparseMatIterator iterator;
2732         CvSparseNode *node;
2733         int minv = INT_MAX;
2734         int maxv = INT_MIN;
2735         CvSparseNode* minNode = 0;
2736         CvSparseNode* maxNode = 0;
2737         const int *_idx_min = 0, *_idx_max = 0;
2738         Cv32suf m;
2739
2740         for( node = cvInitSparseMatIterator( mat, &iterator );
2741              node != 0; node = cvGetNextSparseNode( &iterator ))
2742         {
2743             int value = *(int*)CV_NODE_VAL(mat,node);
2744             value = CV_TOGGLE_FLT(value);
2745             if( value < minv )
2746             {
2747                 minv = value;
2748                 minNode = node;
2749             }
2750
2751             if( value > maxv )
2752             {
2753                 maxv = value;
2754                 maxNode = node;
2755             }
2756         }
2757
2758         if( minNode )
2759         {
2760             _idx_min = CV_NODE_IDX(mat,minNode);
2761             _idx_max = CV_NODE_IDX(mat,maxNode);
2762             m.i = CV_TOGGLE_FLT(minv); minVal = m.f;
2763             m.i = CV_TOGGLE_FLT(maxv); maxVal = m.f;
2764         }
2765         else
2766         {
2767             minVal = maxVal = 0;
2768         }
2769
2770         for(int i = 0; i < dims; i++ )
2771         {
2772             if( idx_min )
2773                 idx_min[i] = _idx_min ? _idx_min[i] : -1;
2774             if( idx_max )
2775                 idx_max[i] = _idx_max ? _idx_max[i] : -1;
2776         }
2777     }
2778
2779     if( value_min )
2780         *value_min = (float)minVal;
2781
2782     if( value_max )
2783         *value_max = (float)maxVal;
2784 }
2785
2786
2787 // Compares two histograms using one of a few methods
2788 CV_IMPL double
2789 cvCompareHist( const CvHistogram* hist1,
2790                const CvHistogram* hist2,
2791                int method )
2792 {
2793     int i;
2794     int size1[CV_MAX_DIM], size2[CV_MAX_DIM], total = 1;
2795
2796     if( !CV_IS_HIST(hist1) || !CV_IS_HIST(hist2) )
2797         CV_Error( CV_StsBadArg, "Invalid histogram header[s]" );
2798
2799     if( CV_IS_SPARSE_MAT(hist1->bins) != CV_IS_SPARSE_MAT(hist2->bins))
2800         CV_Error(CV_StsUnmatchedFormats, "One of histograms is sparse and other is not");
2801
2802     if( !CV_IS_SPARSE_MAT(hist1->bins) )
2803     {
2804         cv::Mat H1 = cv::cvarrToMat(hist1->bins);
2805         cv::Mat H2 = cv::cvarrToMat(hist2->bins);
2806         return cv::compareHist(H1, H2, method);
2807     }
2808
2809     int dims1 = cvGetDims( hist1->bins, size1 );
2810     int dims2 = cvGetDims( hist2->bins, size2 );
2811
2812     if( dims1 != dims2 )
2813         CV_Error( CV_StsUnmatchedSizes,
2814                  "The histograms have different numbers of dimensions" );
2815
2816     for( i = 0; i < dims1; i++ )
2817     {
2818         if( size1[i] != size2[i] )
2819             CV_Error( CV_StsUnmatchedSizes, "The histograms have different sizes" );
2820         total *= size1[i];
2821     }
2822
2823     double result = 0;
2824     CvSparseMat* mat1 = (CvSparseMat*)(hist1->bins);
2825     CvSparseMat* mat2 = (CvSparseMat*)(hist2->bins);
2826     CvSparseMatIterator iterator;
2827     CvSparseNode *node1, *node2;
2828
2829     if( mat1->heap->active_count > mat2->heap->active_count && method != CV_COMP_CHISQR && method != CV_COMP_CHISQR_ALT && method != CV_COMP_KL_DIV )
2830     {
2831         CvSparseMat* t;
2832         CV_SWAP( mat1, mat2, t );
2833     }
2834
2835     if( (method == CV_COMP_CHISQR) || (method == CV_COMP_CHISQR_ALT) )
2836     {
2837         for( node1 = cvInitSparseMatIterator( mat1, &iterator );
2838              node1 != 0; node1 = cvGetNextSparseNode( &iterator ))
2839         {
2840             double v1 = *(float*)CV_NODE_VAL(mat1,node1);
2841             uchar* node2_data = cvPtrND( mat2, CV_NODE_IDX(mat1,node1), 0, 0, &node1->hashval );
2842             double v2 = node2_data ? *(float*)node2_data : 0.f;
2843             double a = v1 - v2;
2844             double b = (method == CV_COMP_CHISQR) ? v1 : v1 + v2;
2845             if( fabs(b) > DBL_EPSILON )
2846                 result += a*a/b;
2847         }
2848     }
2849     else if( method == CV_COMP_CORREL )
2850     {
2851         double s1 = 0, s11 = 0;
2852         double s2 = 0, s22 = 0;
2853         double s12 = 0;
2854         double num, denom2, scale = 1./total;
2855
2856         for( node1 = cvInitSparseMatIterator( mat1, &iterator );
2857              node1 != 0; node1 = cvGetNextSparseNode( &iterator ))
2858         {
2859             double v1 = *(float*)CV_NODE_VAL(mat1,node1);
2860             uchar* node2_data = cvPtrND( mat2, CV_NODE_IDX(mat1,node1),
2861                                         0, 0, &node1->hashval );
2862             if( node2_data )
2863             {
2864                 double v2 = *(float*)node2_data;
2865                 s12 += v1*v2;
2866             }
2867             s1 += v1;
2868             s11 += v1*v1;
2869         }
2870
2871         for( node2 = cvInitSparseMatIterator( mat2, &iterator );
2872              node2 != 0; node2 = cvGetNextSparseNode( &iterator ))
2873         {
2874             double v2 = *(float*)CV_NODE_VAL(mat2,node2);
2875             s2 += v2;
2876             s22 += v2*v2;
2877         }
2878
2879         num = s12 - s1*s2*scale;
2880         denom2 = (s11 - s1*s1*scale)*(s22 - s2*s2*scale);
2881         result = fabs(denom2) > DBL_EPSILON ? num/sqrt(denom2) : 1;
2882     }
2883     else if( method == CV_COMP_INTERSECT )
2884     {
2885         for( node1 = cvInitSparseMatIterator( mat1, &iterator );
2886              node1 != 0; node1 = cvGetNextSparseNode( &iterator ))
2887         {
2888             float v1 = *(float*)CV_NODE_VAL(mat1,node1);
2889             uchar* node2_data = cvPtrND( mat2, CV_NODE_IDX(mat1,node1),
2890                                          0, 0, &node1->hashval );
2891             if( node2_data )
2892             {
2893                 float v2 = *(float*)node2_data;
2894                 if( v1 <= v2 )
2895                     result += v1;
2896                 else
2897                     result += v2;
2898             }
2899         }
2900     }
2901     else if( method == CV_COMP_BHATTACHARYYA )
2902     {
2903         double s1 = 0, s2 = 0;
2904
2905         for( node1 = cvInitSparseMatIterator( mat1, &iterator );
2906              node1 != 0; node1 = cvGetNextSparseNode( &iterator ))
2907         {
2908             double v1 = *(float*)CV_NODE_VAL(mat1,node1);
2909             uchar* node2_data = cvPtrND( mat2, CV_NODE_IDX(mat1,node1),
2910                                          0, 0, &node1->hashval );
2911             s1 += v1;
2912             if( node2_data )
2913             {
2914                 double v2 = *(float*)node2_data;
2915                 result += sqrt(v1 * v2);
2916             }
2917         }
2918
2919         for( node1 = cvInitSparseMatIterator( mat2, &iterator );
2920              node1 != 0; node1 = cvGetNextSparseNode( &iterator ))
2921         {
2922             double v2 = *(float*)CV_NODE_VAL(mat2,node1);
2923             s2 += v2;
2924         }
2925
2926         s1 *= s2;
2927         s1 = fabs(s1) > FLT_EPSILON ? 1./sqrt(s1) : 1.;
2928         result = 1. - result*s1;
2929         result = sqrt(MAX(result,0.));
2930     }
2931     else if( method == CV_COMP_KL_DIV )
2932     {
2933         cv::SparseMat sH1, sH2;
2934         ((const CvSparseMat*)hist1->bins)->copyToSparseMat(sH1);
2935         ((const CvSparseMat*)hist2->bins)->copyToSparseMat(sH2);
2936         result = cv::compareHist( sH1, sH2, CV_COMP_KL_DIV );
2937     }
2938     else
2939         CV_Error( CV_StsBadArg, "Unknown comparison method" );
2940
2941     if( method == CV_COMP_CHISQR_ALT )
2942         result *= 2;
2943
2944     return result;
2945 }
2946
2947 // copies one histogram to another
2948 CV_IMPL void
2949 cvCopyHist( const CvHistogram* src, CvHistogram** _dst )
2950 {
2951     if( !_dst )
2952         CV_Error( CV_StsNullPtr, "Destination double pointer is NULL" );
2953
2954     CvHistogram* dst = *_dst;
2955
2956     if( !CV_IS_HIST(src) || (dst && !CV_IS_HIST(dst)) )
2957         CV_Error( CV_StsBadArg, "Invalid histogram header[s]" );
2958
2959     bool eq = false;
2960     int size1[CV_MAX_DIM];
2961     bool is_sparse = CV_IS_SPARSE_MAT(src->bins);
2962     int dims1 = cvGetDims( src->bins, size1 );
2963
2964     if( dst && (is_sparse == CV_IS_SPARSE_MAT(dst->bins)))
2965     {
2966         int size2[CV_MAX_DIM];
2967         int dims2 = cvGetDims( dst->bins, size2 );
2968
2969         if( dims1 == dims2 )
2970         {
2971             int i;
2972
2973             for( i = 0; i < dims1; i++ )
2974             {
2975                 if( size1[i] != size2[i] )
2976                     break;
2977             }
2978
2979             eq = (i == dims1);
2980         }
2981     }
2982
2983     if( !eq )
2984     {
2985         cvReleaseHist( _dst );
2986         dst = cvCreateHist( dims1, size1, !is_sparse ? CV_HIST_ARRAY : CV_HIST_SPARSE, 0, 0 );
2987         *_dst = dst;
2988     }
2989
2990     if( CV_HIST_HAS_RANGES( src ))
2991     {
2992         float* ranges[CV_MAX_DIM];
2993         float** thresh = 0;
2994
2995         if( CV_IS_UNIFORM_HIST( src ))
2996         {
2997             for( int i = 0; i < dims1; i++ )
2998                 ranges[i] = (float*)src->thresh[i];
2999
3000             thresh = ranges;
3001         }
3002         else
3003         {
3004             thresh = src->thresh2;
3005         }
3006
3007         cvSetHistBinRanges( dst, thresh, CV_IS_UNIFORM_HIST(src));
3008     }
3009
3010     cvCopy( src->bins, dst->bins );
3011 }
3012
3013
3014 // Sets a value range for every histogram bin
3015 CV_IMPL void
3016 cvSetHistBinRanges( CvHistogram* hist, float** ranges, int uniform )
3017 {
3018     int dims, size[CV_MAX_DIM], total = 0;
3019     int i, j;
3020
3021     if( !ranges )
3022         CV_Error( CV_StsNullPtr, "NULL ranges pointer" );
3023
3024     if( !CV_IS_HIST(hist) )
3025         CV_Error( CV_StsBadArg, "Invalid histogram header" );
3026
3027     dims = cvGetDims( hist->bins, size );
3028     for( i = 0; i < dims; i++ )
3029         total += size[i]+1;
3030
3031     if( uniform )
3032     {
3033         for( i = 0; i < dims; i++ )
3034         {
3035             if( !ranges[i] )
3036                 CV_Error( CV_StsNullPtr, "One of <ranges> elements is NULL" );
3037             hist->thresh[i][0] = ranges[i][0];
3038             hist->thresh[i][1] = ranges[i][1];
3039         }
3040
3041         hist->type |= CV_HIST_UNIFORM_FLAG + CV_HIST_RANGES_FLAG;
3042     }
3043     else
3044     {
3045         float* dim_ranges;
3046
3047         if( !hist->thresh2 )
3048         {
3049             hist->thresh2 = (float**)cvAlloc(
3050                         dims*sizeof(hist->thresh2[0])+
3051                         total*sizeof(hist->thresh2[0][0]));
3052         }
3053         dim_ranges = (float*)(hist->thresh2 + dims);
3054
3055         for( i = 0; i < dims; i++ )
3056         {
3057             float val0 = -FLT_MAX;
3058
3059             if( !ranges[i] )
3060                 CV_Error( CV_StsNullPtr, "One of <ranges> elements is NULL" );
3061
3062             for( j = 0; j <= size[i]; j++ )
3063             {
3064                 float val = ranges[i][j];
3065                 if( val <= val0 )
3066                     CV_Error(CV_StsOutOfRange, "Bin ranges should go in ascenting order");
3067                 val0 = dim_ranges[j] = val;
3068             }
3069
3070             hist->thresh2[i] = dim_ranges;
3071             dim_ranges += size[i] + 1;
3072         }
3073
3074         hist->type |= CV_HIST_RANGES_FLAG;
3075         hist->type &= ~CV_HIST_UNIFORM_FLAG;
3076     }
3077 }
3078
3079
3080 CV_IMPL void
3081 cvCalcArrHist( CvArr** img, CvHistogram* hist, int accumulate, const CvArr* mask )
3082 {
3083     if( !CV_IS_HIST(hist))
3084         CV_Error( CV_StsBadArg, "Bad histogram pointer" );
3085
3086     if( !img )
3087         CV_Error( CV_StsNullPtr, "Null double array pointer" );
3088
3089     int size[CV_MAX_DIM];
3090     int i, dims = cvGetDims( hist->bins, size);
3091     bool uniform = CV_IS_UNIFORM_HIST(hist);
3092
3093     std::vector<cv::Mat> images(dims);
3094     for( i = 0; i < dims; i++ )
3095         images[i] = cv::cvarrToMat(img[i]);
3096
3097     cv::Mat _mask;
3098     if( mask )
3099         _mask = cv::cvarrToMat(mask);
3100
3101     const float* uranges[CV_MAX_DIM] = {0};
3102     const float** ranges = 0;
3103
3104     if( hist->type & CV_HIST_RANGES_FLAG )
3105     {
3106         ranges = (const float**)hist->thresh2;
3107         if( uniform )
3108         {
3109             for( i = 0; i < dims; i++ )
3110                 uranges[i] = &hist->thresh[i][0];
3111             ranges = uranges;
3112         }
3113     }
3114
3115     if( !CV_IS_SPARSE_HIST(hist) )
3116     {
3117         cv::Mat H = cv::cvarrToMat(hist->bins);
3118         cv::calcHist( &images[0], (int)images.size(), 0, _mask,
3119                       H, cvGetDims(hist->bins), H.size, ranges, uniform, accumulate != 0 );
3120     }
3121     else
3122     {
3123         CvSparseMat* sparsemat = (CvSparseMat*)hist->bins;
3124
3125         if( !accumulate )
3126             cvZero( hist->bins );
3127         cv::SparseMat sH;
3128         sparsemat->copyToSparseMat(sH);
3129         cv::calcHist( &images[0], (int)images.size(), 0, _mask, sH, sH.dims(),
3130                       sH.dims() > 0 ? sH.hdr->size : 0, ranges, uniform, accumulate != 0, true );
3131
3132         if( accumulate )
3133             cvZero( sparsemat );
3134
3135         cv::SparseMatConstIterator it = sH.begin();
3136         int nz = (int)sH.nzcount();
3137         for( i = 0; i < nz; i++, ++it )
3138             *(float*)cvPtrND(sparsemat, it.node()->idx, 0, -2) = (float)*(const int*)it.ptr;
3139     }
3140 }
3141
3142
3143 CV_IMPL void
3144 cvCalcArrBackProject( CvArr** img, CvArr* dst, const CvHistogram* hist )
3145 {
3146     if( !CV_IS_HIST(hist))
3147         CV_Error( CV_StsBadArg, "Bad histogram pointer" );
3148
3149     if( !img )
3150         CV_Error( CV_StsNullPtr, "Null double array pointer" );
3151
3152     int size[CV_MAX_DIM];
3153     int i, dims = cvGetDims( hist->bins, size );
3154
3155     bool uniform = CV_IS_UNIFORM_HIST(hist);
3156     const float* uranges[CV_MAX_DIM] = {0};
3157     const float** ranges = 0;
3158
3159     if( hist->type & CV_HIST_RANGES_FLAG )
3160     {
3161         ranges = (const float**)hist->thresh2;
3162         if( uniform )
3163         {
3164             for( i = 0; i < dims; i++ )
3165                 uranges[i] = &hist->thresh[i][0];
3166             ranges = uranges;
3167         }
3168     }
3169
3170     std::vector<cv::Mat> images(dims);
3171     for( i = 0; i < dims; i++ )
3172         images[i] = cv::cvarrToMat(img[i]);
3173
3174     cv::Mat _dst = cv::cvarrToMat(dst);
3175
3176     CV_Assert( _dst.size() == images[0].size() && _dst.depth() == images[0].depth() );
3177
3178     if( !CV_IS_SPARSE_HIST(hist) )
3179     {
3180         cv::Mat H = cv::cvarrToMat(hist->bins);
3181         cv::calcBackProject( &images[0], (int)images.size(),
3182                             0, H, _dst, ranges, 1, uniform );
3183     }
3184     else
3185     {
3186         cv::SparseMat sH;
3187         ((const CvSparseMat*)hist->bins)->copyToSparseMat(sH);
3188         cv::calcBackProject( &images[0], (int)images.size(),
3189                              0, sH, _dst, ranges, 1, uniform );
3190     }
3191 }
3192
3193
3194 ////////////////////// B A C K   P R O J E C T   P A T C H /////////////////////////
3195
3196 CV_IMPL void
3197 cvCalcArrBackProjectPatch( CvArr** arr, CvArr* dst, CvSize patch_size, CvHistogram* hist,
3198                            int method, double norm_factor )
3199 {
3200     CvHistogram* model = 0;
3201
3202     IplImage imgstub[CV_MAX_DIM], *img[CV_MAX_DIM];
3203     IplROI roi;
3204     CvMat dststub, *dstmat;
3205     int i, dims;
3206     int x, y;
3207     CvSize size;
3208
3209     if( !CV_IS_HIST(hist))
3210         CV_Error( CV_StsBadArg, "Bad histogram pointer" );
3211
3212     if( !arr )
3213         CV_Error( CV_StsNullPtr, "Null double array pointer" );
3214
3215     if( norm_factor <= 0 )
3216         CV_Error( CV_StsOutOfRange,
3217                   "Bad normalization factor (set it to 1.0 if unsure)" );
3218
3219     if( patch_size.width <= 0 || patch_size.height <= 0 )
3220         CV_Error( CV_StsBadSize, "The patch width and height must be positive" );
3221
3222     dims = cvGetDims( hist->bins );
3223     cvNormalizeHist( hist, norm_factor );
3224
3225     for( i = 0; i < dims; i++ )
3226     {
3227         CvMat stub, *mat;
3228         mat = cvGetMat( arr[i], &stub, 0, 0 );
3229         img[i] = cvGetImage( mat, &imgstub[i] );
3230         img[i]->roi = &roi;
3231     }
3232
3233     dstmat = cvGetMat( dst, &dststub, 0, 0 );
3234     if( CV_MAT_TYPE( dstmat->type ) != CV_32FC1 )
3235         CV_Error( CV_StsUnsupportedFormat, "Resultant image must have 32fC1 type" );
3236
3237     if( dstmat->cols != img[0]->width - patch_size.width + 1 ||
3238         dstmat->rows != img[0]->height - patch_size.height + 1 )
3239         CV_Error( CV_StsUnmatchedSizes,
3240             "The output map must be (W-w+1 x H-h+1), "
3241             "where the input images are (W x H) each and the patch is (w x h)" );
3242
3243     cvCopyHist( hist, &model );
3244
3245     size = cvGetMatSize(dstmat);
3246     roi.coi = 0;
3247     roi.width = patch_size.width;
3248     roi.height = patch_size.height;
3249
3250     for( y = 0; y < size.height; y++ )
3251     {
3252         for( x = 0; x < size.width; x++ )
3253         {
3254             double result;
3255             roi.xOffset = x;
3256             roi.yOffset = y;
3257
3258             cvCalcHist( img, model );
3259             cvNormalizeHist( model, norm_factor );
3260             result = cvCompareHist( model, hist, method );
3261             CV_MAT_ELEM( *dstmat, float, y, x ) = (float)result;
3262         }
3263     }
3264
3265     cvReleaseHist( &model );
3266 }
3267
3268
3269 // Calculates Bayes probabilistic histograms
3270 CV_IMPL void
3271 cvCalcBayesianProb( CvHistogram** src, int count, CvHistogram** dst )
3272 {
3273     int i;
3274
3275     if( !src || !dst )
3276         CV_Error( CV_StsNullPtr, "NULL histogram array pointer" );
3277
3278     if( count < 2 )
3279         CV_Error( CV_StsOutOfRange, "Too small number of histograms" );
3280
3281     for( i = 0; i < count; i++ )
3282     {
3283         if( !CV_IS_HIST(src[i]) || !CV_IS_HIST(dst[i]) )
3284             CV_Error( CV_StsBadArg, "Invalid histogram header" );
3285
3286         if( !CV_IS_MATND(src[i]->bins) || !CV_IS_MATND(dst[i]->bins) )
3287             CV_Error( CV_StsBadArg, "The function supports dense histograms only" );
3288     }
3289
3290     cvZero( dst[0]->bins );
3291     // dst[0] = src[0] + ... + src[count-1]
3292     for( i = 0; i < count; i++ )
3293         cvAdd( src[i]->bins, dst[0]->bins, dst[0]->bins );
3294
3295     cvDiv( 0, dst[0]->bins, dst[0]->bins );
3296
3297     // dst[i] = src[i]*(1/dst[0])
3298     for( i = count - 1; i >= 0; i-- )
3299         cvMul( src[i]->bins, dst[0]->bins, dst[i]->bins );
3300 }
3301
3302
3303 CV_IMPL void
3304 cvCalcProbDensity( const CvHistogram* hist, const CvHistogram* hist_mask,
3305                    CvHistogram* hist_dens, double scale )
3306 {
3307     if( scale <= 0 )
3308         CV_Error( CV_StsOutOfRange, "scale must be positive" );
3309
3310     if( !CV_IS_HIST(hist) || !CV_IS_HIST(hist_mask) || !CV_IS_HIST(hist_dens) )
3311         CV_Error( CV_StsBadArg, "Invalid histogram pointer[s]" );
3312
3313     {
3314         CvArr* arrs[] = { hist->bins, hist_mask->bins, hist_dens->bins };
3315         CvMatND stubs[3];
3316         CvNArrayIterator iterator;
3317
3318         cvInitNArrayIterator( 3, arrs, 0, stubs, &iterator );
3319
3320         if( CV_MAT_TYPE(iterator.hdr[0]->type) != CV_32FC1 )
3321             CV_Error( CV_StsUnsupportedFormat, "All histograms must have 32fC1 type" );
3322
3323         do
3324         {
3325             const float* srcdata = (const float*)(iterator.ptr[0]);
3326             const float* maskdata = (const float*)(iterator.ptr[1]);
3327             float* dstdata = (float*)(iterator.ptr[2]);
3328             int i;
3329
3330             for( i = 0; i < iterator.size.width; i++ )
3331             {
3332                 float s = srcdata[i];
3333                 float m = maskdata[i];
3334                 if( s > FLT_EPSILON )
3335                     if( m <= s )
3336                         dstdata[i] = (float)(m*scale/s);
3337                     else
3338                         dstdata[i] = (float)scale;
3339                 else
3340                     dstdata[i] = (float)0;
3341             }
3342         }
3343         while( cvNextNArraySlice( &iterator ));
3344     }
3345 }
3346
3347 class EqualizeHistCalcHist_Invoker : public cv::ParallelLoopBody
3348 {
3349 public:
3350     enum {HIST_SZ = 256};
3351
3352     EqualizeHistCalcHist_Invoker(cv::Mat& src, int* histogram, cv::Mutex* histogramLock)
3353         : src_(src), globalHistogram_(histogram), histogramLock_(histogramLock)
3354     { }
3355
3356     void operator()( const cv::Range& rowRange ) const
3357     {
3358         int localHistogram[HIST_SZ] = {0, };
3359
3360         const size_t sstep = src_.step;
3361
3362         int width = src_.cols;
3363         int height = rowRange.end - rowRange.start;
3364
3365         if (src_.isContinuous())
3366         {
3367             width *= height;
3368             height = 1;
3369         }
3370
3371         for (const uchar* ptr = src_.ptr<uchar>(rowRange.start); height--; ptr += sstep)
3372         {
3373             int x = 0;
3374             for (; x <= width - 4; x += 4)
3375             {
3376                 int t0 = ptr[x], t1 = ptr[x+1];
3377                 localHistogram[t0]++; localHistogram[t1]++;
3378                 t0 = ptr[x+2]; t1 = ptr[x+3];
3379                 localHistogram[t0]++; localHistogram[t1]++;
3380             }
3381
3382             for (; x < width; ++x)
3383                 localHistogram[ptr[x]]++;
3384         }
3385
3386         cv::AutoLock lock(*histogramLock_);
3387
3388         for( int i = 0; i < HIST_SZ; i++ )
3389             globalHistogram_[i] += localHistogram[i];
3390     }
3391
3392     static bool isWorthParallel( const cv::Mat& src )
3393     {
3394         return ( src.total() >= 640*480 );
3395     }
3396
3397 private:
3398     EqualizeHistCalcHist_Invoker& operator=(const EqualizeHistCalcHist_Invoker&);
3399
3400     cv::Mat& src_;
3401     int* globalHistogram_;
3402     cv::Mutex* histogramLock_;
3403 };
3404
3405 class EqualizeHistLut_Invoker : public cv::ParallelLoopBody
3406 {
3407 public:
3408     EqualizeHistLut_Invoker( cv::Mat& src, cv::Mat& dst, int* lut )
3409         : src_(src),
3410           dst_(dst),
3411           lut_(lut)
3412     { }
3413
3414     void operator()( const cv::Range& rowRange ) const
3415     {
3416         const size_t sstep = src_.step;
3417         const size_t dstep = dst_.step;
3418
3419         int width = src_.cols;
3420         int height = rowRange.end - rowRange.start;
3421         int* lut = lut_;
3422
3423         if (src_.isContinuous() && dst_.isContinuous())
3424         {
3425             width *= height;
3426             height = 1;
3427         }
3428
3429         const uchar* sptr = src_.ptr<uchar>(rowRange.start);
3430         uchar* dptr = dst_.ptr<uchar>(rowRange.start);
3431
3432         for (; height--; sptr += sstep, dptr += dstep)
3433         {
3434             int x = 0;
3435             for (; x <= width - 4; x += 4)
3436             {
3437                 int v0 = sptr[x];
3438                 int v1 = sptr[x+1];
3439                 int x0 = lut[v0];
3440                 int x1 = lut[v1];
3441                 dptr[x] = (uchar)x0;
3442                 dptr[x+1] = (uchar)x1;
3443
3444                 v0 = sptr[x+2];
3445                 v1 = sptr[x+3];
3446                 x0 = lut[v0];
3447                 x1 = lut[v1];
3448                 dptr[x+2] = (uchar)x0;
3449                 dptr[x+3] = (uchar)x1;
3450             }
3451
3452             for (; x < width; ++x)
3453                 dptr[x] = (uchar)lut[sptr[x]];
3454         }
3455     }
3456
3457     static bool isWorthParallel( const cv::Mat& src )
3458     {
3459         return ( src.total() >= 640*480 );
3460     }
3461
3462 private:
3463     EqualizeHistLut_Invoker& operator=(const EqualizeHistLut_Invoker&);
3464
3465     cv::Mat& src_;
3466     cv::Mat& dst_;
3467     int* lut_;
3468 };
3469
3470 CV_IMPL void cvEqualizeHist( const CvArr* srcarr, CvArr* dstarr )
3471 {
3472     cv::equalizeHist(cv::cvarrToMat(srcarr), cv::cvarrToMat(dstarr));
3473 }
3474
3475 #ifdef HAVE_OPENCL
3476
3477 namespace cv {
3478
3479 static bool ocl_equalizeHist(InputArray _src, OutputArray _dst)
3480 {
3481     const ocl::Device & dev = ocl::Device::getDefault();
3482     int compunits = dev.maxComputeUnits();
3483     size_t wgs = dev.maxWorkGroupSize();
3484     Size size = _src.size();
3485     bool use16 = size.width % 16 == 0 && _src.offset() % 16 == 0 && _src.step() % 16 == 0;
3486     int kercn = dev.isAMD() && use16 ? 16 : std::min(4, ocl::predictOptimalVectorWidth(_src));
3487
3488     ocl::Kernel k1("calculate_histogram", ocl::imgproc::histogram_oclsrc,
3489                    format("-D BINS=%d -D HISTS_COUNT=%d -D WGS=%d -D kercn=%d -D T=%s%s",
3490                           BINS, compunits, wgs, kercn,
3491                           kercn == 4 ? "int" : ocl::typeToStr(CV_8UC(kercn)),
3492                           _src.isContinuous() ? " -D HAVE_SRC_CONT" : ""));
3493     if (k1.empty())
3494         return false;
3495
3496     UMat src = _src.getUMat(), ghist(1, BINS * compunits, CV_32SC1);
3497
3498     k1.args(ocl::KernelArg::ReadOnly(src),
3499             ocl::KernelArg::PtrWriteOnly(ghist), (int)src.total());
3500
3501     size_t globalsize = compunits * wgs;
3502     if (!k1.run(1, &globalsize, &wgs, false))
3503         return false;
3504
3505     wgs = std::min<size_t>(ocl::Device::getDefault().maxWorkGroupSize(), BINS);
3506     UMat lut(1, 256, CV_8UC1);
3507     ocl::Kernel k2("calcLUT", ocl::imgproc::histogram_oclsrc,
3508                   format("-D BINS=%d -D HISTS_COUNT=%d -D WGS=%d",
3509                          BINS, compunits, (int)wgs));
3510     k2.args(ocl::KernelArg::PtrWriteOnly(lut),
3511            ocl::KernelArg::PtrReadOnly(ghist), (int)_src.total());
3512
3513     // calculation of LUT
3514     if (!k2.run(1, &wgs, &wgs, false))
3515         return false;
3516
3517     // execute LUT transparently
3518     LUT(_src, lut, _dst);
3519     return true;
3520 }
3521
3522 }
3523
3524 #endif
3525
3526 void cv::equalizeHist( InputArray _src, OutputArray _dst )
3527 {
3528     CV_Assert( _src.type() == CV_8UC1 );
3529
3530     if (_src.empty())
3531         return;
3532
3533     CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
3534                ocl_equalizeHist(_src, _dst))
3535
3536     Mat src = _src.getMat();
3537     _dst.create( src.size(), src.type() );
3538     Mat dst = _dst.getMat();
3539
3540     Mutex histogramLockInstance;
3541
3542     const int hist_sz = EqualizeHistCalcHist_Invoker::HIST_SZ;
3543     int hist[hist_sz] = {0,};
3544     int lut[hist_sz];
3545
3546     EqualizeHistCalcHist_Invoker calcBody(src, hist, &histogramLockInstance);
3547     EqualizeHistLut_Invoker      lutBody(src, dst, lut);
3548     cv::Range heightRange(0, src.rows);
3549
3550     if(EqualizeHistCalcHist_Invoker::isWorthParallel(src))
3551         parallel_for_(heightRange, calcBody);
3552     else
3553         calcBody(heightRange);
3554
3555     int i = 0;
3556     while (!hist[i]) ++i;
3557
3558     int total = (int)src.total();
3559     if (hist[i] == total)
3560     {
3561         dst.setTo(i);
3562         return;
3563     }
3564
3565     float scale = (hist_sz - 1.f)/(total - hist[i]);
3566     int sum = 0;
3567
3568     for (lut[i++] = 0; i < hist_sz; ++i)
3569     {
3570         sum += hist[i];
3571         lut[i] = saturate_cast<uchar>(sum * scale);
3572     }
3573
3574     if(EqualizeHistLut_Invoker::isWorthParallel(src))
3575         parallel_for_(heightRange, lutBody);
3576     else
3577         lutBody(heightRange);
3578 }
3579
3580 // ----------------------------------------------------------------------
3581
3582 /* Implementation of RTTI and Generic Functions for CvHistogram */
3583 #define CV_TYPE_NAME_HIST "opencv-hist"
3584
3585 static int icvIsHist( const void * ptr )
3586 {
3587     return CV_IS_HIST( ((CvHistogram*)ptr) );
3588 }
3589
3590 static CvHistogram * icvCloneHist( const CvHistogram * src )
3591 {
3592     CvHistogram * dst=NULL;
3593     cvCopyHist(src, &dst);
3594     return dst;
3595 }
3596
3597 static void *icvReadHist( CvFileStorage * fs, CvFileNode * node )
3598 {
3599     CvHistogram * h = 0;
3600     int type = 0;
3601     int is_uniform = 0;
3602     int have_ranges = 0;
3603
3604     h = (CvHistogram *)cvAlloc( sizeof(CvHistogram) );
3605
3606     type = cvReadIntByName( fs, node, "type", 0 );
3607     is_uniform = cvReadIntByName( fs, node, "is_uniform", 0 );
3608     have_ranges = cvReadIntByName( fs, node, "have_ranges", 0 );
3609     h->type = CV_HIST_MAGIC_VAL | type |
3610         (is_uniform ? CV_HIST_UNIFORM_FLAG : 0) |
3611         (have_ranges ? CV_HIST_RANGES_FLAG : 0);
3612
3613     if(type == CV_HIST_ARRAY)
3614     {
3615         // read histogram bins
3616         CvMatND* mat = (CvMatND*)cvReadByName( fs, node, "mat" );
3617         int i, sizes[CV_MAX_DIM];
3618
3619         if(!CV_IS_MATND(mat))
3620             CV_Error( CV_StsError, "Expected CvMatND");
3621
3622         for(i=0; i<mat->dims; i++)
3623             sizes[i] = mat->dim[i].size;
3624
3625         cvInitMatNDHeader( &(h->mat), mat->dims, sizes, mat->type, mat->data.ptr );
3626         h->bins = &(h->mat);
3627
3628         // take ownership of refcount pointer as well
3629         h->mat.refcount = mat->refcount;
3630
3631         // increase refcount so freeing temp header doesn't free data
3632         cvIncRefData( mat );
3633
3634         // free temporary header
3635         cvReleaseMatND( &mat );
3636     }
3637     else
3638     {
3639         h->bins = cvReadByName( fs, node, "bins" );
3640         if(!CV_IS_SPARSE_MAT(h->bins)){
3641             CV_Error( CV_StsError, "Unknown Histogram type");
3642         }
3643     }
3644
3645     // read thresholds
3646     if(have_ranges)
3647     {
3648         int i, dims, size[CV_MAX_DIM], total = 0;
3649         CvSeqReader reader;
3650         CvFileNode * thresh_node;
3651
3652         dims = cvGetDims( h->bins, size );
3653         for( i = 0; i < dims; i++ )
3654             total += size[i]+1;
3655
3656         thresh_node = cvGetFileNodeByName( fs, node, "thresh" );
3657         if(!thresh_node)
3658             CV_Error( CV_StsError, "'thresh' node is missing");
3659         cvStartReadRawData( fs, thresh_node, &reader );
3660
3661         if(is_uniform)
3662         {
3663             for(i=0; i<dims; i++)
3664                 cvReadRawDataSlice( fs, &reader, 2, h->thresh[i], "f" );
3665             h->thresh2 = NULL;
3666         }
3667         else
3668         {
3669             float* dim_ranges;
3670             h->thresh2 = (float**)cvAlloc(
3671                 dims*sizeof(h->thresh2[0])+
3672                 total*sizeof(h->thresh2[0][0]));
3673             dim_ranges = (float*)(h->thresh2 + dims);
3674             for(i=0; i < dims; i++)
3675             {
3676                 h->thresh2[i] = dim_ranges;
3677                 cvReadRawDataSlice( fs, &reader, size[i]+1, dim_ranges, "f" );
3678                 dim_ranges += size[i] + 1;
3679             }
3680         }
3681     }
3682
3683     return h;
3684 }
3685
3686 static void icvWriteHist( CvFileStorage* fs, const char* name,
3687                           const void* struct_ptr, CvAttrList /*attributes*/ )
3688 {
3689     const CvHistogram * hist = (const CvHistogram *) struct_ptr;
3690     int sizes[CV_MAX_DIM];
3691     int dims;
3692     int i;
3693     int is_uniform, have_ranges;
3694
3695     cvStartWriteStruct( fs, name, CV_NODE_MAP, CV_TYPE_NAME_HIST );
3696
3697     is_uniform = (CV_IS_UNIFORM_HIST(hist) ? 1 : 0);
3698     have_ranges = (hist->type & CV_HIST_RANGES_FLAG ? 1 : 0);
3699
3700     cvWriteInt( fs, "type", (hist->type & 1) );
3701     cvWriteInt( fs, "is_uniform", is_uniform );
3702     cvWriteInt( fs, "have_ranges", have_ranges );
3703     if(!CV_IS_SPARSE_HIST(hist))
3704         cvWrite( fs, "mat", &(hist->mat) );
3705     else
3706         cvWrite( fs, "bins", hist->bins );
3707
3708     // write thresholds
3709     if(have_ranges){
3710         dims = cvGetDims( hist->bins, sizes );
3711         cvStartWriteStruct( fs, "thresh", CV_NODE_SEQ + CV_NODE_FLOW );
3712         if(is_uniform){
3713             for(i=0; i<dims; i++){
3714                 cvWriteRawData( fs, hist->thresh[i], 2, "f" );
3715             }
3716         }
3717         else{
3718             for(i=0; i<dims; i++){
3719                 cvWriteRawData( fs, hist->thresh2[i], sizes[i]+1, "f" );
3720             }
3721         }
3722         cvEndWriteStruct( fs );
3723     }
3724
3725     cvEndWriteStruct( fs );
3726 }
3727
3728
3729 CvType hist_type( CV_TYPE_NAME_HIST, icvIsHist, (CvReleaseFunc)cvReleaseHist,
3730                   icvReadHist, icvWriteHist, (CvCloneFunc)icvCloneHist );
3731
3732 /* End of file. */