git.maemo.org Git - opencv/blob - cvaux/src/cvhmmobs.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 #include "_cvaux.h"
  42
  43 //*F///////////////////////////////////////////////////////////////////////////////////////
  44 //    Name: icvImgToObs_DCT_8u32f_C1R
  45 //    Purpose: The function takes as input an image and returns the sequnce of observations
  46 //             to be used with an embedded HMM; Each observation is top-left block of DCT
  47 //             coefficient matrix.
  48 //    Context:
  49 //    Parameters: img     - pointer to the original image ROI
  50 //                imgStep - full row width of the image in bytes
  51 //                roi     - width and height of ROI in pixels
  52 //                obs     - pointer to resultant observation vectors
  53 //                dctSize - size of the block for which DCT is calculated
  54 //                obsSize - size of top-left block of DCT coeffs matrix, which is treated
  55 //                          as observation. Each observation vector consists of
  56 //                          obsSize.width * obsSize.height floats.
  57 //                          The following conditions should be satisfied:
  58 //                          0 < objSize.width <= dctSize.width,
  59 //                          0 < objSize.height <= dctSize.height.
  60 //                delta   - dctBlocks are overlapped and this parameter specifies horizontal
  61 //                          and vertical shift.
  62 //    Returns:
  63 //      CV_NO_ERR or error code
  64 //    Notes:
  65 //      The algorithm is following:
  66 //          1. First, number of observation vectors per row and per column are calculated:
  67 //
  68 //             Nx = floor((roi.width - dctSize.width + delta.width)/delta.width);
  69 //             Ny = floor((roi.height - dctSize.height + delta.height)/delta.height);
  70 //
  71 //             So, total number of observation vectors is Nx*Ny, and total size of
  72 //             array obs must be >= Nx*Ny*obsSize.width*obsSize.height*sizeof(float).
  73 //          2. Observation vectors are calculated in the following loop
  74 //               ( actual implementation may be different ), where
  75 //               I[x1:x2,y1:y2] means block of pixels from source image with
  76 //               x1 <= x < x2, y1 <= y < y2,
  77 //               D[x1:x2,y1:y2] means sub matrix of DCT matrix D.
  78 //               O[x,y] means observation vector that corresponds to position
  79 //               (x*delta.width,y*delta.height) in the source image
  80 //               ( all indices are counted from 0 ).
  81 //
  82 //               for( y = 0; y < Ny; y++ )
  83 //               {
  84 //                   for( x = 0; x < Nx; x++ )
  85 //                   {
  86 //                       D = DCT(I[x*delta.width : x*delta.width + dctSize.width,
  87 //                                  y*delta.height : y*delta.height + dctSize.height]);
  88 //                       O[x,y] = D[0:obsSize.width, 0:obsSize.height];
  89 //                   }
  90 //               }
  91 //F*/
  92
  93 /*comment out the following line to make DCT be calculated in floating-point arithmetics*/
  94 //#define _CV_INT_DCT
  95
  96 /* for integer DCT only */
  97 #define DCT_SCALE  15
  98
  99 #ifdef _CV_INT_DCT
 100 typedef int work_t;
 101
 102 #define  DESCALE      CV_DESCALE
 103 #define  SCALE(x)     CV_FLT_TO_FIX((x),DCT_SCALE)
 104 #else
 105 typedef float work_t;
 106
 107 #define  DESCALE(x,n) (float)(x)
 108 #define  SCALE(x)     (float)(x)
 109 #endif
 110
 111 /* calculate dct transform matrix */
 112 static void icvCalcDCTMatrix( work_t * cfs, int n );
 113
 114 #define  MAX_DCT_SIZE  32
 115
 116 static CvStatus CV_STDCALL
 117 icvImgToObs_DCT_8u32f_C1R( uchar * img, int imgStep, CvSize roi,
 118                            float *obs, CvSize dctSize,
 119                            CvSize obsSize, CvSize delta )
 120 {
 121     /* dct transform matrices: horizontal and vertical */
 122     work_t tab_x[MAX_DCT_SIZE * MAX_DCT_SIZE / 2 + 2];
 123     work_t tab_y[MAX_DCT_SIZE * MAX_DCT_SIZE / 2 + 2];
 124
 125     /* temporary buffers for dct */
 126     work_t temp0[MAX_DCT_SIZE * 4];
 127     work_t temp1[MAX_DCT_SIZE * 4];
 128     work_t *buffer = 0;
 129     work_t *buf_limit;
 130
 131     double s;
 132
 133     int y;
 134     int Nx, Ny;
 135
 136     int n1 = dctSize.height, m1 = n1 / 2;
 137     int n2 = dctSize.width, m2 = n2 / 2;
 138
 139     if( !img || !obs )
 140         return CV_NULLPTR_ERR;
 141
 142     if( roi.width <= 0 || roi.height <= 0 )
 143         return CV_BADSIZE_ERR;
 144
 145     if( delta.width <= 0 || delta.height <= 0 )
 146         return CV_BADRANGE_ERR;
 147
 148     if( obsSize.width <= 0 || dctSize.width < obsSize.width ||
 149         obsSize.height <= 0 || dctSize.height < obsSize.height )
 150         return CV_BADRANGE_ERR;
 151
 152     if( dctSize.width > MAX_DCT_SIZE || dctSize.height > MAX_DCT_SIZE )
 153         return CV_BADRANGE_ERR;
 154
 155     Nx = (roi.width - dctSize.width + delta.width) / delta.width;
 156     Ny = (roi.height - dctSize.height + delta.height) / delta.height;
 157
 158     if( Nx <= 0 || Ny <= 0 )
 159         return CV_BADRANGE_ERR;
 160
 161     buffer = (work_t *)cvAlloc( roi.width * obsSize.height * sizeof( buffer[0] ));
 162     if( !buffer )
 163         return CV_OUTOFMEM_ERR;
 164
 165     icvCalcDCTMatrix( tab_x, dctSize.width );
 166     icvCalcDCTMatrix( tab_y, dctSize.height );
 167
 168     buf_limit = buffer + obsSize.height * roi.width;
 169
 170     for( y = 0; y < Ny; y++, img += delta.height * imgStep )
 171     {
 172         int x, i, j, k;
 173         work_t k0 = 0;
 174
 175         /* do transfroms for each column. Calc only first obsSize.height DCT coefficients */
 176         for( x = 0; x < roi.width; x++ )
 177         {
 178             float is = 0;
 179             work_t *buf = buffer + x;
 180             work_t *tab = tab_y + 2;
 181
 182             if( n1 & 1 )
 183             {
 184                 is = img[x + m1 * imgStep];
 185                 k0 = ((work_t) is) * tab[-1];
 186             }
 187
 188             /* first coefficient */
 189             for( j = 0; j < m1; j++ )
 190             {
 191                 float t0 = img[x + j * imgStep];
 192                 float t1 = img[x + (n1 - 1 - j) * imgStep];
 193                 float t2 = t0 + t1;
 194
 195                 t0 -= t1;
 196                 temp0[j] = (work_t) t2;
 197                 is += t2;
 198                 temp1[j] = (work_t) t0;
 199             }
 200
 201             buf[0] = DESCALE( is * tab[-2], PASS1_SHIFT );
 202             if( (buf += roi.width) >= buf_limit )
 203                 continue;
 204
 205             /* other coefficients */
 206             for( ;; )
 207             {
 208                 s = 0;
 209
 210                 for( k = 0; k < m1; k++ )
 211                     s += temp1[k] * tab[k];
 212
 213                 buf[0] = DESCALE( s, PASS1_SHIFT );
 214                 if( (buf += roi.width) >= buf_limit )
 215                     break;
 216
 217                 tab += m1;
 218                 s = 0;
 219
 220                 if( n1 & 1 )
 221                 {
 222                     k0 = -k0;
 223                     s = k0;
 224                 }
 225                 for( k = 0; k < m1; k++ )
 226                     s += temp0[k] * tab[k];
 227
 228                 buf[0] = DESCALE( s, PASS1_SHIFT );
 229                 tab += m1;
 230
 231                 if( (buf += roi.width) >= buf_limit )
 232                     break;
 233             }
 234         }
 235
 236         k0 = 0;
 237
 238         /* do transforms for rows. */
 239         for( x = 0; x + dctSize.width <= roi.width; x += delta.width )
 240         {
 241             for( i = 0; i < obsSize.height; i++ )
 242             {
 243                 work_t *buf = buffer + x + roi.width * i;
 244                 work_t *tab = tab_x + 2;
 245                 float *obs_limit = obs + obsSize.width;
 246
 247                 s = 0;
 248
 249                 if( n2 & 1 )
 250                 {
 251                     s = buf[m2];
 252                     k0 = (work_t) (s * tab[-1]);
 253                 }
 254
 255                 /* first coefficient */
 256                 for( j = 0; j < m2; j++ )
 257                 {
 258                     work_t t0 = buf[j];
 259                     work_t t1 = buf[n2 - 1 - j];
 260                     work_t t2 = t0 + t1;
 261
 262                     t0 -= t1;
 263                     temp0[j] = (work_t) t2;
 264                     s += t2;
 265                     temp1[j] = (work_t) t0;
 266                 }
 267
 268                 *obs++ = (float) DESCALE( s * tab[-2], PASS2_SHIFT );
 269
 270                 if( obs == obs_limit )
 271                     continue;
 272
 273                 /* other coefficients */
 274                 for( ;; )
 275                 {
 276                     s = 0;
 277
 278                     for( k = 0; k < m2; k++ )
 279                         s += temp1[k] * tab[k];
 280
 281                     obs[0] = (float) DESCALE( s, PASS2_SHIFT );
 282                     if( ++obs == obs_limit )
 283                         break;
 284
 285                     tab += m2;
 286
 287                     s = 0;
 288
 289                     if( n2 & 1 )
 290                     {
 291                         k0 = -k0;
 292                         s = k0;
 293                     }
 294                     for( k = 0; k < m2; k++ )
 295                         s += temp0[k] * tab[k];
 296                     obs[0] = (float) DESCALE( s, PASS2_SHIFT );
 297
 298                     tab += m2;
 299                     if( ++obs == obs_limit )
 300                         break;
 301                 }
 302             }
 303         }
 304     }
 305
 306     cvFree( &buffer );
 307     return CV_NO_ERR;
 308 }
 309
 310
 311 static CvStatus CV_STDCALL
 312 icvImgToObs_DCT_32f_C1R( float * img, int imgStep, CvSize roi,
 313                          float *obs, CvSize dctSize,
 314                          CvSize obsSize, CvSize delta )
 315 {
 316     /* dct transform matrices: horizontal and vertical */
 317     work_t tab_x[MAX_DCT_SIZE * MAX_DCT_SIZE / 2 + 2];
 318     work_t tab_y[MAX_DCT_SIZE * MAX_DCT_SIZE / 2 + 2];
 319
 320     /* temporary buffers for dct */
 321     work_t temp0[MAX_DCT_SIZE * 4];
 322     work_t temp1[MAX_DCT_SIZE * 4];
 323     work_t *buffer = 0;
 324     work_t *buf_limit;
 325
 326     double s;
 327
 328     int y;
 329     int Nx, Ny;
 330
 331     int n1 = dctSize.height, m1 = n1 / 2;
 332     int n2 = dctSize.width, m2 = n2 / 2;
 333
 334     if( !img || !obs )
 335         return CV_NULLPTR_ERR;
 336
 337     if( roi.width <= 0 || roi.height <= 0 )
 338         return CV_BADSIZE_ERR;
 339
 340     if( delta.width <= 0 || delta.height <= 0 )
 341         return CV_BADRANGE_ERR;
 342
 343     if( obsSize.width <= 0 || dctSize.width < obsSize.width ||
 344         obsSize.height <= 0 || dctSize.height < obsSize.height )
 345         return CV_BADRANGE_ERR;
 346
 347     if( dctSize.width > MAX_DCT_SIZE || dctSize.height > MAX_DCT_SIZE )
 348         return CV_BADRANGE_ERR;
 349
 350     Nx = (roi.width - dctSize.width + delta.width) / delta.width;
 351     Ny = (roi.height - dctSize.height + delta.height) / delta.height;
 352
 353     if( Nx <= 0 || Ny <= 0 )
 354         return CV_BADRANGE_ERR;
 355
 356     buffer = (work_t *)cvAlloc( roi.width * obsSize.height * sizeof( buffer[0] ));
 357     if( !buffer )
 358         return CV_OUTOFMEM_ERR;
 359
 360     icvCalcDCTMatrix( tab_x, dctSize.width );
 361     icvCalcDCTMatrix( tab_y, dctSize.height );
 362
 363     buf_limit = buffer + obsSize.height * roi.width;
 364
 365     imgStep /= sizeof(img[0]);
 366
 367     for( y = 0; y < Ny; y++, img += delta.height * imgStep )
 368     {
 369         int x, i, j, k;
 370         work_t k0 = 0;
 371
 372         /* do transfroms for each column. Calc only first obsSize.height DCT coefficients */
 373         for( x = 0; x < roi.width; x++ )
 374         {
 375             float is = 0;
 376             work_t *buf = buffer + x;
 377             work_t *tab = tab_y + 2;
 378
 379             if( n1 & 1 )
 380             {
 381                 is = img[x + m1 * imgStep];
 382                 k0 = ((work_t) is) * tab[-1];
 383             }
 384
 385             /* first coefficient */
 386             for( j = 0; j < m1; j++ )
 387             {
 388                 float t0 = img[x + j * imgStep];
 389                 float t1 = img[x + (n1 - 1 - j) * imgStep];
 390                 float t2 = t0 + t1;
 391
 392                 t0 -= t1;
 393                 temp0[j] = (work_t) t2;
 394                 is += t2;
 395                 temp1[j] = (work_t) t0;
 396             }
 397
 398             buf[0] = DESCALE( is * tab[-2], PASS1_SHIFT );
 399             if( (buf += roi.width) >= buf_limit )
 400                 continue;
 401
 402             /* other coefficients */
 403             for( ;; )
 404             {
 405                 s = 0;
 406
 407                 for( k = 0; k < m1; k++ )
 408                     s += temp1[k] * tab[k];
 409
 410                 buf[0] = DESCALE( s, PASS1_SHIFT );
 411                 if( (buf += roi.width) >= buf_limit )
 412                     break;
 413
 414                 tab += m1;
 415                 s = 0;
 416
 417                 if( n1 & 1 )
 418                 {
 419                     k0 = -k0;
 420                     s = k0;
 421                 }
 422                 for( k = 0; k < m1; k++ )
 423                     s += temp0[k] * tab[k];
 424
 425                 buf[0] = DESCALE( s, PASS1_SHIFT );
 426                 tab += m1;
 427
 428                 if( (buf += roi.width) >= buf_limit )
 429                     break;
 430             }
 431         }
 432
 433         k0 = 0;
 434
 435         /* do transforms for rows. */
 436         for( x = 0; x + dctSize.width <= roi.width; x += delta.width )
 437         {
 438             for( i = 0; i < obsSize.height; i++ )
 439             {
 440                 work_t *buf = buffer + x + roi.width * i;
 441                 work_t *tab = tab_x + 2;
 442                 float *obs_limit = obs + obsSize.width;
 443
 444                 s = 0;
 445
 446                 if( n2 & 1 )
 447                 {
 448                     s = buf[m2];
 449                     k0 = (work_t) (s * tab[-1]);
 450                 }
 451
 452                 /* first coefficient */
 453                 for( j = 0; j < m2; j++ )
 454                 {
 455                     work_t t0 = buf[j];
 456                     work_t t1 = buf[n2 - 1 - j];
 457                     work_t t2 = t0 + t1;
 458
 459                     t0 -= t1;
 460                     temp0[j] = (work_t) t2;
 461                     s += t2;
 462                     temp1[j] = (work_t) t0;
 463                 }
 464
 465                 *obs++ = (float) DESCALE( s * tab[-2], PASS2_SHIFT );
 466
 467                 if( obs == obs_limit )
 468                     continue;
 469
 470                 /* other coefficients */
 471                 for( ;; )
 472                 {
 473                     s = 0;
 474
 475                     for( k = 0; k < m2; k++ )
 476                         s += temp1[k] * tab[k];
 477
 478                     obs[0] = (float) DESCALE( s, PASS2_SHIFT );
 479                     if( ++obs == obs_limit )
 480                         break;
 481
 482                     tab += m2;
 483
 484                     s = 0;
 485
 486                     if( n2 & 1 )
 487                     {
 488                         k0 = -k0;
 489                         s = k0;
 490                     }
 491                     for( k = 0; k < m2; k++ )
 492                         s += temp0[k] * tab[k];
 493                     obs[0] = (float) DESCALE( s, PASS2_SHIFT );
 494
 495                     tab += m2;
 496                     if( ++obs == obs_limit )
 497                         break;
 498                 }
 499             }
 500         }
 501     }
 502
 503     cvFree( &buffer );
 504     return CV_NO_ERR;
 505 }
 506
 507
 508 static void
 509 icvCalcDCTMatrix( work_t * cfs, int n )
 510 {
 511     static const double sqrt2 = 1.4142135623730950488016887242097;
 512     static const double pi = 3.1415926535897932384626433832795;
 513
 514     static const double sincos[16 * 2] = {
 515         1.00000000000000000, 0.00000000000000006,
 516         0.70710678118654746, 0.70710678118654757,
 517         0.49999999999999994, 0.86602540378443871,
 518         0.38268343236508978, 0.92387953251128674,
 519         0.30901699437494740, 0.95105651629515353,
 520         0.25881904510252074, 0.96592582628906831,
 521         0.22252093395631439, 0.97492791218182362,
 522         0.19509032201612825, 0.98078528040323043,
 523         0.17364817766693033, 0.98480775301220802,
 524         0.15643446504023087, 0.98768834059513777,
 525         0.14231483827328514, 0.98982144188093268,
 526         0.13052619222005157, 0.99144486137381038,
 527         0.12053668025532305, 0.99270887409805397,
 528         0.11196447610330786, 0.99371220989324260,
 529         0.10452846326765346, 0.99452189536827329,
 530         0.09801714032956060, 0.99518472667219693,
 531     };
 532
 533 #define ROTATE( c, s, dc, ds ) \
 534     {                              \
 535         t = c*dc - s*ds;           \
 536         s = c*ds + s*dc;           \
 537         c = t;                     \
 538     }
 539
 540 #define WRITE2( j, a, b ) \
 541     {                         \
 542         cfs[j]   = SCALE(a);  \
 543         cfs2[j]  = SCALE(b);  \
 544     }
 545
 546     double t, scale = 1. / sqrt( (double)n );
 547     int i, j, m = n / 2;
 548
 549     cfs[0] = SCALE( scale );
 550     scale *= sqrt2;
 551     cfs[1] = SCALE( scale );
 552     cfs += 2 - m;
 553
 554     if( n > 1 )
 555     {
 556         double a0, b0;
 557         double da0, db0;
 558         work_t *cfs2 = cfs + m * n;
 559
 560         if( n <= 16 )
 561         {
 562             da0 = a0 = sincos[2 * n - 1];
 563             db0 = b0 = sincos[2 * n - 2];
 564         }
 565         else
 566         {
 567             t = pi / (2 * n);
 568             da0 = a0 = cos( t );
 569             db0 = b0 = sin( t );
 570         }
 571
 572         /* other rows */
 573         for( i = 1; i <= m; i++ )
 574         {
 575             double a = a0 * scale;
 576             double b = b0 * scale;
 577             double da = a0 * a0 - b0 * b0;
 578             double db = a0 * b0 + a0 * b0;
 579
 580             cfs += m;
 581             cfs2 -= m;
 582
 583             for( j = 0; j < m; j += 2 )
 584             {
 585                 WRITE2( j, a, b );
 586                 ROTATE( a, b, da, db );
 587                 if( j + 1 < m )
 588                 {
 589                     WRITE2( j + 1, a, -b );
 590                     ROTATE( a, b, da, db );
 591                 }
 592             }
 593
 594             ROTATE( a0, b0, da0, db0 );
 595         }
 596     }
 597 #undef ROTATE
 598 #undef WRITE2
 599 }
 600
 601
 602 CV_IMPL void
 603 cvImgToObs_DCT( const void* arr, float *obs, CvSize dctSize,
 604                 CvSize obsSize, CvSize delta )
 605 {
 606     CV_FUNCNAME( "cvImgToObs_DCT" );
 607
 608     __BEGIN__;
 609
 610     CvMat stub, *mat = (CvMat*)arr;
 611
 612     CV_CALL( mat = cvGetMat( arr, &stub ));
 613
 614     switch( CV_MAT_TYPE( mat->type ))
 615     {
 616     case CV_8UC1:
 617         IPPI_CALL( icvImgToObs_DCT_8u32f_C1R( mat->data.ptr, mat->step,
 618                                            cvGetMatSize(mat), obs,
 619                                            dctSize, obsSize, delta ));
 620         break;
 621     case CV_32FC1:
 622         IPPI_CALL( icvImgToObs_DCT_32f_C1R( mat->data.fl, mat->step,
 623                                            cvGetMatSize(mat), obs,
 624                                            dctSize, obsSize, delta ));
 625         break;
 626     default:
 627         CV_ERROR( CV_StsUnsupportedFormat, "" );
 628     }
 629
 630     __END__;
 631 }
 632
 633
 634 /* End of file. */