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