vcs.maemo.org Git - opencv/blob - src/cv/cvhistogram.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 "_cv.h"
  42
  43 namespace cv
  44 {
  45
  46 ////////////////// Helper functions //////////////////////
  47
  48 static const size_t OUT_OF_RANGE = (size_t)1 << (sizeof(size_t)*8 - 2);
  49
  50 static void
  51 calcHistLookupTables_8u( const MatND& hist, const SparseMat& shist,
  52                         const float** ranges,
  53                         const double* uniranges,
  54                         bool uniform, bool issparse, vector<size_t>& _tab )
  55 {
  56     const int low = 0, high = 256;
  57     int i, j, dims = !issparse ? hist.dims : shist.dims();
  58     _tab.resize((high-low)*dims);
  59     size_t* tab = &_tab[0];
  60
  61     if( uniform )
  62     {
  63         for( i = 0; i < dims; i++ )
  64         {
  65             double a = uniranges[i*2];
  66             double b = uniranges[i*2+1];
  67             int sz = !issparse ? hist.size[i] : shist.size(i);
  68             size_t step = !issparse ? hist.step[i] : 1;
  69
  70             for( j = low; j < high; j++ )
  71             {
  72                 int idx = cvFloor(j*a + b);
  73                 size_t written_idx;
  74                 if( (unsigned)idx < (unsigned)sz )
  75                     written_idx = idx*step;
  76                 else
  77                     written_idx = OUT_OF_RANGE;
  78
  79                 tab[i*(high - low) + j - low] = written_idx;
  80             }
  81         }
  82     }
  83     else
  84     {
  85         for( i = 0; i < dims; i++ )
  86         {
  87             int limit = std::min(cvCeil(ranges[i][0]), high);
  88             int idx = -1, sz = !issparse ? hist.size[i] : shist.size(i);
  89             size_t written_idx = OUT_OF_RANGE;
  90             size_t step = !issparse ? hist.step[i] : 1;
  91
  92             for(j = low;;)
  93             {
  94                 for( ; j < limit; j++ )
  95                     tab[i*(high - low) + j - low] = written_idx;
  96
  97                 if( (unsigned)(++idx) < (unsigned)sz )
  98                 {
  99                     limit = std::min(cvCeil(ranges[i][idx+1]), high);
 100                     written_idx = idx*step;
 101                 }
 102                 else
 103                 {
 104                     for( ; j < high; j++ )
 105                         tab[i*(high - low) + j - low] = OUT_OF_RANGE;
 106                     break;
 107                 }
 108             }
 109         }
 110     }
 111 }
 112
 113
 114 static void histPrepareImages( const Mat* images, int nimages, const int* channels,
 115                               const Mat& mask, int dims, const int* histSize,
 116                               const float** ranges, bool uniform,
 117                               vector<uchar*>& ptrs, vector<int>& deltas,
 118                               Size& imsize, vector<double>& uniranges )
 119 {
 120     int i, j, c;
 121     if(!channels)
 122         CV_Assert( nimages == dims );
 123
 124     imsize = images[0].size();
 125     int depth = images[0].depth(), esz1 = (int)images[0].elemSize1();
 126     bool isContinuous = true;
 127
 128     ptrs.resize(dims + 1);
 129     deltas.resize((dims + 1)*2);
 130
 131     for( i = 0; i < dims; i++ )
 132     {
 133         if(!channels)
 134         {
 135             j = i;
 136             c = 0;
 137             CV_Assert( images[j].channels() == 1 );
 138         }
 139         else
 140         {
 141             c = channels[i];
 142             CV_Assert( c >= 0 );
 143             for( j = 0; j < nimages; c -= images[j].channels(), j++ )
 144                 if( c < images[j].channels() )
 145                     break;
 146             CV_Assert( j < nimages );
 147         }
 148
 149         CV_Assert( images[j].size() == imsize && images[j].depth() == depth );
 150         if( !images[j].isContinuous() )
 151             isContinuous = false;
 152         ptrs[i] = images[j].data + c*esz1;
 153         deltas[i*2] = images[j].channels();
 154         deltas[i*2+1] = (int)(images[j].step/esz1 - imsize.width*deltas[i*2]);
 155     }
 156
 157     if( mask.data )
 158     {
 159         CV_Assert( mask.size() == imsize && mask.channels() == 1 );
 160         isContinuous = isContinuous && mask.isContinuous();
 161         ptrs[dims] = mask.data;
 162         deltas[dims*2] = 1;
 163         deltas[dims*2 + 1] = (int)(mask.step/mask.elemSize1());
 164     }
 165
 166     if( isContinuous )
 167     {
 168         imsize.width *= imsize.height;
 169         imsize.height = 1;
 170     }
 171
 172     if( !ranges )
 173     {
 174         CV_Assert( depth == CV_8U );
 175
 176         uniranges.resize( dims*2 );
 177         for( i = 0; i < dims; i++ )
 178         {
 179             uniranges[i*2] = 1;
 180             uniranges[i*2+1] = 0;
 181         }
 182     }
 183     else if( uniform )
 184     {
 185         uniranges.resize( dims*2 );
 186         for( i = 0; i < dims; i++ )
 187         {
 188             CV_Assert( ranges[i] && ranges[i][0] < ranges[i][1] );
 189             double low = ranges[i][0], high = ranges[i][1];
 190             double t = histSize[i]/(high - low);
 191             uniranges[i*2] = t;
 192             uniranges[i*2+1] = -t*low;
 193         }
 194     }
 195     else
 196     {
 197         for( i = 0; i < dims; i++ )
 198         {
 199             size_t j, n = histSize[i];
 200             for( j = 0; j < n; j++ )
 201                 CV_Assert( ranges[i][j] < ranges[i][j+1] );
 202         }
 203     }
 204 }
 205
 206
 207 ////////////////////////////////// C A L C U L A T E    H I S T O G R A M ////////////////////////////////////
 208
 209 template<typename T> static void
 210 calcHist_( vector<uchar*>& _ptrs, const vector<int>& _deltas,
 211            Size imsize, MatND& hist, const float** _ranges,
 212            const double* _uniranges, bool uniform )
 213 {
 214     T** ptrs = (T**)&_ptrs[0];
 215     const int* deltas = &_deltas[0];
 216     uchar* H = hist.data;
 217     int i, x, dims = hist.dims;
 218     const uchar* mask = _ptrs[dims];
 219     int mstep = _deltas[dims*2 + 1];
 220     int size[CV_MAX_DIM];
 221     size_t hstep[CV_MAX_DIM];
 222
 223     for( i = 0; i < dims; i++ )
 224     {
 225         size[i] = hist.size[i];
 226         hstep[i] = hist.step[i];
 227     }
 228
 229     if( uniform )
 230     {
 231         const double* uniranges = &_uniranges[0];
 232
 233         if( dims == 1 )
 234         {
 235             double a = uniranges[0], b = uniranges[1];
 236             int sz = size[0], d0 = deltas[0], step0 = deltas[1];
 237             const T* p0 = (const T*)ptrs[0];
 238
 239             for( ; imsize.height--; p0 += step0, mask += mstep )
 240             {
 241                 if( !mask )
 242                     for( x = 0; x < imsize.width; x++, p0 += d0 )
 243                     {
 244                         int idx = cvFloor(*p0*a + b);
 245                         if( (unsigned)idx < (unsigned)sz )
 246                             ((int*)H)[idx]++;
 247                     }
 248                 else
 249                     for( x = 0; x < imsize.width; x++, p0 += d0 )
 250                         if( mask[x] )
 251                         {
 252                             int idx = cvFloor(*p0*a + b);
 253                             if( (unsigned)idx < (unsigned)sz )
 254                                 ((int*)H)[idx]++;
 255                         }
 256             }
 257         }
 258         else if( dims == 2 )
 259         {
 260             double a0 = uniranges[0], b0 = uniranges[1], a1 = uniranges[2], b1 = uniranges[3];
 261             int sz0 = size[0], sz1 = size[1];
 262             int d0 = deltas[0], step0 = deltas[1],
 263                 d1 = deltas[2], step1 = deltas[3];
 264             size_t hstep0 = hstep[0];
 265             const T* p0 = (const T*)ptrs[0];
 266             const T* p1 = (const T*)ptrs[1];
 267
 268             for( ; imsize.height--; p0 += step0, p1 += step1, mask += mstep )
 269             {
 270                 if( !mask )
 271                     for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
 272                     {
 273                         int idx0 = cvFloor(*p0*a0 + b0);
 274                         int idx1 = cvFloor(*p1*a1 + b1);
 275                         if( (unsigned)idx0 < (unsigned)sz0 && (unsigned)idx1 < (unsigned)sz1 )
 276                             ((int*)(H + hstep0*idx0))[idx1]++;
 277                     }
 278                 else
 279                     for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
 280                         if( mask[x] )
 281                         {
 282                             int idx0 = cvFloor(*p0*a0 + b0);
 283                             int idx1 = cvFloor(*p1*a1 + b1);
 284                             if( (unsigned)idx0 < (unsigned)sz0 && (unsigned)idx1 < (unsigned)sz1 )
 285                                 ((int*)(H + hstep0*idx0))[idx1]++;
 286                         }
 287             }
 288         }
 289         else if( dims == 3 )
 290         {
 291             double a0 = uniranges[0], b0 = uniranges[1],
 292                    a1 = uniranges[2], b1 = uniranges[3],
 293                    a2 = uniranges[4], b2 = uniranges[5];
 294             int sz0 = size[0], sz1 = size[1], sz2 = size[2];
 295             int d0 = deltas[0], step0 = deltas[1],
 296                 d1 = deltas[2], step1 = deltas[3],
 297                 d2 = deltas[4], step2 = deltas[5];
 298             size_t hstep0 = hstep[0], hstep1 = hstep[1];
 299             const T* p0 = (const T*)ptrs[0];
 300             const T* p1 = (const T*)ptrs[1];
 301             const T* p2 = (const T*)ptrs[2];
 302
 303             for( ; imsize.height--; p0 += step0, p1 += step1, p2 += step2, mask += mstep )
 304             {
 305                 if( !mask )
 306                     for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1, p2 += d2 )
 307                     {
 308                         int idx0 = cvFloor(*p0*a0 + b0);
 309                         int idx1 = cvFloor(*p1*a1 + b1);
 310                         int idx2 = cvFloor(*p2*a2 + b2);
 311                         if( (unsigned)idx0 < (unsigned)sz0 &&
 312                             (unsigned)idx1 < (unsigned)sz1 &&
 313                             (unsigned)idx2 < (unsigned)sz2 )
 314                             ((int*)(H + hstep0*idx0 + hstep1*idx1))[idx2]++;
 315                     }
 316                 else
 317                     for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1, p2 += d2 )
 318                         if( mask[x] )
 319                         {
 320                             int idx0 = cvFloor(*p0*a0 + b0);
 321                             int idx1 = cvFloor(*p1*a1 + b1);
 322                             int idx2 = cvFloor(*p2*a2 + b2);
 323                             if( (unsigned)idx0 < (unsigned)sz0 &&
 324                                (unsigned)idx1 < (unsigned)sz1 &&
 325                                (unsigned)idx2 < (unsigned)sz2 )
 326                                 ((int*)(H + hstep0*idx0 + hstep1*idx1))[idx2]++;
 327                         }
 328             }
 329         }
 330         else
 331         {
 332             for( ; imsize.height--; mask += mstep )
 333             {
 334                 if( !mask )
 335                     for( x = 0; x < imsize.width; x++ )
 336                     {
 337                         uchar* Hptr = H;
 338                         for( i = 0; i < dims; i++ )
 339                         {
 340                             int idx = cvFloor(*ptrs[i]*uniranges[i*2] + uniranges[i*2+1]);
 341                             if( (unsigned)idx >= (unsigned)size[i] )
 342                                 break;
 343                             ptrs[i] += deltas[i*2];
 344                             Hptr += idx*hstep[i];
 345                         }
 346
 347                         if( i == dims )
 348                             ++*((int*)Hptr);
 349                         else
 350                             for( ; i < dims; i++ )
 351                                 ptrs[i] += deltas[i*2];
 352                     }
 353                 else
 354                     for( x = 0; x < imsize.width; x++ )
 355                     {
 356                         uchar* Hptr = H;
 357                         i = 0;
 358                         if( mask[x] )
 359                             for( ; i < dims; i++ )
 360                             {
 361                                 int idx = cvFloor(*ptrs[i]*uniranges[i*2] + uniranges[i*2+1]);
 362                                 if( (unsigned)idx >= (unsigned)size[i] )
 363                                     break;
 364                                 ptrs[i] += deltas[i*2];
 365                                 Hptr += idx*hstep[i];
 366                             }
 367
 368                         if( i == dims )
 369                             ++*((int*)Hptr);
 370                         else
 371                             for( ; i < dims; i++ )
 372                                 ptrs[i] += deltas[i*2];
 373                     }
 374                 for( i = 0; i < dims; i++ )
 375                     ptrs[i] += deltas[i*2 + 1];
 376             }
 377         }
 378     }
 379     else
 380     {
 381         // non-uniform histogram
 382         const float* ranges[CV_MAX_DIM];
 383         for( i = 0; i < dims; i++ )
 384             ranges[i] = &_ranges[i][0];
 385
 386         for( ; imsize.height--; mask += mstep )
 387         {
 388             for( x = 0; x < imsize.width; x++ )
 389             {
 390                 uchar* Hptr = H;
 391                 i = 0;
 392
 393                 if( !mask || mask[x] )
 394                     for( ; i < dims; i++ )
 395                     {
 396                         float v = (float)*ptrs[i];
 397                         const float* R = ranges[i];
 398                         int idx = -1, sz = size[i];
 399
 400                         while( v >= R[idx+1] && ++idx < sz )
 401                             ; // nop
 402
 403                         if( (unsigned)idx >= (unsigned)sz )
 404                             break;
 405
 406                         ptrs[i] += deltas[i*2];
 407                         Hptr += idx*hstep[i];
 408                     }
 409
 410                 if( i == dims )
 411                     ++*((int*)Hptr);
 412                 else
 413                     for( ; i < dims; i++ )
 414                         ptrs[i] += deltas[i*2];
 415             }
 416
 417             for( i = 0; i < dims; i++ )
 418                 ptrs[i] += deltas[i*2 + 1];
 419         }
 420     }
 421 }
 422
 423
 424 static void
 425 calcHist_8u( vector<uchar*>& _ptrs, const vector<int>& _deltas,
 426              Size imsize, MatND& hist, const float** _ranges,
 427              const double* _uniranges, bool uniform )
 428 {
 429     uchar** ptrs = &_ptrs[0];
 430     const int* deltas = &_deltas[0];
 431     uchar* H = hist.data;
 432     int i, x, dims = hist.dims;
 433     const uchar* mask = _ptrs[dims];
 434     int mstep = _deltas[dims*2 + 1];
 435     vector<size_t> _tab;
 436
 437     calcHistLookupTables_8u( hist, SparseMat(), _ranges, _uniranges, uniform, false, _tab );
 438     const size_t* tab = &_tab[0];
 439
 440     if( dims == 1 )
 441     {
 442         int d0 = deltas[0], step0 = deltas[1];
 443         int _H[256] = {0};
 444         const uchar* p0 = (const uchar*)ptrs[0];
 445
 446         for( ; imsize.height--; p0 += step0, mask += mstep )
 447         {
 448             if( !mask )
 449             {
 450                 if( d0 == 1 )
 451                 {
 452                     for( x = 0; x <= imsize.width - 4; x += 4 )
 453                     {
 454                         int t0 = p0[x], t1 = p0[x+1];
 455                         _H[t0]++; _H[t1]++;
 456                         t0 = p0[x+2]; t1 = p0[x+3];
 457                         _H[t0]++; _H[t1]++;
 458                     }
 459                     p0 += x;
 460                 }
 461                 else
 462                     for( x = 0; x <= imsize.width - 4; x += 4 )
 463                     {
 464                         int t0 = p0[0], t1 = p0[d0];
 465                         _H[t0]++; _H[t1]++;
 466                         p0 += d0*2;
 467                         t0 = p0[0]; t1 = p0[d0];
 468                         _H[t0]++; _H[t1]++;
 469                         p0 += d0*2;
 470                     }
 471
 472                 for( ; x < imsize.width; x++, p0 += d0 )
 473                     _H[*p0]++;
 474             }
 475             else
 476                 for( x = 0; x < imsize.width; x++, p0 += d0 )
 477                     if( mask[x] )
 478                         _H[*p0]++;
 479         }
 480
 481         for( i = 0; i < 256; i++ )
 482         {
 483             size_t hidx = tab[i];
 484             if( hidx < OUT_OF_RANGE )
 485                 *(int*)(H + hidx) += _H[i];
 486         }
 487     }
 488     else if( dims == 2 )
 489     {
 490         int d0 = deltas[0], step0 = deltas[1],
 491             d1 = deltas[2], step1 = deltas[3];
 492         const uchar* p0 = (const uchar*)ptrs[0];
 493         const uchar* p1 = (const uchar*)ptrs[1];
 494
 495         for( ; imsize.height--; p0 += step0, p1 += step1, mask += mstep )
 496         {
 497             if( !mask )
 498                 for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
 499                 {
 500                     size_t idx = tab[*p0] + tab[*p1 + 256];
 501                     if( idx < OUT_OF_RANGE )
 502                         ++*(int*)(H + idx);
 503                 }
 504             else
 505                 for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
 506                 {
 507                     size_t idx;
 508                     if( mask[x] && (idx = tab[*p0] + tab[*p1 + 256]) < OUT_OF_RANGE )
 509                         ++*(int*)(H + idx);
 510                 }
 511         }
 512     }
 513     else if( dims == 3 )
 514     {
 515         int d0 = deltas[0], step0 = deltas[1],
 516             d1 = deltas[2], step1 = deltas[3],
 517             d2 = deltas[4], step2 = deltas[5];
 518
 519         const uchar* p0 = (const uchar*)ptrs[0];
 520         const uchar* p1 = (const uchar*)ptrs[1];
 521         const uchar* p2 = (const uchar*)ptrs[2];
 522
 523         for( ; imsize.height--; p0 += step0, p1 += step1, p2 += step2, mask += mstep )
 524         {
 525             if( !mask )
 526                 for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1, p2 += d2 )
 527                 {
 528                     size_t idx = tab[*p0] + tab[*p1 + 256] + tab[*p2 + 512];
 529                     if( idx < OUT_OF_RANGE )
 530                         ++*(int*)(H + idx);
 531                 }
 532             else
 533                 for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1, p2 += d2 )
 534                 {
 535                     size_t idx;
 536                     if( mask[x] && (idx = tab[*p0] + tab[*p1 + 256] + tab[*p2 + 512]) < OUT_OF_RANGE )
 537                         ++*(int*)(H + idx);
 538                 }
 539         }
 540     }
 541     else
 542     {
 543         for( ; imsize.height--; mask += mstep )
 544         {
 545             if( !mask )
 546                 for( x = 0; x < imsize.width; x++ )
 547                 {
 548                     uchar* Hptr = H;
 549                     for( i = 0; i < dims; i++ )
 550                     {
 551                         size_t idx = tab[*ptrs[i] + i*256];
 552                         if( idx >= OUT_OF_RANGE )
 553                             break;
 554                         Hptr += idx;
 555                         ptrs[i] += deltas[i*2];
 556                     }
 557
 558                     if( i == dims )
 559                         ++*((int*)Hptr);
 560                     else
 561                         for( ; i < dims; i++ )
 562                             ptrs[i] += deltas[i*2];
 563                 }
 564             else
 565                 for( x = 0; x < imsize.width; x++ )
 566                 {
 567                     uchar* Hptr = H;
 568                     int i = 0;
 569                     if( mask[x] )
 570                         for( ; i < dims; i++ )
 571                         {
 572                             size_t idx = tab[*ptrs[i] + i*256];
 573                             if( idx >= OUT_OF_RANGE )
 574                                 break;
 575                             Hptr += idx;
 576                             ptrs[i] += deltas[i*2];
 577                         }
 578
 579                     if( i == dims )
 580                         ++*((int*)Hptr);
 581                     else
 582                         for( ; i < dims; i++ )
 583                             ptrs[i] += deltas[i*2];
 584                 }
 585             for( i = 0; i < dims; i++ )
 586                 ptrs[i] += deltas[i*2 + 1];
 587         }
 588     }
 589 }
 590
 591
 592 void calcHist( const Mat* images, int nimages, const int* channels,
 593                const Mat& mask, MatND& hist, int dims, const int* histSize,
 594                const float** ranges, bool uniform, bool accumulate )
 595 {
 596     CV_Assert(dims > 0 && histSize);
 597     hist.create(dims, histSize, CV_32F);
 598
 599     MatND ihist = hist;
 600     ihist.flags = (ihist.flags & ~CV_MAT_TYPE_MASK)|CV_32S;
 601
 602     if( !accumulate )
 603         hist = Scalar(0.);
 604     else
 605         hist.convertTo(ihist, CV_32S);
 606
 607     vector<uchar*> ptrs;
 608     vector<int> deltas;
 609     vector<double> uniranges;
 610     Size imsize;
 611
 612     CV_Assert( !mask.data || mask.type() == CV_8UC1 );
 613     histPrepareImages( images, nimages, channels, mask, hist.dims, hist.size, ranges,
 614                        uniform, ptrs, deltas, imsize, uniranges );
 615     const double* _uniranges = uniform ? &uniranges[0] : 0;
 616
 617     int depth = images[0].depth();
 618     if( depth == CV_8U )
 619         calcHist_8u(ptrs, deltas, imsize, ihist, ranges, _uniranges, uniform );
 620     else if( depth == CV_32F )
 621         calcHist_<float>(ptrs, deltas, imsize, ihist, ranges, _uniranges, uniform );
 622
 623     ihist.convertTo(hist, CV_32F);
 624 }
 625
 626
 627 template<typename T> static void
 628 calcSparseHist_( vector<uchar*>& _ptrs, const vector<int>& _deltas,
 629                  Size imsize, SparseMat& hist, const float** _ranges,
 630                  const double* _uniranges, bool uniform )
 631 {
 632     T** ptrs = (T**)&_ptrs[0];
 633     const int* deltas = &_deltas[0];
 634     int i, x, dims = hist.dims();
 635     const uchar* mask = _ptrs[dims];
 636     int mstep = _deltas[dims*2 + 1];
 637     const int* size = hist.hdr->size;
 638     int idx[CV_MAX_DIM];
 639
 640     if( uniform )
 641     {
 642         const double* uniranges = &_uniranges[0];
 643
 644         for( ; imsize.height--; mask += mstep )
 645         {
 646             for( x = 0; x < imsize.width; x++ )
 647             {
 648                 i = 0;
 649                 if( !mask || mask[x] )
 650                     for( ; i < dims; i++ )
 651                     {
 652                         idx[i] = cvFloor(*ptrs[i]*uniranges[i*2] + uniranges[i*2+1]);
 653                         if( (unsigned)idx[i] >= (unsigned)size[i] )
 654                             break;
 655                         ptrs[i] += deltas[i*2];
 656                     }
 657
 658                 if( i == dims )
 659                     ++*(int*)hist.ptr(idx, true);
 660                 else
 661                     for( ; i < dims; i++ )
 662                         ptrs[i] += deltas[i*2];
 663             }
 664             for( i = 0; i < dims; i++ )
 665                 ptrs[i] += deltas[i*2 + 1];
 666         }
 667     }
 668     else
 669     {
 670         // non-uniform histogram
 671         const float* ranges[CV_MAX_DIM];
 672         for( i = 0; i < dims; i++ )
 673             ranges[i] = &_ranges[i][0];
 674
 675         for( ; imsize.height--; mask += mstep )
 676         {
 677             for( x = 0; x < imsize.width; x++ )
 678             {
 679                 i = 0;
 680
 681                 if( !mask || mask[x] )
 682                     for( ; i < dims; i++ )
 683                     {
 684                         float v = (float)*ptrs[i];
 685                         const float* R = ranges[i];
 686                         int j = -1, sz = size[i];
 687
 688                         while( v >= R[j+1] && ++j < sz )
 689                             ; // nop
 690
 691                         if( (unsigned)j >= (unsigned)sz )
 692                             break;
 693                         ptrs[i] += deltas[i*2];
 694                         idx[i] = j;
 695                     }
 696
 697                 if( i == dims )
 698                     ++*(int*)hist.ptr(idx, true);
 699                 else
 700                     for( ; i < dims; i++ )
 701                         ptrs[i] += deltas[i*2];
 702             }
 703
 704             for( i = 0; i < dims; i++ )
 705                 ptrs[i] += deltas[i*2 + 1];
 706         }
 707     }
 708 }
 709
 710
 711 static void
 712 calcSparseHist_8u( vector<uchar*>& _ptrs, const vector<int>& _deltas,
 713                    Size imsize, SparseMat& hist, const float** _ranges,
 714                    const double* _uniranges, bool uniform )
 715 {
 716     uchar** ptrs = (uchar**)&_ptrs[0];
 717     const int* deltas = &_deltas[0];
 718     int i, x, dims = hist.dims();
 719     const uchar* mask = _ptrs[dims];
 720     int mstep = _deltas[dims*2 + 1];
 721     int idx[CV_MAX_DIM];
 722     vector<size_t> _tab;
 723
 724     calcHistLookupTables_8u( MatND(), hist, _ranges, _uniranges, uniform, true, _tab );
 725     const size_t* tab = &_tab[0];
 726
 727     for( ; imsize.height--; mask += mstep )
 728     {
 729         for( x = 0; x < imsize.width; x++ )
 730         {
 731             int i = 0;
 732             if( !mask || mask[x] )
 733                 for( ; i < dims; i++ )
 734                 {
 735                     size_t hidx = tab[*ptrs[i] + i*256];
 736                     if( hidx >= OUT_OF_RANGE )
 737                         break;
 738                     ptrs[i] += deltas[i*2];
 739                     idx[i] = (int)hidx;
 740                 }
 741
 742             if( i == dims )
 743                 ++*(int*)hist.ptr(idx,true);
 744             else
 745                 for( ; i < dims; i++ )
 746                     ptrs[i] += deltas[i*2];
 747         }
 748         for( i = 0; i < dims; i++ )
 749             ptrs[i] += deltas[i*2 + 1];
 750     }
 751 }
 752
 753
 754 static void calcHist( const Mat* images, int nimages, const int* channels,
 755                       const Mat& mask, SparseMat& hist, int dims, const int* histSize,
 756                       const float** ranges, bool uniform, bool accumulate, bool keepInt )
 757 {
 758     SparseMatIterator it;
 759     size_t i, N;
 760
 761     if( !accumulate )
 762         hist.create(dims, histSize, CV_32F);
 763     else
 764         for( i = 0, N = hist.nzcount(); i < N; i++, ++it )
 765         {
 766             int* value = (int*)it.ptr;
 767             *value = cvRound(*(const float*)value);
 768         }
 769
 770     vector<uchar*> ptrs;
 771     vector<int> deltas;
 772     vector<double> uniranges;
 773     Size imsize;
 774
 775     CV_Assert( !mask.data || mask.type() == CV_8UC1 );
 776     histPrepareImages( images, nimages, channels, mask, hist.dims(), hist.hdr->size, ranges,
 777                        uniform, ptrs, deltas, imsize, uniranges );
 778
 779     int depth = images[0].depth();
 780     if( depth == CV_8U )
 781         calcSparseHist_8u(ptrs, deltas, imsize, hist, ranges, &uniranges[0], uniform );
 782     else if( depth == CV_32F )
 783         calcSparseHist_<float>(ptrs, deltas, imsize, hist, ranges, &uniranges[0], uniform );
 784
 785     if( !keepInt )
 786         for( i = 0, N = hist.nzcount(); i < N; i++, ++it )
 787         {
 788             int* value = (int*)it.ptr;
 789             *(float*)value = (float)*value;
 790         }
 791 }
 792
 793
 794 void calcHist( const Mat* images, int nimages, const int* channels,
 795                const Mat& mask, SparseMat& hist, int dims, const int* histSize,
 796                const float** ranges, bool uniform, bool accumulate )
 797 {
 798     calcHist( images, nimages, channels, mask, hist, dims, histSize,
 799               ranges, uniform, accumulate, false );
 800 }
 801
 802
 803 /////////////////////////////////////// B A C K   P R O J E C T ////////////////////////////////////
 804
 805
 806 template<typename T, typename BT> static void
 807 calcBackProj_( vector<uchar*>& _ptrs, const vector<int>& _deltas,
 808                Size imsize, const MatND& hist, const float** _ranges,
 809                const double* _uniranges, float scale, bool uniform )
 810 {
 811     T** ptrs = (T**)&_ptrs[0];
 812     const int* deltas = &_deltas[0];
 813     uchar* H = hist.data;
 814     int i, x, dims = hist.dims;
 815     BT* bproj = (BT*)_ptrs[dims];
 816     int bpstep = _deltas[dims*2 + 1];
 817     int size[CV_MAX_DIM];
 818     size_t hstep[CV_MAX_DIM];
 819
 820     for( i = 0; i < dims; i++ )
 821     {
 822         size[i] = hist.size[i];
 823         hstep[i] = hist.step[i];
 824     }
 825
 826     if( uniform )
 827     {
 828         const double* uniranges = &_uniranges[0];
 829
 830         if( dims == 1 )
 831         {
 832             double a = uniranges[0], b = uniranges[1];
 833             int sz = size[0], d0 = deltas[0], step0 = deltas[1];
 834             const T* p0 = (const T*)ptrs[0];
 835
 836             for( ; imsize.height--; p0 += step0, bproj += bpstep )
 837             {
 838                 for( x = 0; x < imsize.width; x++, p0 += d0 )
 839                 {
 840                     int idx = cvFloor(*p0*a + b);
 841                     bproj[x] = (unsigned)idx < (unsigned)sz ? saturate_cast<BT>(((float*)H)[idx]*scale) : 0;
 842                 }
 843             }
 844         }
 845         else if( dims == 2 )
 846         {
 847             double a0 = uniranges[0], b0 = uniranges[1],
 848                    a1 = uniranges[2], b1 = uniranges[3];
 849             int sz0 = size[0], sz1 = size[1];
 850             int d0 = deltas[0], step0 = deltas[1],
 851                 d1 = deltas[2], step1 = deltas[3];
 852             size_t hstep0 = hstep[0];
 853             const T* p0 = (const T*)ptrs[0];
 854             const T* p1 = (const T*)ptrs[1];
 855
 856             for( ; imsize.height--; p0 += step0, p1 += step1, bproj += bpstep )
 857             {
 858                 for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
 859                 {
 860                     int idx0 = cvFloor(*p0*a0 + b0);
 861                     int idx1 = cvFloor(*p1*a1 + b1);
 862                     bproj[x] = (unsigned)idx0 < (unsigned)sz0 &&
 863                                (unsigned)idx1 < (unsigned)sz1 ?
 864                         saturate_cast<BT>(((float*)(H + hstep0*idx0))[idx1]*scale) : 0;
 865                 }
 866             }
 867         }
 868         else if( dims == 3 )
 869         {
 870             double a0 = uniranges[0], b0 = uniranges[1],
 871                    a1 = uniranges[2], b1 = uniranges[3],
 872                    a2 = uniranges[4], b2 = uniranges[5];
 873             int sz0 = size[0], sz1 = size[1], sz2 = size[2];
 874             int d0 = deltas[0], step0 = deltas[1],
 875                 d1 = deltas[2], step1 = deltas[3],
 876                 d2 = deltas[4], step2 = deltas[5];
 877             size_t hstep0 = hstep[0], hstep1 = hstep[1];
 878             const T* p0 = (const T*)ptrs[0];
 879             const T* p1 = (const T*)ptrs[1];
 880             const T* p2 = (const T*)ptrs[2];
 881
 882             for( ; imsize.height--; p0 += step0, p1 += step1, p2 += step2, bproj += bpstep )
 883             {
 884                 for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1, p2 += d2 )
 885                 {
 886                     int idx0 = cvFloor(*p0*a0 + b0);
 887                     int idx1 = cvFloor(*p1*a1 + b1);
 888                     int idx2 = cvFloor(*p2*a2 + b2);
 889                     bproj[x] = (unsigned)idx0 < (unsigned)sz0 &&
 890                                (unsigned)idx1 < (unsigned)sz1 &&
 891                                (unsigned)idx2 < (unsigned)sz2 ?
 892                         saturate_cast<BT>(((float*)(H + hstep0*idx0 + hstep1*idx1))[idx2]*scale) : 0;
 893                 }
 894             }
 895         }
 896         else
 897         {
 898             for( ; imsize.height--; bproj += bpstep )
 899             {
 900                 for( x = 0; x < imsize.width; x++ )
 901                 {
 902                     uchar* Hptr = H;
 903                     for( i = 0; i < dims; i++ )
 904                     {
 905                         int idx = cvFloor(*ptrs[i]*uniranges[i*2] + uniranges[i*2+1]);
 906                         if( (unsigned)idx >= (unsigned)size[i] )
 907                             break;
 908                         ptrs[i] += deltas[i*2];
 909                         Hptr += idx*hstep[i];
 910                     }
 911
 912                     if( i == dims )
 913                         bproj[x] = saturate_cast<BT>(*(float*)Hptr*scale);
 914                     else
 915                     {
 916                         bproj[x] = 0;
 917                         for( ; i < dims; i++ )
 918                             ptrs[i] += deltas[i*2];
 919                     }
 920                 }
 921                 for( i = 0; i < dims; i++ )
 922                     ptrs[i] += deltas[i*2 + 1];
 923             }
 924         }
 925     }
 926     else
 927     {
 928         // non-uniform histogram
 929         const float* ranges[CV_MAX_DIM];
 930         for( i = 0; i < dims; i++ )
 931             ranges[i] = &_ranges[i][0];
 932
 933         for( ; imsize.height--; bproj += bpstep )
 934         {
 935             for( x = 0; x < imsize.width; x++ )
 936             {
 937                 uchar* Hptr = H;
 938                 for( i = 0; i < dims; i++ )
 939                 {
 940                     float v = (float)*ptrs[i];
 941                     const float* R = ranges[i];
 942                     int idx = -1, sz = size[i];
 943
 944                     while( v >= R[idx+1] && ++idx < sz )
 945                         ; // nop
 946
 947                     if( (unsigned)idx >= (unsigned)sz )
 948                         break;
 949
 950                     ptrs[i] += deltas[i*2];
 951                     Hptr += idx*hstep[i];
 952                 }
 953
 954                 if( i == dims )
 955                     bproj[x] = saturate_cast<BT>(*(float*)Hptr*scale);
 956                 else
 957                 {
 958                     bproj[x] = 0;
 959                     for( ; i < dims; i++ )
 960                         ptrs[i] += deltas[i*2];
 961                 }
 962             }
 963
 964             for( i = 0; i < dims; i++ )
 965                 ptrs[i] += deltas[i*2 + 1];
 966         }
 967     }
 968 }
 969
 970
 971 static void
 972 calcBackProj_8u( vector<uchar*>& _ptrs, const vector<int>& _deltas,
 973                  Size imsize, const MatND& hist, const float** _ranges,
 974                  const double* _uniranges, float scale, bool uniform )
 975 {
 976     uchar** ptrs = &_ptrs[0];
 977     const int* deltas = &_deltas[0];
 978     uchar* H = hist.data;
 979     int i, x, dims = hist.dims;
 980     uchar* bproj = _ptrs[dims];
 981     int bpstep = _deltas[dims*2 + 1];
 982     vector<size_t> _tab;
 983
 984     calcHistLookupTables_8u( hist, SparseMat(), _ranges, _uniranges, uniform, false, _tab );
 985     const size_t* tab = &_tab[0];
 986
 987     if( dims == 1 )
 988     {
 989         int d0 = deltas[0], step0 = deltas[1];
 990         uchar _H[256] = {0};
 991         const uchar* p0 = (const uchar*)ptrs[0];
 992
 993         for( i = 0; i < 256; i++ )
 994         {
 995             size_t hidx = tab[i];
 996             if( hidx < OUT_OF_RANGE )
 997                 _H[i] = saturate_cast<uchar>(*(float*)(H + hidx)*scale);
 998         }
 999
1000         for( ; imsize.height--; p0 += step0, bproj += bpstep )
1001         {
1002             if( d0 == 1 )
1003             {
1004                 for( x = 0; x <= imsize.width - 4; x += 4 )
1005                 {
1006                     uchar t0 = _H[p0[x]], t1 = _H[p0[x+1]];
1007                     bproj[x] = t0; bproj[x+1] = t1;
1008                     t0 = _H[p0[x+2]]; t1 = _H[p0[x+3]];
1009                     bproj[x+2] = t0; bproj[x+3] = t1;
1010                 }
1011                 p0 += x;
1012             }
1013             else
1014                 for( x = 0; x <= imsize.width - 4; x += 4 )
1015                 {
1016                     uchar t0 = _H[p0[0]], t1 = _H[p0[d0]];
1017                     bproj[x] = t0; bproj[x+1] = t1;
1018                     p0 += d0*2;
1019                     t0 = _H[p0[0]]; t1 = _H[p0[d0]];
1020                     bproj[x+2] = t0; bproj[x+3] = t1;
1021                     p0 += d0*2;
1022                 }
1023
1024             for( ; x < imsize.width; x++, p0 += d0 )
1025                 bproj[x] = _H[*p0];
1026         }
1027     }
1028     else if( dims == 2 )
1029     {
1030         int d0 = deltas[0], step0 = deltas[1],
1031             d1 = deltas[2], step1 = deltas[3];
1032         const uchar* p0 = (const uchar*)ptrs[0];
1033         const uchar* p1 = (const uchar*)ptrs[1];
1034
1035         for( ; imsize.height--; p0 += step0, p1 += step1, bproj += bpstep )
1036         {
1037             for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
1038             {
1039                 size_t idx = tab[*p0] + tab[*p1 + 256];
1040                 bproj[x] = idx < OUT_OF_RANGE ? saturate_cast<uchar>(*(float*)(H + idx)*scale) : 0;
1041             }
1042         }
1043     }
1044     else if( dims == 3 )
1045     {
1046         int d0 = deltas[0], step0 = deltas[1],
1047         d1 = deltas[2], step1 = deltas[3],
1048         d2 = deltas[4], step2 = deltas[5];
1049         const uchar* p0 = (const uchar*)ptrs[0];
1050         const uchar* p1 = (const uchar*)ptrs[1];
1051         const uchar* p2 = (const uchar*)ptrs[2];
1052
1053         for( ; imsize.height--; p0 += step0, p1 += step1, p2 += step2, bproj += bpstep )
1054         {
1055             for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1, p2 += d2 )
1056             {
1057                 size_t idx = tab[*p0] + tab[*p1 + 256] + tab[*p2 + 512];
1058                 bproj[x] = idx < OUT_OF_RANGE ? saturate_cast<uchar>(*(float*)(H + idx)*scale) : 0;
1059             }
1060         }
1061     }
1062     else
1063     {
1064         for( ; imsize.height--; bproj += bpstep )
1065         {
1066             for( x = 0; x < imsize.width; x++ )
1067             {
1068                 uchar* Hptr = H;
1069                 for( i = 0; i < dims; i++ )
1070                 {
1071                     size_t idx = tab[*ptrs[i] + i*256];
1072                     if( idx >= OUT_OF_RANGE )
1073                         break;
1074                     ptrs[i] += deltas[i*2];
1075                     Hptr += idx;
1076                 }
1077
1078                 if( i == dims )
1079                     bproj[x] = saturate_cast<uchar>(*(float*)Hptr*scale);
1080                 else
1081                 {
1082                     bproj[x] = 0;
1083                     for( ; i < dims; i++ )
1084                         ptrs[i] += deltas[i*2];
1085                 }
1086             }
1087             for( i = 0; i < dims; i++ )
1088                 ptrs[i] += deltas[i*2 + 1];
1089         }
1090     }
1091 }
1092
1093
1094 void calcBackProject( const Mat* images, int nimages, const int* channels,
1095                       const MatND& hist, Mat& backProject,
1096                       const float** ranges, double scale, bool uniform )
1097 {
1098     vector<uchar*> ptrs;
1099     vector<int> deltas;
1100     vector<double> uniranges;
1101     Size imsize;
1102
1103     CV_Assert( hist.dims > 0 && hist.data );
1104     backProject.create( images[0].size(), images[0].depth() );
1105     histPrepareImages( images, nimages, channels, backProject, hist.dims, hist.size, ranges,
1106                        uniform, ptrs, deltas, imsize, uniranges );
1107     const double* _uniranges = uniform ? &uniranges[0] : 0;
1108
1109     int depth = images[0].depth();
1110     if( depth == CV_8U )
1111         calcBackProj_8u(ptrs, deltas, imsize, hist, ranges, _uniranges, (float)scale, uniform);
1112     else if( depth == CV_32F )
1113         calcBackProj_<float, float>(ptrs, deltas, imsize, hist, ranges, _uniranges, (float)scale, uniform );
1114 }
1115
1116
1117 template<typename T, typename BT> static void
1118 calcSparseBackProj_( vector<uchar*>& _ptrs, const vector<int>& _deltas,
1119                      Size imsize, const SparseMat& hist, const float** _ranges,
1120                      const double* _uniranges, float scale, bool uniform )
1121 {
1122     T** ptrs = (T**)&_ptrs[0];
1123     const int* deltas = &_deltas[0];
1124     int i, x, dims = hist.dims();
1125     BT* bproj = (BT*)_ptrs[dims];
1126     int bpstep = _deltas[dims*2 + 1];
1127     const int* size = hist.hdr->size;
1128     int idx[CV_MAX_DIM];
1129     const SparseMat_<float>& hist_ = (const SparseMat_<float>&)hist;
1130
1131     if( uniform )
1132     {
1133         const double* uniranges = &_uniranges[0];
1134         for( ; imsize.height--; bproj += bpstep )
1135         {
1136             for( x = 0; x < imsize.width; x++ )
1137             {
1138                 for( i = 0; i < dims; i++ )
1139                 {
1140                     idx[i] = cvFloor(*ptrs[i]*uniranges[i*2] + uniranges[i*2+1]);
1141                     if( (unsigned)idx[i] >= (unsigned)size[i] )
1142                         break;
1143                     ptrs[i] += deltas[i*2];
1144                 }
1145
1146                 if( i == dims )
1147                     bproj[x] = saturate_cast<BT>(hist_(idx)*scale);
1148                 else
1149                 {
1150                     bproj[x] = 0;
1151                     for( ; i < dims; i++ )
1152                         ptrs[i] += deltas[i*2];
1153                 }
1154             }
1155             for( i = 0; i < dims; i++ )
1156                 ptrs[i] += deltas[i*2 + 1];
1157         }
1158     }
1159     else
1160     {
1161         // non-uniform histogram
1162         const float* ranges[CV_MAX_DIM];
1163         for( i = 0; i < dims; i++ )
1164             ranges[i] = &_ranges[i][0];
1165
1166         for( ; imsize.height--; bproj += bpstep )
1167         {
1168             for( x = 0; x < imsize.width; x++ )
1169             {
1170                 for( i = 0; i < dims; i++ )
1171                 {
1172                     float v = (float)*ptrs[i];
1173                     const float* R = ranges[i];
1174                     int j = -1, sz = size[i];
1175
1176                     while( v >= R[j+1] && ++j < sz )
1177                         ; // nop
1178
1179                     if( (unsigned)j >= (unsigned)sz )
1180                         break;
1181                     idx[i] = j;
1182                     ptrs[i] += deltas[i*2];
1183                 }
1184
1185                 if( i == dims )
1186                     bproj[x] = saturate_cast<BT>(hist_(idx)*scale);
1187                 else
1188                 {
1189                     bproj[x] = 0;
1190                     for( ; i < dims; i++ )
1191                         ptrs[i] += deltas[i*2];
1192                 }
1193             }
1194
1195             for( i = 0; i < dims; i++ )
1196                 ptrs[i] += deltas[i*2 + 1];
1197         }
1198     }
1199 }
1200
1201
1202 static void
1203 calcSparseBackProj_8u( vector<uchar*>& _ptrs, const vector<int>& _deltas,
1204                        Size imsize, const SparseMat& hist, const float** _ranges,
1205                        const double* _uniranges, float scale, bool uniform )
1206 {
1207     uchar** ptrs = &_ptrs[0];
1208     const int* deltas = &_deltas[0];
1209     int i, x, dims = hist.dims();
1210     uchar* bproj = _ptrs[dims];
1211     int bpstep = _deltas[dims*2 + 1];
1212     vector<size_t> _tab;
1213     int idx[CV_MAX_DIM];
1214
1215     calcHistLookupTables_8u( MatND(), hist, _ranges, _uniranges, uniform, true, _tab );
1216     const size_t* tab = &_tab[0];
1217
1218     for( ; imsize.height--; bproj += bpstep )
1219     {
1220         for( x = 0; x < imsize.width; x++ )
1221         {
1222             for( i = 0; i < dims; i++ )
1223             {
1224                 size_t hidx = tab[*ptrs[i] + i*256];
1225                 if( hidx >= OUT_OF_RANGE )
1226                     break;
1227                 idx[i] = (int)hidx;
1228                 ptrs[i] += deltas[i*2];
1229             }
1230
1231             if( i == dims )
1232                 bproj[x] = saturate_cast<uchar>(hist.value<float>(idx)*scale);
1233             else
1234             {
1235                 bproj[x] = 0;
1236                 for( ; i < dims; i++ )
1237                     ptrs[i] += deltas[i*2];
1238             }
1239         }
1240         for( i = 0; i < dims; i++ )
1241             ptrs[i] += deltas[i*2 + 1];
1242     }
1243 }
1244
1245
1246 void calcBackProject( const Mat* images, int nimages, const int* channels,
1247                       const SparseMat& hist, Mat& backProject,
1248                       const float** ranges, double scale, bool uniform )
1249 {
1250     vector<uchar*> ptrs;
1251     vector<int> deltas;
1252     vector<double> uniranges;
1253     Size imsize;
1254
1255     CV_Assert( hist.dims() > 0 );
1256     backProject.create( images[0].size(), images[0].depth() );
1257     histPrepareImages( images, nimages, channels, backProject,
1258                        hist.dims(), hist.hdr->size, ranges,
1259                        uniform, ptrs, deltas, imsize, uniranges );
1260     const double* _uniranges = uniform ? &uniranges[0] : 0;
1261
1262     int depth = images[0].depth();
1263     if( depth == CV_8U )
1264         calcSparseBackProj_8u(ptrs, deltas, imsize, hist, ranges, _uniranges, (float)scale, uniform);
1265     else if( depth == CV_32F )
1266         calcSparseBackProj_<float, float>(ptrs, deltas, imsize, hist, ranges, _uniranges, (float)scale, uniform );
1267 }
1268
1269
1270 ////////////////// C O M P A R E   H I S T O G R A M S ////////////////////////
1271
1272 double compareHist( const MatND& H1, const MatND& H2, int method )
1273 {
1274     NAryMatNDIterator it(H1, H2);
1275     double result = 0;
1276     int i, len;
1277
1278     CV_Assert( H1.type() == H2.type() && H1.type() == CV_32F );
1279
1280     double s1 = 0, s2 = 0, s11 = 0, s12 = 0, s22 = 0;
1281
1282     CV_Assert( it.planes[0].isContinuous() && it.planes[1].isContinuous() );
1283
1284     for( i = 0; i < it.nplanes; i++, ++it )
1285     {
1286         const float* h1 = (const float*)it.planes[0].data;
1287         const float* h2 = (const float*)it.planes[1].data;
1288         len = it.planes[0].rows*it.planes[0].cols;
1289
1290         if( method == CV_COMP_CHISQR )
1291         {
1292             for( i = 0; i < len; i++ )
1293             {
1294                 double a = h1[i] - h2[i];
1295                 double b = h1[i] + h2[i];
1296                 if( fabs(b) > FLT_EPSILON )
1297                     result += a*a/b;
1298             }
1299         }
1300         else if( method == CV_COMP_CORREL )
1301         {
1302             for( i = 0; i < len; i++ )
1303             {
1304                 double a = h1[i];
1305                 double b = h2[i];
1306
1307                 s12 += a*b;
1308                 s1 += a;
1309                 s11 += a*a;
1310                 s2 += b;
1311                 s22 += b*b;
1312             }
1313         }
1314         else if( method == CV_COMP_INTERSECT )
1315         {
1316             for( i = 0; i < len; i++ )
1317                 result += std::min(h1[i], h2[i]);
1318         }
1319         else if( method == CV_COMP_BHATTACHARYYA )
1320         {
1321             for( i = 0; i < len; i++ )
1322             {
1323                 double a = h1[i];
1324                 double b = h2[i];
1325                 result += std::sqrt(a*b);
1326                 s1 += a;
1327                 s2 += b;
1328             }
1329         }
1330         else
1331             CV_Error( CV_StsBadArg, "Unknown comparison method" );
1332     }
1333
1334     if( method == CV_COMP_CORREL )
1335     {
1336         size_t total = 1;
1337         for( i = 0; i < H1.dims; i++ )
1338             total *= H1.size[i];
1339         double scale = 1./total;
1340         double num = s12 - s1*s2*scale;
1341         double denom2 = (s11 - s1*s1*scale)*(s22 - s2*s2*scale);
1342         result = std::abs(denom2) > DBL_EPSILON ? num/std::sqrt(denom2) : 1.;
1343     }
1344     else if( method == CV_COMP_BHATTACHARYYA )
1345     {
1346         s1 *= s2;
1347         s1 = fabs(s1) > FLT_EPSILON ? 1./std::sqrt(s1) : 1.;
1348         result = std::sqrt(std::max(1. - result*s1, 0.));
1349     }
1350
1351     return result;
1352 }
1353
1354
1355 double compareHist( const SparseMat& H1, const SparseMat& H2, int method )
1356 {
1357     double result = 0;
1358     int i, dims = H1.dims();
1359
1360     CV_Assert( dims > 0 && dims == H2.dims() && H1.type() == H2.type() && H1.type() == CV_32F );
1361     for( i = 0; i < dims; i++ )
1362         CV_Assert( H1.size(i) == H2.size(i) );
1363
1364     const SparseMat *PH1 = &H1, *PH2 = &H2;
1365     if( PH1->nzcount() > PH2->nzcount() )
1366         std::swap(PH1, PH2);
1367
1368     SparseMatConstIterator it = PH1->begin();
1369     int N1 = (int)PH1->nzcount(), N2 = (int)PH2->nzcount();
1370
1371     if( method == CV_COMP_CHISQR )
1372     {
1373         for( i = 0; i < N1; i++, ++it )
1374         {
1375             float v1 = it.value<float>();
1376             const SparseMat::Node* node = it.node();
1377             float v2 = PH2->value<float>(node->idx, (size_t*)&node->hashval);
1378             if( !v2 )
1379                 result += v1;
1380             else
1381             {
1382                 double a = v1 - v2;
1383                 double b = v1 + v2;
1384                 if( fabs(b) > FLT_EPSILON )
1385                     result += a*a/b;
1386             }
1387         }
1388
1389         it = PH2->begin();
1390         for( i = 0; i < N2; i++, ++it )
1391         {
1392             float v2 = it.value<float>();
1393             const SparseMat::Node* node = it.node();
1394             if( !PH1->find<float>(node->idx, (size_t*)&node->hashval) )
1395                 result += v2;
1396         }
1397     }
1398     else if( method == CV_COMP_CORREL )
1399     {
1400         double s1 = 0, s2 = 0, s11 = 0, s12 = 0, s22 = 0;
1401
1402         for( i = 0; i < N1; i++, ++it )
1403         {
1404             double v1 = it.value<float>();
1405             const SparseMat::Node* node = it.node();
1406             s12 += v1*PH2->value<float>(node->idx, (size_t*)&node->hashval);
1407             s1 += v1;
1408             s11 += v1*v1;
1409         }
1410
1411         it = PH2->begin();
1412         for( i = 0; i < N2; i++, ++it )
1413         {
1414             double v2 = it.value<float>();
1415             s2 += v2;
1416             s22 += v2*v2;
1417         }
1418
1419         size_t total = 1;
1420         for( i = 0; i < H1.dims(); i++ )
1421             total *= H1.size(i);
1422         double scale = 1./total;
1423         double num = s12 - s1*s2*scale;
1424         double denom2 = (s11 - s1*s1*scale)*(s22 - s2*s2*scale);
1425         result = std::abs(denom2) > DBL_EPSILON ? num/std::sqrt(denom2) : 1.;
1426     }
1427     else if( method == CV_COMP_INTERSECT )
1428     {
1429         for( i = 0; i < N1; i++, ++it )
1430         {
1431             float v1 = it.value<float>();
1432             const SparseMat::Node* node = it.node();
1433             float v2 = PH2->value<float>(node->idx, (size_t*)&node->hashval);
1434             if( v2 )
1435                 result += std::min(v1, v2);
1436         }
1437     }
1438     else if( method == CV_COMP_BHATTACHARYYA )
1439     {
1440         double s1 = 0, s2 = 0;
1441
1442         for( i = 0; i < N1; i++, ++it )
1443         {
1444             double v1 = it.value<float>();
1445             const SparseMat::Node* node = it.node();
1446             double v2 = PH2->value<float>(node->idx, (size_t*)&node->hashval);
1447             result += std::sqrt(v1*v2);
1448             s1 += v1;
1449         }
1450
1451         it = PH2->begin();
1452         for( i = 0; i < N2; i++, ++it )
1453             s2 += it.value<float>();
1454
1455         s1 *= s2;
1456         s1 = fabs(s1) > FLT_EPSILON ? 1./std::sqrt(s1) : 1.;
1457         result = std::sqrt(std::max(1. - result*s1, 0.));
1458     }
1459     else
1460         CV_Error( CV_StsBadArg, "Unknown comparison method" );
1461
1462     return result;
1463 }
1464
1465 }
1466
1467
1468 const int CV_HIST_DEFAULT_TYPE = CV_32F;
1469
1470 /* Creates new histogram */
1471 CvHistogram *
1472 cvCreateHist( int dims, int *sizes, CvHistType type, float** ranges, int uniform )
1473 {
1474     CvHistogram *hist = 0;
1475
1476     if( (unsigned)dims > CV_MAX_DIM )
1477         CV_Error( CV_BadOrder, "Number of dimensions is out of range" );
1478
1479     if( !sizes )
1480         CV_Error( CV_HeaderIsNull, "Null <sizes> pointer" );
1481
1482     hist = (CvHistogram *)cvAlloc( sizeof( CvHistogram ));
1483     hist->type = CV_HIST_MAGIC_VAL + ((int)type & 1);
1484         if (uniform) hist->type|= CV_HIST_UNIFORM_FLAG;
1485     hist->thresh2 = 0;
1486     hist->bins = 0;
1487     if( type == CV_HIST_ARRAY )
1488     {
1489         hist->bins = cvInitMatNDHeader( &hist->mat, dims, sizes,
1490                                         CV_HIST_DEFAULT_TYPE );
1491         cvCreateData( hist->bins );
1492     }
1493     else if( type == CV_HIST_SPARSE )
1494         hist->bins = cvCreateSparseMat( dims, sizes, CV_HIST_DEFAULT_TYPE );
1495     else
1496         CV_Error( CV_StsBadArg, "Invalid histogram type" );
1497
1498     if( ranges )
1499         cvSetHistBinRanges( hist, ranges, uniform );
1500
1501     return hist;
1502 }
1503
1504
1505 /* Creates histogram wrapping header for given array */
1506 CV_IMPL CvHistogram*
1507 cvMakeHistHeaderForArray( int dims, int *sizes, CvHistogram *hist,
1508                           float *data, float **ranges, int uniform )
1509 {
1510     if( !hist )
1511         CV_Error( CV_StsNullPtr, "Null histogram header pointer" );
1512
1513     if( !data )
1514         CV_Error( CV_StsNullPtr, "Null data pointer" );
1515
1516     hist->thresh2 = 0;
1517     hist->type = CV_HIST_MAGIC_VAL;
1518     hist->bins = cvInitMatNDHeader( &hist->mat, dims, sizes, CV_HIST_DEFAULT_TYPE, data );
1519
1520     if( ranges )
1521     {
1522         if( !uniform )
1523             CV_Error( CV_StsBadArg, "Only uniform bin ranges can be used here "
1524                                     "(to avoid memory allocation)" );
1525         cvSetHistBinRanges( hist, ranges, uniform );
1526     }
1527
1528     return hist;
1529 }
1530
1531
1532 CV_IMPL void
1533 cvReleaseHist( CvHistogram **hist )
1534 {
1535     if( !hist )
1536         CV_Error( CV_StsNullPtr, "" );
1537
1538     if( *hist )
1539     {
1540         CvHistogram* temp = *hist;
1541
1542         if( !CV_IS_HIST(temp))
1543             CV_Error( CV_StsBadArg, "Invalid histogram header" );
1544         *hist = 0;
1545
1546         if( CV_IS_SPARSE_HIST( temp ))
1547             cvReleaseSparseMat( (CvSparseMat**)&temp->bins );
1548         else
1549         {
1550             cvReleaseData( temp->bins );
1551             temp->bins = 0;
1552         }
1553
1554         if( temp->thresh2 )
1555             cvFree( &temp->thresh2 );
1556         cvFree( &temp );
1557     }
1558 }
1559
1560 CV_IMPL void
1561 cvClearHist( CvHistogram *hist )
1562 {
1563     if( !CV_IS_HIST(hist) )
1564         CV_Error( CV_StsBadArg, "Invalid histogram header" );
1565     cvZero( hist->bins );
1566 }
1567
1568
1569 // Clears histogram bins that are below than threshold
1570 CV_IMPL void
1571 cvThreshHist( CvHistogram* hist, double thresh )
1572 {
1573     if( !CV_IS_HIST(hist) )
1574         CV_Error( CV_StsBadArg, "Invalid histogram header" );
1575
1576     if( !CV_IS_SPARSE_MAT(hist->bins) )
1577     {
1578         CvMat mat;
1579         cvGetMat( hist->bins, &mat, 0, 1 );
1580         cvThreshold( &mat, &mat, thresh, 0, CV_THRESH_TOZERO );
1581     }
1582     else
1583     {
1584         CvSparseMat* mat = (CvSparseMat*)hist->bins;
1585         CvSparseMatIterator iterator;
1586         CvSparseNode *node;
1587
1588         for( node = cvInitSparseMatIterator( mat, &iterator );
1589              node != 0; node = cvGetNextSparseNode( &iterator ))
1590         {
1591             float* val = (float*)CV_NODE_VAL( mat, node );
1592             if( *val <= thresh )
1593                 *val = 0;
1594         }
1595     }
1596 }
1597
1598
1599 // Normalizes histogram (make sum of the histogram bins == factor)
1600 CV_IMPL void
1601 cvNormalizeHist( CvHistogram* hist, double factor )
1602 {
1603     double sum = 0;
1604
1605     if( !CV_IS_HIST(hist) )
1606         CV_Error( CV_StsBadArg, "Invalid histogram header" );
1607
1608     if( !CV_IS_SPARSE_HIST(hist) )
1609     {
1610         CvMat mat;
1611         cvGetMat( hist->bins, &mat, 0, 1 );
1612         sum = cvSum( &mat ).val[0];
1613         if( fabs(sum) < DBL_EPSILON )
1614             sum = 1;
1615         cvScale( &mat, &mat, factor/sum, 0 );
1616     }
1617     else
1618     {
1619         CvSparseMat* mat = (CvSparseMat*)hist->bins;
1620         CvSparseMatIterator iterator;
1621         CvSparseNode *node;
1622         float scale;
1623
1624         for( node = cvInitSparseMatIterator( mat, &iterator );
1625              node != 0; node = cvGetNextSparseNode( &iterator ))
1626         {
1627             sum += *(float*)CV_NODE_VAL(mat,node);
1628         }
1629
1630         if( fabs(sum) < DBL_EPSILON )
1631             sum = 1;
1632         scale = (float)(factor/sum);
1633
1634         for( node = cvInitSparseMatIterator( mat, &iterator );
1635              node != 0; node = cvGetNextSparseNode( &iterator ))
1636         {
1637             *(float*)CV_NODE_VAL(mat,node) *= scale;
1638         }
1639     }
1640 }
1641
1642
1643 // Retrieves histogram global min, max and their positions
1644 CV_IMPL void
1645 cvGetMinMaxHistValue( const CvHistogram* hist,
1646                       float *value_min, float* value_max,
1647                       int* idx_min, int* idx_max )
1648 {
1649     double minVal, maxVal;
1650     int i, dims, size[CV_MAX_DIM];
1651
1652     if( !CV_IS_HIST(hist) )
1653         CV_Error( CV_StsBadArg, "Invalid histogram header" );
1654
1655     dims = cvGetDims( hist->bins, size );
1656
1657     if( !CV_IS_SPARSE_HIST(hist) )
1658     {
1659         CvMat mat;
1660         CvPoint minPt, maxPt;
1661
1662         cvGetMat( hist->bins, &mat, 0, 1 );
1663         cvMinMaxLoc( &mat, &minVal, &maxVal, &minPt, &maxPt );
1664
1665         if( dims == 1 )
1666         {
1667             if( idx_min )
1668                 *idx_min = minPt.y + minPt.x;
1669             if( idx_max )
1670                 *idx_max = maxPt.y + maxPt.x;
1671         }
1672         else if( dims == 2 )
1673         {
1674             if( idx_min )
1675                 idx_min[0] = minPt.y, idx_min[1] = minPt.x;
1676             if( idx_max )
1677                 idx_max[0] = maxPt.y, idx_max[1] = maxPt.x;
1678         }
1679         else if( idx_min || idx_max )
1680         {
1681             int imin = minPt.y*mat.cols + minPt.x;
1682             int imax = maxPt.y*mat.cols + maxPt.x;
1683             int i;
1684
1685             for( i = dims - 1; i >= 0; i-- )
1686             {
1687                 if( idx_min )
1688                 {
1689                     int t = imin / size[i];
1690                     idx_min[i] = imin - t*size[i];
1691                     imin = t;
1692                 }
1693
1694                 if( idx_max )
1695                 {
1696                     int t = imax / size[i];
1697                     idx_max[i] = imax - t*size[i];
1698                     imax = t;
1699                 }
1700             }
1701         }
1702     }
1703     else
1704     {
1705         CvSparseMat* mat = (CvSparseMat*)hist->bins;
1706         CvSparseMatIterator iterator;
1707         CvSparseNode *node;
1708         int minv = INT_MAX;
1709         int maxv = INT_MIN;
1710         CvSparseNode* minNode = 0;
1711         CvSparseNode* maxNode = 0;
1712         const int *_idx_min = 0, *_idx_max = 0;
1713         Cv32suf m;
1714
1715         for( node = cvInitSparseMatIterator( mat, &iterator );
1716              node != 0; node = cvGetNextSparseNode( &iterator ))
1717         {
1718             int value = *(int*)CV_NODE_VAL(mat,node);
1719             value = CV_TOGGLE_FLT(value);
1720             if( value < minv )
1721             {
1722                 minv = value;
1723                 minNode = node;
1724             }
1725
1726             if( value > maxv )
1727             {
1728                 maxv = value;
1729                 maxNode = node;
1730             }
1731         }
1732
1733         if( minNode )
1734         {
1735             _idx_min = CV_NODE_IDX(mat,minNode);
1736             _idx_max = CV_NODE_IDX(mat,maxNode);
1737             m.i = CV_TOGGLE_FLT(minv); minVal = m.f;
1738             m.i = CV_TOGGLE_FLT(maxv); maxVal = m.f;
1739         }
1740         else
1741         {
1742             minVal = maxVal = 0;
1743         }
1744
1745         for( i = 0; i < dims; i++ )
1746         {
1747             if( idx_min )
1748                 idx_min[i] = _idx_min ? _idx_min[i] : -1;
1749             if( idx_max )
1750                 idx_max[i] = _idx_max ? _idx_max[i] : -1;
1751         }
1752     }
1753
1754     if( value_min )
1755         *value_min = (float)minVal;
1756
1757     if( value_max )
1758         *value_max = (float)maxVal;
1759 }
1760
1761
1762 // Compares two histograms using one of a few methods
1763 CV_IMPL double
1764 cvCompareHist( const CvHistogram* hist1,
1765                const CvHistogram* hist2,
1766                int method )
1767 {
1768     int i;
1769     int size1[CV_MAX_DIM], size2[CV_MAX_DIM], total = 1;
1770
1771     if( !CV_IS_HIST(hist1) || !CV_IS_HIST(hist2) )
1772         CV_Error( CV_StsBadArg, "Invalid histogram header[s]" );
1773
1774     if( CV_IS_SPARSE_MAT(hist1->bins) != CV_IS_SPARSE_MAT(hist2->bins))
1775         CV_Error(CV_StsUnmatchedFormats, "One of histograms is sparse and other is not");
1776
1777     if( !CV_IS_SPARSE_MAT(hist1->bins) )
1778     {
1779         cv::MatND H1((const CvMatND*)hist1->bins), H2((const CvMatND*)hist2->bins);
1780         return cv::compareHist(H1, H2, method);
1781     }
1782
1783     int dims1 = cvGetDims( hist1->bins, size1 );
1784     int dims2 = cvGetDims( hist2->bins, size2 );
1785
1786     if( dims1 != dims2 )
1787         CV_Error( CV_StsUnmatchedSizes,
1788                  "The histograms have different numbers of dimensions" );
1789
1790     for( i = 0; i < dims1; i++ )
1791     {
1792         if( size1[i] != size2[i] )
1793             CV_Error( CV_StsUnmatchedSizes, "The histograms have different sizes" );
1794         total *= size1[i];
1795     }
1796
1797     double result = 0;
1798     CvSparseMat* mat1 = (CvSparseMat*)(hist1->bins);
1799     CvSparseMat* mat2 = (CvSparseMat*)(hist2->bins);
1800     CvSparseMatIterator iterator;
1801     CvSparseNode *node1, *node2;
1802
1803     if( mat1->heap->active_count > mat2->heap->active_count )
1804     {
1805         CvSparseMat* t;
1806         CV_SWAP( mat1, mat2, t );
1807     }
1808
1809     if( method == CV_COMP_CHISQR )
1810     {
1811         for( node1 = cvInitSparseMatIterator( mat1, &iterator );
1812              node1 != 0; node1 = cvGetNextSparseNode( &iterator ))
1813         {
1814             double v1 = *(float*)CV_NODE_VAL(mat1,node1);
1815             uchar* node2_data = cvPtrND( mat2, CV_NODE_IDX(mat1,node1), 0, 0, &node1->hashval );
1816             if( !node2_data )
1817                 result += v1;
1818             else
1819             {
1820                 double v2 = *(float*)node2_data;
1821                 double a = v1 - v2;
1822                 double b = v1 + v2;
1823                 if( fabs(b) > DBL_EPSILON )
1824                     result += a*a/b;
1825             }
1826         }
1827
1828         for( node2 = cvInitSparseMatIterator( mat2, &iterator );
1829              node2 != 0; node2 = cvGetNextSparseNode( &iterator ))
1830         {
1831             double v2 = *(float*)CV_NODE_VAL(mat2,node2);
1832             if( !cvPtrND( mat1, CV_NODE_IDX(mat2,node2), 0, 0, &node2->hashval ))
1833                 result += v2;
1834         }
1835     }
1836     else if( method == CV_COMP_CORREL )
1837     {
1838         double s1 = 0, s11 = 0;
1839         double s2 = 0, s22 = 0;
1840         double s12 = 0;
1841         double num, denom2, scale = 1./total;
1842
1843         for( node1 = cvInitSparseMatIterator( mat1, &iterator );
1844              node1 != 0; node1 = cvGetNextSparseNode( &iterator ))
1845         {
1846             double v1 = *(float*)CV_NODE_VAL(mat1,node1);
1847             uchar* node2_data = cvPtrND( mat2, CV_NODE_IDX(mat1,node1),
1848                                         0, 0, &node1->hashval );
1849             if( node2_data )
1850             {
1851                 double v2 = *(float*)node2_data;
1852                 s12 += v1*v2;
1853             }
1854             s1 += v1;
1855             s11 += v1*v1;
1856         }
1857
1858         for( node2 = cvInitSparseMatIterator( mat2, &iterator );
1859              node2 != 0; node2 = cvGetNextSparseNode( &iterator ))
1860         {
1861             double v2 = *(float*)CV_NODE_VAL(mat2,node2);
1862             s2 += v2;
1863             s22 += v2*v2;
1864         }
1865
1866         num = s12 - s1*s2*scale;
1867         denom2 = (s11 - s1*s1*scale)*(s22 - s2*s2*scale);
1868         result = fabs(denom2) > DBL_EPSILON ? num/sqrt(denom2) : 1;
1869     }
1870     else if( method == CV_COMP_INTERSECT )
1871     {
1872         for( node1 = cvInitSparseMatIterator( mat1, &iterator );
1873              node1 != 0; node1 = cvGetNextSparseNode( &iterator ))
1874         {
1875             float v1 = *(float*)CV_NODE_VAL(mat1,node1);
1876             uchar* node2_data = cvPtrND( mat2, CV_NODE_IDX(mat1,node1),
1877                                          0, 0, &node1->hashval );
1878             if( node2_data )
1879             {
1880                 float v2 = *(float*)node2_data;
1881                 if( v1 <= v2 )
1882                     result += v1;
1883                 else
1884                     result += v2;
1885             }
1886         }
1887     }
1888     else if( method == CV_COMP_BHATTACHARYYA )
1889     {
1890         double s1 = 0, s2 = 0;
1891
1892         for( node1 = cvInitSparseMatIterator( mat1, &iterator );
1893              node1 != 0; node1 = cvGetNextSparseNode( &iterator ))
1894         {
1895             double v1 = *(float*)CV_NODE_VAL(mat1,node1);
1896             uchar* node2_data = cvPtrND( mat2, CV_NODE_IDX(mat1,node1),
1897                                          0, 0, &node1->hashval );
1898             s1 += v1;
1899             if( node2_data )
1900             {
1901                 double v2 = *(float*)node2_data;
1902                 result += sqrt(v1 * v2);
1903             }
1904         }
1905
1906         for( node1 = cvInitSparseMatIterator( mat2, &iterator );
1907              node1 != 0; node1 = cvGetNextSparseNode( &iterator ))
1908         {
1909             double v2 = *(float*)CV_NODE_VAL(mat2,node1);
1910             s2 += v2;
1911         }
1912
1913         s1 *= s2;
1914         s1 = fabs(s1) > FLT_EPSILON ? 1./sqrt(s1) : 1.;
1915         result = 1. - result*s1;
1916         result = sqrt(MAX(result,0.));
1917     }
1918     else
1919         CV_Error( CV_StsBadArg, "Unknown comparison method" );
1920
1921     return result;
1922 }
1923
1924 // copies one histogram to another
1925 CV_IMPL void
1926 cvCopyHist( const CvHistogram* src, CvHistogram** _dst )
1927 {
1928     int eq = 0;
1929     int is_sparse;
1930     int i, dims1, dims2;
1931     int size1[CV_MAX_DIM], size2[CV_MAX_DIM], total = 1;
1932     float* ranges[CV_MAX_DIM];
1933     float** thresh = 0;
1934     CvHistogram* dst;
1935
1936     if( !_dst )
1937         CV_Error( CV_StsNullPtr, "Destination double pointer is NULL" );
1938
1939     dst = *_dst;
1940
1941     if( !CV_IS_HIST(src) || (dst && !CV_IS_HIST(dst)) )
1942         CV_Error( CV_StsBadArg, "Invalid histogram header[s]" );
1943
1944     is_sparse = CV_IS_SPARSE_MAT(src->bins);
1945     dims1 = cvGetDims( src->bins, size1 );
1946     for( i = 0; i < dims1; i++ )
1947         total *= size1[i];
1948
1949     if( dst && is_sparse == CV_IS_SPARSE_MAT(dst->bins))
1950     {
1951         dims2 = cvGetDims( dst->bins, size2 );
1952
1953         if( dims1 == dims2 )
1954         {
1955             for( i = 0; i < dims1; i++ )
1956                 if( size1[i] != size2[i] )
1957                     break;
1958         }
1959
1960         eq = i == dims1;
1961     }
1962
1963     if( !eq )
1964     {
1965         cvReleaseHist( _dst );
1966         dst = cvCreateHist( dims1, size1, !is_sparse ? CV_HIST_ARRAY : CV_HIST_SPARSE, 0, 0 );
1967         *_dst = dst;
1968     }
1969
1970     if( CV_HIST_HAS_RANGES( src ))
1971     {
1972         if( CV_IS_UNIFORM_HIST( src ))
1973         {
1974             for( i = 0; i < dims1; i++ )
1975                 ranges[i] = (float*)src->thresh[i];
1976             thresh = ranges;
1977         }
1978         else
1979             thresh = src->thresh2;
1980         cvSetHistBinRanges( dst, thresh, CV_IS_UNIFORM_HIST(src));
1981     }
1982
1983     cvCopy( src->bins, dst->bins );
1984 }
1985
1986
1987 // Sets a value range for every histogram bin
1988 CV_IMPL void
1989 cvSetHistBinRanges( CvHistogram* hist, float** ranges, int uniform )
1990 {
1991     int dims, size[CV_MAX_DIM], total = 0;
1992     int i, j;
1993
1994     if( !ranges )
1995         CV_Error( CV_StsNullPtr, "NULL ranges pointer" );
1996
1997     if( !CV_IS_HIST(hist) )
1998         CV_Error( CV_StsBadArg, "Invalid histogram header" );
1999
2000     dims = cvGetDims( hist->bins, size );
2001     for( i = 0; i < dims; i++ )
2002         total += size[i]+1;
2003
2004     if( uniform )
2005     {
2006         for( i = 0; i < dims; i++ )
2007         {
2008             if( !ranges[i] )
2009                 CV_Error( CV_StsNullPtr, "One of <ranges> elements is NULL" );
2010             hist->thresh[i][0] = ranges[i][0];
2011             hist->thresh[i][1] = ranges[i][1];
2012         }
2013
2014         hist->type |= CV_HIST_UNIFORM_FLAG + CV_HIST_RANGES_FLAG;
2015     }
2016     else
2017     {
2018         float* dim_ranges;
2019
2020         if( !hist->thresh2 )
2021         {
2022             hist->thresh2 = (float**)cvAlloc(
2023                         dims*sizeof(hist->thresh2[0])+
2024                         total*sizeof(hist->thresh2[0][0]));
2025         }
2026         dim_ranges = (float*)(hist->thresh2 + dims);
2027
2028         for( i = 0; i < dims; i++ )
2029         {
2030             float val0 = -FLT_MAX;
2031
2032             if( !ranges[i] )
2033                 CV_Error( CV_StsNullPtr, "One of <ranges> elements is NULL" );
2034
2035             for( j = 0; j <= size[i]; j++ )
2036             {
2037                 float val = ranges[i][j];
2038                 if( val <= val0 )
2039                     CV_Error(CV_StsOutOfRange, "Bin ranges should go in ascenting order");
2040                 val0 = dim_ranges[j] = val;
2041             }
2042
2043             hist->thresh2[i] = dim_ranges;
2044             dim_ranges += size[i] + 1;
2045         }
2046
2047         hist->type |= CV_HIST_RANGES_FLAG;
2048         hist->type &= ~CV_HIST_UNIFORM_FLAG;
2049     }
2050 }
2051
2052
2053 CV_IMPL void
2054 cvCalcArrHist( CvArr** img, CvHistogram* hist, int accumulate, const CvArr* mask )
2055 {
2056     if( !CV_IS_HIST(hist))
2057         CV_Error( CV_StsBadArg, "Bad histogram pointer" );
2058
2059     if( !img )
2060         CV_Error( CV_StsNullPtr, "Null double array pointer" );
2061
2062     int size[CV_MAX_DIM];
2063     int i, dims = cvGetDims( hist->bins, size);
2064     bool uniform = CV_IS_UNIFORM_HIST(hist);
2065
2066     cv::vector<cv::Mat> images(dims);
2067     for( i = 0; i < dims; i++ )
2068         images[i] = cv::cvarrToMat(img[i]);
2069
2070     cv::Mat _mask;
2071     if( mask )
2072         _mask = cv::cvarrToMat(mask);
2073
2074     const float* uranges[CV_MAX_DIM] = {0};
2075     const float** ranges = 0;
2076
2077     if( hist->type & CV_HIST_RANGES_FLAG )
2078     {
2079         ranges = (const float**)hist->thresh2;
2080         if( uniform )
2081         {
2082             for( i = 0; i < dims; i++ )
2083                 uranges[i] = &hist->thresh[i][0];
2084             ranges = uranges;
2085         }
2086     }
2087
2088     if( !CV_IS_SPARSE_HIST(hist) )
2089     {
2090         cv::MatND H((const CvMatND*)hist->bins);
2091         cv::calcHist( &images[0], (int)images.size(), 0, _mask,
2092                       H, H.dims, H.size, ranges, uniform, accumulate != 0 );
2093     }
2094     else
2095     {
2096         CvSparseMat* sparsemat = (CvSparseMat*)hist->bins;
2097
2098         if( !accumulate )
2099             cvZero( hist->bins );
2100         cv::SparseMat sH(sparsemat);
2101         cv::calcHist( &images[0], (int)images.size(), 0, _mask, sH, sH.dims(),
2102                       sH.dims() > 0 ? sH.hdr->size : 0, ranges, uniform, accumulate != 0, true );
2103
2104         if( accumulate )
2105             cvZero( sparsemat );
2106
2107         cv::SparseMatConstIterator it = sH.begin();
2108         int nz = (int)sH.nzcount();
2109         for( i = 0; i < nz; i++, ++it )
2110             *(float*)cvPtrND(sparsemat, it.node()->idx, 0, -2) = (float)*(const int*)it.ptr;
2111     }
2112 }
2113
2114
2115 CV_IMPL void
2116 cvCalcArrBackProject( CvArr** img, CvArr* dst, const CvHistogram* hist )
2117 {
2118     if( !CV_IS_HIST(hist))
2119         CV_Error( CV_StsBadArg, "Bad histogram pointer" );
2120
2121     if( !img )
2122         CV_Error( CV_StsNullPtr, "Null double array pointer" );
2123
2124     int size[CV_MAX_DIM];
2125     int i, dims = cvGetDims( hist->bins, size );
2126
2127     bool uniform = CV_IS_UNIFORM_HIST(hist);
2128     const float* uranges[CV_MAX_DIM] = {0};
2129     const float** ranges = 0;
2130
2131     if( hist->type & CV_HIST_RANGES_FLAG )
2132     {
2133         ranges = (const float**)hist->thresh2;
2134         if( uniform )
2135         {
2136             for( i = 0; i < dims; i++ )
2137                 uranges[i] = &hist->thresh[i][0];
2138             ranges = uranges;
2139         }
2140     }
2141
2142     cv::vector<cv::Mat> images(dims);
2143     for( i = 0; i < dims; i++ )
2144         images[i] = cv::cvarrToMat(img[i]);
2145
2146     cv::Mat _dst = cv::cvarrToMat(dst);
2147
2148     CV_Assert( _dst.size() == images[0].size() && _dst.depth() == images[0].depth() );
2149
2150     if( !CV_IS_SPARSE_HIST(hist) )
2151     {
2152         cv::MatND H((const CvMatND*)hist->bins);
2153         cv::calcBackProject( &images[0], (int)images.size(),
2154                             0, H, _dst, ranges, 1, uniform );
2155     }
2156     else
2157     {
2158         cv::SparseMat sH((const CvSparseMat*)hist->bins);
2159         cv::calcBackProject( &images[0], (int)images.size(),
2160                              0, sH, _dst, ranges, 1, uniform );
2161     }
2162 }
2163
2164
2165 ////////////////////// B A C K   P R O J E C T   P A T C H /////////////////////////
2166
2167 CV_IMPL void
2168 cvCalcArrBackProjectPatch( CvArr** arr, CvArr* dst, CvSize patch_size, CvHistogram* hist,
2169                            int method, double norm_factor )
2170 {
2171     CvHistogram* model = 0;
2172
2173     IplImage imgstub[CV_MAX_DIM], *img[CV_MAX_DIM];
2174     IplROI roi;
2175     CvMat dststub, *dstmat;
2176     int i, dims;
2177     int x, y;
2178     CvSize size;
2179
2180     if( !CV_IS_HIST(hist))
2181         CV_Error( CV_StsBadArg, "Bad histogram pointer" );
2182
2183     if( !arr )
2184         CV_Error( CV_StsNullPtr, "Null double array pointer" );
2185
2186     if( norm_factor <= 0 )
2187         CV_Error( CV_StsOutOfRange,
2188                   "Bad normalization factor (set it to 1.0 if unsure)" );
2189
2190     if( patch_size.width <= 0 || patch_size.height <= 0 )
2191         CV_Error( CV_StsBadSize, "The patch width and height must be positive" );
2192
2193     dims = cvGetDims( hist->bins );
2194     cvNormalizeHist( hist, norm_factor );
2195
2196     for( i = 0; i < dims; i++ )
2197     {
2198         CvMat stub, *mat;
2199         mat = cvGetMat( arr[i], &stub, 0, 0 );
2200         img[i] = cvGetImage( mat, &imgstub[i] );
2201         img[i]->roi = &roi;
2202     }
2203
2204     dstmat = cvGetMat( dst, &dststub, 0, 0 );
2205     if( CV_MAT_TYPE( dstmat->type ) != CV_32FC1 )
2206         CV_Error( CV_StsUnsupportedFormat, "Resultant image must have 32fC1 type" );
2207
2208     if( dstmat->cols != img[0]->width - patch_size.width + 1 ||
2209         dstmat->rows != img[0]->height - patch_size.height + 1 )
2210         CV_Error( CV_StsUnmatchedSizes,
2211             "The output map must be (W-w+1 x H-h+1), "
2212             "where the input images are (W x H) each and the patch is (w x h)" );
2213
2214     cvCopyHist( hist, &model );
2215
2216     size = cvGetMatSize(dstmat);
2217     roi.coi = 0;
2218     roi.width = patch_size.width;
2219     roi.height = patch_size.height;
2220
2221     for( y = 0; y < size.height; y++ )
2222     {
2223         for( x = 0; x < size.width; x++ )
2224         {
2225             double result;
2226             roi.xOffset = x;
2227             roi.yOffset = y;
2228
2229             cvCalcHist( img, model );
2230             cvNormalizeHist( model, norm_factor );
2231             result = cvCompareHist( model, hist, method );
2232             CV_MAT_ELEM( *dstmat, float, y, x ) = (float)result;
2233         }
2234     }
2235
2236     cvReleaseHist( &model );
2237 }
2238
2239
2240 // Calculates Bayes probabilistic histograms
2241 CV_IMPL void
2242 cvCalcBayesianProb( CvHistogram** src, int count, CvHistogram** dst )
2243 {
2244     int i;
2245
2246     if( !src || !dst )
2247         CV_Error( CV_StsNullPtr, "NULL histogram array pointer" );
2248
2249     if( count < 2 )
2250         CV_Error( CV_StsOutOfRange, "Too small number of histograms" );
2251
2252     for( i = 0; i < count; i++ )
2253     {
2254         if( !CV_IS_HIST(src[i]) || !CV_IS_HIST(dst[i]) )
2255             CV_Error( CV_StsBadArg, "Invalid histogram header" );
2256
2257         if( !CV_IS_MATND(src[i]->bins) || !CV_IS_MATND(dst[i]->bins) )
2258             CV_Error( CV_StsBadArg, "The function supports dense histograms only" );
2259     }
2260
2261     cvZero( dst[0]->bins );
2262     // dst[0] = src[0] + ... + src[count-1]
2263     for( i = 0; i < count; i++ )
2264         cvAdd( src[i]->bins, dst[0]->bins, dst[0]->bins );
2265
2266     cvDiv( 0, dst[0]->bins, dst[0]->bins );
2267
2268     // dst[i] = src[i]*(1/dst[0])
2269     for( i = count - 1; i >= 0; i-- )
2270         cvMul( src[i]->bins, dst[0]->bins, dst[i]->bins );
2271 }
2272
2273
2274 CV_IMPL void
2275 cvCalcProbDensity( const CvHistogram* hist, const CvHistogram* hist_mask,
2276                    CvHistogram* hist_dens, double scale )
2277 {
2278     if( scale <= 0 )
2279         CV_Error( CV_StsOutOfRange, "scale must be positive" );
2280
2281     if( !CV_IS_HIST(hist) || !CV_IS_HIST(hist_mask) || !CV_IS_HIST(hist_dens) )
2282         CV_Error( CV_StsBadArg, "Invalid histogram pointer[s]" );
2283
2284     {
2285         CvArr* arrs[] = { hist->bins, hist_mask->bins, hist_dens->bins };
2286         CvMatND stubs[3];
2287         CvNArrayIterator iterator;
2288
2289         cvInitNArrayIterator( 3, arrs, 0, stubs, &iterator );
2290
2291         if( CV_MAT_TYPE(iterator.hdr[0]->type) != CV_32FC1 )
2292             CV_Error( CV_StsUnsupportedFormat, "All histograms must have 32fC1 type" );
2293
2294         do
2295         {
2296             const float* srcdata = (const float*)(iterator.ptr[0]);
2297             const float* maskdata = (const float*)(iterator.ptr[1]);
2298             float* dstdata = (float*)(iterator.ptr[2]);
2299             int i;
2300
2301             for( i = 0; i < iterator.size.width; i++ )
2302             {
2303                 float s = srcdata[i];
2304                 float m = maskdata[i];
2305                 if( s > FLT_EPSILON )
2306                     if( m <= s )
2307                         dstdata[i] = (float)(m*scale/s);
2308                     else
2309                         dstdata[i] = (float)scale;
2310                 else
2311                     dstdata[i] = (float)0;
2312             }
2313         }
2314         while( cvNextNArraySlice( &iterator ));
2315     }
2316 }
2317
2318
2319 CV_IMPL void cvEqualizeHist( const CvArr* srcarr, CvArr* dstarr )
2320 {
2321     CvMat sstub, *src = cvGetMat(srcarr, &sstub);
2322     CvMat dstub, *dst = cvGetMat(dstarr, &dstub);
2323
2324     CV_Assert( CV_ARE_SIZES_EQ(src, dst) && CV_ARE_TYPES_EQ(src, dst) &&
2325                CV_MAT_TYPE(src->type) == CV_8UC1 );
2326     CvSize size = cvGetMatSize(src);
2327     if( CV_IS_MAT_CONT(src->type & dst->type) )
2328     {
2329         size.width *= size.height;
2330         size.height = 1;
2331     }
2332     int x, y;
2333     const int hist_sz = 256;
2334     int hist[hist_sz];
2335     memset(hist, 0, sizeof(hist));
2336
2337     for( y = 0; y < size.height; y++ )
2338     {
2339         const uchar* sptr = src->data.ptr + src->step*y;
2340         for( x = 0; x < size.width; x++ )
2341             hist[sptr[x]]++;
2342     }
2343
2344     float scale = 255.f/(size.width*size.height);
2345     int sum = 0;
2346     uchar lut[hist_sz+1];
2347
2348     for( int i = 0; i < hist_sz; i++ )
2349     {
2350         sum += hist[i];
2351         int val = cvRound(sum*scale);
2352         lut[i] = CV_CAST_8U(val);
2353     }
2354
2355     lut[0] = 0;
2356     for( y = 0; y < size.height; y++ )
2357     {
2358         const uchar* sptr = src->data.ptr + src->step*y;
2359         uchar* dptr = dst->data.ptr + dst->step*y;
2360         for( x = 0; x < size.width; x++ )
2361             dptr[x] = lut[sptr[x]];
2362     }
2363 }
2364
2365
2366 void cv::equalizeHist( const Mat& src, Mat& dst )
2367 {
2368     dst.create( src.size(), src.type() );
2369     CvMat _src = src, _dst = dst;
2370     cvEqualizeHist( &_src, &_dst );
2371 }
2372
2373 /* Implementation of RTTI and Generic Functions for CvHistogram */
2374 #define CV_TYPE_NAME_HIST "opencv-hist"
2375
2376 static int icvIsHist( const void * ptr )
2377 {
2378     return CV_IS_HIST( ((CvHistogram*)ptr) );
2379 }
2380
2381 static CvHistogram * icvCloneHist( const CvHistogram * src )
2382 {
2383     CvHistogram * dst=NULL;
2384     cvCopyHist(src, &dst);
2385     return dst;
2386 }
2387
2388 static void *icvReadHist( CvFileStorage * fs, CvFileNode * node )
2389 {
2390     CvHistogram * h = 0;
2391     int type = 0;
2392     int is_uniform = 0;
2393     int have_ranges = 0;
2394
2395     h = (CvHistogram *)cvAlloc( sizeof(CvHistogram) );
2396
2397     type = cvReadIntByName( fs, node, "type", 0 );
2398     is_uniform = cvReadIntByName( fs, node, "is_uniform", 0 );
2399     have_ranges = cvReadIntByName( fs, node, "have_ranges", 0 );
2400     h->type = CV_HIST_MAGIC_VAL | type |
2401         (is_uniform ? CV_HIST_UNIFORM_FLAG : 0) |
2402         (have_ranges ? CV_HIST_RANGES_FLAG : 0);
2403
2404     if(type == CV_HIST_ARRAY)
2405     {
2406         // read histogram bins
2407         CvMatND* mat = (CvMatND*)cvReadByName( fs, node, "mat" );
2408         int i, sizes[CV_MAX_DIM];
2409
2410         if(!CV_IS_MATND(mat))
2411             CV_Error( CV_StsError, "Expected CvMatND");
2412
2413         for(i=0; i<mat->dims; i++)
2414             sizes[i] = mat->dim[i].size;
2415
2416         cvInitMatNDHeader( &(h->mat), mat->dims, sizes, mat->type, mat->data.ptr );
2417         h->bins = &(h->mat);
2418
2419         // take ownership of refcount pointer as well
2420         h->mat.refcount = mat->refcount;
2421
2422         // increase refcount so freeing temp header doesn't free data
2423         cvIncRefData( mat );
2424
2425         // free temporary header
2426         cvReleaseMatND( &mat );
2427     }
2428     else
2429     {
2430         h->bins = cvReadByName( fs, node, "bins" );
2431         if(!CV_IS_SPARSE_MAT(h->bins)){
2432             CV_Error( CV_StsError, "Unknown Histogram type");
2433         }
2434     }
2435
2436     // read thresholds
2437     if(have_ranges)
2438     {
2439         int i, dims, size[CV_MAX_DIM], total = 0;
2440         CvSeqReader reader;
2441         CvFileNode * thresh_node;
2442
2443         dims = cvGetDims( h->bins, size );
2444         for( i = 0; i < dims; i++ )
2445             total += size[i]+1;
2446
2447         thresh_node = cvGetFileNodeByName( fs, node, "thresh" );
2448         if(!thresh_node)
2449             CV_Error( CV_StsError, "'thresh' node is missing");
2450         cvStartReadRawData( fs, thresh_node, &reader );
2451
2452         if(is_uniform)
2453         {
2454             for(i=0; i<dims; i++)
2455                 cvReadRawDataSlice( fs, &reader, 2, h->thresh[i], "f" );
2456             h->thresh2 = NULL;
2457         }
2458         else
2459         {
2460             float* dim_ranges;
2461             h->thresh2 = (float**)cvAlloc(
2462                 dims*sizeof(h->thresh2[0])+
2463                 total*sizeof(h->thresh2[0][0]));
2464             dim_ranges = (float*)(h->thresh2 + dims);
2465             for(i=0; i < dims; i++)
2466             {
2467                 h->thresh2[i] = dim_ranges;
2468                 cvReadRawDataSlice( fs, &reader, size[i]+1, dim_ranges, "f" );
2469                 dim_ranges += size[i] + 1;
2470             }
2471         }
2472     }
2473
2474     return h;
2475 }
2476
2477 static void icvWriteHist( CvFileStorage* fs, const char* name,
2478                           const void* struct_ptr, CvAttrList /*attributes*/ )
2479 {
2480     const CvHistogram * hist = (const CvHistogram *) struct_ptr;
2481     int sizes[CV_MAX_DIM];
2482     int dims;
2483     int i;
2484     int is_uniform, have_ranges;
2485
2486     cvStartWriteStruct( fs, name, CV_NODE_MAP, CV_TYPE_NAME_HIST );
2487
2488     is_uniform = (CV_IS_UNIFORM_HIST(hist) ? 1 : 0);
2489     have_ranges = (hist->type & CV_HIST_RANGES_FLAG ? 1 : 0);
2490
2491     cvWriteInt( fs, "type", (hist->type & 1) );
2492     cvWriteInt( fs, "is_uniform", is_uniform );
2493     cvWriteInt( fs, "have_ranges", have_ranges );
2494     if(!CV_IS_SPARSE_HIST(hist))
2495         cvWrite( fs, "mat", &(hist->mat) );
2496     else
2497         cvWrite( fs, "bins", hist->bins );
2498
2499     // write thresholds
2500     if(have_ranges){
2501         dims = cvGetDims( hist->bins, sizes );
2502         cvStartWriteStruct( fs, "thresh", CV_NODE_SEQ + CV_NODE_FLOW );
2503         if(is_uniform){
2504             for(i=0; i<dims; i++){
2505                 cvWriteRawData( fs, hist->thresh[i], 2, "f" );
2506             }
2507         }
2508         else{
2509             for(i=0; i<dims; i++){
2510                 cvWriteRawData( fs, hist->thresh2[i], sizes[i]+1, "f" );
2511             }
2512         }
2513         cvEndWriteStruct( fs );
2514     }
2515
2516     cvEndWriteStruct( fs );
2517 }
2518
2519
2520 CvType hist_type( CV_TYPE_NAME_HIST, icvIsHist, (CvReleaseFunc)cvReleaseHist,
2521                   icvReadHist, icvWriteHist, (CvCloneFunc)icvCloneHist );
2522
2523 /* End of file. */
2524