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