Update to 2.0.0 tree from current Fremantle build
[opencv] / src / ml / ml_inner_functions.cpp
1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
5 //  By downloading, copying, installing or using the software you agree to this license.
6 //  If you do not agree to this license, do not download, install,
7 //  copy or use the software.
8 //
9 //
10 //                        Intel License Agreement
11 //
12 // Copyright (C) 2000, Intel Corporation, all rights reserved.
13 // Third party copyrights are property of their respective owners.
14 //
15 // Redistribution and use in source and binary forms, with or without modification,
16 // are permitted provided that the following conditions are met:
17 //
18 //   * Redistribution's of source code must retain the above copyright notice,
19 //     this list of conditions and the following disclaimer.
20 //
21 //   * Redistribution's in binary form must reproduce the above copyright notice,
22 //     this list of conditions and the following disclaimer in the documentation
23 //     and/or other materials provided with the distribution.
24 //
25 //   * The name of Intel Corporation may not be used to endorse or promote products
26 //     derived from this software without specific prior written permission.
27 //
28 // This software is provided by the copyright holders and contributors "as is" and
29 // any express or implied warranties, including, but not limited to, the implied
30 // warranties of merchantability and fitness for a particular purpose are disclaimed.
31 // In no event shall the Intel Corporation or contributors be liable for any direct,
32 // indirect, incidental, special, exemplary, or consequential damages
33 // (including, but not limited to, procurement of substitute goods or services;
34 // loss of use, data, or profits; or business interruption) however caused
35 // and on any theory of liability, whether in contract, strict liability,
36 // or tort (including negligence or otherwise) arising in any way out of
37 // the use of this software, even if advised of the possibility of such damage.
38 //
39 //M*/
40
41 #include "_ml.h"
42
43
44 CvStatModel::CvStatModel()
45 {
46     default_model_name = "my_stat_model";
47 }
48
49
50 CvStatModel::~CvStatModel()
51 {
52     clear();
53 }
54
55
56 void CvStatModel::clear()
57 {
58 }
59
60
61 void CvStatModel::save( const char* filename, const char* name ) const
62 {
63     CvFileStorage* fs = 0;
64
65     CV_FUNCNAME( "CvStatModel::save" );
66
67     __BEGIN__;
68
69     CV_CALL( fs = cvOpenFileStorage( filename, 0, CV_STORAGE_WRITE ));
70     if( !fs )
71         CV_ERROR( CV_StsError, "Could not open the file storage. Check the path and permissions" );
72
73     write( fs, name ? name : default_model_name );
74
75     __END__;
76
77     cvReleaseFileStorage( &fs );
78 }
79
80
81 void CvStatModel::load( const char* filename, const char* name )
82 {
83     CvFileStorage* fs = 0;
84
85     CV_FUNCNAME( "CvStatModel::load" );
86
87     __BEGIN__;
88
89     CvFileNode* model_node = 0;
90
91     CV_CALL( fs = cvOpenFileStorage( filename, 0, CV_STORAGE_READ ));
92     if( !fs )
93         EXIT;
94
95     if( name )
96         model_node = cvGetFileNodeByName( fs, 0, name );
97     else
98     {
99         CvFileNode* root = cvGetRootFileNode( fs );
100         if( root->data.seq->total > 0 )
101             model_node = (CvFileNode*)cvGetSeqElem( root->data.seq, 0 );
102     }
103
104     read( fs, model_node );
105
106     __END__;
107
108     cvReleaseFileStorage( &fs );
109 }
110
111
112 void CvStatModel::write( CvFileStorage*, const char* ) const
113 {
114     OPENCV_ERROR( CV_StsNotImplemented, "CvStatModel::write", "" );
115 }
116
117
118 void CvStatModel::read( CvFileStorage*, CvFileNode* )
119 {
120     OPENCV_ERROR( CV_StsNotImplemented, "CvStatModel::read", "" );
121 }
122
123
124 /* Calculates upper triangular matrix S, where A is a symmetrical matrix A=S'*S */
125 CV_IMPL void cvChol( CvMat* A, CvMat* S )
126 {
127     int dim = A->rows;
128
129     int i, j, k;
130     float sum;
131
132     for( i = 0; i < dim; i++ )
133     {
134         for( j = 0; j < i; j++ )
135             CV_MAT_ELEM(*S, float, i, j) = 0;
136
137         sum = 0;
138         for( k = 0; k < i; k++ )
139             sum += CV_MAT_ELEM(*S, float, k, i) * CV_MAT_ELEM(*S, float, k, i);
140
141         CV_MAT_ELEM(*S, float, i, i) = (float)sqrt(CV_MAT_ELEM(*A, float, i, i) - sum);
142
143         for( j = i + 1; j < dim; j++ )
144         {
145             sum = 0;
146             for( k = 0; k < i; k++ )
147                 sum += CV_MAT_ELEM(*S, float, k, i) * CV_MAT_ELEM(*S, float, k, j);
148
149             CV_MAT_ELEM(*S, float, i, j) =
150                 (CV_MAT_ELEM(*A, float, i, j) - sum) / CV_MAT_ELEM(*S, float, i, i);
151
152         }
153     }
154 }
155
156 /* Generates <sample> from multivariate normal distribution, where <mean> - is an
157    average row vector, <cov> - symmetric covariation matrix */
158 CV_IMPL void cvRandMVNormal( CvMat* mean, CvMat* cov, CvMat* sample, CvRNG* rng )
159 {
160     int dim = sample->cols;
161     int amount = sample->rows;
162
163     CvRNG state = rng ? *rng : cvRNG( cvGetTickCount() );
164     cvRandArr(&state, sample, CV_RAND_NORMAL, cvScalarAll(0), cvScalarAll(1) );
165
166     CvMat* utmat = cvCreateMat(dim, dim, sample->type);
167     CvMat* vect = cvCreateMatHeader(1, dim, sample->type);
168
169     cvChol(cov, utmat);
170
171     int i;
172     for( i = 0; i < amount; i++ )
173     {
174         cvGetRow(sample, vect, i);
175         cvMatMulAdd(vect, utmat, mean, vect);
176     }
177
178     cvReleaseMat(&vect);
179     cvReleaseMat(&utmat);
180 }
181
182
183 /* Generates <sample> of <amount> points from a discrete variate xi,
184    where Pr{xi = k} == probs[k], 0 < k < len - 1. */
185 CV_IMPL void cvRandSeries( float probs[], int len, int sample[], int amount )
186 {
187     CvMat* univals = cvCreateMat(1, amount, CV_32FC1);
188     float* knots = (float*)cvAlloc( len * sizeof(float) );
189
190     int i, j;
191
192     CvRNG state = cvRNG(-1);
193     cvRandArr(&state, univals, CV_RAND_UNI, cvScalarAll(0), cvScalarAll(1) );
194
195     knots[0] = probs[0];
196     for( i = 1; i < len; i++ )
197         knots[i] = knots[i - 1] + probs[i];
198
199     for( i = 0; i < amount; i++ )
200         for( j = 0; j < len; j++ )
201         {
202             if ( CV_MAT_ELEM(*univals, float, 0, i) <= knots[j] )
203             {
204                 sample[i] = j;
205                 break;
206             }
207         }
208
209     cvFree(&knots);
210 }
211
212 /* Generates <sample> from gaussian mixture distribution */
213 CV_IMPL void cvRandGaussMixture( CvMat* means[],
214                                  CvMat* covs[],
215                                  float weights[],
216                                  int clsnum,
217                                  CvMat* sample,
218                                  CvMat* sampClasses )
219 {
220     int dim = sample->cols;
221     int amount = sample->rows;
222
223     int i, clss;
224
225     int* sample_clsnum = (int*)cvAlloc( amount * sizeof(int) );
226     CvMat** utmats = (CvMat**)cvAlloc( clsnum * sizeof(CvMat*) );
227     CvMat* vect = cvCreateMatHeader(1, dim, CV_32FC1);
228
229     CvMat* classes;
230     if( sampClasses )
231         classes = sampClasses;
232     else
233         classes = cvCreateMat(1, amount, CV_32FC1);
234
235     CvRNG state = cvRNG(-1);
236     cvRandArr(&state, sample, CV_RAND_NORMAL, cvScalarAll(0), cvScalarAll(1));
237
238     cvRandSeries(weights, clsnum, sample_clsnum, amount);
239
240     for( i = 0; i < clsnum; i++ )
241     {
242         utmats[i] = cvCreateMat(dim, dim, CV_32FC1);
243         cvChol(covs[i], utmats[i]);
244     }
245
246     for( i = 0; i < amount; i++ )
247     {
248         CV_MAT_ELEM(*classes, float, 0, i) = (float)sample_clsnum[i];
249         cvGetRow(sample, vect, i);
250         clss = sample_clsnum[i];
251         cvMatMulAdd(vect, utmats[clss], means[clss], vect);
252     }
253
254     if( !sampClasses )
255         cvReleaseMat(&classes);
256     for( i = 0; i < clsnum; i++ )
257         cvReleaseMat(&utmats[i]);
258     cvFree(&utmats);
259     cvFree(&sample_clsnum);
260     cvReleaseMat(&vect);
261 }
262
263
264 CvMat* icvGenerateRandomClusterCenters ( int seed, const CvMat* data,
265                                          int num_of_clusters, CvMat* _centers )
266 {
267     CvMat* centers = _centers;
268
269     CV_FUNCNAME("icvGenerateRandomClusterCenters");
270     __BEGIN__;
271
272     CvRNG rng;
273     CvMat data_comp, centers_comp;
274     CvPoint minLoc, maxLoc; // Not used, just for function "cvMinMaxLoc"
275     double minVal, maxVal;
276     int i;
277     int dim = data ? data->cols : 0;
278
279     if( ICV_IS_MAT_OF_TYPE(data, CV_32FC1) )
280     {
281         if( _centers && !ICV_IS_MAT_OF_TYPE (_centers, CV_32FC1) )
282         {
283             CV_ERROR(CV_StsBadArg,"");
284         }
285         else if( !_centers )
286             CV_CALL(centers = cvCreateMat (num_of_clusters, dim, CV_32FC1));
287     }
288     else if( ICV_IS_MAT_OF_TYPE(data, CV_64FC1) )
289     {
290         if( _centers && !ICV_IS_MAT_OF_TYPE (_centers, CV_64FC1) )
291         {
292             CV_ERROR(CV_StsBadArg,"");
293         }
294         else if( !_centers )
295             CV_CALL(centers = cvCreateMat (num_of_clusters, dim, CV_64FC1));
296     }
297     else
298         CV_ERROR (CV_StsBadArg,"");
299
300     if( num_of_clusters < 1 )
301         CV_ERROR (CV_StsBadArg,"");
302
303     rng = cvRNG(seed);
304     for (i = 0; i < dim; i++)
305     {
306         CV_CALL(cvGetCol (data, &data_comp, i));
307         CV_CALL(cvMinMaxLoc (&data_comp, &minVal, &maxVal, &minLoc, &maxLoc));
308         CV_CALL(cvGetCol (centers, &centers_comp, i));
309         CV_CALL(cvRandArr (&rng, &centers_comp, CV_RAND_UNI, cvScalarAll(minVal), cvScalarAll(maxVal)));
310     }
311
312     __END__;
313
314     if( (cvGetErrStatus () < 0) || (centers != _centers) )
315         cvReleaseMat (&centers);
316
317     return _centers ? _centers : centers;
318 } // end of icvGenerateRandomClusterCenters
319
320 // By S. Dilman - begin -
321
322 #define ICV_RAND_MAX    4294967296 // == 2^32
323
324 CV_IMPL void cvRandRoundUni (CvMat* center,
325                              float radius_small,
326                              float radius_large,
327                              CvMat* desired_matrix,
328                              CvRNG* rng_state_ptr)
329 {
330     float rad, norm, coefficient;
331     int dim, size, i, j;
332     CvMat *cov, sample;
333     CvRNG rng_local;
334
335     CV_FUNCNAME("cvRandRoundUni");
336     __BEGIN__
337
338     rng_local = *rng_state_ptr;
339
340     CV_ASSERT ((radius_small >= 0) &&
341                (radius_large > 0) &&
342                (radius_small <= radius_large));
343     CV_ASSERT (center && desired_matrix && rng_state_ptr);
344     CV_ASSERT (center->rows == 1);
345     CV_ASSERT (center->cols == desired_matrix->cols);
346
347     dim = desired_matrix->cols;
348     size = desired_matrix->rows;
349     cov = cvCreateMat (dim, dim, CV_32FC1);
350     cvSetIdentity (cov);
351     cvRandMVNormal (center, cov, desired_matrix, &rng_local);
352
353     for (i = 0; i < size; i++)
354     {
355         rad = (float)(cvRandReal(&rng_local)*(radius_large - radius_small) + radius_small);
356         cvGetRow (desired_matrix, &sample, i);
357         norm = (float) cvNorm (&sample, 0, CV_L2);
358         coefficient = rad / norm;
359         for (j = 0; j < dim; j++)
360              CV_MAT_ELEM (sample, float, 0, j) *= coefficient;
361     }
362
363     __END__
364
365 }
366
367 // By S. Dilman - end -
368
369 static int CV_CDECL
370 icvCmpIntegers( const void* a, const void* b )
371 {
372     return *(const int*)a - *(const int*)b;
373 }
374
375
376 static int CV_CDECL
377 icvCmpIntegersPtr( const void* _a, const void* _b )
378 {
379     int a = **(const int**)_a;
380     int b = **(const int**)_b;
381     return (a < b ? -1 : 0)|(a > b);
382 }
383
384
385 static int icvCmpSparseVecElems( const void* a, const void* b )
386 {
387     return ((CvSparseVecElem32f*)a)->idx - ((CvSparseVecElem32f*)b)->idx;
388 }
389
390
391 CvMat*
392 cvPreprocessIndexArray( const CvMat* idx_arr, int data_arr_size, bool check_for_duplicates )
393 {
394     CvMat* idx = 0;
395
396     CV_FUNCNAME( "cvPreprocessIndexArray" );
397
398     __BEGIN__;
399
400     int i, idx_total, idx_selected = 0, step, type, prev = INT_MIN, is_sorted = 1;
401     uchar* srcb = 0;
402     int* srci = 0;
403     int* dsti;
404
405     if( !CV_IS_MAT(idx_arr) )
406         CV_ERROR( CV_StsBadArg, "Invalid index array" );
407
408     if( idx_arr->rows != 1 && idx_arr->cols != 1 )
409         CV_ERROR( CV_StsBadSize, "the index array must be 1-dimensional" );
410
411     idx_total = idx_arr->rows + idx_arr->cols - 1;
412     srcb = idx_arr->data.ptr;
413     srci = idx_arr->data.i;
414
415     type = CV_MAT_TYPE(idx_arr->type);
416     step = CV_IS_MAT_CONT(idx_arr->type) ? 1 : idx_arr->step/CV_ELEM_SIZE(type);
417
418     switch( type )
419     {
420     case CV_8UC1:
421     case CV_8SC1:
422         // idx_arr is array of 1's and 0's -
423         // i.e. it is a mask of the selected components
424         if( idx_total != data_arr_size )
425             CV_ERROR( CV_StsUnmatchedSizes,
426             "Component mask should contain as many elements as the total number of input variables" );
427
428         for( i = 0; i < idx_total; i++ )
429             idx_selected += srcb[i*step] != 0;
430
431         if( idx_selected == 0 )
432             CV_ERROR( CV_StsOutOfRange, "No components/input_variables is selected!" );
433
434         if( idx_selected == idx_total )
435             EXIT;
436         break;
437     case CV_32SC1:
438         // idx_arr is array of integer indices of selected components
439         if( idx_total > data_arr_size )
440             CV_ERROR( CV_StsOutOfRange,
441             "index array may not contain more elements than the total number of input variables" );
442         idx_selected = idx_total;
443         // check if sorted already
444         for( i = 0; i < idx_total; i++ )
445         {
446             int val = srci[i*step];
447             if( val >= prev )
448             {
449                 is_sorted = 0;
450                 break;
451             }
452             prev = val;
453         }
454         break;
455     default:
456         CV_ERROR( CV_StsUnsupportedFormat, "Unsupported index array data type "
457                                            "(it should be 8uC1, 8sC1 or 32sC1)" );
458     }
459
460     CV_CALL( idx = cvCreateMat( 1, idx_selected, CV_32SC1 ));
461     dsti = idx->data.i;
462
463     if( type < CV_32SC1 )
464     {
465         for( i = 0; i < idx_total; i++ )
466             if( srcb[i*step] )
467                 *dsti++ = i;
468     }
469     else
470     {
471         for( i = 0; i < idx_total; i++ )
472             dsti[i] = srci[i*step];
473
474         if( !is_sorted )
475             qsort( dsti, idx_total, sizeof(dsti[0]), icvCmpIntegers );
476
477         if( dsti[0] < 0 || dsti[idx_total-1] >= data_arr_size )
478             CV_ERROR( CV_StsOutOfRange, "the index array elements are out of range" );
479
480         if( check_for_duplicates )
481         {
482             for( i = 1; i < idx_total; i++ )
483                 if( dsti[i] <= dsti[i-1] )
484                     CV_ERROR( CV_StsBadArg, "There are duplicated index array elements" );
485         }
486     }
487
488     __END__;
489
490     if( cvGetErrStatus() < 0 )
491         cvReleaseMat( &idx );
492
493     return idx;
494 }
495
496
497 CvMat*
498 cvPreprocessVarType( const CvMat* var_type, const CvMat* var_idx,
499                      int var_count, int* response_type )
500 {
501     CvMat* out_var_type = 0;
502     CV_FUNCNAME( "cvPreprocessVarType" );
503
504     if( response_type )
505         *response_type = -1;
506
507     __BEGIN__;
508
509     int i, tm_size, tm_step;
510     //int* map = 0;
511     const uchar* src;
512     uchar* dst;
513
514     if( !CV_IS_MAT(var_type) )
515         CV_ERROR( var_type ? CV_StsBadArg : CV_StsNullPtr, "Invalid or absent var_type array" );
516
517     if( var_type->rows != 1 && var_type->cols != 1 )
518         CV_ERROR( CV_StsBadSize, "var_type array must be 1-dimensional" );
519
520     if( !CV_IS_MASK_ARR(var_type))
521         CV_ERROR( CV_StsUnsupportedFormat, "type mask must be 8uC1 or 8sC1 array" );
522
523     tm_size = var_type->rows + var_type->cols - 1;
524     tm_step = var_type->rows == 1 ? 1 : var_type->step/CV_ELEM_SIZE(var_type->type);
525
526     if( /*tm_size != var_count &&*/ tm_size != var_count + 1 )
527         CV_ERROR( CV_StsBadArg,
528         "type mask must be of <input var count> + 1 size" );
529
530     if( response_type && tm_size > var_count )
531         *response_type = var_type->data.ptr[var_count*tm_step] != 0;
532
533     if( var_idx )
534     {
535         if( !CV_IS_MAT(var_idx) || CV_MAT_TYPE(var_idx->type) != CV_32SC1 ||
536             (var_idx->rows != 1 && var_idx->cols != 1) || !CV_IS_MAT_CONT(var_idx->type) )
537             CV_ERROR( CV_StsBadArg, "var index array should be continuous 1-dimensional integer vector" );
538         if( var_idx->rows + var_idx->cols - 1 > var_count )
539             CV_ERROR( CV_StsBadSize, "var index array is too large" );
540         //map = var_idx->data.i;
541         var_count = var_idx->rows + var_idx->cols - 1;
542     }
543
544     CV_CALL( out_var_type = cvCreateMat( 1, var_count, CV_8UC1 ));
545     src = var_type->data.ptr;
546     dst = out_var_type->data.ptr;
547
548     for( i = 0; i < var_count; i++ )
549     {
550         //int idx = map ? map[i] : i;
551         assert( (unsigned)/*idx*/i < (unsigned)tm_size );
552         dst[i] = (uchar)(src[/*idx*/i*tm_step] != 0);
553     }
554
555     __END__;
556
557     return out_var_type;
558 }
559
560
561 CvMat*
562 cvPreprocessOrderedResponses( const CvMat* responses, const CvMat* sample_idx, int sample_all )
563 {
564     CvMat* out_responses = 0;
565
566     CV_FUNCNAME( "cvPreprocessOrderedResponses" );
567
568     __BEGIN__;
569
570     int i, r_type, r_step;
571     const int* map = 0;
572     float* dst;
573     int sample_count = sample_all;
574
575     if( !CV_IS_MAT(responses) )
576         CV_ERROR( CV_StsBadArg, "Invalid response array" );
577
578     if( responses->rows != 1 && responses->cols != 1 )
579         CV_ERROR( CV_StsBadSize, "Response array must be 1-dimensional" );
580
581     if( responses->rows + responses->cols - 1 != sample_count )
582         CV_ERROR( CV_StsUnmatchedSizes,
583         "Response array must contain as many elements as the total number of samples" );
584
585     r_type = CV_MAT_TYPE(responses->type);
586     if( r_type != CV_32FC1 && r_type != CV_32SC1 )
587         CV_ERROR( CV_StsUnsupportedFormat, "Unsupported response type" );
588
589     r_step = responses->step ? responses->step / CV_ELEM_SIZE(responses->type) : 1;
590
591     if( r_type == CV_32FC1 && CV_IS_MAT_CONT(responses->type) && !sample_idx )
592     {
593         out_responses = (CvMat*)responses;
594         EXIT;
595     }
596
597     if( sample_idx )
598     {
599         if( !CV_IS_MAT(sample_idx) || CV_MAT_TYPE(sample_idx->type) != CV_32SC1 ||
600             (sample_idx->rows != 1 && sample_idx->cols != 1) || !CV_IS_MAT_CONT(sample_idx->type) )
601             CV_ERROR( CV_StsBadArg, "sample index array should be continuous 1-dimensional integer vector" );
602         if( sample_idx->rows + sample_idx->cols - 1 > sample_count )
603             CV_ERROR( CV_StsBadSize, "sample index array is too large" );
604         map = sample_idx->data.i;
605         sample_count = sample_idx->rows + sample_idx->cols - 1;
606     }
607
608     CV_CALL( out_responses = cvCreateMat( 1, sample_count, CV_32FC1 ));
609
610     dst = out_responses->data.fl;
611     if( r_type == CV_32FC1 )
612     {
613         const float* src = responses->data.fl;
614         for( i = 0; i < sample_count; i++ )
615         {
616             int idx = map ? map[i] : i;
617             assert( (unsigned)idx < (unsigned)sample_all );
618             dst[i] = src[idx*r_step];
619         }
620     }
621     else
622     {
623         const int* src = responses->data.i;
624         for( i = 0; i < sample_count; i++ )
625         {
626             int idx = map ? map[i] : i;
627             assert( (unsigned)idx < (unsigned)sample_all );
628             dst[i] = (float)src[idx*r_step];
629         }
630     }
631
632     __END__;
633
634     return out_responses;
635 }
636
637 CvMat*
638 cvPreprocessCategoricalResponses( const CvMat* responses,
639     const CvMat* sample_idx, int sample_all,
640     CvMat** out_response_map, CvMat** class_counts )
641 {
642     CvMat* out_responses = 0;
643     int** response_ptr = 0;
644
645     CV_FUNCNAME( "cvPreprocessCategoricalResponses" );
646
647     if( out_response_map )
648         *out_response_map = 0;
649
650     if( class_counts )
651         *class_counts = 0;
652
653     __BEGIN__;
654
655     int i, r_type, r_step;
656     int cls_count = 1, prev_cls, prev_i;
657     const int* map = 0;
658     const int* srci;
659     const float* srcfl;
660     int* dst;
661     int* cls_map;
662     int* cls_counts = 0;
663     int sample_count = sample_all;
664
665     if( !CV_IS_MAT(responses) )
666         CV_ERROR( CV_StsBadArg, "Invalid response array" );
667
668     if( responses->rows != 1 && responses->cols != 1 )
669         CV_ERROR( CV_StsBadSize, "Response array must be 1-dimensional" );
670
671     if( responses->rows + responses->cols - 1 != sample_count )
672         CV_ERROR( CV_StsUnmatchedSizes,
673         "Response array must contain as many elements as the total number of samples" );
674
675     r_type = CV_MAT_TYPE(responses->type);
676     if( r_type != CV_32FC1 && r_type != CV_32SC1 )
677         CV_ERROR( CV_StsUnsupportedFormat, "Unsupported response type" );
678
679     r_step = responses->step ? responses->step / CV_ELEM_SIZE(responses->type) : 1;
680
681     if( sample_idx )
682     {
683         if( !CV_IS_MAT(sample_idx) || CV_MAT_TYPE(sample_idx->type) != CV_32SC1 ||
684             (sample_idx->rows != 1 && sample_idx->cols != 1) || !CV_IS_MAT_CONT(sample_idx->type) )
685             CV_ERROR( CV_StsBadArg, "sample index array should be continuous 1-dimensional integer vector" );
686         if( sample_idx->rows + sample_idx->cols - 1 > sample_count )
687             CV_ERROR( CV_StsBadSize, "sample index array is too large" );
688         map = sample_idx->data.i;
689         sample_count = sample_idx->rows + sample_idx->cols - 1;
690     }
691
692     CV_CALL( out_responses = cvCreateMat( 1, sample_count, CV_32SC1 ));
693
694     if( !out_response_map )
695         CV_ERROR( CV_StsNullPtr, "out_response_map pointer is NULL" );
696
697     CV_CALL( response_ptr = (int**)cvAlloc( sample_count*sizeof(response_ptr[0])));
698
699     srci = responses->data.i;
700     srcfl = responses->data.fl;
701     dst = out_responses->data.i;
702
703     for( i = 0; i < sample_count; i++ )
704     {
705         int idx = map ? map[i] : i;
706         assert( (unsigned)idx < (unsigned)sample_all );
707         if( r_type == CV_32SC1 )
708             dst[i] = srci[idx*r_step];
709         else
710         {
711             float rf = srcfl[idx*r_step];
712             int ri = cvRound(rf);
713             if( ri != rf )
714             {
715                 char buf[100];
716                 sprintf( buf, "response #%d is not integral", idx );
717                 CV_ERROR( CV_StsBadArg, buf );
718             }
719             dst[i] = ri;
720         }
721         response_ptr[i] = dst + i;
722     }
723
724     qsort( response_ptr, sample_count, sizeof(int*), icvCmpIntegersPtr );
725
726     // count the classes
727     for( i = 1; i < sample_count; i++ )
728         cls_count += *response_ptr[i] != *response_ptr[i-1];
729
730     if( cls_count < 2 )
731         CV_ERROR( CV_StsBadArg, "There is only a single class" );
732
733     CV_CALL( *out_response_map = cvCreateMat( 1, cls_count, CV_32SC1 ));
734
735     if( class_counts )
736     {
737         CV_CALL( *class_counts = cvCreateMat( 1, cls_count, CV_32SC1 ));
738         cls_counts = (*class_counts)->data.i;
739     }
740
741     // compact the class indices and build the map
742     prev_cls = ~*response_ptr[0];
743     cls_count = -1;
744     cls_map = (*out_response_map)->data.i;
745
746     for( i = 0, prev_i = -1; i < sample_count; i++ )
747     {
748         int cur_cls = *response_ptr[i];
749         if( cur_cls != prev_cls )
750         {
751             if( cls_counts && cls_count >= 0 )
752                 cls_counts[cls_count] = i - prev_i;
753             cls_map[++cls_count] = prev_cls = cur_cls;
754             prev_i = i;
755         }
756         *response_ptr[i] = cls_count;
757     }
758
759     if( cls_counts )
760         cls_counts[cls_count] = i - prev_i;
761
762     __END__;
763
764     cvFree( &response_ptr );
765
766     return out_responses;
767 }
768
769
770 const float**
771 cvGetTrainSamples( const CvMat* train_data, int tflag,
772                    const CvMat* var_idx, const CvMat* sample_idx,
773                    int* _var_count, int* _sample_count,
774                    bool always_copy_data )
775 {
776     float** samples = 0;
777
778     CV_FUNCNAME( "cvGetTrainSamples" );
779
780     __BEGIN__;
781
782     int i, j, var_count, sample_count, s_step, v_step;
783     bool copy_data;
784     const float* data;
785     const int *s_idx, *v_idx;
786
787     if( !CV_IS_MAT(train_data) )
788         CV_ERROR( CV_StsBadArg, "Invalid or NULL training data matrix" );
789
790     var_count = var_idx ? var_idx->cols + var_idx->rows - 1 :
791                 tflag == CV_ROW_SAMPLE ? train_data->cols : train_data->rows;
792     sample_count = sample_idx ? sample_idx->cols + sample_idx->rows - 1 :
793                    tflag == CV_ROW_SAMPLE ? train_data->rows : train_data->cols;
794
795     if( _var_count )
796         *_var_count = var_count;
797
798     if( _sample_count )
799         *_sample_count = sample_count;
800
801     copy_data = tflag != CV_ROW_SAMPLE || var_idx || always_copy_data;
802
803     CV_CALL( samples = (float**)cvAlloc(sample_count*sizeof(samples[0]) +
804                 (copy_data ? 1 : 0)*var_count*sample_count*sizeof(samples[0][0])) );
805     data = train_data->data.fl;
806     s_step = train_data->step / sizeof(samples[0][0]);
807     v_step = 1;
808     s_idx = sample_idx ? sample_idx->data.i : 0;
809     v_idx = var_idx ? var_idx->data.i : 0;
810
811     if( !copy_data )
812     {
813         for( i = 0; i < sample_count; i++ )
814             samples[i] = (float*)(data + (s_idx ? s_idx[i] : i)*s_step);
815     }
816     else
817     {
818         samples[0] = (float*)(samples + sample_count);
819         if( tflag != CV_ROW_SAMPLE )
820             CV_SWAP( s_step, v_step, i );
821
822         for( i = 0; i < sample_count; i++ )
823         {
824             float* dst = samples[i] = samples[0] + i*var_count;
825             const float* src = data + (s_idx ? s_idx[i] : i)*s_step;
826
827             if( !v_idx )
828                 for( j = 0; j < var_count; j++ )
829                     dst[j] = src[j*v_step];
830             else
831                 for( j = 0; j < var_count; j++ )
832                     dst[j] = src[v_idx[j]*v_step];
833         }
834     }
835
836     __END__;
837
838     return (const float**)samples;
839 }
840
841
842 void
843 cvCheckTrainData( const CvMat* train_data, int tflag,
844                   const CvMat* missing_mask,
845                   int* var_all, int* sample_all )
846 {
847     CV_FUNCNAME( "cvCheckTrainData" );
848
849     if( var_all )
850         *var_all = 0;
851
852     if( sample_all )
853         *sample_all = 0;
854
855     __BEGIN__;
856
857     // check parameter types and sizes
858     if( !CV_IS_MAT(train_data) || CV_MAT_TYPE(train_data->type) != CV_32FC1 )
859         CV_ERROR( CV_StsBadArg, "train data must be floating-point matrix" );
860
861     if( missing_mask )
862     {
863         if( !CV_IS_MAT(missing_mask) || !CV_IS_MASK_ARR(missing_mask) ||
864             !CV_ARE_SIZES_EQ(train_data, missing_mask) )
865             CV_ERROR( CV_StsBadArg,
866             "missing value mask must be 8-bit matrix of the same size as training data" );
867     }
868
869     if( tflag != CV_ROW_SAMPLE && tflag != CV_COL_SAMPLE )
870         CV_ERROR( CV_StsBadArg,
871         "Unknown training data layout (must be CV_ROW_SAMPLE or CV_COL_SAMPLE)" );
872
873     if( var_all )
874         *var_all = tflag == CV_ROW_SAMPLE ? train_data->cols : train_data->rows;
875
876     if( sample_all )
877         *sample_all = tflag == CV_ROW_SAMPLE ? train_data->rows : train_data->cols;
878
879     __END__;
880 }
881
882
883 int
884 cvPrepareTrainData( const char* /*funcname*/,
885                     const CvMat* train_data, int tflag,
886                     const CvMat* responses, int response_type,
887                     const CvMat* var_idx,
888                     const CvMat* sample_idx,
889                     bool always_copy_data,
890                     const float*** out_train_samples,
891                     int* _sample_count,
892                     int* _var_count,
893                     int* _var_all,
894                     CvMat** out_responses,
895                     CvMat** out_response_map,
896                     CvMat** out_var_idx,
897                     CvMat** out_sample_idx )
898 {
899     int ok = 0;
900     CvMat* _var_idx = 0;
901     CvMat* _sample_idx = 0;
902     CvMat* _responses = 0;
903     int sample_all = 0, sample_count = 0, var_all = 0, var_count = 0;
904
905     CV_FUNCNAME( "cvPrepareTrainData" );
906
907     // step 0. clear all the output pointers to ensure we do not try
908     // to call free() with uninitialized pointers
909     if( out_responses )
910         *out_responses = 0;
911
912     if( out_response_map )
913         *out_response_map = 0;
914
915     if( out_var_idx )
916         *out_var_idx = 0;
917
918     if( out_sample_idx )
919         *out_sample_idx = 0;
920
921     if( out_train_samples )
922         *out_train_samples = 0;
923
924     if( _sample_count )
925         *_sample_count = 0;
926
927     if( _var_count )
928         *_var_count = 0;
929
930     if( _var_all )
931         *_var_all = 0;
932
933     __BEGIN__;
934
935     if( !out_train_samples )
936         CV_ERROR( CV_StsBadArg, "output pointer to train samples is NULL" );
937
938     CV_CALL( cvCheckTrainData( train_data, tflag, 0, &var_all, &sample_all ));
939
940     if( sample_idx )
941         CV_CALL( _sample_idx = cvPreprocessIndexArray( sample_idx, sample_all ));
942     if( var_idx )
943         CV_CALL( _var_idx = cvPreprocessIndexArray( var_idx, var_all ));
944
945     if( responses )
946     {
947         if( !out_responses )
948             CV_ERROR( CV_StsNullPtr, "output response pointer is NULL" );
949
950         if( response_type == CV_VAR_NUMERICAL )
951         {
952             CV_CALL( _responses = cvPreprocessOrderedResponses( responses,
953                                                 _sample_idx, sample_all ));
954         }
955         else
956         {
957             CV_CALL( _responses = cvPreprocessCategoricalResponses( responses,
958                                 _sample_idx, sample_all, out_response_map, 0 ));
959         }
960     }
961
962     CV_CALL( *out_train_samples =
963                 cvGetTrainSamples( train_data, tflag, _var_idx, _sample_idx,
964                                    &var_count, &sample_count, always_copy_data ));
965
966     ok = 1;
967
968     __END__;
969
970     if( ok )
971     {
972         if( out_responses )
973             *out_responses = _responses, _responses = 0;
974
975         if( out_var_idx )
976             *out_var_idx = _var_idx, _var_idx = 0;
977
978         if( out_sample_idx )
979             *out_sample_idx = _sample_idx, _sample_idx = 0;
980
981         if( _sample_count )
982             *_sample_count = sample_count;
983
984         if( _var_count )
985             *_var_count = var_count;
986
987         if( _var_all )
988             *_var_all = var_all;
989     }
990     else
991     {
992         if( out_response_map )
993             cvReleaseMat( out_response_map );
994         cvFree( out_train_samples );
995     }
996
997     if( _responses != responses )
998         cvReleaseMat( &_responses );
999     cvReleaseMat( &_var_idx );
1000     cvReleaseMat( &_sample_idx );
1001
1002     return ok;
1003 }
1004
1005
1006 typedef struct CvSampleResponsePair
1007 {
1008     const float* sample;
1009     const uchar* mask;
1010     int response;
1011     int index;
1012 }
1013 CvSampleResponsePair;
1014
1015
1016 static int
1017 CV_CDECL icvCmpSampleResponsePairs( const void* a, const void* b )
1018 {
1019     int ra = ((const CvSampleResponsePair*)a)->response;
1020     int rb = ((const CvSampleResponsePair*)b)->response;
1021     int ia = ((const CvSampleResponsePair*)a)->index;
1022     int ib = ((const CvSampleResponsePair*)b)->index;
1023
1024     return ra < rb ? -1 : ra > rb ? 1 : ia - ib;
1025     //return (ra > rb ? -1 : 0)|(ra < rb);
1026 }
1027
1028
1029 void
1030 cvSortSamplesByClasses( const float** samples, const CvMat* classes,
1031                         int* class_ranges, const uchar** mask )
1032 {
1033     CvSampleResponsePair* pairs = 0;
1034     CV_FUNCNAME( "cvSortSamplesByClasses" );
1035
1036     __BEGIN__;
1037
1038     int i, k = 0, sample_count;
1039
1040     if( !samples || !classes || !class_ranges )
1041         CV_ERROR( CV_StsNullPtr, "INTERNAL ERROR: some of the args are NULL pointers" );
1042
1043     if( classes->rows != 1 || CV_MAT_TYPE(classes->type) != CV_32SC1 )
1044         CV_ERROR( CV_StsBadArg, "classes array must be a single row of integers" );
1045
1046     sample_count = classes->cols;
1047     CV_CALL( pairs = (CvSampleResponsePair*)cvAlloc( (sample_count+1)*sizeof(pairs[0])));
1048
1049     for( i = 0; i < sample_count; i++ )
1050     {
1051         pairs[i].sample = samples[i];
1052         pairs[i].mask = (mask) ? (mask[i]) : 0;
1053         pairs[i].response = classes->data.i[i];
1054         pairs[i].index = i;
1055         assert( classes->data.i[i] >= 0 );
1056     }
1057
1058     qsort( pairs, sample_count, sizeof(pairs[0]), icvCmpSampleResponsePairs );
1059     pairs[sample_count].response = -1;
1060     class_ranges[0] = 0;
1061
1062     for( i = 0; i < sample_count; i++ )
1063     {
1064         samples[i] = pairs[i].sample;
1065         if (mask)
1066             mask[i] = pairs[i].mask;
1067         classes->data.i[i] = pairs[i].response;
1068
1069         if( pairs[i].response != pairs[i+1].response )
1070             class_ranges[++k] = i+1;
1071     }
1072
1073     __END__;
1074
1075     cvFree( &pairs );
1076 }
1077
1078
1079 void
1080 cvPreparePredictData( const CvArr* _sample, int dims_all,
1081                       const CvMat* comp_idx, int class_count,
1082                       const CvMat* prob, float** _row_sample,
1083                       int as_sparse )
1084 {
1085     float* row_sample = 0;
1086     int* inverse_comp_idx = 0;
1087
1088     CV_FUNCNAME( "cvPreparePredictData" );
1089
1090     __BEGIN__;
1091
1092     const CvMat* sample = (const CvMat*)_sample;
1093     float* sample_data;
1094     int sample_step;
1095     int is_sparse = CV_IS_SPARSE_MAT(sample);
1096     int d, sizes[CV_MAX_DIM];
1097     int i, dims_selected;
1098     int vec_size;
1099
1100     if( !is_sparse && !CV_IS_MAT(sample) )
1101         CV_ERROR( !sample ? CV_StsNullPtr : CV_StsBadArg, "The sample is not a valid vector" );
1102
1103     if( cvGetElemType( sample ) != CV_32FC1 )
1104         CV_ERROR( CV_StsUnsupportedFormat, "Input sample must have 32fC1 type" );
1105
1106     CV_CALL( d = cvGetDims( sample, sizes ));
1107
1108     if( !((is_sparse && d == 1) || (!is_sparse && d == 2 && (sample->rows == 1 || sample->cols == 1))) )
1109         CV_ERROR( CV_StsBadSize, "Input sample must be 1-dimensional vector" );
1110
1111     if( d == 1 )
1112         sizes[1] = 1;
1113
1114     if( sizes[0] + sizes[1] - 1 != dims_all )
1115         CV_ERROR( CV_StsUnmatchedSizes,
1116         "The sample size is different from what has been used for training" );
1117
1118     if( !_row_sample )
1119         CV_ERROR( CV_StsNullPtr, "INTERNAL ERROR: The row_sample pointer is NULL" );
1120
1121     if( comp_idx && (!CV_IS_MAT(comp_idx) || comp_idx->rows != 1 ||
1122         CV_MAT_TYPE(comp_idx->type) != CV_32SC1) )
1123         CV_ERROR( CV_StsBadArg, "INTERNAL ERROR: invalid comp_idx" );
1124
1125     dims_selected = comp_idx ? comp_idx->cols : dims_all;
1126
1127     if( prob )
1128     {
1129         if( !CV_IS_MAT(prob) )
1130             CV_ERROR( CV_StsBadArg, "The output matrix of probabilities is invalid" );
1131
1132         if( (prob->rows != 1 && prob->cols != 1) ||
1133             (CV_MAT_TYPE(prob->type) != CV_32FC1 &&
1134             CV_MAT_TYPE(prob->type) != CV_64FC1) )
1135             CV_ERROR( CV_StsBadSize,
1136             "The matrix of probabilities must be 1-dimensional vector of 32fC1 type" );
1137
1138         if( prob->rows + prob->cols - 1 != class_count )
1139             CV_ERROR( CV_StsUnmatchedSizes,
1140             "The vector of probabilities must contain as many elements as "
1141             "the number of classes in the training set" );
1142     }
1143
1144     vec_size = !as_sparse ? dims_selected*sizeof(row_sample[0]) :
1145                 (dims_selected + 1)*sizeof(CvSparseVecElem32f);
1146
1147     if( CV_IS_MAT(sample) )
1148     {
1149         sample_data = sample->data.fl;
1150         sample_step = sample->step / sizeof(row_sample[0]);
1151
1152         if( !comp_idx && CV_IS_MAT_CONT(sample->type) && !as_sparse )
1153             *_row_sample = sample_data;
1154         else
1155         {
1156             CV_CALL( row_sample = (float*)cvAlloc( vec_size ));
1157
1158             if( !comp_idx )
1159                 for( i = 0; i < dims_selected; i++ )
1160                     row_sample[i] = sample_data[sample_step*i];
1161             else
1162             {
1163                 int* comp = comp_idx->data.i;
1164                 if( !sample_step )
1165                     for( i = 0; i < dims_selected; i++ )
1166                         row_sample[i] = sample_data[comp[i]];
1167                 else
1168                     for( i = 0; i < dims_selected; i++ )
1169                         row_sample[i] = sample_data[sample_step*comp[i]];
1170             }
1171
1172             *_row_sample = row_sample;
1173         }
1174
1175         if( as_sparse )
1176         {
1177             const float* src = (const float*)row_sample;
1178             CvSparseVecElem32f* dst = (CvSparseVecElem32f*)row_sample;
1179
1180             dst[dims_selected].idx = -1;
1181             for( i = dims_selected - 1; i >= 0; i-- )
1182             {
1183                 dst[i].idx = i;
1184                 dst[i].val = src[i];
1185             }
1186         }
1187     }
1188     else
1189     {
1190         CvSparseNode* node;
1191         CvSparseMatIterator mat_iterator;
1192         const CvSparseMat* sparse = (const CvSparseMat*)sample;
1193         assert( is_sparse );
1194
1195         node = cvInitSparseMatIterator( sparse, &mat_iterator );
1196         CV_CALL( row_sample = (float*)cvAlloc( vec_size ));
1197
1198         if( comp_idx )
1199         {
1200             CV_CALL( inverse_comp_idx = (int*)cvAlloc( dims_all*sizeof(int) ));
1201             memset( inverse_comp_idx, -1, dims_all*sizeof(int) );
1202             for( i = 0; i < dims_selected; i++ )
1203                 inverse_comp_idx[comp_idx->data.i[i]] = i;
1204         }
1205
1206         if( !as_sparse )
1207         {
1208             memset( row_sample, 0, vec_size );
1209
1210             for( ; node != 0; node = cvGetNextSparseNode(&mat_iterator) )
1211             {
1212                 int idx = *CV_NODE_IDX( sparse, node );
1213                 if( inverse_comp_idx )
1214                 {
1215                     idx = inverse_comp_idx[idx];
1216                     if( idx < 0 )
1217                         continue;
1218                 }
1219                 row_sample[idx] = *(float*)CV_NODE_VAL( sparse, node );
1220             }
1221         }
1222         else
1223         {
1224             CvSparseVecElem32f* ptr = (CvSparseVecElem32f*)row_sample;
1225
1226             for( ; node != 0; node = cvGetNextSparseNode(&mat_iterator) )
1227             {
1228                 int idx = *CV_NODE_IDX( sparse, node );
1229                 if( inverse_comp_idx )
1230                 {
1231                     idx = inverse_comp_idx[idx];
1232                     if( idx < 0 )
1233                         continue;
1234                 }
1235                 ptr->idx = idx;
1236                 ptr->val = *(float*)CV_NODE_VAL( sparse, node );
1237                 ptr++;
1238             }
1239
1240             qsort( row_sample, ptr - (CvSparseVecElem32f*)row_sample,
1241                    sizeof(ptr[0]), icvCmpSparseVecElems );
1242             ptr->idx = -1;
1243         }
1244
1245         *_row_sample = row_sample;
1246     }
1247
1248     __END__;
1249
1250     if( inverse_comp_idx )
1251         cvFree( &inverse_comp_idx );
1252
1253     if( cvGetErrStatus() < 0 && _row_sample )
1254     {
1255         cvFree( &row_sample );
1256         *_row_sample = 0;
1257     }
1258 }
1259
1260
1261 static void
1262 icvConvertDataToSparse( const uchar* src, int src_step, int src_type,
1263                         uchar* dst, int dst_step, int dst_type,
1264                         CvSize size, int* idx )
1265 {
1266     CV_FUNCNAME( "icvConvertDataToSparse" );
1267
1268     __BEGIN__;
1269
1270     int i, j;
1271     src_type = CV_MAT_TYPE(src_type);
1272     dst_type = CV_MAT_TYPE(dst_type);
1273
1274     if( CV_MAT_CN(src_type) != 1 || CV_MAT_CN(dst_type) != 1 )
1275         CV_ERROR( CV_StsUnsupportedFormat, "The function supports only single-channel arrays" );
1276
1277     if( src_step == 0 )
1278         src_step = CV_ELEM_SIZE(src_type);
1279
1280     if( dst_step == 0 )
1281         dst_step = CV_ELEM_SIZE(dst_type);
1282
1283     // if there is no "idx" and if both arrays are continuous,
1284     // do the whole processing (copying or conversion) in a single loop
1285     if( !idx && CV_ELEM_SIZE(src_type)*size.width == src_step &&
1286         CV_ELEM_SIZE(dst_type)*size.width == dst_step )
1287     {
1288         size.width *= size.height;
1289         size.height = 1;
1290     }
1291
1292     if( src_type == dst_type )
1293     {
1294         int full_width = CV_ELEM_SIZE(dst_type)*size.width;
1295
1296         if( full_width == sizeof(int) ) // another common case: copy int's or float's
1297             for( i = 0; i < size.height; i++, src += src_step )
1298                 *(int*)(dst + dst_step*(idx ? idx[i] : i)) = *(int*)src;
1299         else
1300             for( i = 0; i < size.height; i++, src += src_step )
1301                 memcpy( dst + dst_step*(idx ? idx[i] : i), src, full_width );
1302     }
1303     else if( src_type == CV_32SC1 && (dst_type == CV_32FC1 || dst_type == CV_64FC1) )
1304         for( i = 0; i < size.height; i++, src += src_step )
1305         {
1306             uchar* _dst = dst + dst_step*(idx ? idx[i] : i);
1307             if( dst_type == CV_32FC1 )
1308                 for( j = 0; j < size.width; j++ )
1309                     ((float*)_dst)[j] = (float)((int*)src)[j];
1310             else
1311                 for( j = 0; j < size.width; j++ )
1312                     ((double*)_dst)[j] = ((int*)src)[j];
1313         }
1314     else if( (src_type == CV_32FC1 || src_type == CV_64FC1) && dst_type == CV_32SC1 )
1315         for( i = 0; i < size.height; i++, src += src_step )
1316         {
1317             uchar* _dst = dst + dst_step*(idx ? idx[i] : i);
1318             if( src_type == CV_32FC1 )
1319                 for( j = 0; j < size.width; j++ )
1320                     ((int*)_dst)[j] = cvRound(((float*)src)[j]);
1321             else
1322                 for( j = 0; j < size.width; j++ )
1323                     ((int*)_dst)[j] = cvRound(((double*)src)[j]);
1324         }
1325     else if( (src_type == CV_32FC1 && dst_type == CV_64FC1) ||
1326              (src_type == CV_64FC1 && dst_type == CV_32FC1) )
1327         for( i = 0; i < size.height; i++, src += src_step )
1328         {
1329             uchar* _dst = dst + dst_step*(idx ? idx[i] : i);
1330             if( src_type == CV_32FC1 )
1331                 for( j = 0; j < size.width; j++ )
1332                     ((double*)_dst)[j] = ((float*)src)[j];
1333             else
1334                 for( j = 0; j < size.width; j++ )
1335                     ((float*)_dst)[j] = (float)((double*)src)[j];
1336         }
1337     else
1338         CV_ERROR( CV_StsUnsupportedFormat, "Unsupported combination of input and output vectors" );
1339
1340     __END__;
1341 }
1342
1343
1344 void
1345 cvWritebackLabels( const CvMat* labels, CvMat* dst_labels,
1346                    const CvMat* centers, CvMat* dst_centers,
1347                    const CvMat* probs, CvMat* dst_probs,
1348                    const CvMat* sample_idx, int samples_all,
1349                    const CvMat* comp_idx, int dims_all )
1350 {
1351     CV_FUNCNAME( "cvWritebackLabels" );
1352
1353     __BEGIN__;
1354
1355     int samples_selected = samples_all, dims_selected = dims_all;
1356
1357     if( dst_labels && !CV_IS_MAT(dst_labels) )
1358         CV_ERROR( CV_StsBadArg, "Array of output labels is not a valid matrix" );
1359
1360     if( dst_centers )
1361         if( !ICV_IS_MAT_OF_TYPE(dst_centers, CV_32FC1) &&
1362             !ICV_IS_MAT_OF_TYPE(dst_centers, CV_64FC1) )
1363             CV_ERROR( CV_StsBadArg, "Array of cluster centers is not a valid matrix" );
1364
1365     if( dst_probs && !CV_IS_MAT(dst_probs) )
1366         CV_ERROR( CV_StsBadArg, "Probability matrix is not valid" );
1367
1368     if( sample_idx )
1369     {
1370         CV_ASSERT( sample_idx->rows == 1 && CV_MAT_TYPE(sample_idx->type) == CV_32SC1 );
1371         samples_selected = sample_idx->cols;
1372     }
1373
1374     if( comp_idx )
1375     {
1376         CV_ASSERT( comp_idx->rows == 1 && CV_MAT_TYPE(comp_idx->type) == CV_32SC1 );
1377         dims_selected = comp_idx->cols;
1378     }
1379
1380     if( dst_labels && (!labels || labels->data.ptr != dst_labels->data.ptr) )
1381     {
1382         if( !labels )
1383             CV_ERROR( CV_StsNullPtr, "NULL labels" );
1384
1385         CV_ASSERT( labels->rows == 1 );
1386
1387         if( dst_labels->rows != 1 && dst_labels->cols != 1 )
1388             CV_ERROR( CV_StsBadSize, "Array of output labels should be 1d vector" );
1389
1390         if( dst_labels->rows + dst_labels->cols - 1 != samples_all )
1391             CV_ERROR( CV_StsUnmatchedSizes,
1392             "Size of vector of output labels is not equal to the total number of input samples" );
1393
1394         CV_ASSERT( labels->cols == samples_selected );
1395
1396         CV_CALL( icvConvertDataToSparse( labels->data.ptr, labels->step, labels->type,
1397                         dst_labels->data.ptr, dst_labels->step, dst_labels->type,
1398                         cvSize( 1, samples_selected ), sample_idx ? sample_idx->data.i : 0 ));
1399     }
1400
1401     if( dst_centers && (!centers || centers->data.ptr != dst_centers->data.ptr) )
1402     {
1403         int i;
1404
1405         if( !centers )
1406             CV_ERROR( CV_StsNullPtr, "NULL centers" );
1407
1408         if( centers->rows != dst_centers->rows )
1409             CV_ERROR( CV_StsUnmatchedSizes, "Invalid number of rows in matrix of output centers" );
1410
1411         if( dst_centers->cols != dims_all )
1412             CV_ERROR( CV_StsUnmatchedSizes,
1413             "Number of columns in matrix of output centers is "
1414             "not equal to the total number of components in the input samples" );
1415
1416         CV_ASSERT( centers->cols == dims_selected );
1417
1418         for( i = 0; i < centers->rows; i++ )
1419             CV_CALL( icvConvertDataToSparse( centers->data.ptr + i*centers->step, 0, centers->type,
1420                         dst_centers->data.ptr + i*dst_centers->step, 0, dst_centers->type,
1421                         cvSize( 1, dims_selected ), comp_idx ? comp_idx->data.i : 0 ));
1422     }
1423
1424     if( dst_probs && (!probs || probs->data.ptr != dst_probs->data.ptr) )
1425     {
1426         if( !probs )
1427             CV_ERROR( CV_StsNullPtr, "NULL probs" );
1428
1429         if( probs->cols != dst_probs->cols )
1430             CV_ERROR( CV_StsUnmatchedSizes, "Invalid number of columns in output probability matrix" );
1431
1432         if( dst_probs->rows != samples_all )
1433             CV_ERROR( CV_StsUnmatchedSizes,
1434             "Number of rows in output probability matrix is "
1435             "not equal to the total number of input samples" );
1436
1437         CV_ASSERT( probs->rows == samples_selected );
1438
1439         CV_CALL( icvConvertDataToSparse( probs->data.ptr, probs->step, probs->type,
1440                         dst_probs->data.ptr, dst_probs->step, dst_probs->type,
1441                         cvSize( probs->cols, samples_selected ),
1442                         sample_idx ? sample_idx->data.i : 0 ));
1443     }
1444
1445     __END__;
1446 }
1447
1448 #if 0
1449 CV_IMPL void
1450 cvStatModelMultiPredict( const CvStatModel* stat_model,
1451                          const CvArr* predict_input,
1452                          int flags, CvMat* predict_output,
1453                          CvMat* probs, const CvMat* sample_idx )
1454 {
1455     CvMemStorage* storage = 0;
1456     CvMat* sample_idx_buffer = 0;
1457     CvSparseMat** sparse_rows = 0;
1458     int samples_selected = 0;
1459
1460     CV_FUNCNAME( "cvStatModelMultiPredict" );
1461
1462     __BEGIN__;
1463
1464     int i;
1465     int predict_output_step = 1, sample_idx_step = 1;
1466     int type;
1467     int d, sizes[CV_MAX_DIM];
1468     int tflag = flags == CV_COL_SAMPLE;
1469     int samples_all, dims_all;
1470     int is_sparse = CV_IS_SPARSE_MAT(predict_input);
1471     CvMat predict_input_part;
1472     CvArr* sample = &predict_input_part;
1473     CvMat probs_part;
1474     CvMat* probs1 = probs ? &probs_part : 0;
1475
1476     if( !CV_IS_STAT_MODEL(stat_model) )
1477         CV_ERROR( !stat_model ? CV_StsNullPtr : CV_StsBadArg, "Invalid statistical model" );
1478
1479     if( !stat_model->predict )
1480         CV_ERROR( CV_StsNotImplemented, "There is no \"predict\" method" );
1481
1482     if( !predict_input || !predict_output )
1483         CV_ERROR( CV_StsNullPtr, "NULL input or output matrices" );
1484
1485     if( !is_sparse && !CV_IS_MAT(predict_input) )
1486         CV_ERROR( CV_StsBadArg, "predict_input should be a matrix or a sparse matrix" );
1487
1488     if( !CV_IS_MAT(predict_output) )
1489         CV_ERROR( CV_StsBadArg, "predict_output should be a matrix" );
1490
1491     type = cvGetElemType( predict_input );
1492     if( type != CV_32FC1 ||
1493         (CV_MAT_TYPE(predict_output->type) != CV_32FC1 &&
1494          CV_MAT_TYPE(predict_output->type) != CV_32SC1 ))
1495          CV_ERROR( CV_StsUnsupportedFormat, "The input or output matrix has unsupported format" );
1496
1497     CV_CALL( d = cvGetDims( predict_input, sizes ));
1498     if( d > 2 )
1499         CV_ERROR( CV_StsBadSize, "The input matrix should be 1- or 2-dimensional" );
1500
1501     if( !tflag )
1502     {
1503         samples_all = samples_selected = sizes[0];
1504         dims_all = sizes[1];
1505     }
1506     else
1507     {
1508         samples_all = samples_selected = sizes[1];
1509         dims_all = sizes[0];
1510     }
1511
1512     if( sample_idx )
1513     {
1514         if( !CV_IS_MAT(sample_idx) )
1515             CV_ERROR( CV_StsBadArg, "Invalid sample_idx matrix" );
1516
1517         if( sample_idx->cols != 1 && sample_idx->rows != 1 )
1518             CV_ERROR( CV_StsBadSize, "sample_idx must be 1-dimensional matrix" );
1519
1520         samples_selected = sample_idx->rows + sample_idx->cols - 1;
1521
1522         if( CV_MAT_TYPE(sample_idx->type) == CV_32SC1 )
1523         {
1524             if( samples_selected > samples_all )
1525                 CV_ERROR( CV_StsBadSize, "sample_idx is too large vector" );
1526         }
1527         else if( samples_selected != samples_all )
1528             CV_ERROR( CV_StsUnmatchedSizes, "sample_idx has incorrect size" );
1529
1530         sample_idx_step = sample_idx->step ?
1531             sample_idx->step / CV_ELEM_SIZE(sample_idx->type) : 1;
1532     }
1533
1534     if( predict_output->rows != 1 && predict_output->cols != 1 )
1535         CV_ERROR( CV_StsBadSize, "predict_output should be a 1-dimensional matrix" );
1536
1537     if( predict_output->rows + predict_output->cols - 1 != samples_all )
1538         CV_ERROR( CV_StsUnmatchedSizes, "predict_output and predict_input have uncoordinated sizes" );
1539
1540     predict_output_step = predict_output->step ?
1541         predict_output->step / CV_ELEM_SIZE(predict_output->type) : 1;
1542
1543     if( probs )
1544     {
1545         if( !CV_IS_MAT(probs) )
1546             CV_ERROR( CV_StsBadArg, "Invalid matrix of probabilities" );
1547
1548         if( probs->rows != samples_all )
1549             CV_ERROR( CV_StsUnmatchedSizes,
1550             "matrix of probabilities must have as many rows as the total number of samples" );
1551
1552         if( CV_MAT_TYPE(probs->type) != CV_32FC1 )
1553             CV_ERROR( CV_StsUnsupportedFormat, "matrix of probabilities must have 32fC1 type" );
1554     }
1555
1556     if( is_sparse )
1557     {
1558         CvSparseNode* node;
1559         CvSparseMatIterator mat_iterator;
1560         CvSparseMat* sparse = (CvSparseMat*)predict_input;
1561
1562         if( sample_idx && CV_MAT_TYPE(sample_idx->type) == CV_32SC1 )
1563         {
1564             CV_CALL( sample_idx_buffer = cvCreateMat( 1, samples_all, CV_8UC1 ));
1565             cvZero( sample_idx_buffer );
1566             for( i = 0; i < samples_selected; i++ )
1567                 sample_idx_buffer->data.ptr[sample_idx->data.i[i*sample_idx_step]] = 1;
1568             samples_selected = samples_all;
1569             sample_idx = sample_idx_buffer;
1570             sample_idx_step = 1;
1571         }
1572
1573         CV_CALL( sparse_rows = (CvSparseMat**)cvAlloc( samples_selected*sizeof(sparse_rows[0])));
1574         for( i = 0; i < samples_selected; i++ )
1575         {
1576             if( sample_idx && sample_idx->data.ptr[i*sample_idx_step] == 0 )
1577                 continue;
1578             CV_CALL( sparse_rows[i] = cvCreateSparseMat( 1, &dims_all, type ));
1579             if( !storage )
1580                 storage = sparse_rows[i]->heap->storage;
1581             else
1582             {
1583                 // hack: to decrease memory footprint, make all the sparse matrices
1584                 // reside in the same storage
1585                 int elem_size = sparse_rows[i]->heap->elem_size;
1586                 cvReleaseMemStorage( &sparse_rows[i]->heap->storage );
1587                 sparse_rows[i]->heap = cvCreateSet( 0, sizeof(CvSet), elem_size, storage );
1588             }
1589         }
1590
1591         // put each row (or column) of predict_input into separate sparse matrix.
1592         node = cvInitSparseMatIterator( sparse, &mat_iterator );
1593         for( ; node != 0; node = cvGetNextSparseNode( &mat_iterator ))
1594         {
1595             int* idx = CV_NODE_IDX( sparse, node );
1596             int idx0 = idx[tflag ^ 1];
1597             int idx1 = idx[tflag];
1598
1599             if( sample_idx && sample_idx->data.ptr[idx0*sample_idx_step] == 0 )
1600                 continue;
1601
1602             assert( sparse_rows[idx0] != 0 );
1603             *(float*)cvPtrND( sparse, &idx1, 0, 1, 0 ) = *(float*)CV_NODE_VAL( sparse, node );
1604         }
1605     }
1606
1607     for( i = 0; i < samples_selected; i++ )
1608     {
1609         int idx = i;
1610         float response;
1611
1612         if( sample_idx )
1613         {
1614             if( CV_MAT_TYPE(sample_idx->type) == CV_32SC1 )
1615             {
1616                 idx = sample_idx->data.i[i*sample_idx_step];
1617                 if( (unsigned)idx >= (unsigned)samples_all )
1618                     CV_ERROR( CV_StsOutOfRange, "Some of sample_idx elements are out of range" );
1619             }
1620             else if( CV_MAT_TYPE(sample_idx->type) == CV_8UC1 &&
1621                      sample_idx->data.ptr[i*sample_idx_step] == 0 )
1622                 continue;
1623         }
1624
1625         if( !is_sparse )
1626         {
1627             if( !tflag )
1628                 cvGetRow( predict_input, &predict_input_part, idx );
1629             else
1630             {
1631                 cvGetCol( predict_input, &predict_input_part, idx );
1632             }
1633         }
1634         else
1635             sample = sparse_rows[idx];
1636
1637         if( probs )
1638             cvGetRow( probs, probs1, idx );
1639
1640         CV_CALL( response = stat_model->predict( stat_model, (CvMat*)sample, probs1 ));
1641
1642         if( CV_MAT_TYPE(predict_output->type) == CV_32FC1 )
1643             predict_output->data.fl[idx*predict_output_step] = response;
1644         else
1645         {
1646             CV_ASSERT( cvRound(response) == response );
1647             predict_output->data.i[idx*predict_output_step] = cvRound(response);
1648         }
1649     }
1650
1651     __END__;
1652
1653     if( sparse_rows )
1654     {
1655         int i;
1656         for( i = 0; i < samples_selected; i++ )
1657             if( sparse_rows[i] )
1658             {
1659                 sparse_rows[i]->heap->storage = 0;
1660                 cvReleaseSparseMat( &sparse_rows[i] );
1661             }
1662         cvFree( &sparse_rows );
1663     }
1664
1665     cvReleaseMat( &sample_idx_buffer );
1666     cvReleaseMemStorage( &storage );
1667 }
1668 #endif
1669
1670 // By P. Yarykin - begin -
1671
1672 void cvCombineResponseMaps (CvMat*  _responses,
1673                       const CvMat*  old_response_map,
1674                             CvMat*  new_response_map,
1675                             CvMat** out_response_map)
1676 {
1677     int** old_data = NULL;
1678     int** new_data = NULL;
1679
1680         CV_FUNCNAME ("cvCombineResponseMaps");
1681         __BEGIN__
1682
1683     int i,j;
1684     int old_n, new_n, out_n;
1685     int samples, free_response;
1686     int* first;
1687     int* responses;
1688     int* out_data;
1689
1690     if( out_response_map )
1691         *out_response_map = 0;
1692
1693 // Check input data.
1694     if ((!ICV_IS_MAT_OF_TYPE (_responses, CV_32SC1)) ||
1695         (!ICV_IS_MAT_OF_TYPE (old_response_map, CV_32SC1)) ||
1696         (!ICV_IS_MAT_OF_TYPE (new_response_map, CV_32SC1)))
1697     {
1698         CV_ERROR (CV_StsBadArg, "Some of input arguments is not the CvMat")
1699     }
1700
1701 // Prepare sorted responses.
1702     first = new_response_map->data.i;
1703     new_n = new_response_map->cols;
1704     CV_CALL (new_data = (int**)cvAlloc (new_n * sizeof (new_data[0])));
1705     for (i = 0; i < new_n; i++)
1706         new_data[i] = first + i;
1707     qsort (new_data, new_n, sizeof(int*), icvCmpIntegersPtr);
1708
1709     first = old_response_map->data.i;
1710     old_n = old_response_map->cols;
1711     CV_CALL (old_data = (int**)cvAlloc (old_n * sizeof (old_data[0])));
1712     for (i = 0; i < old_n; i++)
1713         old_data[i] = first + i;
1714     qsort (old_data, old_n, sizeof(int*), icvCmpIntegersPtr);
1715
1716 // Count the number of different responses.
1717     for (i = 0, j = 0, out_n = 0; i < old_n && j < new_n; out_n++)
1718     {
1719         if (*old_data[i] == *new_data[j])
1720         {
1721             i++;
1722             j++;
1723         }
1724         else if (*old_data[i] < *new_data[j])
1725             i++;
1726         else
1727             j++;
1728     }
1729     out_n += old_n - i + new_n - j;
1730
1731 // Create and fill the result response maps.
1732     CV_CALL (*out_response_map = cvCreateMat (1, out_n, CV_32SC1));
1733     out_data = (*out_response_map)->data.i;
1734     memcpy (out_data, first, old_n * sizeof (int));
1735
1736     free_response = old_n;
1737     for (i = 0, j = 0; i < old_n && j < new_n; )
1738     {
1739         if (*old_data[i] == *new_data[j])
1740         {
1741             *new_data[j] = (int)(old_data[i] - first);
1742             i++;
1743             j++;
1744         }
1745         else if (*old_data[i] < *new_data[j])
1746             i++;
1747         else
1748         {
1749             out_data[free_response] = *new_data[j];
1750             *new_data[j] = free_response++;
1751             j++;
1752         }
1753     }
1754     for (; j < new_n; j++)
1755     {
1756         out_data[free_response] = *new_data[j];
1757         *new_data[j] = free_response++;
1758     }
1759     CV_ASSERT (free_response == out_n);
1760
1761 // Change <responses> according to out response map.
1762     samples = _responses->cols + _responses->rows - 1;
1763     responses = _responses->data.i;
1764     first = new_response_map->data.i;
1765     for (i = 0; i < samples; i++)
1766     {
1767         responses[i] = first[responses[i]];
1768     }
1769
1770         __END__
1771
1772     cvFree(&old_data);
1773     cvFree(&new_data);
1774
1775 }
1776
1777
1778 int icvGetNumberOfCluster( double* prob_vector, int num_of_clusters, float r,
1779                            float outlier_thresh, int normalize_probs )
1780 {
1781     int max_prob_loc = 0;
1782
1783     CV_FUNCNAME("icvGetNumberOfCluster");
1784     __BEGIN__;
1785
1786     double prob, maxprob, sum;
1787     int i;
1788
1789     CV_ASSERT(prob_vector);
1790     CV_ASSERT(num_of_clusters >= 0);
1791
1792     maxprob = prob_vector[0];
1793     max_prob_loc = 0;
1794     sum = maxprob;
1795     for( i = 1; i < num_of_clusters; i++ )
1796     {
1797         prob = prob_vector[i];
1798         sum += prob;
1799         if( prob > maxprob )
1800         {
1801             max_prob_loc = i;
1802             maxprob = prob;
1803         }
1804     }
1805     if( normalize_probs && fabs(sum - 1.) > FLT_EPSILON )
1806     {
1807         for( i = 0; i < num_of_clusters; i++ )
1808             prob_vector[i] /= sum;
1809     }
1810     if( fabs(r - 1.) > FLT_EPSILON && fabs(sum - 1.) < outlier_thresh )
1811         max_prob_loc = -1;
1812
1813     __END__;
1814
1815     return max_prob_loc;
1816
1817 } // End of icvGetNumberOfCluster
1818
1819
1820 void icvFindClusterLabels( const CvMat* probs, float outlier_thresh, float r,
1821                           const CvMat* labels )
1822 {
1823     CvMat* counts = 0;
1824
1825     CV_FUNCNAME("icvFindClusterLabels");
1826     __BEGIN__;
1827
1828     int nclusters, nsamples;
1829     int i, j;
1830     double* probs_data;
1831
1832     CV_ASSERT( ICV_IS_MAT_OF_TYPE(probs, CV_64FC1) );
1833     CV_ASSERT( ICV_IS_MAT_OF_TYPE(labels, CV_32SC1) );
1834
1835     nclusters = probs->cols;
1836     nsamples  = probs->rows;
1837     CV_ASSERT( nsamples == labels->cols );
1838
1839     CV_CALL( counts = cvCreateMat( 1, nclusters + 1, CV_32SC1 ) );
1840     CV_CALL( cvSetZero( counts ));
1841     for( i = 0; i < nsamples; i++ )
1842     {
1843         labels->data.i[i] = icvGetNumberOfCluster( probs->data.db + i*probs->cols,
1844             nclusters, r, outlier_thresh, 1 );
1845         counts->data.i[labels->data.i[i] + 1]++;
1846     }
1847     CV_ASSERT((int)cvSum(counts).val[0] == nsamples);
1848     // Filling empty clusters with the vector, that has the maximal probability
1849     for( j = 0; j < nclusters; j++ ) // outliers are ignored
1850     {
1851         int maxprob_loc = -1;
1852         double maxprob = 0;
1853
1854         if( counts->data.i[j+1] ) // j-th class is not empty
1855             continue;
1856         // look for the presentative, which is not lonely in it's cluster
1857         // and that has a maximal probability among all these vectors
1858         probs_data = probs->data.db;
1859         for( i = 0; i < nsamples; i++, probs_data++ )
1860         {
1861             int label = labels->data.i[i];
1862             double prob;
1863             if( counts->data.i[label+1] == 0 ||
1864                 (counts->data.i[label+1] <= 1 && label != -1) )
1865                 continue;
1866             prob = *probs_data;
1867             if( prob >= maxprob )
1868             {
1869                 maxprob = prob;
1870                 maxprob_loc = i;
1871             }
1872         }
1873         // maxprob_loc == 0 <=> number of vectors less then number of clusters
1874         CV_ASSERT( maxprob_loc >= 0 );
1875         counts->data.i[labels->data.i[maxprob_loc] + 1]--;
1876         labels->data.i[maxprob_loc] = j;
1877         counts->data.i[j + 1]++;
1878     }
1879
1880     __END__;
1881
1882     cvReleaseMat( &counts );
1883 } // End of icvFindClusterLabels
1884
1885 /* End of file */