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