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