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 */
44CvHistogram *
45cvCreateHist( 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 */
91CV_IMPL CvHistogram*
92cvMakeHistHeaderForArray( 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
134CV_IMPL void
135cvReleaseHist( 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
170CV_IMPL void
171cvClearHist( 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
187CV_IMPL void
188cvThreshHist( 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)
223CV_IMPL void
224cvNormalizeHist( 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
272CV_IMPL void
273cvGetMinMaxHistValue( 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
398CV_IMPL double
399cvCompareHist( 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
650CV_IMPL void
651cvCopyHist( 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
720CV_IMPL void
721cvSetHistBinRanges( 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
793static CvStatus
794icvCalcHistLookupTables8u( 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
870static 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
1123static 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
1401CV_IMPL void
1402cvCalcArrHist( 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
1552static 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
1814static 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
2047CV_IMPL void
2048cvCalcArrBackProject( 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
2140CV_IMPL void
2141cvCalcArrBackProjectPatch( 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( !arr )
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
2221CV_IMPL void
2222cvCalcBayesianProb( 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
2260CV_IMPL void
2261cvCalcProbDensity( 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
2311CV_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
2356static int icvIsHist( const void * ptr ){
2357	return CV_IS_HIST( ((CvHistogram*)ptr) );
2358}
2359
2360static CvHistogram * icvCloneHist( const CvHistogram * src ){
2361	CvHistogram * dst=NULL;
2362	cvCopyHist(src, &dst);
2363	return dst;
2364}
2365
2366static 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
2459static 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
2509CvType hist_type( CV_TYPE_NAME_HIST, icvIsHist, (CvReleaseFunc)cvReleaseHist,
2510                  icvReadHist, icvWriteHist, (CvCloneFunc)icvCloneHist );
2511
2512/* End of file. */
2513
2514