git.maemo.org Git - opencv/blob - tests/cxts/cxts_math.cpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 //
  10 //                        Intel License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
  14 // Third party copyrights are property of their respective owners.
  15 //
  16 // Redistribution and use in source and binary forms, with or without modification,
  17 // are permitted provided that the following conditions are met:
  18 //
  19 //   * Redistribution's of source code must retain the above copyright notice,
  20 //     this list of conditions and the following disclaimer.
  21 //
  22 //   * Redistribution's in binary form must reproduce the above copyright notice,
  23 //     this list of conditions and the following disclaimer in the documentation
  24 //     and/or other materials provided with the distribution.
  25 //
  26 //   * The name of Intel Corporation may not be used to endorse or promote products
  27 //     derived from this software without specific prior written permission.
  28 //
  29 // This software is provided by the copyright holders and contributors "as is" and
  30 // any express or implied warranties, including, but not limited to, the implied
  31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  32 // In no event shall the Intel Corporation or contributors be liable for any direct,
  33 // indirect, incidental, special, exemplary, or consequential damages
  34 // (including, but not limited to, procurement of substitute goods or services;
  35 // loss of use, data, or profits; or business interruption) however caused
  36 // and on any theory of liability, whether in contract, strict liability,
  37 // or tort (including negligence or otherwise) arising in any way out of
  38 // the use of this software, even if advised of the possibility of such damage.
  39 //
  40 //M*/
  41
  42 #include "_cxts.h"
  43
  44 /****************************************************************************************\
  45 *                                   Utility Functions                                    *
  46 \****************************************************************************************/
  47
  48 const char* cvTsGetTypeName( int type )
  49 {
  50     static const char* type_names[] = { "8u", "8s", "16u", "16s", "32s", "32f", "64f", "ptr" };
  51     return type_names[CV_MAT_DEPTH(type)];
  52 }
  53
  54
  55 int cvTsTypeByName( const char* name )
  56 {
  57     int i;
  58     for( i = 0; i < CV_DEPTH_MAX; i++ )
  59         if( strcmp(name, cvTsGetTypeName(i)) == 0 )
  60             return i;
  61     return -1;
  62 }
  63
  64
  65 void cvTsRandUni( CvRNG* rng, CvMat* a, CvScalar param0, CvScalar param1 )
  66 {
  67     int i, j, k, cn, ncols;
  68     CvScalar scale = param0;
  69     CvScalar delta = param1;
  70     double C = 1./(65536.*65536.);
  71
  72     cn = CV_MAT_CN(a->type);
  73     ncols = a->cols*cn;
  74
  75     for( k = 0; k < 4; k++ )
  76     {
  77         double s = scale.val[k] - delta.val[k];
  78         if( s >= 0 )
  79             scale.val[k] = s;
  80         else
  81         {
  82             delta.val[k] = scale.val[k];
  83             scale.val[k] = -s;
  84         }
  85         scale.val[k] *= C;
  86     }
  87
  88     for( i = 0; i < a->rows; i++ )
  89     {
  90         uchar* data = a->data.ptr + i*a->step;
  91
  92         switch( CV_MAT_DEPTH(a->type) )
  93         {
  94         case CV_8U:
  95             for( j = 0; j < ncols; j += cn )
  96                 for( k = 0; k < cn; k++ )
  97                 {
  98                     int val = cvFloor( cvTsRandInt(rng)*scale.val[k] + delta.val[k] );
  99                     ((uchar*)data)[j + k] = CV_CAST_8U(val);
 100                 }
 101             break;
 102         case CV_8S:
 103             for( j = 0; j < ncols; j += cn )
 104                 for( k = 0; k < cn; k++ )
 105                 {
 106                     int val = cvFloor( cvTsRandInt(rng)*scale.val[k] + delta.val[k] );
 107                     ((schar*)data)[j + k] = CV_CAST_8S(val);
 108                 }
 109             break;
 110         case CV_16U:
 111             for( j = 0; j < ncols; j += cn )
 112                 for( k = 0; k < cn; k++ )
 113                 {
 114                     int val = cvFloor( cvTsRandInt(rng)*scale.val[k] + delta.val[k] );
 115                     ((ushort*)data)[j + k] = CV_CAST_16U(val);
 116                 }
 117             break;
 118         case CV_16S:
 119             for( j = 0; j < ncols; j += cn )
 120                 for( k = 0; k < cn; k++ )
 121                 {
 122                     int val = cvFloor( cvTsRandInt(rng)*scale.val[k] + delta.val[k] );
 123                     ((short*)data)[j + k] = CV_CAST_16S(val);
 124                 }
 125             break;
 126         case CV_32S:
 127             for( j = 0; j < ncols; j += cn )
 128                 for( k = 0; k < cn; k++ )
 129                 {
 130                     int val = cvFloor( cvTsRandInt(rng)*scale.val[k] + delta.val[k] );
 131                     ((int*)data)[j + k] = val;
 132                 }
 133             break;
 134         case CV_32F:
 135             for( j = 0; j < ncols; j += cn )
 136                 for( k = 0; k < cn; k++ )
 137                 {
 138                     double val = cvTsRandInt(rng)*scale.val[k] + delta.val[k];
 139                     ((float*)data)[j + k] = (float)val;
 140                 }
 141             break;
 142         case CV_64F:
 143             for( j = 0; j < ncols; j += cn )
 144                 for( k = 0; k < cn; k++ )
 145                 {
 146                     double val = cvTsRandInt(rng);
 147                     val = (val + cvTsRandInt(rng)*C)*scale.val[k] + delta.val[k];
 148                     ((double*)data)[j + k] = val;
 149                 }
 150             break;
 151         default:
 152             assert(0);
 153             return;
 154         }
 155     }
 156 }
 157
 158
 159 void cvTsZero( CvMat* c, const CvMat* mask )
 160 {
 161     int i, j, elem_size = CV_ELEM_SIZE(c->type), width = c->cols;
 162
 163     for( i = 0; i < c->rows; i++ )
 164     {
 165         if( !mask )
 166             memset( c->data.ptr + i*c->step, 0, width*elem_size );
 167         else
 168         {
 169             const uchar* mrow = mask->data.ptr + mask->step*i;
 170             uchar* cptr = c->data.ptr + c->step*i;
 171             for( j = 0; j < width; j++, cptr += elem_size )
 172                 if( mrow[j] )
 173                     memset( cptr, 0, elem_size );
 174         }
 175     }
 176 }
 177
 178
 179 // initializes scaled identity matrix
 180 void cvTsSetIdentity( CvMat* c, CvScalar diag_value )
 181 {
 182     int i, width;
 183     cvTsZero( c );
 184     width = MIN(c->rows, c->cols);
 185     for( i = 0; i < width; i++ )
 186         cvSet2D( c, i, i, diag_value );
 187 }
 188
 189
 190 // copies selected region of one array to another array
 191 void cvTsCopy( const CvMat* a, CvMat* b, const CvMat* mask )
 192 {
 193     int i = 0, j = 0, k;
 194     int el_size, ncols;
 195
 196     el_size = CV_ELEM_SIZE(a->type);
 197     ncols = a->cols;
 198     if( mask )
 199     {
 200         assert( CV_ARE_SIZES_EQ(a,mask) &&
 201             (CV_MAT_TYPE(mask->type) == CV_8UC1 ||
 202              CV_MAT_TYPE(mask->type) == CV_8SC1 ));
 203     }
 204
 205     assert( CV_ARE_TYPES_EQ(a,b) && CV_ARE_SIZES_EQ(a,b) );
 206
 207     if( !mask )
 208         ncols *= el_size;
 209
 210     for( i = 0; i < a->rows; i++ )
 211     {
 212         uchar* a_data = a->data.ptr + a->step*i;
 213         uchar* b_data = b->data.ptr + b->step*i;
 214
 215         if( !mask )
 216             memcpy( b_data, a_data, ncols );
 217         else
 218         {
 219             uchar* m_data = mask->data.ptr + mask->step*i;
 220
 221             for( j = 0; j < ncols; j++, b_data += el_size, a_data += el_size )
 222             {
 223                 if( m_data[j] )
 224                 {
 225                     for( k = 0; k < el_size; k++ )
 226                         b_data[k] = a_data[k];
 227                 }
 228             }
 229         }
 230     }
 231 }
 232
 233
 234 void cvTsConvert( const CvMat* a, CvMat* b )
 235 {
 236     int i, j, ncols = b->cols*CV_MAT_CN(b->type);
 237     double* buf = 0;
 238
 239     assert( CV_ARE_SIZES_EQ(a,b) && CV_ARE_CNS_EQ(a,b) );
 240     buf = (double*)cvStackAlloc(ncols*sizeof(buf[0]));
 241
 242     for( i = 0; i < b->rows; i++ )
 243     {
 244         uchar* a_data = a->data.ptr + i*a->step;
 245         uchar* b_data = b->data.ptr + i*b->step;
 246
 247         switch( CV_MAT_DEPTH(a->type) )
 248         {
 249         case CV_8U:
 250             for( j = 0; j < ncols; j++ )
 251                 buf[j] = ((uchar*)a_data)[j];
 252             break;
 253         case CV_8S:
 254             for( j = 0; j < ncols; j++ )
 255                 buf[j] = ((schar*)a_data)[j];
 256             break;
 257         case CV_16U:
 258             for( j = 0; j < ncols; j++ )
 259                 buf[j] = ((ushort*)a_data)[j];
 260             break;
 261         case CV_16S:
 262             for( j = 0; j < ncols; j++ )
 263                 buf[j] = ((short*)a_data)[j];
 264             break;
 265         case CV_32S:
 266             for( j = 0; j < ncols; j++ )
 267                 buf[j] = ((int*)a_data)[j];
 268             break;
 269         case CV_32F:
 270             for( j = 0; j < ncols; j++ )
 271                 buf[j] = ((float*)a_data)[j];
 272             break;
 273         case CV_64F:
 274             for( j = 0; j < ncols; j++ )
 275                 buf[j] = ((double*)a_data)[j];
 276             break;
 277         default:
 278             assert(0);
 279             return;
 280         }
 281
 282         switch( CV_MAT_DEPTH(b->type) )
 283         {
 284         case CV_8U:
 285             for( j = 0; j < ncols; j++ )
 286             {
 287                 int val = cvRound(buf[j]);
 288                 ((uchar*)b_data)[j] = CV_CAST_8U(val);
 289             }
 290             break;
 291         case CV_8S:
 292             for( j = 0; j < ncols; j++ )
 293             {
 294                 int val = cvRound(buf[j]);
 295                 ((schar*)b_data)[j] = CV_CAST_8S(val);
 296             }
 297             break;
 298         case CV_16U:
 299             for( j = 0; j < ncols; j++ )
 300             {
 301                 int val = cvRound(buf[j]);
 302                 ((ushort*)b_data)[j] = CV_CAST_16U(val);
 303             }
 304             break;
 305         case CV_16S:
 306             for( j = 0; j < ncols; j++ )
 307             {
 308                 int val = cvRound(buf[j]);
 309                 ((short*)b_data)[j] = CV_CAST_16S(val);
 310             }
 311             break;
 312         case CV_32S:
 313             for( j = 0; j < ncols; j++ )
 314             {
 315                 int val = cvRound(buf[j]);
 316                 ((int*)b_data)[j] = CV_CAST_32S(val);
 317             }
 318             break;
 319         case CV_32F:
 320             for( j = 0; j < ncols; j++ )
 321                 ((float*)b_data)[j] = CV_CAST_32F(buf[j]);
 322             break;
 323         case CV_64F:
 324             for( j = 0; j < ncols; j++ )
 325                 ((double*)b_data)[j] = CV_CAST_64F(buf[j]);
 326             break;
 327         default:
 328             assert(0);
 329         }
 330     }
 331 }
 332
 333
 334 // extracts a single channel from a multi-channel array
 335 void cvTsExtract( const CvMat* a, CvMat* b, int coi )
 336 {
 337     int i = 0, j = 0, k;
 338     int el_size, el_size1, ncols;
 339
 340     el_size = CV_ELEM_SIZE(a->type);
 341     el_size1 = CV_ELEM_SIZE(b->type);
 342     ncols = a->cols;
 343
 344     assert( CV_ARE_DEPTHS_EQ(a,b) && CV_ARE_SIZES_EQ(a,b) &&
 345             (unsigned)coi < (unsigned)CV_MAT_CN(a->type) &&
 346             CV_MAT_CN(b->type) == 1 );
 347
 348     for( i = 0; i < a->rows; i++ )
 349     {
 350         uchar* a_data = a->data.ptr + a->step*i;
 351         uchar* b_data = b->data.ptr + b->step*i;
 352         a_data += el_size1*coi;
 353         for( j = 0; j < ncols; j++, b_data += el_size1, a_data += el_size )
 354         {
 355             for( k = 0; k < el_size1; k++ )
 356                 b_data[k] = a_data[k];
 357         }
 358     }
 359 }
 360
 361 // replaces a single channel in a multi-channel array
 362 void cvTsInsert( const CvMat* a, CvMat* b, int coi )
 363 {
 364     int i = 0, j = 0, k;
 365     int el_size, el_size1, ncols;
 366
 367     el_size = CV_ELEM_SIZE(b->type);
 368     el_size1 = CV_ELEM_SIZE(a->type);
 369     ncols = a->cols;
 370
 371     assert( CV_ARE_DEPTHS_EQ(a,b) && CV_ARE_SIZES_EQ(a,b) &&
 372             (unsigned)coi < (unsigned)CV_MAT_CN(b->type) &&
 373             CV_MAT_CN(a->type) == 1 );
 374
 375     for( i = 0; i < a->rows; i++ )
 376     {
 377         uchar* a_data = a->data.ptr + a->step*i;
 378         uchar* b_data = b->data.ptr + b->step*i;
 379         b_data += el_size1*coi;
 380         for( j = 0; j < ncols; j++, b_data += el_size, a_data += el_size1 )
 381         {
 382             for( k = 0; k < el_size1; k++ )
 383                 b_data[k] = a_data[k];
 384         }
 385     }
 386 }
 387
 388
 389 // c = alpha*a + beta*b + gamma
 390 void cvTsAdd( const CvMat* a, CvScalar alpha, const CvMat* b, CvScalar beta,
 391               CvScalar gamma, CvMat* c, int calc_abs )
 392 {
 393     int i, j, k, cn, ncols;
 394     double* buf = 0;
 395     double* alpha_buf = 0;
 396     double* beta_buf = 0;
 397     double* gamma_buf = 0;
 398
 399     if( !c )
 400     {
 401         assert(0);
 402         return;
 403     }
 404
 405     cn = CV_MAT_CN(c->type);
 406     ncols = c->cols;
 407
 408     if( !a )
 409     {
 410         a = b;
 411         alpha = beta;
 412         b = 0;
 413         beta = cvScalar(0);
 414     }
 415
 416     if( a )
 417     {
 418         assert( CV_ARE_SIZES_EQ(a,c) && CV_MAT_CN(a->type) == cn );
 419         buf = (double*)malloc( a->cols * cn * sizeof(buf[0]) );
 420         alpha_buf = (double*)malloc( a->cols * cn * sizeof(alpha_buf[0]) );
 421     }
 422
 423     if( b )
 424     {
 425         assert( CV_ARE_SIZES_EQ(b,c) && CV_MAT_CN(b->type) == cn );
 426         beta_buf = (double*)malloc( b->cols * cn * sizeof(beta_buf[0]) );
 427     }
 428
 429     ncols *= cn;
 430     gamma_buf = (double*)malloc( ncols * sizeof(gamma_buf[0]) );
 431     if( !buf )
 432         buf = gamma_buf;
 433
 434     if( !a && !b && calc_abs )
 435     {
 436         for( k = 0; k < cn; k++ )
 437             gamma.val[k] = fabs(gamma.val[k]);
 438     }
 439
 440     for( i = 0; i < 1 + (a != 0) + (b != 0); i++ )
 441     {
 442         double* scalar_buf = i == 0 ? gamma_buf : i == 1 ? alpha_buf : beta_buf;
 443         CvScalar scalar = i == 0 ? gamma : i == 1 ? alpha : beta;
 444         for( j = 0; j < ncols; j += cn )
 445             for( k = 0; k < cn; k++ )
 446                 scalar_buf[j + k] = scalar.val[k];
 447     }
 448
 449     for( i = 0; i < c->rows; i++ )
 450     {
 451         uchar* c_data = c->data.ptr + i*c->step;
 452
 453         if( a )
 454         {
 455             uchar* a_data = a->data.ptr + i*a->step;
 456
 457             switch( CV_MAT_DEPTH(a->type) )
 458             {
 459             case CV_8U:
 460                 for( j = 0; j < ncols; j++ )
 461                     buf[j] = ((uchar*)a_data)[j]*alpha_buf[j] + gamma_buf[j];
 462                 break;
 463             case CV_8S:
 464                 for( j = 0; j < ncols; j++ )
 465                     buf[j] = ((schar*)a_data)[j]*alpha_buf[j] + gamma_buf[j];
 466                 break;
 467             case CV_16U:
 468                 for( j = 0; j < ncols; j++ )
 469                     buf[j] = ((ushort*)a_data)[j]*alpha_buf[j] + gamma_buf[j];
 470                 break;
 471             case CV_16S:
 472                 for( j = 0; j < ncols; j++ )
 473                     buf[j] = ((short*)a_data)[j]*alpha_buf[j] + gamma_buf[j];
 474                 break;
 475             case CV_32S:
 476                 for( j = 0; j < ncols; j++ )
 477                     buf[j] = ((int*)a_data)[j]*alpha_buf[j] + gamma_buf[j];
 478                 break;
 479             case CV_32F:
 480                 for( j = 0; j < ncols; j++ )
 481                     buf[j] = ((float*)a_data)[j]*alpha_buf[j] + gamma_buf[j];
 482                 break;
 483             case CV_64F:
 484                 for( j = 0; j < ncols; j++ )
 485                     buf[j] = ((double*)a_data)[j]*alpha_buf[j] + gamma_buf[j];
 486                 break;
 487             default:
 488                 assert(0);
 489                 return;
 490             }
 491         }
 492
 493         if( b )
 494         {
 495             uchar* b_data = b->data.ptr + i*b->step;
 496
 497             switch( CV_MAT_DEPTH(b->type) )
 498             {
 499             case CV_8U:
 500                 for( j = 0; j < ncols; j++ )
 501                     buf[j] += ((uchar*)b_data)[j]*beta_buf[j];
 502                 break;
 503             case CV_8S:
 504                 for( j = 0; j < ncols; j++ )
 505                     buf[j] += ((schar*)b_data)[j]*beta_buf[j];
 506                 break;
 507             case CV_16U:
 508                 for( j = 0; j < ncols; j++ )
 509                     buf[j] += ((ushort*)b_data)[j]*beta_buf[j];
 510                 break;
 511             case CV_16S:
 512                 for( j = 0; j < ncols; j++ )
 513                     buf[j] += ((short*)b_data)[j]*beta_buf[j];
 514                 break;
 515             case CV_32S:
 516                 for( j = 0; j < ncols; j++ )
 517                     buf[j] += ((int*)b_data)[j]*beta_buf[j];
 518                 break;
 519             case CV_32F:
 520                 for( j = 0; j < ncols; j++ )
 521                     buf[j] += ((float*)b_data)[j]*beta_buf[j];
 522                 break;
 523             case CV_64F:
 524                 for( j = 0; j < ncols; j++ )
 525                     buf[j] += ((double*)b_data)[j]*beta_buf[j];
 526                 break;
 527             default:
 528                 assert(0);
 529                 return;
 530             }
 531         }
 532
 533         if( a || b )
 534         {
 535             if( calc_abs )
 536             {
 537                 for( j = 0; j < ncols; j++ )
 538                     buf[j] = fabs(buf[j]);
 539             }
 540         }
 541         else if( i > 0 )
 542         {
 543             memcpy( c_data, c_data - c->step, c->cols*CV_ELEM_SIZE(c->type) );
 544             continue;
 545         }
 546
 547         switch( CV_MAT_DEPTH(c->type) )
 548         {
 549         case CV_8U:
 550             for( j = 0; j < ncols; j++ )
 551             {
 552                 int val = cvRound(buf[j]);
 553                 ((uchar*)c_data)[j] = CV_CAST_8U(val);
 554             }
 555             break;
 556         case CV_8S:
 557             for( j = 0; j < ncols; j++ )
 558             {
 559                 int val = cvRound(buf[j]);
 560                 ((schar*)c_data)[j] = CV_CAST_8S(val);
 561             }
 562             break;
 563         case CV_16U:
 564             for( j = 0; j < ncols; j++ )
 565             {
 566                 int val = cvRound(buf[j]);
 567                 ((ushort*)c_data)[j] = CV_CAST_16U(val);
 568             }
 569             break;
 570         case CV_16S:
 571             for( j = 0; j < ncols; j++ )
 572             {
 573                 int val = cvRound(buf[j]);
 574                 ((short*)c_data)[j] = CV_CAST_16S(val);
 575             }
 576             break;
 577         case CV_32S:
 578             for( j = 0; j < ncols; j++ )
 579             {
 580                 int val = cvRound(buf[j]);
 581                 ((int*)c_data)[j] = CV_CAST_32S(val);
 582             }
 583             break;
 584         case CV_32F:
 585             for( j = 0; j < ncols; j++ )
 586                 ((float*)c_data)[j] = CV_CAST_32F(buf[j]);
 587             break;
 588         case CV_64F:
 589             for( j = 0; j < ncols; j++ )
 590                 ((double*)c_data)[j] = CV_CAST_64F(buf[j]);
 591             break;
 592         default:
 593             assert(0);
 594         }
 595     }
 596
 597     if( buf && buf != gamma_buf )
 598         free( buf );
 599     if( gamma_buf )
 600         free( gamma_buf );
 601     if( alpha_buf )
 602         free( alpha_buf );
 603     if( beta_buf )
 604         free( beta_buf );
 605 }
 606
 607
 608 // c = a*b*alpha
 609 void cvTsMul( const CvMat* a, const CvMat* b, CvScalar alpha, CvMat* c )
 610 {
 611     int i, j, k, cn, ncols;
 612     double* buf = 0;
 613     double* alpha_buf = 0;
 614
 615     if( !a || !b || !c )
 616     {
 617         assert(0);
 618         return;
 619     }
 620
 621     assert( CV_ARE_SIZES_EQ(a,c) && CV_ARE_SIZES_EQ(b,c) &&
 622             CV_ARE_TYPES_EQ(a,b) && CV_ARE_CNS_EQ(a,c) );
 623
 624     cn = CV_MAT_CN(c->type);
 625     ncols = c->cols * cn;
 626     alpha_buf = (double*)malloc( ncols * sizeof(alpha_buf[0]) );
 627     buf = (double*)malloc( ncols * sizeof(buf[0]) );
 628
 629     for( j = 0; j < ncols; j += cn )
 630         for( k = 0; k < cn; k++ )
 631             alpha_buf[j + k] = alpha.val[k];
 632
 633     for( i = 0; i < c->rows; i++ )
 634     {
 635         uchar* c_data = c->data.ptr + i*c->step;
 636         uchar* a_data = a->data.ptr + i*a->step;
 637         uchar* b_data = b->data.ptr + i*b->step;
 638
 639         switch( CV_MAT_DEPTH(a->type) )
 640         {
 641         case CV_8U:
 642             for( j = 0; j < ncols; j++ )
 643                 buf[j] = (alpha_buf[j]*((uchar*)a_data)[j])*((uchar*)b_data)[j];
 644             break;
 645         case CV_8S:
 646             for( j = 0; j < ncols; j++ )
 647                 buf[j] = (alpha_buf[j]*((schar*)a_data)[j])*((schar*)b_data)[j];
 648             break;
 649         case CV_16U:
 650             for( j = 0; j < ncols; j++ )
 651                 buf[j] = (alpha_buf[j]*((ushort*)a_data)[j])*((ushort*)b_data)[j];
 652             break;
 653         case CV_16S:
 654             for( j = 0; j < ncols; j++ )
 655                 buf[j] = (alpha_buf[j]*((short*)a_data)[j])*((short*)b_data)[j];
 656             break;
 657         case CV_32S:
 658             for( j = 0; j < ncols; j++ )
 659                 buf[j] = (alpha_buf[j]*((int*)a_data)[j])*((int*)b_data)[j];
 660             break;
 661         case CV_32F:
 662             for( j = 0; j < ncols; j++ )
 663                 buf[j] = (alpha_buf[j]*((float*)a_data)[j])*((float*)b_data)[j];
 664             break;
 665         case CV_64F:
 666             for( j = 0; j < ncols; j++ )
 667                 buf[j] = (alpha_buf[j]*((double*)a_data)[j])*((double*)b_data)[j];
 668             break;
 669         default:
 670             assert(0);
 671             return;
 672         }
 673
 674         switch( CV_MAT_DEPTH(c->type) )
 675         {
 676         case CV_8U:
 677             for( j = 0; j < ncols; j++ )
 678             {
 679                 int val = cvRound(buf[j]);
 680                 ((uchar*)c_data)[j] = CV_CAST_8U(val);
 681             }
 682             break;
 683         case CV_8S:
 684             for( j = 0; j < ncols; j++ )
 685             {
 686                 int val = cvRound(buf[j]);
 687                 ((schar*)c_data)[j] = CV_CAST_8S(val);
 688             }
 689             break;
 690         case CV_16U:
 691             for( j = 0; j < ncols; j++ )
 692             {
 693                 int val = cvRound(buf[j]);
 694                 ((ushort*)c_data)[j] = CV_CAST_16U(val);
 695             }
 696             break;
 697         case CV_16S:
 698             for( j = 0; j < ncols; j++ )
 699             {
 700                 int val = cvRound(buf[j]);
 701                 ((short*)c_data)[j] = CV_CAST_16S(val);
 702             }
 703             break;
 704         case CV_32S:
 705             for( j = 0; j < ncols; j++ )
 706             {
 707                 int val = cvRound(buf[j]);
 708                 ((int*)c_data)[j] = CV_CAST_32S(val);
 709             }
 710             break;
 711         case CV_32F:
 712             for( j = 0; j < ncols; j++ )
 713                 ((float*)c_data)[j] = CV_CAST_32F(buf[j]);
 714             break;
 715         case CV_64F:
 716             for( j = 0; j < ncols; j++ )
 717                 ((double*)c_data)[j] = CV_CAST_64F(buf[j]);
 718             break;
 719         default:
 720             assert(0);
 721         }
 722     }
 723
 724     if( alpha_buf )
 725         free( alpha_buf );
 726
 727     if( buf )
 728         free( buf );
 729 }
 730
 731
 732 // c = a*alpha/b
 733 void cvTsDiv( const CvMat* a, const CvMat* b, CvScalar alpha, CvMat* c )
 734 {
 735     int i, j, k, cn, ncols;
 736     double* buf = 0;
 737     double* alpha_buf = 0;
 738
 739     if( !b || !c )
 740     {
 741         assert(0);
 742         return;
 743     }
 744
 745     if( a )
 746     {
 747         assert( CV_ARE_SIZES_EQ(a,c) &&
 748                 CV_ARE_TYPES_EQ(a,b) && CV_ARE_CNS_EQ(a,c) );
 749     }
 750
 751     assert( CV_ARE_SIZES_EQ(b,c) && CV_ARE_CNS_EQ(b,c) );
 752
 753     cn = CV_MAT_CN(c->type);
 754     ncols = c->cols * cn;
 755     alpha_buf = (double*)malloc( ncols * sizeof(alpha_buf[0]) );
 756     buf = (double*)malloc( ncols * sizeof(buf[0]) );
 757
 758     for( j = 0; j < ncols; j += cn )
 759         for( k = 0; k < cn; k++ )
 760             alpha_buf[j + k] = alpha.val[k];
 761
 762     for( i = 0; i < c->rows; i++ )
 763     {
 764         uchar* c_data = c->data.ptr + i*c->step;
 765         uchar* a_data = a ? a->data.ptr + i*a->step : 0;
 766         uchar* b_data = b->data.ptr + i*b->step;
 767
 768         switch( CV_MAT_DEPTH(b->type) )
 769         {
 770         case CV_8U:
 771             for( j = 0; j < ncols; j++ )
 772             {
 773                 int denom = ((uchar*)b_data)[j];
 774                 int num = a_data ? ((uchar*)a_data)[j] : 1;
 775                 buf[j] = !denom ? 0 : (alpha_buf[j]*num/denom);
 776             }
 777             break;
 778         case CV_8S:
 779             for( j = 0; j < ncols; j++ )
 780             {
 781                 int denom = ((schar*)b_data)[j];
 782                 int num = a_data ? ((schar*)a_data)[j] : 1;
 783                 buf[j] = !denom ? 0 : (alpha_buf[j]*num/denom);
 784             }
 785             break;
 786         case CV_16U:
 787             for( j = 0; j < ncols; j++ )
 788             {
 789                 int denom = ((ushort*)b_data)[j];
 790                 int num = a_data ? ((ushort*)a_data)[j] : 1;
 791                 buf[j] = !denom ? 0 : (alpha_buf[j]*num/denom);
 792             }
 793             break;
 794         case CV_16S:
 795             for( j = 0; j < ncols; j++ )
 796             {
 797                 int denom = ((short*)b_data)[j];
 798                 int num = a_data ? ((short*)a_data)[j] : 1;
 799                 buf[j] = !denom ? 0 : (alpha_buf[j]*num/denom);
 800             }
 801             break;
 802         case CV_32S:
 803             for( j = 0; j < ncols; j++ )
 804             {
 805                 int denom = ((int*)b_data)[j];
 806                 int num = a_data ? ((int*)a_data)[j] : 1;
 807                 buf[j] = !denom ? 0 : (alpha_buf[j]*num/denom);
 808             }
 809             break;
 810         case CV_32F:
 811             for( j = 0; j < ncols; j++ )
 812             {
 813                 double denom = ((float*)b_data)[j];
 814                 double num = a_data ? ((float*)a_data)[j] : 1;
 815                 buf[j] = !denom ? 0 : (alpha_buf[j]*num/denom);
 816             }
 817             break;
 818         case CV_64F:
 819             for( j = 0; j < ncols; j++ )
 820             {
 821                 double denom = ((double*)b_data)[j];
 822                 double num = a_data ? ((double*)a_data)[j] : 1;
 823                 buf[j] = !denom ? 0 : (alpha_buf[j]*num/denom);
 824             }
 825             break;
 826         default:
 827             assert(0);
 828             return;
 829         }
 830
 831         switch( CV_MAT_DEPTH(c->type) )
 832         {
 833         case CV_8U:
 834             for( j = 0; j < ncols; j++ )
 835             {
 836                 int val = cvRound(buf[j]);
 837                 ((uchar*)c_data)[j] = CV_CAST_8U(val);
 838             }
 839             break;
 840         case CV_8S:
 841             for( j = 0; j < ncols; j++ )
 842             {
 843                 int val = cvRound(buf[j]);
 844                 ((schar*)c_data)[j] = CV_CAST_8S(val);
 845             }
 846             break;
 847         case CV_16U:
 848             for( j = 0; j < ncols; j++ )
 849             {
 850                 int val = cvRound(buf[j]);
 851                 ((ushort*)c_data)[j] = CV_CAST_16U(val);
 852             }
 853             break;
 854         case CV_16S:
 855             for( j = 0; j < ncols; j++ )
 856             {
 857                 int val = cvRound(buf[j]);
 858                 ((short*)c_data)[j] = CV_CAST_16S(val);
 859             }
 860             break;
 861         case CV_32S:
 862             for( j = 0; j < ncols; j++ )
 863             {
 864                 int val = cvRound(buf[j]);
 865                 ((int*)c_data)[j] = CV_CAST_32S(val);
 866             }
 867             break;
 868         case CV_32F:
 869             for( j = 0; j < ncols; j++ )
 870                 ((float*)c_data)[j] = CV_CAST_32F(buf[j]);
 871             break;
 872         case CV_64F:
 873             for( j = 0; j < ncols; j++ )
 874                 ((double*)c_data)[j] = CV_CAST_64F(buf[j]);
 875             break;
 876         default:
 877             assert(0);
 878         }
 879     }
 880
 881     if( alpha_buf )
 882         free( alpha_buf );
 883
 884     if( buf )
 885         free( buf );
 886 }
 887
 888
 889 // c = min(a,b) or c = max(a,b)
 890 void cvTsMinMax( const CvMat* a, const CvMat* b, CvMat* c, int op_type )
 891 {
 892     int i, j, ncols;
 893     int calc_max = op_type == CV_TS_MAX;
 894
 895     if( !a || !b || !c )
 896     {
 897         assert(0);
 898         return;
 899     }
 900
 901     assert( CV_ARE_SIZES_EQ(a,c) && CV_ARE_TYPES_EQ(a,c) &&
 902             CV_ARE_SIZES_EQ(b,c) && CV_ARE_TYPES_EQ(b,c) &&
 903             CV_MAT_CN(a->type) == 1 );
 904     ncols = c->cols;
 905
 906     for( i = 0; i < c->rows; i++ )
 907     {
 908         uchar* c_data = c->data.ptr + i*c->step;
 909         uchar* a_data = a->data.ptr + i*a->step;
 910         uchar* b_data = b->data.ptr + i*b->step;
 911
 912         switch( CV_MAT_DEPTH(a->type) )
 913         {
 914         case CV_8U:
 915             for( j = 0; j < ncols; j++ )
 916             {
 917                 int aj = ((uchar*)a_data)[j];
 918                 int bj = ((uchar*)b_data)[j];
 919                 ((uchar*)c_data)[j] = (uchar)(calc_max ? MAX(aj, bj) : MIN(aj,bj));
 920             }
 921             break;
 922         case CV_8S:
 923             for( j = 0; j < ncols; j++ )
 924             {
 925                 int aj = ((schar*)a_data)[j];
 926                 int bj = ((schar*)b_data)[j];
 927                 ((schar*)c_data)[j] = (schar)(calc_max ? MAX(aj, bj) : MIN(aj,bj));
 928             }
 929             break;
 930         case CV_16U:
 931             for( j = 0; j < ncols; j++ )
 932             {
 933                 int aj = ((ushort*)a_data)[j];
 934                 int bj = ((ushort*)b_data)[j];
 935                 ((ushort*)c_data)[j] = (ushort)(calc_max ? MAX(aj, bj) : MIN(aj,bj));
 936             }
 937             break;
 938         case CV_16S:
 939             for( j = 0; j < ncols; j++ )
 940             {
 941                 int aj = ((short*)a_data)[j];
 942                 int bj = ((short*)b_data)[j];
 943                 ((short*)c_data)[j] = (short)(calc_max ? MAX(aj, bj) : MIN(aj,bj));
 944             }
 945             break;
 946         case CV_32S:
 947             for( j = 0; j < ncols; j++ )
 948             {
 949                 int aj = ((int*)a_data)[j];
 950                 int bj = ((int*)b_data)[j];
 951                 ((int*)c_data)[j] = calc_max ? MAX(aj, bj) : MIN(aj,bj);
 952             }
 953             break;
 954         case CV_32F:
 955             for( j = 0; j < ncols; j++ )
 956             {
 957                 float aj = ((float*)a_data)[j];
 958                 float bj = ((float*)b_data)[j];
 959                 ((float*)c_data)[j] = calc_max ? MAX(aj, bj) : MIN(aj,bj);
 960             }
 961             break;
 962         case CV_64F:
 963             for( j = 0; j < ncols; j++ )
 964             {
 965                 double aj = ((double*)a_data)[j];
 966                 double bj = ((double*)b_data)[j];
 967                 ((double*)c_data)[j] = calc_max ? MAX(aj, bj) : MIN(aj,bj);
 968             }
 969             break;
 970         default:
 971             assert(0);
 972             return;
 973         }
 974     }
 975 }
 976
 977 // c = min(a,b) or c = max(a,b)
 978 void cvTsMinMaxS( const CvMat* a, double s, CvMat* c, int op_type )
 979 {
 980     int i, j, ncols;
 981     int calc_max = op_type == CV_TS_MAX;
 982     float fs = (float)s;
 983     int is = cvRound(s);
 984
 985     if( !a || !c )
 986     {
 987         assert(0);
 988         return;
 989     }
 990
 991     assert( CV_ARE_SIZES_EQ(a,c) && CV_ARE_TYPES_EQ(a,c) &&
 992             CV_MAT_CN(a->type) == 1 );
 993     ncols = c->cols;
 994
 995     switch( CV_MAT_DEPTH(a->type) )
 996     {
 997     case CV_8U:
 998         is = CV_CAST_8U(is);
 999         break;
1000     case CV_8S:
1001         is = CV_CAST_8S(is);
1002         break;
1003     case CV_16U:
1004         is = CV_CAST_16U(is);
1005         break;
1006     case CV_16S:
1007         is = CV_CAST_16S(is);
1008         break;
1009     default:
1010         ;
1011     }
1012
1013     for( i = 0; i < c->rows; i++ )
1014     {
1015         uchar* c_data = c->data.ptr + i*c->step;
1016         uchar* a_data = a->data.ptr + i*a->step;
1017
1018         switch( CV_MAT_DEPTH(a->type) )
1019         {
1020         case CV_8U:
1021             for( j = 0; j < ncols; j++ )
1022             {
1023                 int aj = ((uchar*)a_data)[j];
1024                 ((uchar*)c_data)[j] = (uchar)(calc_max ? MAX(aj, is) : MIN(aj, is));
1025             }
1026             break;
1027         case CV_8S:
1028             for( j = 0; j < ncols; j++ )
1029             {
1030                 int aj = ((schar*)a_data)[j];
1031                 ((schar*)c_data)[j] = (schar)(calc_max ? MAX(aj, is) : MIN(aj, is));
1032             }
1033             break;
1034         case CV_16U:
1035             for( j = 0; j < ncols; j++ )
1036             {
1037                 int aj = ((ushort*)a_data)[j];
1038                 ((ushort*)c_data)[j] = (ushort)(calc_max ? MAX(aj, is) : MIN(aj, is));
1039             }
1040             break;
1041         case CV_16S:
1042             for( j = 0; j < ncols; j++ )
1043             {
1044                 int aj = ((short*)a_data)[j];
1045                 ((short*)c_data)[j] = (short)(calc_max ? MAX(aj, is) : MIN(aj, is));
1046             }
1047             break;
1048         case CV_32S:
1049             for( j = 0; j < ncols; j++ )
1050             {
1051                 int aj = ((int*)a_data)[j];
1052                 ((int*)c_data)[j] = calc_max ? MAX(aj, is) : MIN(aj, is);
1053             }
1054             break;
1055         case CV_32F:
1056             for( j = 0; j < ncols; j++ )
1057             {
1058                 float aj = ((float*)a_data)[j];
1059                 ((float*)c_data)[j] = calc_max ? MAX(aj, fs) : MIN(aj, fs);
1060             }
1061             break;
1062         case CV_64F:
1063             for( j = 0; j < ncols; j++ )
1064             {
1065                 double aj = ((double*)a_data)[j];
1066                 ((double*)c_data)[j] = calc_max ? MAX(aj, s) : MIN(aj, s);
1067             }
1068             break;
1069         default:
1070             assert(0);
1071             return;
1072         }
1073     }
1074 }
1075
1076 // checks that the array does not have NaNs and/or Infs and all the elements are
1077 // within [min_val,max_val). idx is the index of the first "bad" element.
1078 int cvTsCheck( const CvMat* a, double min_val, double max_val, CvPoint* idx )
1079 {
1080     int i = 0, j = 0;
1081     int cn, ncols;
1082     int imin = 0, imax = 0;
1083     cn = CV_MAT_CN(a->type);
1084     ncols = a->cols*cn;
1085
1086     if( CV_MAT_DEPTH(a->type) <= CV_32S )
1087     {
1088         imin = cvCeil(min_val);
1089         imax = cvFloor(max_val);
1090     }
1091
1092     for( i = 0; i < a->rows; i++ )
1093     {
1094         uchar* data = a->data.ptr + a->step*i;
1095
1096         switch( CV_MAT_DEPTH(a->type) )
1097         {
1098         case CV_8U:
1099             for( j = 0; j < ncols; j++ )
1100             {
1101                 int val = ((uchar*)data)[j];
1102                 if( val < imin || imax < val )
1103                     goto _exit_;
1104             }
1105             break;
1106         case CV_8S:
1107             for( j = 0; j < ncols; j++ )
1108             {
1109                 int val = ((schar*)data)[j];
1110                 if( val < imin || imax < val )
1111                     goto _exit_;
1112             }
1113             break;
1114         case CV_16U:
1115             for( j = 0; j < ncols; j++ )
1116             {
1117                 int val = ((ushort*)data)[j];
1118                 if( val < imin || imax < val )
1119                     goto _exit_;
1120             }
1121             break;
1122         case CV_16S:
1123             for( j = 0; j < ncols; j++ )
1124             {
1125                 int val = ((short*)data)[j];
1126                 if( val < imin || imax < val )
1127                     goto _exit_;
1128             }
1129             break;
1130         case CV_32S:
1131             for( j = 0; j < ncols; j++ )
1132             {
1133                 int val = ((int*)data)[j];
1134                 if( val < imin || imax < val )
1135                     goto _exit_;
1136             }
1137             break;
1138         case CV_32F:
1139             for( j = 0; j < ncols; j++ )
1140             {
1141                 double val = ((float*)data)[j];
1142                 if( cvIsNaN(val) || cvIsInf(val) || val < min_val || max_val < val )
1143                     goto _exit_;
1144             }
1145             break;
1146         case CV_64F:
1147             for( j = 0; j < ncols; j++ )
1148             {
1149                 double val = ((double*)data)[j];
1150                 if( cvIsNaN(val) || cvIsInf(val) || val < min_val || max_val < val )
1151                     goto _exit_;
1152             }
1153             break;
1154         default:
1155             assert(0);
1156             return -1;
1157         }
1158     }
1159
1160     return 0;
1161 _exit_:
1162
1163     idx->x = j;
1164     idx->y = i;
1165     return -1;
1166 }
1167
1168 // compares two arrays. max_diff is the maximum actual difference,
1169 // success_err_level is maximum allowed difference, idx is the index of the first
1170 // element for which difference is >success_err_level
1171 // (or index of element with the maximum difference)
1172 int cvTsCmpEps( const CvMat* check_arr, const CvMat* etalon, double* _max_diff,
1173                 double success_err_level, CvPoint* idx, bool element_wise_relative_error )
1174 {
1175     int i = 0, j = 0;
1176     int cn, ncols;
1177     double maxdiff = 0;
1178     double maxval = 0;
1179     int imaxdiff = 0;
1180     int ilevel = 0;
1181     int result = -1;
1182
1183     cn = CV_MAT_CN(check_arr->type);
1184     ncols = check_arr->cols*cn;
1185
1186     *idx = cvPoint(0,0);
1187
1188     assert( CV_ARE_TYPES_EQ(check_arr,etalon) && CV_ARE_SIZES_EQ(check_arr,etalon) );
1189
1190     if( CV_MAT_DEPTH(check_arr->type) < CV_32S )
1191         ilevel = cvFloor(success_err_level);
1192
1193     if( CV_MAT_DEPTH(check_arr->type) >= CV_32F && !element_wise_relative_error )
1194     {
1195         double maxval0 = 1.;
1196         maxval = cvTsNorm( etalon, 0, CV_C, 0 );
1197         maxval = MAX(maxval, maxval0);
1198     }
1199
1200     for( i = 0; i < check_arr->rows; i++ )
1201     {
1202         uchar* a_data = check_arr->data.ptr + check_arr->step*i;
1203         uchar* b_data = etalon->data.ptr + etalon->step*i;
1204
1205         switch( CV_MAT_DEPTH(check_arr->type) )
1206         {
1207         case CV_8U:
1208             for( j = 0; j < ncols; j++ )
1209             {
1210                 int val = abs(((uchar*)a_data)[j] - ((uchar*)b_data)[j]);
1211                 if( val > imaxdiff )
1212                 {
1213                     imaxdiff = val;
1214                     *idx = cvPoint(j,i);
1215                     if( val > ilevel )
1216                         goto _exit_;
1217                 }
1218             }
1219             break;
1220         case CV_8S:
1221             for( j = 0; j < ncols; j++ )
1222             {
1223                 int val = abs(((schar*)a_data)[j] - ((schar*)b_data)[j]);
1224                 if( val > imaxdiff )
1225                 {
1226                     imaxdiff = val;
1227                     *idx = cvPoint(j,i);
1228                     if( val > ilevel )
1229                         goto _exit_;
1230                 }
1231             }
1232             break;
1233         case CV_16U:
1234             for( j = 0; j < ncols; j++ )
1235             {
1236                 int val = abs(((ushort*)a_data)[j] - ((ushort*)b_data)[j]);
1237                 if( val > imaxdiff )
1238                 {
1239                     imaxdiff = val;
1240                     *idx = cvPoint(j,i);
1241                     if( val > ilevel )
1242                         goto _exit_;
1243                 }
1244             }
1245             break;
1246         case CV_16S:
1247             for( j = 0; j < ncols; j++ )
1248             {
1249                 int val = abs(((short*)a_data)[j] - ((short*)b_data)[j]);
1250                 if( val > imaxdiff )
1251                 {
1252                     imaxdiff = val;
1253                     *idx = cvPoint(j,i);
1254                     if( val > ilevel )
1255                         goto _exit_;
1256                 }
1257             }
1258             break;
1259         case CV_32S:
1260             for( j = 0; j < ncols; j++ )
1261             {
1262                 double val = fabs((double)((int*)a_data)[j] - (double)((int*)b_data)[j]);
1263                 if( val > maxdiff )
1264                 {
1265                     maxdiff = val;
1266                     *idx = cvPoint(j,i);
1267                     if( val > success_err_level )
1268                         goto _exit_;
1269                 }
1270             }
1271             break;
1272         case CV_32F:
1273             for( j = 0; j < ncols; j++ )
1274             {
1275                 double a_val = ((float*)a_data)[j];
1276                 double b_val = ((float*)b_data)[j];
1277                 double threshold;
1278                 if( cvIsNaN(a_val) || cvIsInf(a_val) )
1279                 {
1280                     result = -2;
1281                     *idx = cvPoint(j,i);
1282                     goto _exit_;
1283                 }
1284                 if( cvIsNaN(b_val) || cvIsInf(b_val) )
1285                 {
1286                     result = -3;
1287                     *idx = cvPoint(j,i);
1288                     goto _exit_;
1289                 }
1290                 a_val = fabs(a_val - b_val);
1291                 threshold = element_wise_relative_error ? fabs(b_val) + 1 : maxval;
1292                 if( a_val > threshold*success_err_level )
1293                 {
1294                     maxdiff = a_val/threshold;
1295                     *idx = cvPoint(j,i);
1296                     goto _exit_;
1297                 }
1298             }
1299             break;
1300         case CV_64F:
1301             for( j = 0; j < ncols; j++ )
1302             {
1303                 double a_val = ((double*)a_data)[j];
1304                 double b_val = ((double*)b_data)[j];
1305                 double threshold;
1306                 if( cvIsNaN(a_val) || cvIsInf(a_val) )
1307                 {
1308                     result = -2;
1309                     *idx = cvPoint(j,i);
1310                     goto _exit_;
1311                 }
1312                 if( cvIsNaN(b_val) || cvIsInf(b_val) )
1313                 {
1314                     result = -3;
1315                     *idx = cvPoint(j,i);
1316                     goto _exit_;
1317                 }
1318                 a_val = fabs(a_val - b_val);
1319                 threshold = element_wise_relative_error ? fabs(b_val)+1 : maxval;
1320                 if( a_val > threshold*success_err_level )
1321                 {
1322                     maxdiff = a_val/threshold;
1323                     *idx = cvPoint(j,i);
1324                     goto _exit_;
1325                 }
1326             }
1327             break;
1328         default:
1329             assert(0);
1330             return -1;
1331         }
1332     }
1333
1334     result = 0;
1335 _exit_:
1336
1337     if( CV_MAT_DEPTH(check_arr->type) < CV_32S )
1338         maxdiff = imaxdiff;
1339
1340     if( result < -1 )
1341         maxdiff = exp(1000.);
1342     *_max_diff = maxdiff;
1343     return result;
1344 }
1345
1346
1347 int cvTsCmpEps2( CvTS* ts, const CvArr* _a, const CvArr* _b, double success_err_level,
1348                  bool element_wise_relative_error, const char* desc )
1349 {
1350     char msg[100];
1351     double diff = 0;
1352     CvMat astub, bstub, *a, *b;
1353     CvPoint idx = {0,0};
1354     int code;
1355
1356     a = cvGetMat( _a, &astub );
1357     b = cvGetMat( _b, &bstub );
1358     code = cvTsCmpEps( a, b, &diff, success_err_level, &idx,
1359         element_wise_relative_error );
1360
1361     switch( code )
1362     {
1363     case -1:
1364         sprintf( msg, "%s: Too big difference (=%g)", desc, diff );
1365         code = CvTS::FAIL_BAD_ACCURACY;
1366         break;
1367     case -2:
1368         sprintf( msg, "%s: Invalid output", desc );
1369         code = CvTS::FAIL_INVALID_OUTPUT;
1370         break;
1371     case -3:
1372         sprintf( msg, "%s: Invalid reference output", desc );
1373         code = CvTS::FAIL_INVALID_OUTPUT;
1374         break;
1375     default:
1376         ;
1377     }
1378
1379     if( code < 0 )
1380     {
1381         if( a->rows == 1 && a->cols == 1 )
1382         {
1383             assert( idx.x == 0 && idx.y == 0 );
1384             ts->printf( CvTS::LOG, "%s\n", msg );
1385         }
1386         else if( a->rows == 1 || a->cols == 1 )
1387         {
1388             assert( idx.x == 0 || idx.y == 0 );
1389             ts->printf( CvTS::LOG, "%s at element %d\n", msg, idx.x + idx.y );
1390         }
1391         else
1392             ts->printf( CvTS::LOG, "%s at (%d,%d)\n", msg, idx.x, idx.y );
1393     }
1394
1395     return code;
1396 }
1397
1398
1399 int cvTsCmpEps2_64f( CvTS* ts, const double* val, const double* ref_val, int len,
1400                      double eps, const char* param_name )
1401 {
1402     CvMat _val = cvMat( 1, len, CV_64F, (void*)val );
1403     CvMat _ref_val = cvMat( 1, len, CV_64F, (void*)ref_val );
1404
1405     return cvTsCmpEps2( ts, &_val, &_ref_val, eps, true, param_name );
1406 }
1407
1408 // compares two arrays. the result is 8s image that takes values -1, 0, 1
1409 void cvTsCmp( const CvMat* a, const CvMat* b, CvMat* result, int cmp_op )
1410 {
1411     int i = 0, j = 0, ncols;
1412     ncols = a->cols;
1413
1414     assert( CV_ARE_TYPES_EQ(a,b) && CV_ARE_SIZES_EQ(a,b) && CV_MAT_CN(a->type) == 1 );
1415     assert( CV_ARE_SIZES_EQ(a,result) &&
1416             (CV_MAT_TYPE(result->type) == CV_8UC1 ||
1417              CV_MAT_TYPE(result->type) == CV_8SC1 ));
1418
1419     for( i = 0; i < a->rows; i++ )
1420     {
1421         uchar* a_data = a->data.ptr + a->step*i;
1422         uchar* b_data = b->data.ptr + b->step*i;
1423         schar* r_data = (schar*)(result->data.ptr + result->step*i);
1424
1425         switch( CV_MAT_DEPTH(a->type) )
1426         {
1427         case CV_8U:
1428             for( j = 0; j < ncols; j++ )
1429             {
1430                 int a_val = ((uchar*)a_data)[j];
1431                 int b_val = ((uchar*)b_data)[j];
1432                 r_data[j] = (schar)CV_CMP(a_val,b_val);
1433             }
1434             break;
1435         case CV_8S:
1436             for( j = 0; j < ncols; j++ )
1437             {
1438                 int a_val = ((schar*)a_data)[j];
1439                 int b_val = ((schar*)b_data)[j];
1440                 r_data[j] = (schar)CV_CMP(a_val,b_val);
1441             }
1442             break;
1443         case CV_16U:
1444             for( j = 0; j < ncols; j++ )
1445             {
1446                 int a_val = ((ushort*)a_data)[j];
1447                 int b_val = ((ushort*)b_data)[j];
1448                 r_data[j] = (schar)CV_CMP(a_val,b_val);
1449             }
1450             break;
1451         case CV_16S:
1452             for( j = 0; j < ncols; j++ )
1453             {
1454                 int a_val = ((short*)a_data)[j];
1455                 int b_val = ((short*)b_data)[j];
1456                 r_data[j] = (schar)CV_CMP(a_val,b_val);
1457             }
1458             break;
1459         case CV_32S:
1460             for( j = 0; j < ncols; j++ )
1461             {
1462                 int a_val = ((int*)a_data)[j];
1463                 int b_val = ((int*)b_data)[j];
1464                 r_data[j] = (schar)CV_CMP(a_val,b_val);
1465             }
1466             break;
1467         case CV_32F:
1468             for( j = 0; j < ncols; j++ )
1469             {
1470                 float a_val = ((float*)a_data)[j];
1471                 float b_val = ((float*)b_data)[j];
1472                 r_data[j] = (schar)CV_CMP(a_val,b_val);
1473             }
1474             break;
1475         case CV_64F:
1476             for( j = 0; j < ncols; j++ )
1477             {
1478                 double a_val = ((double*)a_data)[j];
1479                 double b_val = ((double*)b_data)[j];
1480                 r_data[j] = (schar)CV_CMP(a_val,b_val);
1481             }
1482             break;
1483         default:
1484             assert(0);
1485         }
1486
1487         switch( cmp_op )
1488         {
1489         case CV_CMP_EQ:
1490             for( j = 0; j < ncols; j++ )
1491                 r_data[j] = (schar)(r_data[j] == 0 ? -1 : 0);
1492             break;
1493         case CV_CMP_NE:
1494             for( j = 0; j < ncols; j++ )
1495                 r_data[j] = (schar)(r_data[j] != 0 ? -1 : 0);
1496             break;
1497         case CV_CMP_LT:
1498             for( j = 0; j < ncols; j++ )
1499                 r_data[j] = (schar)(r_data[j] < 0 ? -1 : 0);
1500             break;
1501         case CV_CMP_LE:
1502             for( j = 0; j < ncols; j++ )
1503                 r_data[j] = (schar)(r_data[j] <= 0 ? -1 : 0);
1504             break;
1505         case CV_CMP_GE:
1506             for( j = 0; j < ncols; j++ )
1507                 r_data[j] = (schar)(r_data[j] >= 0 ? -1 : 0);
1508             break;
1509         case CV_CMP_GT:
1510             for( j = 0; j < ncols; j++ )
1511                 r_data[j] = (schar)(r_data[j] > 0 ? -1 : 0);
1512             break;
1513         default:
1514             ;
1515         }
1516     }
1517 }
1518
1519 // compares two arrays. the result is 8s image that takes values -1, 0, 1
1520 void cvTsCmpS( const CvMat* a, double fval, CvMat* result, int cmp_op )
1521 {
1522     int i = 0, j = 0;
1523     int ncols, ival = 0;
1524     ncols = a->cols;
1525
1526     if( CV_MAT_DEPTH(a->type) <= CV_32S )
1527         ival = cvRound(fval);
1528
1529     assert( CV_MAT_CN(a->type) == 1 );
1530     assert( CV_ARE_SIZES_EQ(a,result) &&
1531             (CV_MAT_TYPE(result->type) == CV_8UC1 ||
1532              CV_MAT_TYPE(result->type) == CV_8SC1 ));
1533
1534     for( i = 0; i < a->rows; i++ )
1535     {
1536         uchar* a_data = a->data.ptr + a->step*i;
1537         schar* r_data = (schar*)(result->data.ptr + result->step*i);
1538
1539         switch( CV_MAT_DEPTH(a->type) )
1540         {
1541         case CV_8U:
1542             for( j = 0; j < ncols; j++ )
1543             {
1544                 int a_val = ((uchar*)a_data)[j];
1545                 r_data[j] = (schar)CV_CMP(a_val,ival);
1546             }
1547             break;
1548         case CV_8S:
1549             for( j = 0; j < ncols; j++ )
1550             {
1551                 int a_val = ((schar*)a_data)[j];
1552                 r_data[j] = (schar)CV_CMP(a_val,ival);
1553             }
1554             break;
1555         case CV_16U:
1556             for( j = 0; j < ncols; j++ )
1557             {
1558                 int a_val = ((ushort*)a_data)[j];
1559                 r_data[j] = (schar)CV_CMP(a_val,ival);
1560             }
1561             break;
1562         case CV_16S:
1563             for( j = 0; j < ncols; j++ )
1564             {
1565                 int a_val = ((short*)a_data)[j];
1566                 r_data[j] = (schar)CV_CMP(a_val,ival);
1567             }
1568             break;
1569         case CV_32S:
1570             for( j = 0; j < ncols; j++ )
1571             {
1572                 int a_val = ((int*)a_data)[j];
1573                 r_data[j] = (schar)CV_CMP(a_val,ival);
1574             }
1575             break;
1576         case CV_32F:
1577             for( j = 0; j < ncols; j++ )
1578             {
1579                 float a_val = ((float*)a_data)[j];
1580                 r_data[j] = (schar)CV_CMP(a_val,fval);
1581             }
1582             break;
1583         case CV_64F:
1584             for( j = 0; j < ncols; j++ )
1585             {
1586                 double a_val = ((double*)a_data)[j];
1587                 r_data[j] = (schar)CV_CMP(a_val,fval);
1588             }
1589             break;
1590         default:
1591             assert(0);
1592         }
1593
1594         switch( cmp_op )
1595         {
1596         case CV_CMP_EQ:
1597             for( j = 0; j < ncols; j++ )
1598                 r_data[j] = (schar)(r_data[j] == 0 ? -1 : 0);
1599             break;
1600         case CV_CMP_NE:
1601             for( j = 0; j < ncols; j++ )
1602                 r_data[j] = (schar)(r_data[j] != 0 ? -1 : 0);
1603             break;
1604         case CV_CMP_LT:
1605             for( j = 0; j < ncols; j++ )
1606                 r_data[j] = (schar)(r_data[j] < 0 ? -1 : 0);
1607             break;
1608         case CV_CMP_LE:
1609             for( j = 0; j < ncols; j++ )
1610                 r_data[j] = (schar)(r_data[j] <= 0 ? -1 : 0);
1611             break;
1612         case CV_CMP_GE:
1613             for( j = 0; j < ncols; j++ )
1614                 r_data[j] = (schar)(r_data[j] >= 0 ? -1 : 0);
1615             break;
1616         case CV_CMP_GT:
1617             for( j = 0; j < ncols; j++ )
1618                 r_data[j] = (schar)(r_data[j] > 0 ? -1 : 0);
1619             break;
1620         default:
1621             ;
1622         }
1623     }
1624 }
1625
1626
1627 // calculates norm of a matrix
1628 double cvTsNorm( const CvMat* arr, const CvMat* mask, int norm_type, int coi )
1629 {
1630     int i = 0, j = 0, k;
1631     int depth, cn0, cn, ncols, el_size1;
1632     int inorm = 0;
1633     double fnorm = 0;
1634     void* buffer = 0;
1635     uchar* zerobuf = 0;
1636
1637     cn0 = cn = CV_MAT_CN(arr->type);
1638     ncols = arr->cols*cn;
1639     depth = CV_MAT_DEPTH(arr->type);
1640     el_size1 = CV_ELEM_SIZE(depth);
1641     zerobuf = (uchar*)cvStackAlloc(el_size1*cn);
1642     memset( zerobuf, 0, el_size1*cn);
1643
1644     if( mask )
1645     {
1646         assert( CV_ARE_SIZES_EQ( arr, mask ) && CV_IS_MASK_ARR(mask) );
1647         buffer = cvStackAlloc( el_size1*ncols );
1648     }
1649
1650     if( coi == 0 )
1651         cn = 1;
1652
1653     for( i = 0; i < arr->rows; i++ )
1654     {
1655         const uchar* data = arr->data.ptr + arr->step*i + (coi - (coi != 0))*el_size1;
1656
1657         if( mask )
1658         {
1659             const uchar* mdata = mask->data.ptr + mask->step*i;
1660             switch( depth )
1661             {
1662             case CV_8U:
1663             case CV_8S:
1664                 for( j = 0; j < ncols; j += cn0 )
1665                 {
1666                     const uchar* src = *mdata++ ? (uchar*)data + j : zerobuf;
1667                     for( k = 0; k < cn0; k++ )
1668                         ((uchar*)buffer)[j+k] = src[k];
1669                 }
1670                 break;
1671             case CV_16U:
1672             case CV_16S:
1673                 for( j = 0; j < ncols; j += cn0 )
1674                 {
1675                     const short* src = *mdata++ ? (short*)data + j : (short*)zerobuf;
1676                     for( k = 0; k < cn0; k++ )
1677                         ((short*)buffer)[j+k] = src[k];
1678                 }
1679                 break;
1680             case CV_32S:
1681             case CV_32F:
1682                 for( j = 0; j < ncols; j += cn0 )
1683                 {
1684                     const int* src = *mdata++ ? (int*)data + j : (int*)zerobuf;
1685                     for( k = 0; k < cn0; k++ )
1686                         ((int*)buffer)[j+k] = src[k];
1687                 }
1688                 break;
1689             case CV_64F:
1690                 for( j = 0; j < ncols; j += cn0 )
1691                 {
1692                     const double* src = *mdata++ ? (double*)data + j : (double*)zerobuf;
1693                     for( k = 0; k < cn0; k++ )
1694                         ((double*)buffer)[j+k] = src[k];
1695                 }
1696                 break;
1697             default:
1698                 assert(0);
1699                 return -1;
1700             }
1701             data = (const uchar*)buffer;
1702         }
1703
1704         switch( depth )
1705         {
1706         case CV_8U:
1707             if( norm_type == CV_C )
1708             {
1709                 for( j = 0; j < ncols; j += cn )
1710                 {
1711                     int val = ((const uchar*)data)[j];
1712                     inorm = MAX( inorm, val );
1713                 }
1714             }
1715             else if( norm_type == CV_L1 )
1716             {
1717                 inorm = 0;
1718                 for( j = 0; j < ncols; j += cn )
1719                 {
1720                     int val = ((const uchar*)data)[j];
1721                     inorm += val;
1722                 }
1723                 fnorm += inorm;
1724             }
1725             else
1726             {
1727                 inorm = 0;
1728                 for( j = 0; j < ncols; j += cn )
1729                 {
1730                     int val = ((const uchar*)data)[j];
1731                     inorm += val*val;
1732                 }
1733                 fnorm += inorm;
1734             }
1735             break;
1736         case CV_8S:
1737             if( norm_type == CV_C )
1738             {
1739                 for( j = 0; j < ncols; j += cn )
1740                 {
1741                     int val = abs(((const schar*)data)[j]);
1742                     inorm = MAX( inorm, val );
1743                 }
1744             }
1745             else if( norm_type == CV_L1 )
1746             {
1747                 inorm = 0;
1748                 for( j = 0; j < ncols; j += cn )
1749                 {
1750                     int val = abs(((const schar*)data)[j]);
1751                     inorm += val;
1752                 }
1753                 fnorm += inorm;
1754             }
1755             else
1756             {
1757                 inorm = 0;
1758                 for( j = 0; j < ncols; j += cn )
1759                 {
1760                     int val = ((const schar*)data)[j];
1761                     inorm += val*val;
1762                 }
1763                 fnorm += inorm;
1764             }
1765             break;
1766         case CV_16U:
1767             if( norm_type == CV_C )
1768             {
1769                 for( j = 0; j < ncols; j += cn )
1770                 {
1771                     int val = ((const ushort*)data)[j];
1772                     inorm = MAX( inorm, val );
1773                 }
1774             }
1775             else if( norm_type == CV_L1 )
1776             {
1777                 inorm = 0;
1778                 for( j = 0; j < ncols; j += cn )
1779                 {
1780                     int val = ((const ushort*)data)[j];
1781                     inorm += val;
1782                 }
1783                 fnorm += inorm;
1784             }
1785             else
1786             {
1787                 for( j = 0; j < ncols; j += cn )
1788                 {
1789                     double val = ((const ushort*)data)[j];
1790                     fnorm += val*val;
1791                 }
1792             }
1793             break;
1794         case CV_16S:
1795             if( norm_type == CV_C )
1796             {
1797                 for( j = 0; j < ncols; j += cn )
1798                 {
1799                     int val = abs(((const short*)data)[j]);
1800                     inorm = MAX( inorm, val );
1801                 }
1802             }
1803             else if( norm_type == CV_L1 )
1804             {
1805                 inorm = 0;
1806                 for( j = 0; j < ncols; j += cn )
1807                 {
1808                     int val = abs(((const short*)data)[j]);
1809                     inorm += val;
1810                 }
1811                 fnorm += inorm;
1812             }
1813             else
1814             {
1815                 for( j = 0; j < ncols; j += cn )
1816                 {
1817                     double val = ((const short*)data)[j];
1818                     fnorm += val*val;
1819                 }
1820             }
1821             break;
1822         case CV_32S:
1823             if( norm_type == CV_C )
1824             {
1825                 for( j = 0; j < ncols; j += cn )
1826                 {
1827                     int val = abs(((const int*)data)[j]);
1828                     inorm = MAX( inorm, val );
1829                 }
1830             }
1831             else if( norm_type == CV_L1 )
1832             {
1833                 for( j = 0; j < ncols; j += cn )
1834                 {
1835                     double val = fabs((double)((const int*)data)[j]);
1836                     fnorm += val;
1837                 }
1838             }
1839             else
1840             {
1841                 for( j = 0; j < ncols; j += cn )
1842                 {
1843                     double val = ((const int*)data)[j];
1844                     fnorm += val*val;
1845                 }
1846             }
1847             break;
1848         case CV_32F:
1849             if( norm_type == CV_C )
1850             {
1851                 for( j = 0; j < ncols; j += cn )
1852                 {
1853                     double val = fabs((double)((const float*)data)[j]);
1854                     fnorm = MAX( fnorm, val );
1855                 }
1856             }
1857             else if( norm_type == CV_L1 )
1858             {
1859                 for( j = 0; j < ncols; j += cn )
1860                 {
1861                     double val = fabs((double)((const float*)data)[j]);
1862                     fnorm += val;
1863                 }
1864             }
1865             else
1866             {
1867                 for( j = 0; j < ncols; j += cn )
1868                 {
1869                     double val = ((const float*)data)[j];
1870                     fnorm += val*val;
1871                 }
1872             }
1873             break;
1874         case CV_64F:
1875             if( norm_type == CV_C )
1876             {
1877                 for( j = 0; j < ncols; j += cn )
1878                 {
1879                     double val = fabs(((const double*)data)[j]);
1880                     fnorm = MAX( fnorm, val );
1881                 }
1882             }
1883             else if( norm_type == CV_L1 )
1884             {
1885                 for( j = 0; j < ncols; j += cn )
1886                 {
1887                     double val = fabs(((const double*)data)[j]);
1888                     fnorm += val;
1889                 }
1890             }
1891             else
1892             {
1893                 for( j = 0; j < ncols; j += cn )
1894                 {
1895                     double val = ((const double*)data)[j];
1896                     fnorm += val*val;
1897                 }
1898             }
1899             break;
1900         default:
1901             assert(0);
1902             return -1;
1903         }
1904     }
1905
1906     if( norm_type == CV_L2 )
1907         fnorm = sqrt( fnorm );
1908     else if( depth < CV_32F && norm_type == CV_C )
1909         fnorm = inorm;
1910
1911     return fnorm;
1912 }
1913
1914
1915 // retrieves mean, standard deviation and the number of nonzero mask pixels
1916 int cvTsMeanStdDevNonZero( const CvMat* arr, const CvMat* mask,
1917                            CvScalar* _mean, CvScalar* _stddev, int coi )
1918 {
1919     int i = 0, j = 0, k;
1920     int depth, cn0, cn, cols, ncols, el_size1;
1921     CvScalar sum = cvScalar(0), sqsum = cvScalar(0);
1922     double inv_area;
1923     int isum[4], isqsum[4];
1924     int nonzero = 0;
1925     uchar* maskbuf = 0;
1926
1927     cn0 = cn = CV_MAT_CN(arr->type);
1928     cols = arr->cols;
1929     ncols = arr->cols*cn;
1930     depth = CV_MAT_DEPTH(arr->type);
1931     el_size1 = CV_ELEM_SIZE(depth);
1932     if( mask )
1933     {
1934         assert( CV_ARE_SIZES_EQ( arr, mask ) && CV_IS_MASK_ARR(mask) );
1935     }
1936     else
1937     {
1938         maskbuf = (uchar*)cvStackAlloc( cols );
1939         memset( maskbuf, 1, cols );
1940         nonzero = cols*arr->rows;
1941     }
1942
1943     if( coi != 0 )
1944         cn = 1;
1945
1946     for( i = 0; i < arr->rows; i++ )
1947     {
1948         const uchar* data = arr->data.ptr + arr->step*i + (coi - (coi != 0))*el_size1;
1949         const uchar* mdata;
1950
1951         if( mask )
1952         {
1953             mdata = mask->data.ptr + mask->step*i;
1954             for( j = 0; j < cols; j++ )
1955                 nonzero += mdata[j] != 0;
1956         }
1957         else
1958         {
1959             mdata = maskbuf;
1960         }
1961
1962         // if only a number of pixels in the mask is needed, skip the rest of the loop body
1963         if( !_mean && !_stddev )
1964             continue;
1965
1966         switch( depth )
1967         {
1968         case CV_8U:
1969             for( k = 0; k < cn; k++ )
1970                 isum[k] = isqsum[k] = 0;
1971             for( j = 0; j < ncols; j += cn0 )
1972                 if( *mdata++ )
1973                 {
1974                     for( k = 0; k < cn; k++ )
1975                     {
1976                         int val = ((const uchar*)data)[j+k];
1977                         isum[k] += val;
1978                         isqsum[k] += val*val;
1979                     }
1980                 }
1981             for( k = 0; k < cn; k++ )
1982             {
1983                 sum.val[k] += isum[k];
1984                 sqsum.val[k] += isqsum[k];
1985             }
1986             break;
1987         case CV_8S:
1988             for( k = 0; k < cn; k++ )
1989                 isum[k] = isqsum[k] = 0;
1990             for( j = 0; j < ncols; j += cn0 )
1991                 if( *mdata++ )
1992                 {
1993                     for( k = 0; k < cn; k++ )
1994                     {
1995                         int val = ((const schar*)data)[j+k];
1996                         isum[k] += val;
1997                         isqsum[k] += val*val;
1998                     }
1999                 }
2000             for( k = 0; k < cn; k++ )
2001             {
2002                 sum.val[k] += isum[k];
2003                 sqsum.val[k] += isqsum[k];
2004             }
2005             break;
2006         case CV_16U:
2007             for( k = 0; k < cn; k++ )
2008                 isum[k] = 0;
2009             for( j = 0; j < ncols; j += cn0 )
2010                 if( *mdata++ )
2011                 {
2012                     for( k = 0; k < cn; k++ )
2013                     {
2014                         int val = ((const ushort*)data)[j+k];
2015                         isum[k] += val;
2016                         sqsum.val[k] += ((double)val)*val;
2017                     }
2018                 }
2019             for( k = 0; k < cn; k++ )
2020                 sum.val[k] += isum[k];
2021             break;
2022         case CV_16S:
2023             for( k = 0; k < cn; k++ )
2024                 isum[k] = 0;
2025             for( j = 0; j < ncols; j += cn0 )
2026                 if( *mdata++ )
2027                 {
2028                     for( k = 0; k < cn; k++ )
2029                     {
2030                         int val = ((const short*)data)[j+k];
2031                         isum[k] += val;
2032                         sqsum.val[k] += ((double)val)*val;
2033                     }
2034                 }
2035             for( k = 0; k < cn; k++ )
2036                 sum.val[k] += isum[k];
2037             break;
2038         case CV_32S:
2039             for( j = 0; j < ncols; j += cn0 )
2040                 if( *mdata++ )
2041                 {
2042                     for( k = 0; k < cn; k++ )
2043                     {
2044                         double val = ((const int*)data)[j+k];
2045                         sum.val[k] += val;
2046                         sqsum.val[k] += val*val;
2047                     }
2048                 }
2049             break;
2050         case CV_32F:
2051             for( j = 0; j < ncols; j += cn0 )
2052                 if( *mdata++ )
2053                 {
2054                     for( k = 0; k < cn; k++ )
2055                     {
2056                         double val = ((const float*)data)[j+k];
2057                         sum.val[k] += val;
2058                         sqsum.val[k] += val*val;
2059                     }
2060                 }
2061             break;
2062         case CV_64F:
2063             for( j = 0; j < ncols; j += cn0 )
2064                 if( *mdata++ )
2065                 {
2066                     for( k = 0; k < cn; k++ )
2067                     {
2068                         double val = ((const double*)data)[j+k];
2069                         sum.val[k] += val;
2070                         sqsum.val[k] += val*val;
2071                     }
2072                 }
2073             break;
2074         default:
2075             assert(0);
2076             return -1;
2077         }
2078     }
2079
2080     inv_area = nonzero ? 1./nonzero : 0.;
2081     for( k = 0; k < cn; k++ )
2082     {
2083         sum.val[k] *= inv_area;
2084         double t = sqsum.val[k]*inv_area - sum.val[k]*sum.val[k];
2085         sqsum.val[k] = sqrt(MAX(t, 0));
2086     }
2087     if( _mean )
2088         *_mean = sum;
2089     if( _stddev )
2090         *_stddev = sqsum;
2091     return nonzero;
2092 }
2093
2094 // retrieves global extremums and their positions
2095 void cvTsMinMaxLoc( const CvMat* arr, const CvMat* mask,
2096                     double* _minval, double* _maxval,
2097                     CvPoint* _minidx, CvPoint* _maxidx, int coi )
2098 {
2099     int i = 0, j = 0;
2100     int depth, cn, cols, ncols, el_size1;
2101     CvPoint minidx = {-1,-1}, maxidx = {-1,-1};
2102     uchar* maskbuf = 0;
2103     int iminval = INT_MAX, imaxval = INT_MIN;
2104     double minval = DBL_MAX, maxval = -minval;
2105
2106     cn = CV_MAT_CN(arr->type);
2107     cols = arr->cols;
2108     ncols = arr->cols*cn;
2109     depth = CV_MAT_DEPTH(arr->type);
2110     el_size1 = CV_ELEM_SIZE(depth);
2111
2112     if( mask )
2113     {
2114         assert( CV_ARE_SIZES_EQ( arr, mask ) && CV_IS_MASK_ARR(mask) );
2115     }
2116     else
2117     {
2118         maskbuf = (uchar*)cvStackAlloc( cols );
2119         memset( maskbuf, 1, cols );
2120     }
2121
2122     if( coi == 0 && cn > 1 )
2123     {
2124         assert(0);
2125         return;
2126     }
2127
2128     for( i = 0; i < arr->rows; i++ )
2129     {
2130         const uchar* data = arr->data.ptr + arr->step*i + (coi - (coi != 0))*el_size1;
2131         const uchar* mdata = mask ? mask->data.ptr + mask->step*i : maskbuf;
2132
2133         switch( depth )
2134         {
2135         case CV_8U:
2136             for( j = 0; j < ncols; j += cn, mdata++ )
2137             {
2138                 int val = ((const uchar*)data)[j];
2139                 if( val < iminval && *mdata )
2140                 {
2141                     iminval = val;
2142                     minidx = cvPoint(j,i);
2143                 }
2144                 if( val > imaxval && *mdata )
2145                 {
2146                     imaxval = val;
2147                     maxidx = cvPoint(j,i);
2148                 }
2149             }
2150             break;
2151         case CV_8S:
2152             for( j = 0; j < ncols; j += cn, mdata++ )
2153             {
2154                 int val = ((const schar*)data)[j];
2155                 if( val < iminval && *mdata )
2156                 {
2157                     iminval = val;
2158                     minidx = cvPoint(j,i);
2159                 }
2160                 if( val > imaxval && *mdata )
2161                 {
2162                     imaxval = val;
2163                     maxidx = cvPoint(j,i);
2164                 }
2165             }
2166             break;
2167         case CV_16U:
2168             for( j = 0; j < ncols; j += cn, mdata++ )
2169             {
2170                 int val = ((const ushort*)data)[j];
2171                 if( val < iminval && *mdata )
2172                 {
2173                     iminval = val;
2174                     minidx = cvPoint(j,i);
2175                 }
2176                 if( val > imaxval && *mdata )
2177                 {
2178                     imaxval = val;
2179                     maxidx = cvPoint(j,i);
2180                 }
2181             }
2182             break;
2183         case CV_16S:
2184             for( j = 0; j < ncols; j += cn, mdata++ )
2185             {
2186                 int val = ((const short*)data)[j];
2187                 if( val < iminval && *mdata )
2188                 {
2189                     iminval = val;
2190                     minidx = cvPoint(j,i);
2191                 }
2192                 if( val > imaxval && *mdata )
2193                 {
2194                     imaxval = val;
2195                     maxidx = cvPoint(j,i);
2196                 }
2197             }
2198             break;
2199         case CV_32S:
2200             for( j = 0; j < ncols; j += cn, mdata++ )
2201             {
2202                 int val = ((const int*)data)[j];
2203                 if( val < iminval && *mdata )
2204                 {
2205                     iminval = val;
2206                     minidx = cvPoint(j,i);
2207                 }
2208                 if( val > imaxval && *mdata )
2209                 {
2210                     imaxval = val;
2211                     maxidx = cvPoint(j,i);
2212                 }
2213             }
2214             break;
2215         case CV_32F:
2216             for( j = 0; j < ncols; j += cn, mdata++ )
2217             {
2218                 float val = ((const float*)data)[j];
2219                 if( val < minval && *mdata )
2220                 {
2221                     minval = val;
2222                     minidx = cvPoint(j,i);
2223                 }
2224                 if( val > maxval && *mdata )
2225                 {
2226                     maxval = val;
2227                     maxidx = cvPoint(j,i);
2228                 }
2229             }
2230             break;
2231         case CV_64F:
2232             for( j = 0; j < ncols; j += cn, mdata++ )
2233             {
2234                 double val = ((const double*)data)[j];
2235                 if( val < minval && *mdata )
2236                 {
2237                     minval = val;
2238                     minidx = cvPoint(j,i);
2239                 }
2240                 if( val > maxval && *mdata )
2241                 {
2242                     maxval = val;
2243                     maxidx = cvPoint(j,i);
2244                 }
2245             }
2246             break;
2247         default:
2248             assert(0);
2249             return;
2250         }
2251     }
2252
2253     if( minidx.x < 0 )
2254         minval = maxval = 0;
2255     else
2256     {
2257         if( depth < CV_32F )
2258             minval = iminval, maxval = imaxval;
2259         minidx.x /= cn;
2260         maxidx.x /= cn;
2261     }
2262
2263     if( _minval )
2264         *_minval = minval;
2265
2266     if( _maxval )
2267         *_maxval = maxval;
2268
2269     if( _minidx )
2270         *_minidx = minidx;
2271
2272     if( _maxidx )
2273         *_maxidx = maxidx;
2274 }
2275
2276
2277 void cvTsLogic( const CvMat* a, const CvMat* b, CvMat* c, int logic_op )
2278 {
2279     int i = 0, j = 0, ncols;
2280     ncols = a->cols*CV_ELEM_SIZE(a->type);
2281
2282     assert( CV_ARE_TYPES_EQ(a,b) && CV_ARE_SIZES_EQ(a,b) );
2283     assert( CV_ARE_TYPES_EQ(a,c) && CV_ARE_SIZES_EQ(a,c) );
2284
2285     for( i = 0; i < a->rows; i++ )
2286     {
2287         uchar* a_data = a->data.ptr + a->step*i;
2288         uchar* b_data = b->data.ptr + b->step*i;
2289         uchar* c_data = c->data.ptr + c->step*i;
2290
2291         switch( logic_op )
2292         {
2293         case CV_TS_LOGIC_AND:
2294             for( j = 0; j < ncols; j++ )
2295                 c_data[j] = (uchar)(a_data[j] & b_data[j]);
2296             break;
2297         case CV_TS_LOGIC_OR:
2298             for( j = 0; j < ncols; j++ )
2299                 c_data[j] = (uchar)(a_data[j] | b_data[j]);
2300             break;
2301         case CV_TS_LOGIC_XOR:
2302             for( j = 0; j < ncols; j++ )
2303                 c_data[j] = (uchar)(a_data[j] ^ b_data[j]);
2304             break;
2305         default:
2306             assert(0);
2307             return;
2308         }
2309     }
2310 }
2311
2312 void cvTsLogicS( const CvMat* a, CvScalar s, CvMat* c, int logic_op )
2313 {
2314     int i = 0, j = 0, k;
2315     int cn, ncols, elem_size;
2316     uchar* b_data;
2317     union
2318     {
2319         uchar ptr[4];
2320         schar c[4];
2321         short s[4];
2322         ushort w[4];
2323         int i[4];
2324         float f[4];
2325         double d[4];
2326     } buf;
2327     cn = CV_MAT_CN(a->type);
2328     elem_size = CV_ELEM_SIZE(a->type);
2329     ncols = a->cols * elem_size;
2330     b_data = (uchar*)malloc( ncols );
2331
2332     assert( CV_ARE_TYPES_EQ(a,c) && CV_ARE_SIZES_EQ(a,c) );
2333
2334     if( logic_op == CV_TS_LOGIC_NOT )
2335     {
2336         memset( b_data, -1, ncols );
2337         logic_op = CV_TS_LOGIC_XOR;
2338     }
2339     else
2340     {
2341         switch( CV_MAT_DEPTH(a->type) )
2342         {
2343         case CV_8U:
2344             for( k = 0; k < cn; k++ )
2345             {
2346                 int val = cvRound(s.val[k]);
2347                 buf.ptr[k] = CV_CAST_8U(val);
2348             }
2349             break;
2350         case CV_8S:
2351             for( k = 0; k < cn; k++ )
2352             {
2353                 int val = cvRound(s.val[k]);
2354                 buf.c[k] = CV_CAST_8S(val);
2355             }
2356             break;
2357         case CV_16U:
2358             for( k = 0; k < cn; k++ )
2359             {
2360                 int val = cvRound(s.val[k]);
2361                 buf.w[k] = CV_CAST_16U(val);
2362             }
2363             break;
2364         case CV_16S:
2365             for( k = 0; k < cn; k++ )
2366             {
2367                 int val = cvRound(s.val[k]);
2368                 buf.s[k] = CV_CAST_16S(val);
2369             }
2370             break;
2371         case CV_32S:
2372             for( k = 0; k < cn; k++ )
2373             {
2374                 int val = cvRound(s.val[k]);
2375                 buf.i[k] = CV_CAST_32S(val);
2376             }
2377             break;
2378         case CV_32F:
2379             for( k = 0; k < cn; k++ )
2380             {
2381                 double val = s.val[k];
2382                 buf.f[k] = CV_CAST_32F(val);
2383             }
2384             break;
2385         case CV_64F:
2386             for( k = 0; k < cn; k++ )
2387             {
2388                 double val = s.val[k];
2389                 buf.d[k] = CV_CAST_64F(val);
2390             }
2391             break;
2392         default:
2393             assert(0);
2394             return;
2395         }
2396
2397         for( j = 0; j < ncols; j += elem_size )
2398             memcpy( b_data + j, buf.ptr, elem_size );
2399     }
2400
2401     for( i = 0; i < a->rows; i++ )
2402     {
2403         uchar* a_data = a->data.ptr + a->step*i;
2404         uchar* c_data = c->data.ptr + c->step*i;
2405
2406         switch( logic_op )
2407         {
2408         case CV_TS_LOGIC_AND:
2409             for( j = 0; j < ncols; j++ )
2410                 c_data[j] = (uchar)(a_data[j] & b_data[j]);
2411             break;
2412         case CV_TS_LOGIC_OR:
2413             for( j = 0; j < ncols; j++ )
2414                 c_data[j] = (uchar)(a_data[j] | b_data[j]);
2415             break;
2416         case CV_TS_LOGIC_XOR:
2417             for( j = 0; j < ncols; j++ )
2418                 c_data[j] = (uchar)(a_data[j] ^ b_data[j]);
2419             break;
2420         default:
2421             assert(0);
2422             return;
2423         }
2424     }
2425
2426     if( b_data )
2427         free( b_data );
2428 }
2429
2430
2431 void cvTsGEMM( const CvMat* a, const CvMat* b, double alpha,
2432                const CvMat* c, double beta, CvMat* d, int flags )
2433 {
2434     int i, j, k;
2435     int a_rows, a_cols, b_rows, b_cols;
2436     int c_rows, c_cols, d_rows, d_cols;
2437     int cn, el_size;
2438     int a_step, a_delta, b_step, b_delta;
2439     int c_step, c_delta, d_step;
2440
2441     a_rows = a->rows; a_cols = a->cols;
2442     cn = CV_MAT_CN(a->type);
2443     el_size = CV_ELEM_SIZE(a->type & ~CV_MAT_CN_MASK);
2444     a_step = a->step / el_size; a_delta = cn;
2445     d_rows = d->rows; d_cols = d->cols;
2446     b_rows = b->rows; b_cols = b->cols;
2447     b_step = b->step / el_size; b_delta = cn;
2448     c_rows = c ? c->rows : 0; c_cols = c ? c->cols : 0;
2449     c_step = c ? c->step / el_size : 0; c_delta = c ? cn : 0;
2450     d_step = d->step / el_size;
2451
2452     assert( CV_ARE_TYPES_EQ(a,b) && CV_ARE_TYPES_EQ(a,d) );
2453     assert( CV_MAT_CN(a->type) <= 2 );
2454
2455     if( flags & CV_TS_GEMM_A_T )
2456     {
2457         CV_SWAP( a_rows, a_cols, i );
2458         CV_SWAP( a_step, a_delta, i );
2459     }
2460
2461     if( flags & CV_TS_GEMM_B_T )
2462     {
2463         CV_SWAP( b_rows, b_cols, i );
2464         CV_SWAP( b_step, b_delta, i );
2465     }
2466
2467     if( flags & CV_TS_GEMM_C_T )
2468     {
2469         CV_SWAP( c_rows, c_cols, i );
2470         CV_SWAP( c_step, c_delta, i );
2471     }
2472
2473     assert( a_rows == d_rows && a_cols == b_rows && b_cols == d_cols );
2474     assert( a->data.ptr != d->data.ptr && b->data.ptr != d->data.ptr );
2475
2476     if( c )
2477     {
2478         assert( CV_ARE_TYPES_EQ(a,c) && c_rows == d_rows && c_cols == d_cols );
2479         assert( c->data.ptr != d->data.ptr || (flags & CV_TS_GEMM_C_T) == 0 );
2480     }
2481
2482     if( CV_MAT_DEPTH(a->type) == CV_32F )
2483     {
2484         float* a_data0 = a->data.fl;
2485         float* b_data0 = b->data.fl;
2486         float* c_data0 = c ? c->data.fl : 0;
2487         float* d_data = d->data.fl;
2488
2489         for( i = 0; i < d_rows; i++, d_data += d_step, c_data0 += c_step, a_data0 += a_step )
2490         {
2491             for( j = 0; j < d_cols; j++ )
2492             {
2493                 float* a_data = a_data0;
2494                 float* b_data = b_data0 + j*b_delta;
2495                 float* c_data = c_data0 + j*c_delta;
2496
2497                 if( cn == 1 )
2498                 {
2499                     double s = 0;
2500                     for( k = 0; k < a_cols; k++ )
2501                     {
2502                         s += ((double)a_data[0])*b_data[0];
2503                         a_data += a_delta;
2504                         b_data += b_step;
2505                     }
2506                     d_data[j] = (float)(s*alpha + (c_data ? c_data[0]*beta : 0));
2507                 }
2508                 else
2509                 {
2510                     double s_re = 0, s_im = 0;
2511
2512                     for( k = 0; k < a_cols; k++ )
2513                     {
2514                         s_re += ((double)a_data[0])*b_data[0] - ((double)a_data[1])*b_data[1];
2515                         s_im += ((double)a_data[0])*b_data[1] + ((double)a_data[1])*b_data[0];
2516                         a_data += a_delta;
2517                         b_data += b_step;
2518                     }
2519
2520                     s_re *= alpha;
2521                     s_im *= alpha;
2522
2523                     if( c_data )
2524                     {
2525                         s_re += c_data[0]*beta;
2526                         s_im += c_data[1]*beta;
2527                     }
2528
2529                     d_data[j*2] = (float)s_re;
2530                     d_data[j*2+1] = (float)s_im;
2531                 }
2532             }
2533         }
2534     }
2535     else if( CV_MAT_DEPTH(a->type) == CV_64F )
2536     {
2537         double* a_data0 = a->data.db;
2538         double* b_data0 = b->data.db;
2539         double* c_data0 = c ? c->data.db : 0;
2540         double* d_data = d->data.db;
2541
2542         for( i = 0; i < d_rows; i++, d_data += d_step, c_data0 += c_step, a_data0 += a_step )
2543         {
2544             for( j = 0; j < d_cols; j++ )
2545             {
2546                 double* a_data = a_data0;
2547                 double* b_data = b_data0 + j*b_delta;
2548                 double* c_data = c_data0 + j*c_delta;
2549
2550                 if( cn == 1 )
2551                 {
2552                     double s = 0;
2553                     for( k = 0; k < a_cols; k++ )
2554                     {
2555                         s += a_data[0]*b_data[0];
2556                         a_data += a_delta;
2557                         b_data += b_step;
2558                     }
2559                     d_data[j] = s*alpha + (c_data ? c_data[0]*beta : 0);
2560                 }
2561                 else
2562                 {
2563                     double s_re = 0, s_im = 0;
2564
2565                     for( k = 0; k < a_cols; k++ )
2566                     {
2567                         s_re += a_data[0]*b_data[0] - a_data[1]*b_data[1];
2568                         s_im += a_data[0]*b_data[1] + a_data[1]*b_data[0];
2569                         a_data += a_delta;
2570                         b_data += b_step;
2571                     }
2572                     s_re *= alpha;
2573                     s_im *= alpha;
2574
2575                     if( c_data )
2576                     {
2577                         s_re += c_data[0]*beta;
2578                         s_im += c_data[1]*beta;
2579                     }
2580
2581                     d_data[j*2] = s_re;
2582                     d_data[j*2+1] = s_im;
2583                 }
2584             }
2585         }
2586     }
2587     else
2588     {
2589         assert(0);
2590     }
2591 }
2592
2593
2594 CvMat* cvTsSelect( const CvMat* a, CvMat* header, CvRect rect )
2595 {
2596     CvMat h;
2597     int el_size;
2598
2599     h = cvMat( rect.height, rect.width, a->type );
2600     el_size = CV_ELEM_SIZE(a->type);
2601
2602     h.data.ptr = a->data.ptr + rect.y*a->step + rect.x*el_size;
2603     h.step = rect.height > 1 ? a->step : 0;
2604     h.type &= ~CV_MAT_CONT_FLAG;
2605     if( rect.height == 1 || h.step == h.cols*el_size )
2606         h.type |= CV_MAT_CONT_FLAG;
2607     *header = h;
2608     return header;
2609 }
2610
2611
2612 double cvTsMinVal( int type )
2613 {
2614     switch( CV_MAT_DEPTH(type) )
2615     {
2616     case CV_8U:
2617         return 0;
2618     case CV_8S:
2619         return -128;
2620     case CV_16U:
2621         return 0;
2622     case CV_16S:
2623         return -32768;
2624     case CV_32S:
2625         return -1000000;
2626     case CV_32F:
2627         return -1000.;
2628     case CV_64F:
2629         return -1000.;
2630     }
2631     return log(-1.);
2632 }
2633
2634
2635 double cvTsMaxVal( int type )
2636 {
2637     switch( CV_MAT_DEPTH(type) )
2638     {
2639     case CV_8U:
2640         return 256;
2641     case CV_8S:
2642         return 128;
2643     case CV_16U:
2644         return 65536;
2645     case CV_16S:
2646         return 32768;
2647     case CV_32S:
2648         return 1000000;
2649     case CV_32F:
2650         return 1000.;
2651     case CV_64F:
2652         return 1000.;
2653     }
2654     return log(-1.);
2655 }
2656
2657
2658 void cvTsPrepareToFilter( const CvMat* a, CvMat* b, CvPoint ofs,
2659                           int border_mode, CvScalar fill_val )
2660 {
2661     int i, j, dir;
2662     CvMat temp, temp2;
2663
2664     assert( 0 <= ofs.x && ofs.x <= b->cols - a->cols &&
2665             0 <= ofs.y && ofs.y <= b->rows - a->rows );
2666
2667     cvTsSelect( b, &temp, cvRect( ofs.x, ofs.y, a->cols, a->rows ));
2668     cvTsCopy( a, &temp, 0 );
2669
2670     assert( border_mode == CV_TS_BORDER_FILL ||
2671             border_mode == CV_TS_BORDER_REPLICATE ||
2672             border_mode == CV_TS_BORDER_REFLECT );
2673
2674     if( ofs.y > 0 )
2675     {
2676         if( border_mode == CV_TS_BORDER_FILL )
2677         {
2678             cvTsSelect( b, &temp, cvRect( ofs.x, 0, a->cols, ofs.y ));
2679             cvTsAdd( 0, cvScalar(0), 0, cvScalar(0), fill_val, &temp, 0 );
2680         }
2681         else if( border_mode == CV_TS_BORDER_REPLICATE || a->rows == 1 )
2682         {
2683             cvTsSelect( b, &temp, cvRect( ofs.x, ofs.y, a->cols, 1 ));
2684             for( i = ofs.y-1; i >= 0; i-- )
2685             {
2686                 cvTsSelect( b, &temp2, cvRect( ofs.x, i, a->cols, 1 ));
2687                 cvTsCopy( &temp, &temp2, 0 );
2688             }
2689         }
2690         else if( border_mode == CV_TS_BORDER_REFLECT )
2691         {
2692             j = 1; dir = 1;
2693             for( i = ofs.y-1; i >= 0; i-- )
2694             {
2695                 cvTsSelect( b, &temp, cvRect( ofs.x, ofs.y+j, a->cols, 1 ));
2696                 cvTsSelect( b, &temp2, cvRect( ofs.x, i, a->cols, 1 ));
2697                 cvTsCopy( &temp, &temp2, 0 );
2698                 if( (unsigned)(j + dir) >= (unsigned)a->rows )
2699                     dir = -dir;
2700                 j += dir;
2701             }
2702         }
2703     }
2704
2705     ofs.y += a->rows;
2706     if( ofs.y < b->rows )
2707     {
2708         if( border_mode == CV_TS_BORDER_FILL )
2709         {
2710             cvTsSelect( b, &temp, cvRect( ofs.x, ofs.y, a->cols, b->rows - ofs.y ));
2711             cvTsAdd( 0, cvScalar(0), 0, cvScalar(0), fill_val, &temp, 0 );
2712         }
2713         else if( border_mode == CV_TS_BORDER_REPLICATE || a->rows == 1 )
2714         {
2715             cvTsSelect( b, &temp, cvRect( ofs.x, ofs.y - 1, a->cols, 1 ));
2716             for( i = ofs.y; i < b->rows; i++ )
2717             {
2718                 cvTsSelect( b, &temp2, cvRect( ofs.x, i, a->cols, 1 ));
2719                 cvTsCopy( &temp, &temp2, 0 );
2720             }
2721         }
2722         else
2723         {
2724             j = a->rows - 2; dir = -1;
2725             for( i = ofs.y; i < b->rows; i++ )
2726             {
2727                 cvTsSelect( b, &temp, cvRect( ofs.x, ofs.y-a->rows+j, a->cols, 1 ));
2728                 cvTsSelect( b, &temp2, cvRect( ofs.x, i, a->cols, 1 ));
2729                 cvTsCopy( &temp, &temp2, 0 );
2730                 if( (unsigned)(j + dir) >= (unsigned)a->rows )
2731                     dir = -dir;
2732                 j += dir;
2733             }
2734         }
2735     }
2736
2737     if( ofs.x > 0 )
2738     {
2739         if( border_mode == CV_TS_BORDER_FILL )
2740         {
2741             cvTsSelect( b, &temp, cvRect( 0, 0, ofs.x, b->rows ));
2742             cvTsAdd( 0, cvScalar(0), 0, cvScalar(0), fill_val, &temp, 0 );
2743         }
2744         else if( border_mode == CV_TS_BORDER_REPLICATE || a->cols == 1 )
2745         {
2746             cvTsSelect( b, &temp, cvRect( ofs.x, 0, 1, b->rows ));
2747             for( i = ofs.x-1; i >= 0; i-- )
2748             {
2749                 cvTsSelect( b, &temp2, cvRect( i, 0, 1, b->rows ));
2750                 cvTsCopy( &temp, &temp2, 0 );
2751             }
2752         }
2753         else if( border_mode == CV_TS_BORDER_REFLECT )
2754         {
2755             j = 1; dir = 1;
2756             for( i = ofs.x-1; i >= 0; i-- )
2757             {
2758                 cvTsSelect( b, &temp, cvRect( ofs.x+j, 0, 1, b->rows ));
2759                 cvTsSelect( b, &temp2, cvRect( i, 0, 1, b->rows ));
2760                 cvTsCopy( &temp, &temp2, 0 );
2761                 if( (unsigned)(j + dir) >= (unsigned)a->cols )
2762                     dir = -dir;
2763                 j += dir;
2764             }
2765         }
2766     }
2767
2768     ofs.x += a->cols;
2769     if( ofs.x < b->cols )
2770     {
2771         if( border_mode == CV_TS_BORDER_FILL )
2772         {
2773             cvTsSelect( b, &temp, cvRect( ofs.x, 0, b->cols - ofs.x, b->rows ));
2774             cvTsAdd( 0, cvScalar(0), 0, cvScalar(0), fill_val, &temp, 0 );
2775         }
2776         else if( border_mode == CV_TS_BORDER_REPLICATE || a->cols == 1 )
2777         {
2778             cvTsSelect( b, &temp, cvRect( ofs.x-1, 0, 1, b->rows ));
2779             for( i = ofs.x; i < b->cols; i++ )
2780             {
2781                 cvTsSelect( b, &temp2, cvRect( i, 0, 1, b->rows ));
2782                 cvTsCopy( &temp, &temp2, 0 );
2783             }
2784         }
2785         else if( border_mode == CV_TS_BORDER_REFLECT )
2786         {
2787             j = a->cols - 2; dir = -1;
2788             for( i = ofs.x; i < b->cols; i++ )
2789             {
2790                 cvTsSelect( b, &temp, cvRect( ofs.x-a->cols+j, 0, 1, b->rows ));
2791                 cvTsSelect( b, &temp2, cvRect( i, 0, 1, b->rows ));
2792                 cvTsCopy( &temp, &temp2, 0 );
2793                 if( (unsigned)(j + dir) >= (unsigned)a->cols )
2794                     dir = -dir;
2795                 j += dir;
2796             }
2797         }
2798     }
2799 }
2800
2801
2802 void cvTsConvolve2D( const CvMat* a, CvMat* b, const CvMat* kernel, CvPoint anchor )
2803 {
2804     int i, j, k;
2805     int cn, ncols, a_step;
2806     int ker_size = kernel->rows*kernel->cols;
2807     int* offset = (int*)malloc( ker_size*sizeof(offset[0]));
2808     float* k_data = (float*)malloc( ker_size*sizeof(k_data[0]));
2809     int all_same = 1;
2810     float first = kernel->data.fl[0];
2811     uchar *a_data, *b_data;
2812
2813     cn = CV_MAT_CN(a->type);
2814     ncols = b->cols*cn;
2815     a_step = a->step / CV_ELEM_SIZE(a->type & ~CV_MAT_CN_MASK);
2816
2817     assert( a->cols == b->cols + kernel->cols - 1 &&
2818             a->rows == b->rows + kernel->rows - 1 && CV_ARE_TYPES_EQ( a, b ) );
2819     assert( CV_MAT_TYPE(kernel->type) == CV_32FC1 );
2820     assert( 0 <= anchor.x && anchor.x < kernel->cols &&
2821             0 <= anchor.y && anchor.y < kernel->rows );
2822
2823     for( i = 0, k = 0; i < kernel->rows; i++ )
2824         for( j = 0; j < kernel->cols; j++ )
2825         {
2826             float f = ((float*)(kernel->data.ptr + kernel->step*i))[j];
2827             if( f )
2828             {
2829                 k_data[k] = f;
2830                 offset[k++] = (i - anchor.y)*a_step + (j - anchor.x)*cn;
2831             }
2832             if( f != first )
2833                 all_same = 0;
2834         }
2835
2836     ker_size = k;
2837     a_data = a->data.ptr + a->step*anchor.y + CV_ELEM_SIZE(a->type)*anchor.x;
2838     b_data = b->data.ptr;
2839
2840     for( i = 0; i < b->rows; i++, a_data += a->step, b_data += b->step )
2841     {
2842         switch( CV_MAT_DEPTH(a->type) )
2843         {
2844         case CV_8U:
2845             for( j = 0; j < ncols; j++ )
2846             {
2847                 double s = 0;
2848                 int val;
2849                 for( k = 0; k < ker_size; k++ )
2850                     s += ((uchar*)a_data)[j+offset[k]]*k_data[k];
2851                 val = cvRound(s);
2852                 ((uchar*)b_data)[j] = CV_CAST_8U(val);
2853             }
2854             break;
2855         case CV_8S:
2856             for( j = 0; j < ncols; j++ )
2857             {
2858                 double s = 0;
2859                 int val;
2860                 for( k = 0; k < ker_size; k++ )
2861                     s += ((schar*)a_data)[j+offset[k]]*k_data[k];
2862                 val = cvRound(s);
2863                 ((schar*)b_data)[j] = CV_CAST_8S(val);
2864             }
2865             break;
2866         case CV_16U:
2867             for( j = 0; j < ncols; j++ )
2868             {
2869                 double s = 0;
2870                 int val;
2871                 for( k = 0; k < ker_size; k++ )
2872                     s += ((ushort*)a_data)[j+offset[k]]*k_data[k];
2873                 val = cvRound(s);
2874                 ((ushort*)b_data)[j] = CV_CAST_16U(val);
2875             }
2876             break;
2877         case CV_16S:
2878             for( j = 0; j < ncols; j++ )
2879             {
2880                 double s = 0;
2881                 int val;
2882                 for( k = 0; k < ker_size; k++ )
2883                     s += ((short*)a_data)[j+offset[k]]*k_data[k];
2884                 val = cvRound(s);
2885                 ((short*)b_data)[j] = CV_CAST_16S(val);
2886             }
2887             break;
2888         case CV_32S:
2889             for( j = 0; j < ncols; j++ )
2890             {
2891                 double s = 0;
2892                 for( k = 0; k < ker_size; k++ )
2893                     s += ((int*)a_data)[j+offset[k]]*k_data[k];
2894                 ((int*)b_data)[j] = cvRound(s);
2895             }
2896             break;
2897         case CV_32F:
2898             if( !all_same )
2899             {
2900                 for( j = 0; j < ncols; j++ )
2901                 {
2902                     double s = 0;
2903                     for( k = 0; k < ker_size; k++ )
2904                         s += (double)((float*)a_data)[j+offset[k]]*k_data[k];
2905                     ((float*)b_data)[j] = (float)s;
2906                 }
2907             }
2908             else
2909             {
2910                 // special branch to speedup feature selection and blur tests
2911                 for( j = 0; j < ncols; j++ )
2912                 {
2913                     double s = 0;
2914                     for( k = 0; k < ker_size; k++ )
2915                         s += (double)((float*)a_data)[j+offset[k]];
2916                     ((float*)b_data)[j] = (float)(s*first);
2917                 }
2918             }
2919             break;
2920         case CV_64F:
2921             for( j = 0; j < ncols; j++ )
2922             {
2923                 double s = 0;
2924                 for( k = 0; k < ker_size; k++ )
2925                     s += ((double*)a_data)[j+offset[k]]*k_data[k];
2926                 ((double*)b_data)[j] = (double)s;
2927             }
2928             break;
2929         default:
2930             assert(0);
2931         }
2932     }
2933
2934     free( offset );
2935     free( k_data );
2936 }
2937
2938
2939 void cvTsMinMaxFilter( const CvMat* a, CvMat* b, const IplConvKernel* kernel, int op_type )
2940 {
2941     int i, j, k;
2942     int cn, ncols, a_step;
2943     int ker_size = kernel->nRows*kernel->nCols;
2944     int* offset = (int*)malloc( ker_size*sizeof(offset[0]));
2945     int calc_max = op_type == CV_TS_MAX;
2946     uchar *a_data, *b_data;
2947
2948     cn = CV_MAT_CN(a->type);
2949     ncols = b->cols*cn;
2950     a_step = a->step / CV_ELEM_SIZE(a->type & ~CV_MAT_CN_MASK);
2951
2952     assert( a->cols == b->cols + kernel->nCols - 1 &&
2953             a->rows == b->rows + kernel->nRows - 1 && CV_ARE_TYPES_EQ( a, b ) );
2954     assert( 0 <= kernel->anchorX && kernel->anchorX < kernel->nCols &&
2955             0 <= kernel->anchorY && kernel->anchorY < kernel->nRows );
2956
2957     for( i = 0, k = 0; i < kernel->nRows; i++ )
2958         for( j = 0; j < kernel->nCols; j++ )
2959         {
2960             if( !kernel->values || kernel->values[i*kernel->nCols + j] )
2961                 offset[k++] = (i - kernel->anchorY)*a_step + (j - kernel->anchorX)*cn;
2962         }
2963
2964     if( k == 0 )
2965         offset[k++] = 0;
2966
2967     ker_size = k;
2968
2969     a_data = a->data.ptr + kernel->anchorY*a->step + kernel->anchorX*CV_ELEM_SIZE(a->type);
2970     b_data = b->data.ptr;
2971
2972     for( i = 0; i < b->rows; i++, a_data += a->step, b_data += b->step )
2973     {
2974         switch( CV_MAT_DEPTH(a->type) )
2975         {
2976         case CV_8U:
2977             for( j = 0; j < ncols; j++ )
2978             {
2979                 int m = ((uchar*)a_data)[j+offset[0]];
2980                 for( k = 1; k < ker_size; k++ )
2981                 {
2982                     int v = ((uchar*)a_data)[j+offset[k]];
2983                     if( calc_max )
2984                     {
2985                         if( m < v )
2986                             m = v;
2987                     }
2988                     else if( m > v )
2989                         m = v;
2990                 }
2991                 ((uchar*)b_data)[j] = (uchar)m;
2992             }
2993             break;
2994         case CV_16U:
2995             for( j = 0; j < ncols; j++ )
2996             {
2997                 int m = ((ushort*)a_data)[j+offset[0]];
2998                 for( k = 1; k < ker_size; k++ )
2999                 {
3000                     int v = ((ushort*)a_data)[j+offset[k]];
3001                     if( calc_max )
3002                     {
3003                         if( m < v )
3004                             m = v;
3005                     }
3006                     else if( m > v )
3007                         m = v;
3008                 }
3009                 ((ushort*)b_data)[j] = (ushort)m;
3010             }
3011             break;
3012         case CV_16S:
3013             for( j = 0; j < ncols; j++ )
3014             {
3015                 int m = ((short*)a_data)[j+offset[0]];
3016                 for( k = 1; k < ker_size; k++ )
3017                 {
3018                     int v = ((short*)a_data)[j+offset[k]];
3019                     if( calc_max )
3020                     {
3021                         if( m < v )
3022                             m = v;
3023                     }
3024                     else if( m > v )
3025                         m = v;
3026                 }
3027                 ((short*)b_data)[j] = (short)m;
3028             }
3029             break;
3030         case CV_32S:
3031             for( j = 0; j < ncols; j++ )
3032             {
3033                 int m = ((int*)a_data)[j+offset[0]];
3034                 for( k = 1; k < ker_size; k++ )
3035                 {
3036                     int v = ((int*)a_data)[j+offset[k]];
3037                     if( calc_max )
3038                     {
3039                         if( m < v )
3040                             m = v;
3041                     }
3042                     else if( m > v )
3043                         m = v;
3044                 }
3045                 ((int*)b_data)[j] = m;
3046             }
3047             break;
3048         case CV_32F:
3049             for( j = 0; j < ncols; j++ )
3050             {
3051                 float m = ((float*)a_data)[j+offset[0]];
3052                 for( k = 1; k < ker_size; k++ )
3053                 {
3054                     float v = ((float*)a_data)[j+offset[k]];
3055                     if( calc_max )
3056                     {
3057                         if( m < v )
3058                             m = v;
3059                     }
3060                     else if( m > v )
3061                         m = v;
3062                 }
3063                 ((float*)b_data)[j] = (float)m;
3064             }
3065             break;
3066         case CV_64F:
3067             for( j = 0; j < ncols; j++ )
3068             {
3069                 double m = ((double*)a_data)[j+offset[0]];
3070                 for( k = 1; k < ker_size; k++ )
3071                 {
3072                     double v = ((double*)a_data)[j+offset[k]];
3073                     if( calc_max )
3074                     {
3075                         if( m < v )
3076                             m = v;
3077                     }
3078                     else if( m > v )
3079                         m = v;
3080                 }
3081                 ((double*)b_data)[j] = (double)m;
3082             }
3083             break;
3084         default:
3085             assert(0);
3086         }
3087     }
3088
3089     free( offset );
3090 }
3091
3092
3093 double cvTsCrossCorr( const CvMat* a, const CvMat* b )
3094 {
3095     int i, j;
3096     int cn, ncols;
3097     double s = 0;
3098
3099     cn = CV_MAT_CN(a->type);
3100     ncols = a->cols*cn;
3101
3102     assert( CV_ARE_SIZES_EQ( a, b ) && CV_ARE_TYPES_EQ( a, b ) );
3103     for( i = 0; i < a->rows; i++ )
3104     {
3105         uchar* a_data = a->data.ptr + a->step*i;
3106         uchar* b_data = b->data.ptr + b->step*i;
3107
3108         switch( CV_MAT_DEPTH(a->type) )
3109         {
3110         case CV_8U:
3111             for( j = 0; j < ncols; j++ )
3112                 s += ((uchar*)a_data)[j]*((uchar*)b_data)[j];
3113             break;
3114         case CV_8S:
3115             for( j = 0; j < ncols; j++ )
3116                 s += ((schar*)a_data)[j]*((schar*)b_data)[j];
3117             break;
3118         case CV_16U:
3119             for( j = 0; j < ncols; j++ )
3120                 s += (double)((ushort*)a_data)[j]*((ushort*)b_data)[j];
3121             break;
3122         case CV_16S:
3123             for( j = 0; j < ncols; j++ )
3124                 s += ((short*)a_data)[j]*((short*)b_data)[j];
3125             break;
3126         case CV_32S:
3127             for( j = 0; j < ncols; j++ )
3128                 s += ((double)((int*)a_data)[j])*((int*)b_data)[j];
3129             break;
3130         case CV_32F:
3131             for( j = 0; j < ncols; j++ )
3132                 s += ((double)((float*)a_data)[j])*((float*)b_data)[j];
3133             break;
3134         case CV_64F:
3135             for( j = 0; j < ncols; j++ )
3136                 s += ((double*)a_data)[j]*((double*)b_data)[j];
3137             break;
3138         default:
3139             assert(0);
3140             return log(-1.);
3141         }
3142     }
3143
3144     return s;
3145 }
3146
3147
3148 void cvTsTransform( const CvMat* a, CvMat* b, const CvMat* transmat, const CvMat* shift )
3149 {
3150     int i, j, k, cols, dst_cols;
3151     int cn, dst_cn, depth, mat_depth, shiftstep;
3152     double mat[20], *buf, *dst_buf;
3153
3154     cn = CV_MAT_CN(a->type);
3155     dst_cn = CV_MAT_CN(b->type);
3156     depth = CV_MAT_DEPTH(a->type);
3157     mat_depth = CV_MAT_DEPTH(transmat->type);
3158     cols = transmat->cols;
3159
3160     // prepare cn x (cn + 1) transform matrix
3161     if( mat_depth == CV_32F )
3162     {
3163         shiftstep = shift && shift->rows > 1 ? shift->step/sizeof(float) : 1;
3164         for( i = 0; i < transmat->rows; i++ )
3165         {
3166             mat[i*(cn+1) + cn] = 0.;
3167             for( j = 0; j < cols; j++ )
3168                 mat[i*(cn+1) + j] = ((float*)(transmat->data.ptr + transmat->step*i))[j];
3169             if( shift )
3170                 mat[i*(cn+1) + cn] = shift->data.fl[i*shiftstep];
3171         }
3172     }
3173     else
3174     {
3175         assert( mat_depth == CV_64F );
3176
3177         shiftstep = shift && shift->rows > 1 ? shift->step/sizeof(double) : 1;
3178         for( i = 0; i < transmat->rows; i++ )
3179         {
3180             mat[i*(cn+1) + cn] = 0.;
3181             for( j = 0; j < cols; j++ )
3182                 mat[i*(cn+1) + j] = ((double*)(transmat->data.ptr + transmat->step*i))[j];
3183             if( shift )
3184                 mat[i*(cn+1) + cn] = shift->data.db[i*shiftstep];
3185         }
3186     }
3187
3188     // transform data
3189     cols = a->cols * cn;
3190     dst_cols = a->cols * dst_cn;
3191     buf = (double*)cvStackAlloc( cols * sizeof(double) );
3192     dst_buf = (double*)cvStackAlloc( dst_cols * sizeof(double) );
3193
3194     for( i = 0; i < a->rows; i++ )
3195     {
3196         uchar* src = a->data.ptr + i*a->step;
3197         uchar* dst = b->data.ptr + i*b->step;
3198         double* _dst = dst_buf;
3199
3200         switch( depth )
3201         {
3202         case CV_8U:
3203             for( j = 0; j < cols; j++ )
3204                 buf[j] = ((uchar*)src)[j];
3205             break;
3206         case CV_16U:
3207             for( j = 0; j < cols; j++ )
3208                 buf[j] = ((ushort*)src)[j];
3209             break;
3210         case CV_16S:
3211             for( j = 0; j < cols; j++ )
3212                 buf[j] = ((short*)src)[j];
3213             break;
3214         case CV_32S:
3215             for( j = 0; j < cols; j++ )
3216                 buf[j] = ((int*)src)[j];
3217             break;
3218         case CV_32F:
3219             for( j = 0; j < cols; j++ )
3220                 buf[j] = ((float*)src)[j];
3221             break;
3222         case CV_64F:
3223             for( j = 0; j < cols; j++ )
3224                 buf[j] = ((double*)src)[j];
3225             break;
3226         default:
3227             assert(0);
3228         }
3229
3230         switch( cn )
3231         {
3232         case 1:
3233             for( j = 0; j < cols; j++, _dst += dst_cn )
3234                 for( k = 0; k < dst_cn; k++ )
3235                     _dst[k] = buf[j]*mat[2*k] + mat[2*k+1];
3236             break;
3237         case 2:
3238             for( j = 0; j < cols; j += 2, _dst += dst_cn )
3239                 for( k = 0; k < dst_cn; k++ )
3240                     _dst[k] = buf[j]*mat[3*k] + buf[j+1]*mat[3*k+1] + mat[3*k+2];
3241             break;
3242         case 3:
3243             for( j = 0; j < cols; j += 3, _dst += dst_cn )
3244                 for( k = 0; k < dst_cn; k++ )
3245                     _dst[k] = buf[j]*mat[4*k] + buf[j+1]*mat[4*k+1] +
3246                               buf[j+2]*mat[4*k+2] + mat[4*k+3];
3247             break;
3248         case 4:
3249             for( j = 0; j < cols; j += 4, _dst += dst_cn )
3250                 for( k = 0; k < dst_cn; k++ )
3251                     _dst[k] = buf[j]*mat[5*k] + buf[j+1]*mat[5*k+1] +
3252                         buf[j+2]*mat[5*k+2] + buf[j+3]*mat[5*k+3] + mat[5*k+4];
3253             break;
3254         default:
3255             assert(0);
3256         }
3257
3258         switch( depth )
3259         {
3260         case CV_8U:
3261             for( j = 0; j < dst_cols; j++ )
3262             {
3263                 int val = cvRound(dst_buf[j]);
3264                 ((uchar*)dst)[j] = CV_CAST_8U(val);
3265             }
3266             break;
3267         case CV_16U:
3268             for( j = 0; j < dst_cols; j++ )
3269             {
3270                 int val = cvRound(dst_buf[j]);
3271                 ((ushort*)dst)[j] = CV_CAST_16U(val);
3272             }
3273             break;
3274         case CV_16S:
3275             for( j = 0; j < dst_cols; j++ )
3276             {
3277                 int val = cvRound(dst_buf[j]);
3278                 ((short*)dst)[j] = CV_CAST_16S(val);
3279             }
3280             break;
3281         case CV_32S:
3282             for( j = 0; j < dst_cols; j++ )
3283                 ((int*)dst)[j] = cvRound(dst_buf[j]);
3284             break;
3285         case CV_32F:
3286             for( j = 0; j < dst_cols; j++ )
3287                 ((float*)dst)[j] = (float)dst_buf[j];
3288             break;
3289         case CV_64F:
3290             for( j = 0; j < dst_cols; j++ )
3291                 ((double*)dst)[j] = dst_buf[j];
3292             break;
3293         default:
3294             assert(0);
3295         }
3296     }
3297 }
3298
3299
3300 CvMat* cvTsTranspose( const CvMat* a, CvMat* b )
3301 {
3302     int i, j, k, rows, cols, elem_size;
3303     uchar *a_data, *b_data;
3304     int a_step, b_step;
3305
3306     elem_size = CV_ELEM_SIZE(a->type);
3307     rows = a->rows;
3308     cols = a->cols;
3309
3310     assert( a->rows == b->cols && a->cols == b->rows && CV_ARE_TYPES_EQ(a,b) );
3311     a_data = a->data.ptr;
3312     a_step = a->step;
3313     b_data = b->data.ptr;
3314     b_step = b->step;
3315
3316     if( rows == cols )
3317     {
3318         for( i = 0; i < rows; i++ )
3319         {
3320             for( j = 0; j <= i; j++ )
3321             {
3322                 uchar* a_ij = a_data + a_step*i + elem_size*j;
3323                 uchar* a_ji = a_data + a_step*j + elem_size*i;
3324                 uchar* b_ij = b_data + b_step*i + elem_size*j;
3325                 uchar* b_ji = b_data + b_step*j + elem_size*i;
3326                 for( k = 0; k < elem_size; k++ )
3327                 {
3328                     uchar t0 = a_ij[k];
3329                     uchar t1 = a_ji[k];
3330                     b_ji[k] = t0;
3331                     b_ij[k] = t1;
3332                 }
3333             }
3334         }
3335     }
3336     else
3337     {
3338         for( i = 0; i < cols; i++ )
3339         {
3340             for( j = 0; j < rows; j++ )
3341             {
3342                 uchar* a_ji = a_data + a_step*j + elem_size*i;
3343                 uchar* b_ij = b_data + b_step*i + elem_size*j;
3344                 for( k = 0; k < elem_size; k++ )
3345                     b_ij[k] = a_ji[k];
3346             }
3347         }
3348     }
3349
3350     return b;
3351 }
3352
3353
3354 void cvTsFlip( const CvMat* a, CvMat* b, int flip_type )
3355 {
3356     int i, j, k, rows, cols, elem_size;
3357     uchar *a_data, *b_data;
3358     int a_step, b_step;
3359
3360     elem_size = CV_ELEM_SIZE(a->type);
3361     rows = a->rows;
3362     cols = a->cols*elem_size;
3363
3364     assert( CV_ARE_SIZES_EQ(a,b) && CV_ARE_TYPES_EQ(a,b) && a->data.ptr != b->data.ptr );
3365     a_data = a->data.ptr;
3366     a_step = a->step;
3367     b_data = b->data.ptr;
3368     b_step = b->step;
3369
3370     if( flip_type <= 0 )
3371     {
3372         a_data += a_step*(rows-1);
3373         a_step = -a_step;
3374     }
3375
3376     for( i = 0; i < rows; i++ )
3377     {
3378         if( flip_type == 0 )
3379             memcpy( b_data, a_data, cols );
3380         else
3381         {
3382             for( j = 0; j < cols; j += elem_size )
3383                 for( k = 0; k < elem_size; k++ )
3384                     b_data[j+k] = a_data[cols - elem_size - j + k];
3385         }
3386         a_data += a_step;
3387         b_data += b_step;
3388     }
3389 }
3390
3391
3392 void  cvTsPatchZeros( CvMat* mat, double level )
3393 {
3394     int i, j, ncols = mat->cols * CV_MAT_CN(mat->type);
3395
3396     for( i = 0; i < mat->rows; i++ )
3397     {
3398         switch( CV_MAT_DEPTH(mat->type) )
3399         {
3400         case CV_32F:
3401             {
3402             float* data = (float*)(mat->data.ptr + i*mat->step);
3403             for( j = 0; j < ncols; j++ )
3404                 if( fabs(data[j]) < level )
3405                     data[j] += 1;
3406             }
3407             break;
3408         case CV_64F:
3409             {
3410             double* data = (double*)(mat->data.ptr + i*mat->step);
3411             for( j = 0; j < ncols; j++ )
3412                 if( fabs(data[j]) < level )
3413                     data[j] += 1;
3414             }
3415             break;
3416         default:
3417             assert(0);
3418             return;
3419         }
3420     }
3421 }
3422
3423
3424 /* End of file. */