git.maemo.org Git - opencv/blob - src/cvaux/cveigenobjects.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 "_cvaux.h"
  43
  44 CvStatus CV_STDCALL
  45 icvJacobiEigens_32f(float *A, float *V, float *E, int n, float eps)
  46 {
  47     int i, j, k, ind;
  48     float *AA = A, *VV = V;
  49     double Amax, anorm = 0, ax;
  50
  51     if( A == NULL || V == NULL || E == NULL )
  52         return CV_NULLPTR_ERR;
  53     if( n <= 0 )
  54         return CV_BADSIZE_ERR;
  55     if( eps < 1.0e-7f )
  56         eps = 1.0e-7f;
  57
  58     /*-------- Prepare --------*/
  59     for( i = 0; i < n; i++, VV += n, AA += n )
  60     {
  61         for( j = 0; j < i; j++ )
  62         {
  63             double Am = AA[j];
  64
  65             anorm += Am * Am;
  66         }
  67         for( j = 0; j < n; j++ )
  68             VV[j] = 0.f;
  69         VV[i] = 1.f;
  70     }
  71
  72     anorm = sqrt( anorm + anorm );
  73     ax = anorm * eps / n;
  74     Amax = anorm;
  75
  76     while( Amax > ax )
  77     {
  78         Amax /= n;
  79         do                      /* while (ind) */
  80         {
  81             int p, q;
  82             float *V1 = V, *A1 = A;
  83
  84             ind = 0;
  85             for( p = 0; p < n - 1; p++, A1 += n, V1 += n )
  86             {
  87                 float *A2 = A + n * (p + 1), *V2 = V + n * (p + 1);
  88
  89                 for( q = p + 1; q < n; q++, A2 += n, V2 += n )
  90                 {
  91                     double x, y, c, s, c2, s2, a;
  92                     float *A3, Apq = A1[q], App, Aqq, Aip, Aiq, Vpi, Vqi;
  93
  94                     if( fabs( Apq ) < Amax )
  95                         continue;
  96
  97                     ind = 1;
  98
  99                     /*---- Calculation of rotation angle's sine & cosine ----*/
 100                     App = A1[p];
 101                     Aqq = A2[q];
 102                     y = 5.0e-1 * (App - Aqq);
 103                     x = -Apq / sqrt( (double)Apq * Apq + (double)y * y );
 104                     if( y < 0.0 )
 105                         x = -x;
 106                     s = x / sqrt( 2.0 * (1.0 + sqrt( 1.0 - (double)x * x )));
 107                     s2 = s * s;
 108                     c = sqrt( 1.0 - s2 );
 109                     c2 = c * c;
 110                     a = 2.0 * Apq * c * s;
 111
 112                     /*---- Apq annulation ----*/
 113                     A3 = A;
 114                     for( i = 0; i < p; i++, A3 += n )
 115                     {
 116                         Aip = A3[p];
 117                         Aiq = A3[q];
 118                         Vpi = V1[i];
 119                         Vqi = V2[i];
 120                         A3[p] = (float) (Aip * c - Aiq * s);
 121                         A3[q] = (float) (Aiq * c + Aip * s);
 122                         V1[i] = (float) (Vpi * c - Vqi * s);
 123                         V2[i] = (float) (Vqi * c + Vpi * s);
 124                     }
 125                     for( ; i < q; i++, A3 += n )
 126                     {
 127                         Aip = A1[i];
 128                         Aiq = A3[q];
 129                         Vpi = V1[i];
 130                         Vqi = V2[i];
 131                         A1[i] = (float) (Aip * c - Aiq * s);
 132                         A3[q] = (float) (Aiq * c + Aip * s);
 133                         V1[i] = (float) (Vpi * c - Vqi * s);
 134                         V2[i] = (float) (Vqi * c + Vpi * s);
 135                     }
 136                     for( ; i < n; i++ )
 137                     {
 138                         Aip = A1[i];
 139                         Aiq = A2[i];
 140                         Vpi = V1[i];
 141                         Vqi = V2[i];
 142                         A1[i] = (float) (Aip * c - Aiq * s);
 143                         A2[i] = (float) (Aiq * c + Aip * s);
 144                         V1[i] = (float) (Vpi * c - Vqi * s);
 145                         V2[i] = (float) (Vqi * c + Vpi * s);
 146                     }
 147                     A1[p] = (float) (App * c2 + Aqq * s2 - a);
 148                     A2[q] = (float) (App * s2 + Aqq * c2 + a);
 149                     A1[q] = A2[p] = 0.0f;
 150                 }               /*q */
 151             }                   /*p */
 152         }
 153         while( ind );
 154         Amax /= n;
 155     }                           /* while ( Amax > ax ) */
 156
 157     for( i = 0, k = 0; i < n; i++, k += n + 1 )
 158         E[i] = A[k];
 159     /*printf(" M = %d\n", M); */
 160
 161     /* -------- ordering -------- */
 162     for( i = 0; i < n; i++ )
 163     {
 164         int m = i;
 165         float Em = (float) fabs( E[i] );
 166
 167         for( j = i + 1; j < n; j++ )
 168         {
 169             float Ej = (float) fabs( E[j] );
 170
 171             m = (Em < Ej) ? j : m;
 172             Em = (Em < Ej) ? Ej : Em;
 173         }
 174         if( m != i )
 175         {
 176             int l;
 177             float b = E[i];
 178
 179             E[i] = E[m];
 180             E[m] = b;
 181             for( j = 0, k = i * n, l = m * n; j < n; j++, k++, l++ )
 182             {
 183                 b = V[k];
 184                 V[k] = V[l];
 185                 V[l] = b;
 186             }
 187         }
 188     }
 189
 190     return CV_NO_ERR;
 191 }
 192
 193 /*F///////////////////////////////////////////////////////////////////////////////////////
 194 //    Name: icvCalcCovarMatrixEx_8u32fR
 195 //    Purpose: The function calculates a covariance matrix for a group of input objects
 196 //             (images, vectors, etc.). ROI supported.
 197 //    Context:
 198 //    Parameters:  nObjects    - number of source objects
 199 //                 objects     - array of pointers to ROIs of the source objects
 200 //                 imgStep     - full width of each source object row in bytes
 201 //                 avg         - pointer to averaged object
 202 //                 avgStep     - full width of averaged object row in bytes
 203 //                 size        - ROI size of each source and averaged objects
 204 //                 covarMatrix - covariance matrix (output parameter; must be allocated
 205 //                               before call)
 206 //
 207 //    Returns: CV_NO_ERR or error code
 208 //
 209 //    Notes:
 210 //F*/
 211 static CvStatus  CV_STDCALL
 212 icvCalcCovarMatrixEx_8u32fR( int nObjects, void *input, int objStep1,
 213                              int ioFlags, int ioBufSize, uchar* buffer,
 214                              void *userData, float *avg, int avgStep,
 215                              CvSize size, float *covarMatrix )
 216 {
 217     int objStep = objStep1;
 218
 219     /* ---- TEST OF PARAMETERS ---- */
 220
 221     if( nObjects < 2 )
 222         return CV_BADFACTOR_ERR;
 223     if( ioFlags < 0 || ioFlags > 3 )
 224         return CV_BADFACTOR_ERR;
 225     if( ioFlags && ioBufSize < 1024 )
 226         return CV_BADFACTOR_ERR;
 227     if( ioFlags && buffer == NULL )
 228         return CV_NULLPTR_ERR;
 229     if( input == NULL || avg == NULL || covarMatrix == NULL )
 230         return CV_NULLPTR_ERR;
 231     if( size.width > objStep || 4 * size.width > avgStep || size.height < 1 )
 232         return CV_BADSIZE_ERR;
 233
 234     avgStep /= 4;
 235
 236     if( ioFlags & CV_EIGOBJ_INPUT_CALLBACK )    /* ==== USE INPUT CALLBACK ==== */
 237     {
 238         int nio, ngr, igr, n = size.width * size.height, mm = 0;
 239         CvCallback read_callback = ((CvInput *) & input)->callback;
 240         uchar *buffer2;
 241
 242         objStep = n;
 243         nio = ioBufSize / n;    /* number of objects in buffer */
 244         ngr = nObjects / nio;   /* number of io groups */
 245         if( nObjects % nio )
 246             mm = 1;
 247         ngr += mm;
 248
 249         buffer2 = (uchar *)cvAlloc( sizeof( uchar ) * n );
 250         if( buffer2 == NULL )
 251             return CV_OUTOFMEM_ERR;
 252
 253         for( igr = 0; igr < ngr; igr++ )
 254         {
 255             int k, l;
 256             int io, jo, imin = igr * nio, imax = imin + nio;
 257             uchar *bu1 = buffer, *bu2;
 258
 259             if( imax > nObjects )
 260                 imax = nObjects;
 261
 262             /* read igr group */
 263             for( io = imin; io < imax; io++, bu1 += n )
 264             {
 265                 CvStatus r;
 266
 267                 r = (CvStatus)read_callback( io, (void *) bu1, userData );
 268                 if( r )
 269                     return r;
 270             }
 271
 272             /* diagonal square calc */
 273             bu1 = buffer;
 274             for( io = imin; io < imax; io++, bu1 += n )
 275             {
 276                 bu2 = bu1;
 277                 for( jo = io; jo < imax; jo++, bu2 += n )
 278                 {
 279                     float w = 0.f;
 280                     float *fu = avg;
 281                     int ij = 0;
 282
 283                     for( k = 0; k < size.height; k++, fu += avgStep )
 284                         for( l = 0; l < size.width; l++, ij++ )
 285                         {
 286                             float f = fu[l], u1 = bu1[ij], u2 = bu2[ij];
 287
 288                             w += (u1 - f) * (u2 - f);
 289                         }
 290                     covarMatrix[io * nObjects + jo] = covarMatrix[jo * nObjects + io] = w;
 291                 }
 292             }
 293
 294             /* non-diagonal elements calc */
 295             for( jo = imax; jo < nObjects; jo++ )
 296             {
 297                 CvStatus r;
 298
 299                 bu1 = buffer;
 300                 bu2 = buffer2;
 301
 302                 /* read jo object */
 303                 r = (CvStatus)read_callback( jo, (void *) bu2, userData );
 304                 if( r )
 305                     return r;
 306
 307                 for( io = imin; io < imax; io++, bu1 += n )
 308                 {
 309                     float w = 0.f;
 310                     float *fu = avg;
 311                     int ij = 0;
 312
 313                     for( k = 0; k < size.height; k++, fu += avgStep )
 314                     {
 315                         for( l = 0; l < size.width - 3; l += 4, ij += 4 )
 316                         {
 317                             float f = fu[l];
 318                             uchar u1 = bu1[ij];
 319                             uchar u2 = bu2[ij];
 320
 321                             w += (u1 - f) * (u2 - f);
 322                             f = fu[l + 1];
 323                             u1 = bu1[ij + 1];
 324                             u2 = bu2[ij + 1];
 325                             w += (u1 - f) * (u2 - f);
 326                             f = fu[l + 2];
 327                             u1 = bu1[ij + 2];
 328                             u2 = bu2[ij + 2];
 329                             w += (u1 - f) * (u2 - f);
 330                             f = fu[l + 3];
 331                             u1 = bu1[ij + 3];
 332                             u2 = bu2[ij + 3];
 333                             w += (u1 - f) * (u2 - f);
 334                         }
 335                         for( ; l < size.width; l++, ij++ )
 336                         {
 337                             float f = fu[l], u1 = bu1[ij], u2 = bu2[ij];
 338
 339                             w += (u1 - f) * (u2 - f);
 340                         }
 341                     }
 342                     covarMatrix[io * nObjects + jo] = covarMatrix[jo * nObjects + io] = w;
 343                 }
 344             }
 345         }                       /* igr */
 346
 347         cvFree( &buffer2 );
 348     }                           /* if() */
 349
 350     else
 351         /* ==== NOT USE INPUT CALLBACK ==== */
 352     {
 353         int i, j;
 354         uchar **objects = (uchar **) (((CvInput *) & input)->data);
 355
 356         for( i = 0; i < nObjects; i++ )
 357         {
 358             uchar *bu = objects[i];
 359
 360             for( j = i; j < nObjects; j++ )
 361             {
 362                 int k, l;
 363                 float w = 0.f;
 364                 float *a = avg;
 365                 uchar *bu1 = bu;
 366                 uchar *bu2 = objects[j];
 367
 368                 for( k = 0; k < size.height;
 369                      k++, bu1 += objStep, bu2 += objStep, a += avgStep )
 370                 {
 371                     for( l = 0; l < size.width - 3; l += 4 )
 372                     {
 373                         float f = a[l];
 374                         uchar u1 = bu1[l];
 375                         uchar u2 = bu2[l];
 376
 377                         w += (u1 - f) * (u2 - f);
 378                         f = a[l + 1];
 379                         u1 = bu1[l + 1];
 380                         u2 = bu2[l + 1];
 381                         w += (u1 - f) * (u2 - f);
 382                         f = a[l + 2];
 383                         u1 = bu1[l + 2];
 384                         u2 = bu2[l + 2];
 385                         w += (u1 - f) * (u2 - f);
 386                         f = a[l + 3];
 387                         u1 = bu1[l + 3];
 388                         u2 = bu2[l + 3];
 389                         w += (u1 - f) * (u2 - f);
 390                     }
 391                     for( ; l < size.width; l++ )
 392                     {
 393                         float f = a[l];
 394                         uchar u1 = bu1[l];
 395                         uchar u2 = bu2[l];
 396
 397                         w += (u1 - f) * (u2 - f);
 398                     }
 399                 }
 400
 401                 covarMatrix[i * nObjects + j] = covarMatrix[j * nObjects + i] = w;
 402             }
 403         }
 404     }                           /* else */
 405
 406     return CV_NO_ERR;
 407 }
 408
 409 /*======================== end of icvCalcCovarMatrixEx_8u32fR ===========================*/
 410
 411
 412 static int
 413 icvDefaultBufferSize( void )
 414 {
 415     return 10 * 1024 * 1024;
 416 }
 417
 418 /*F///////////////////////////////////////////////////////////////////////////////////////
 419 //    Name: icvCalcEigenObjects_8u32fR
 420 //    Purpose: The function calculates an orthonormal eigen basis and a mean (averaged)
 421 //             object for a group of input objects (images, vectors, etc.). ROI supported.
 422 //    Context:
 423 //    Parameters: nObjects  - number of source objects
 424 //                input     - pointer either to array of pointers to input objects
 425 //                            or to read callback function (depending on ioFlags)
 426 //                imgStep   - full width of each source object row in bytes
 427 //                output    - pointer either to array of pointers to output eigen objects
 428 //                            or to write callback function (depending on ioFlags)
 429 //                eigStep   - full width of each eigenobject row in bytes
 430 //                size      - ROI size of each source object
 431 //                ioFlags   - input/output flags (see Notes)
 432 //                ioBufSize - input/output buffer size
 433 //                userData  - pointer to the structure which contains all necessary
 434 //                            data for the callback functions
 435 //                calcLimit - determines the calculation finish conditions
 436 //                avg       - pointer to averaged object (has the same size as ROI)
 437 //                avgStep   - full width of averaged object row in bytes
 438 //                eigVals   - pointer to corresponding eigenvalues (array of <nObjects>
 439 //                            elements in descending order)
 440 //
 441 //    Returns: CV_NO_ERR or error code
 442 //
 443 //    Notes: 1. input/output data (that is, input objects and eigen ones) may either
 444 //              be allocated in the RAM or be read from/written to the HDD (or any
 445 //              other device) by read/write callback functions. It depends on the
 446 //              value of ioFlags paramater, which may be the following:
 447 //                  CV_EIGOBJ_NO_CALLBACK, or 0;
 448 //                  CV_EIGOBJ_INPUT_CALLBACK;
 449 //                  CV_EIGOBJ_OUTPUT_CALLBACK;
 450 //                  CV_EIGOBJ_BOTH_CALLBACK, or
 451 //                            CV_EIGOBJ_INPUT_CALLBACK | CV_EIGOBJ_OUTPUT_CALLBACK.
 452 //              The callback functions as well as the user data structure must be
 453 //              developed by the user.
 454 //
 455 //           2. If ioBufSize = 0, or it's too large, the function dermines buffer size
 456 //              itself.
 457 //
 458 //           3. Depending on calcLimit parameter, calculations are finished either if
 459 //              eigenfaces number comes up to certain value or the relation of the
 460 //              current eigenvalue and the largest one comes down to certain value
 461 //              (or any of the above conditions takes place). The calcLimit->type value
 462 //              must be CV_TERMCRIT_NUMB, CV_TERMCRIT_EPS or
 463 //              CV_TERMCRIT_NUMB | CV_TERMCRIT_EPS. The function returns the real
 464 //              values calcLimit->max_iter and calcLimit->epsilon.
 465 //
 466 //           4. eigVals may be equal to NULL (if you don't need eigen values in further).
 467 //
 468 //F*/
 469 static CvStatus CV_STDCALL
 470 icvCalcEigenObjects_8u32fR( int nObjects, void* input, int objStep,
 471                             void* output, int eigStep, CvSize size,
 472                             int  ioFlags, int ioBufSize, void* userData,
 473                             CvTermCriteria* calcLimit, float* avg,
 474                             int    avgStep, float  *eigVals )
 475 {
 476     int i, j, n, iev = 0, m1 = nObjects - 1, objStep1 = objStep, eigStep1 = eigStep / 4;
 477     CvSize objSize, eigSize, avgSize;
 478     float *c = 0;
 479     float *ev = 0;
 480     float *bf = 0;
 481     uchar *buf = 0;
 482     void *buffer = 0;
 483     float m = 1.0f / (float) nObjects;
 484     CvStatus r;
 485
 486     if( m1 > calcLimit->max_iter && calcLimit->type != CV_TERMCRIT_EPS )
 487         m1 = calcLimit->max_iter;
 488
 489     /* ---- TEST OF PARAMETERS ---- */
 490
 491     if( nObjects < 2 )
 492         return CV_BADFACTOR_ERR;
 493     if( ioFlags < 0 || ioFlags > 3 )
 494         return CV_BADFACTOR_ERR;
 495     if( input == NULL || output == NULL || avg == NULL )
 496         return CV_NULLPTR_ERR;
 497     if( size.width > objStep || 4 * size.width > eigStep ||
 498         4 * size.width > avgStep || size.height < 1 )
 499         return CV_BADSIZE_ERR;
 500     if( !(ioFlags & CV_EIGOBJ_INPUT_CALLBACK) )
 501         for( i = 0; i < nObjects; i++ )
 502             if( ((uchar **) input)[i] == NULL )
 503                 return CV_NULLPTR_ERR;
 504     if( !(ioFlags & CV_EIGOBJ_OUTPUT_CALLBACK) )
 505         for( i = 0; i < m1; i++ )
 506             if( ((float **) output)[i] == NULL )
 507                 return CV_NULLPTR_ERR;
 508
 509     avgStep /= 4;
 510     eigStep /= 4;
 511
 512     if( objStep == size.width && eigStep == size.width && avgStep == size.width )
 513     {
 514         size.width *= size.height;
 515         size.height = 1;
 516         objStep = objStep1 = eigStep = eigStep1 = avgStep = size.width;
 517     }
 518     objSize = eigSize = avgSize = size;
 519
 520     if( ioFlags & CV_EIGOBJ_INPUT_CALLBACK )
 521     {
 522         objSize.width *= objSize.height;
 523         objSize.height = 1;
 524         objStep = objSize.width;
 525         objStep1 = size.width;
 526     }
 527
 528     if( ioFlags & CV_EIGOBJ_OUTPUT_CALLBACK )
 529     {
 530         eigSize.width *= eigSize.height;
 531         eigSize.height = 1;
 532         eigStep = eigSize.width;
 533         eigStep1 = size.width;
 534     }
 535
 536     n = objSize.height * objSize.width * (ioFlags & CV_EIGOBJ_INPUT_CALLBACK) +
 537         2 * eigSize.height * eigSize.width * (ioFlags & CV_EIGOBJ_OUTPUT_CALLBACK);
 538
 539     /* Buffer size determination */
 540     if( ioFlags )
 541     {
 542         int size = icvDefaultBufferSize();
 543         ioBufSize = MIN( size, n );
 544     }
 545
 546     /* memory allocation (if necesseay) */
 547
 548     if( ioFlags & CV_EIGOBJ_INPUT_CALLBACK )
 549     {
 550         buf = (uchar *) cvAlloc( sizeof( uchar ) * objSize.width );
 551         if( buf == NULL )
 552             return CV_OUTOFMEM_ERR;
 553     }
 554
 555     if( ioFlags )
 556     {
 557         buffer = (void *) cvAlloc( ioBufSize );
 558         if( buffer == NULL )
 559         {
 560             if( buf )
 561                 cvFree( &buf );
 562             return CV_OUTOFMEM_ERR;
 563         }
 564     }
 565
 566     /* Calculation of averaged object */
 567     bf = avg;
 568     for( i = 0; i < avgSize.height; i++, bf += avgStep )
 569         for( j = 0; j < avgSize.width; j++ )
 570             bf[j] = 0.f;
 571
 572     for( i = 0; i < nObjects; i++ )
 573     {
 574         int k, l;
 575         uchar *bu = (ioFlags & CV_EIGOBJ_INPUT_CALLBACK) ? buf : ((uchar **) input)[i];
 576
 577         if( ioFlags & CV_EIGOBJ_INPUT_CALLBACK )
 578         {
 579             CvCallback read_callback = ((CvInput *) & input)->callback;
 580
 581             r = (CvStatus)read_callback( i, (void *) buf, userData );
 582             if( r )
 583             {
 584                 if( buffer )
 585                     cvFree( &buffer );
 586                 if( buf )
 587                     cvFree( &buf );
 588                 return r;
 589             }
 590         }
 591
 592         bf = avg;
 593         for( k = 0; k < avgSize.height; k++, bf += avgStep, bu += objStep1 )
 594             for( l = 0; l < avgSize.width; l++ )
 595                 bf[l] += bu[l];
 596     }
 597
 598     bf = avg;
 599     for( i = 0; i < avgSize.height; i++, bf += avgStep )
 600         for( j = 0; j < avgSize.width; j++ )
 601             bf[j] *= m;
 602
 603     /* Calculation of covariance matrix */
 604     c = (float *) cvAlloc( sizeof( float ) * nObjects * nObjects );
 605
 606     if( c == NULL )
 607     {
 608         if( buffer )
 609             cvFree( &buffer );
 610         if( buf )
 611             cvFree( &buf );
 612         return CV_OUTOFMEM_ERR;
 613     }
 614
 615     r = icvCalcCovarMatrixEx_8u32fR( nObjects, input, objStep1, ioFlags, ioBufSize,
 616                                       (uchar *) buffer, userData, avg, 4 * avgStep, size, c );
 617     if( r )
 618     {
 619         cvFree( &c );
 620         if( buffer )
 621             cvFree( &buffer );
 622         if( buf )
 623             cvFree( &buf );
 624         return r;
 625     }
 626
 627     /* Calculation of eigenvalues & eigenvectors */
 628     ev = (float *) cvAlloc( sizeof( float ) * nObjects * nObjects );
 629
 630     if( ev == NULL )
 631     {
 632         cvFree( &c );
 633         if( buffer )
 634             cvFree( &buffer );
 635         if( buf )
 636             cvFree( &buf );
 637         return CV_OUTOFMEM_ERR;
 638     }
 639
 640     if( eigVals == NULL )
 641     {
 642         eigVals = (float *) cvAlloc( sizeof( float ) * nObjects );
 643
 644         if( eigVals == NULL )
 645         {
 646             cvFree( &c );
 647             cvFree( &ev );
 648             if( buffer )
 649                 cvFree( &buffer );
 650             if( buf )
 651                 cvFree( &buf );
 652             return CV_OUTOFMEM_ERR;
 653         }
 654         iev = 1;
 655     }
 656
 657     r = icvJacobiEigens_32f( c, ev, eigVals, nObjects, 0.0f );
 658     cvFree( &c );
 659     if( r )
 660     {
 661         cvFree( &ev );
 662         if( buffer )
 663             cvFree( &buffer );
 664         if( buf )
 665             cvFree( &buf );
 666         if( iev )
 667             cvFree( &eigVals );
 668         return r;
 669     }
 670
 671     /* Eigen objects number determination */
 672     if( calcLimit->type != CV_TERMCRIT_NUMBER )
 673     {
 674         for( i = 0; i < m1; i++ )
 675             if( fabs( eigVals[i] / eigVals[0] ) < calcLimit->epsilon )
 676                 break;
 677         m1 = calcLimit->max_iter = i;
 678     }
 679     else
 680         m1 = calcLimit->max_iter;
 681     calcLimit->epsilon = (float) fabs( eigVals[m1 - 1] / eigVals[0] );
 682
 683     for( i = 0; i < m1; i++ )
 684         eigVals[i] = (float) (1.0 / sqrt( (double)eigVals[i] ));
 685
 686     /* ----------------- Calculation of eigenobjects ----------------------- */
 687     if( ioFlags & CV_EIGOBJ_OUTPUT_CALLBACK )
 688     {
 689         int nio, ngr, igr;
 690
 691         nio = ioBufSize / (4 * eigSize.width);  /* number of eigen objects in buffer */
 692         ngr = m1 / nio;         /* number of io groups */
 693         if( nObjects % nio )
 694             ngr += 1;
 695
 696         for( igr = 0; igr < ngr; igr++ )
 697         {
 698             int i, io, ie, imin = igr * nio, imax = imin + nio;
 699
 700             if( imax > m1 )
 701                 imax = m1;
 702
 703             for( i = 0; i < eigSize.width * (imax - imin); i++ )
 704                 ((float *) buffer)[i] = 0.f;
 705
 706             for( io = 0; io < nObjects; io++ )
 707             {
 708                 uchar *bu = ioFlags & CV_EIGOBJ_INPUT_CALLBACK ? buf : ((uchar **) input)[io];
 709
 710                 if( ioFlags & CV_EIGOBJ_INPUT_CALLBACK )
 711                 {
 712                     CvCallback read_callback = ((CvInput *) & input)->callback;
 713
 714                     r = (CvStatus)read_callback( io, (void *) buf, userData );
 715                     if( r )
 716                     {
 717                         cvFree( &ev );
 718                         if( iev )
 719                             cvFree( &eigVals );
 720                         if( buffer )
 721                             cvFree( &buffer );
 722                         if( buf )
 723                             cvFree( &buf );
 724                         return r;
 725                     }
 726                 }
 727
 728                 for( ie = imin; ie < imax; ie++ )
 729                 {
 730                     int k, l;
 731                     uchar *bv = bu;
 732                     float e = ev[ie * nObjects + io] * eigVals[ie];
 733                     float *be = ((float *) buffer) + ((ie - imin) * eigStep);
 734
 735                     bf = avg;
 736                     for( k = 0; k < size.height; k++, bv += objStep1,
 737                          bf += avgStep, be += eigStep1 )
 738                     {
 739                         for( l = 0; l < size.width - 3; l += 4 )
 740                         {
 741                             float f = bf[l];
 742                             uchar v = bv[l];
 743
 744                             be[l] += e * (v - f);
 745                             f = bf[l + 1];
 746                             v = bv[l + 1];
 747                             be[l + 1] += e * (v - f);
 748                             f = bf[l + 2];
 749                             v = bv[l + 2];
 750                             be[l + 2] += e * (v - f);
 751                             f = bf[l + 3];
 752                             v = bv[l + 3];
 753                             be[l + 3] += e * (v - f);
 754                         }
 755                         for( ; l < size.width; l++ )
 756                             be[l] += e * (bv[l] - bf[l]);
 757                     }
 758                 }
 759             }                   /* io */
 760
 761             for( ie = imin; ie < imax; ie++ )   /* calculated eigen objects writting */
 762             {
 763                 CvCallback write_callback = ((CvInput *) & output)->callback;
 764                 float *be = ((float *) buffer) + ((ie - imin) * eigStep);
 765
 766                 r = (CvStatus)write_callback( ie, (void *) be, userData );
 767                 if( r )
 768                 {
 769                     cvFree( &ev );
 770                     if( iev )
 771                         cvFree( &eigVals );
 772                     if( buffer )
 773                         cvFree( &buffer );
 774                     if( buf )
 775                         cvFree( &buf );
 776                     return r;
 777                 }
 778             }
 779         }                       /* igr */
 780     }
 781
 782     else
 783     {
 784         int k, p, l;
 785
 786         for( i = 0; i < m1; i++ )       /* e.o. annulation */
 787         {
 788             float *be = ((float **) output)[i];
 789
 790             for( p = 0; p < eigSize.height; p++, be += eigStep )
 791                 for( l = 0; l < eigSize.width; l++ )
 792                     be[l] = 0.0f;
 793         }
 794
 795         for( k = 0; k < nObjects; k++ )
 796         {
 797             uchar *bv = (ioFlags & CV_EIGOBJ_INPUT_CALLBACK) ? buf : ((uchar **) input)[k];
 798
 799             if( ioFlags & CV_EIGOBJ_INPUT_CALLBACK )
 800             {
 801                 CvCallback read_callback = ((CvInput *) & input)->callback;
 802
 803                 r = (CvStatus)read_callback( k, (void *) buf, userData );
 804                 if( r )
 805                 {
 806                     cvFree( &ev );
 807                     if( iev )
 808                         cvFree( &eigVals );
 809                     if( buffer )
 810                         cvFree( &buffer );
 811                     if( buf )
 812                         cvFree( &buf );
 813                     return r;
 814                 }
 815             }
 816
 817             for( i = 0; i < m1; i++ )
 818             {
 819                 float v = eigVals[i] * ev[i * nObjects + k];
 820                 float *be = ((float **) output)[i];
 821                 uchar *bu = bv;
 822
 823                 bf = avg;
 824
 825                 for( p = 0; p < size.height; p++, bu += objStep1,
 826                      bf += avgStep, be += eigStep1 )
 827                 {
 828                     for( l = 0; l < size.width - 3; l += 4 )
 829                     {
 830                         float f = bf[l];
 831                         uchar u = bu[l];
 832
 833                         be[l] += v * (u - f);
 834                         f = bf[l + 1];
 835                         u = bu[l + 1];
 836                         be[l + 1] += v * (u - f);
 837                         f = bf[l + 2];
 838                         u = bu[l + 2];
 839                         be[l + 2] += v * (u - f);
 840                         f = bf[l + 3];
 841                         u = bu[l + 3];
 842                         be[l + 3] += v * (u - f);
 843                     }
 844                     for( ; l < size.width; l++ )
 845                         be[l] += v * (bu[l] - bf[l]);
 846                 }
 847             }                   /* i */
 848         }                       /* k */
 849     }                           /* else */
 850
 851     cvFree( &ev );
 852     if( iev )
 853         cvFree( &eigVals );
 854     else
 855         for( i = 0; i < m1; i++ )
 856             eigVals[i] = 1.f / (eigVals[i] * eigVals[i]);
 857     if( buffer )
 858         cvFree( &buffer );
 859     if( buf )
 860         cvFree( &buf );
 861     return CV_NO_ERR;
 862 }
 863
 864 /* --- End of icvCalcEigenObjects_8u32fR --- */
 865
 866 /*F///////////////////////////////////////////////////////////////////////////////////////
 867 //    Name: icvCalcDecompCoeff_8u32fR
 868 //    Purpose: The function calculates one decomposition coefficient of input object
 869 //             using previously calculated eigen object and the mean (averaged) object
 870 //    Context:
 871 //    Parameters:  obj     - input object
 872 //                 objStep - its step (in bytes)
 873 //                 eigObj  - pointer to eigen object
 874 //                 eigStep - its step (in bytes)
 875 //                 avg     - pointer to averaged object
 876 //                 avgStep - its step (in bytes)
 877 //                 size    - ROI size of each source object
 878 //
 879 //    Returns: decomposition coefficient value or large negative value (if error)
 880 //
 881 //    Notes:
 882 //F*/
 883 static float CV_STDCALL
 884 icvCalcDecompCoeff_8u32fR( uchar* obj, int objStep,
 885                            float *eigObj, int eigStep,
 886                            float *avg, int avgStep, CvSize size )
 887 {
 888     int i, k;
 889     float w = 0.0f;
 890
 891     if( size.width > objStep || 4 * size.width > eigStep
 892         || 4 * size.width > avgStep || size.height < 1 )
 893         return -1.0e30f;
 894     if( obj == NULL || eigObj == NULL || avg == NULL )
 895         return -1.0e30f;
 896
 897     eigStep /= 4;
 898     avgStep /= 4;
 899
 900     if( size.width == objStep && size.width == eigStep && size.width == avgStep )
 901     {
 902         size.width *= size.height;
 903         size.height = 1;
 904         objStep = eigStep = avgStep = size.width;
 905     }
 906
 907     for( i = 0; i < size.height; i++, obj += objStep, eigObj += eigStep, avg += avgStep )
 908     {
 909         for( k = 0; k < size.width - 4; k += 4 )
 910         {
 911             float o = (float) obj[k];
 912             float e = eigObj[k];
 913             float a = avg[k];
 914
 915             w += e * (o - a);
 916             o = (float) obj[k + 1];
 917             e = eigObj[k + 1];
 918             a = avg[k + 1];
 919             w += e * (o - a);
 920             o = (float) obj[k + 2];
 921             e = eigObj[k + 2];
 922             a = avg[k + 2];
 923             w += e * (o - a);
 924             o = (float) obj[k + 3];
 925             e = eigObj[k + 3];
 926             a = avg[k + 3];
 927             w += e * (o - a);
 928         }
 929         for( ; k < size.width; k++ )
 930             w += eigObj[k] * ((float) obj[k] - avg[k]);
 931     }
 932
 933     return w;
 934 }
 935
 936 /*F///////////////////////////////////////////////////////////////////////////////////////
 937 //    Names: icvEigenDecomposite_8u32fR
 938 //    Purpose: The function calculates all decomposition coefficients for input object
 939 //             using previously calculated eigen objects basis and the mean (averaged)
 940 //             object
 941 //    Context:
 942 //    Parameters:  obj         - input object
 943 //                 objStep     - its step (in bytes)
 944 //                 nEigObjs    - number of eigen objects
 945 //                 eigInput    - pointer either to array of pointers to eigen objects
 946 //                               or to read callback function (depending on ioFlags)
 947 //                 eigStep     - eigen objects step (in bytes)
 948 //                 ioFlags     - input/output flags
 949 //                 iserData    - pointer to the structure which contains all necessary
 950 //                               data for the callback function
 951 //                 avg         - pointer to averaged object
 952 //                 avgStep     - its step (in bytes)
 953 //                 size        - ROI size of each source object
 954 //                 coeffs      - calculated coefficients (output data)
 955 //
 956 //    Returns: icv status
 957 //
 958 //    Notes:   see notes for icvCalcEigenObjects_8u32fR function
 959 //F*/
 960 static CvStatus CV_STDCALL
 961 icvEigenDecomposite_8u32fR( uchar * obj, int objStep, int nEigObjs,
 962                             void *eigInput, int eigStep, int ioFlags,
 963                             void *userData, float *avg, int avgStep,
 964                             CvSize size, float *coeffs )
 965 {
 966     int i;
 967
 968     if( nEigObjs < 2 )
 969         return CV_BADFACTOR_ERR;
 970     if( ioFlags < 0 || ioFlags > 1 )
 971         return CV_BADFACTOR_ERR;
 972     if( size.width > objStep || 4 * size.width > eigStep ||
 973         4 * size.width > avgStep || size.height < 1 )
 974         return CV_BADSIZE_ERR;
 975     if( obj == NULL || eigInput == NULL || coeffs == NULL || avg == NULL )
 976         return CV_NULLPTR_ERR;
 977     if( !ioFlags )
 978         for( i = 0; i < nEigObjs; i++ )
 979             if( ((uchar **) eigInput)[i] == NULL )
 980                 return CV_NULLPTR_ERR;
 981
 982     if( ioFlags )               /* callback */
 983
 984     {
 985         float *buffer;
 986         CvCallback read_callback = ((CvInput *) & eigInput)->callback;
 987
 988         eigStep = 4 * size.width;
 989
 990         /* memory allocation */
 991         buffer = (float *) cvAlloc( sizeof( float ) * size.width * size.height );
 992
 993         if( buffer == NULL )
 994             return CV_OUTOFMEM_ERR;
 995
 996         for( i = 0; i < nEigObjs; i++ )
 997         {
 998             float w;
 999             CvStatus r = (CvStatus)read_callback( i, (void *) buffer, userData );
1000
1001             if( r )
1002             {
1003                 cvFree( &buffer );
1004                 return r;
1005             }
1006             w = icvCalcDecompCoeff_8u32fR( obj, objStep, buffer,
1007                                             eigStep, avg, avgStep, size );
1008             if( w < -1.0e29f )
1009             {
1010                 cvFree( &buffer );
1011                 return CV_NOTDEFINED_ERR;
1012             }
1013             coeffs[i] = w;
1014         }
1015         cvFree( &buffer );
1016     }
1017
1018     else
1019         /* no callback */
1020         for( i = 0; i < nEigObjs; i++ )
1021         {
1022             float w = icvCalcDecompCoeff_8u32fR( obj, objStep, ((float **) eigInput)[i],
1023                                                   eigStep, avg, avgStep, size );
1024
1025             if( w < -1.0e29f )
1026                 return CV_NOTDEFINED_ERR;
1027             coeffs[i] = w;
1028         }
1029
1030     return CV_NO_ERR;
1031 }
1032
1033
1034 /*F///////////////////////////////////////////////////////////////////////////////////////
1035 //    Names: icvEigenProjection_8u32fR
1036 //    Purpose: The function calculates object projection to the eigen sub-space (restores
1037 //             an object) using previously calculated eigen objects basis, mean (averaged)
1038 //             object and decomposition coefficients of the restored object
1039 //    Context:
1040 //    Parameters:  nEigObjs - Number of eigen objects
1041 //                 eigens   - Array of pointers to eigen objects
1042 //                 eigStep  - Eigen objects step (in bytes)
1043 //                 coeffs   - Previously calculated decomposition coefficients
1044 //                 avg      - Pointer to averaged object
1045 //                 avgStep  - Its step (in bytes)
1046 //                 rest     - Pointer to restored object
1047 //                 restStep - Its step (in bytes)
1048 //                 size     - ROI size of each object
1049 //
1050 //    Returns: CV status
1051 //
1052 //    Notes:
1053 //F*/
1054 static CvStatus CV_STDCALL
1055 icvEigenProjection_8u32fR( int nEigObjs, void *eigInput, int eigStep,
1056                            int ioFlags, void *userData, float *coeffs,
1057                            float *avg, int avgStep, uchar * rest,
1058                            int restStep, CvSize size )
1059 {
1060     int i, j, k;
1061     float *buf;
1062     float *buffer = NULL;
1063     float *b;
1064     CvCallback read_callback = ((CvInput *) & eigInput)->callback;
1065
1066     if( size.width > avgStep || 4 * size.width > eigStep || size.height < 1 )
1067         return CV_BADSIZE_ERR;
1068     if( rest == NULL || eigInput == NULL || avg == NULL || coeffs == NULL )
1069         return CV_NULLPTR_ERR;
1070     if( ioFlags < 0 || ioFlags > 1 )
1071         return CV_BADFACTOR_ERR;
1072     if( !ioFlags )
1073         for( i = 0; i < nEigObjs; i++ )
1074             if( ((uchar **) eigInput)[i] == NULL )
1075                 return CV_NULLPTR_ERR;
1076     eigStep /= 4;
1077     avgStep /= 4;
1078
1079     if( size.width == restStep && size.width == eigStep && size.width == avgStep )
1080     {
1081         size.width *= size.height;
1082         size.height = 1;
1083         restStep = eigStep = avgStep = size.width;
1084     }
1085
1086     buf = (float *) cvAlloc( sizeof( float ) * size.width * size.height );
1087
1088     if( buf == NULL )
1089         return CV_OUTOFMEM_ERR;
1090     b = buf;
1091     for( i = 0; i < size.height; i++, avg += avgStep, b += size.width )
1092         for( j = 0; j < size.width; j++ )
1093             b[j] = avg[j];
1094
1095     if( ioFlags )
1096     {
1097         buffer = (float *) cvAlloc( sizeof( float ) * size.width * size.height );
1098
1099         if( buffer == NULL )
1100         {
1101             cvFree( &buf );
1102             return CV_OUTOFMEM_ERR;
1103         }
1104         eigStep = size.width;
1105     }
1106
1107     for( k = 0; k < nEigObjs; k++ )
1108     {
1109         float *e = ioFlags ? buffer : ((float **) eigInput)[k];
1110         float c = coeffs[k];
1111
1112         if( ioFlags )           /* read eigen object */
1113         {
1114             CvStatus r = (CvStatus)read_callback( k, (void *) buffer, userData );
1115
1116             if( r )
1117             {
1118                 cvFree( &buf );
1119                 cvFree( &buffer );
1120                 return r;
1121             }
1122         }
1123
1124         b = buf;
1125         for( i = 0; i < size.height; i++, e += eigStep, b += size.width )
1126         {
1127             for( j = 0; j < size.width - 3; j += 4 )
1128             {
1129                 float b0 = c * e[j];
1130                 float b1 = c * e[j + 1];
1131                 float b2 = c * e[j + 2];
1132                 float b3 = c * e[j + 3];
1133
1134                 b[j] += b0;
1135                 b[j + 1] += b1;
1136                 b[j + 2] += b2;
1137                 b[j + 3] += b3;
1138             }
1139             for( ; j < size.width; j++ )
1140                 b[j] += c * e[j];
1141         }
1142     }
1143
1144     b = buf;
1145     for( i = 0; i < size.height; i++, avg += avgStep, b += size.width, rest += restStep )
1146         for( j = 0; j < size.width; j++ )
1147         {
1148             int w = cvRound( b[j] );
1149
1150             w = !(w & ~255) ? w : w < 0 ? 0 : 255;
1151             rest[j] = (uchar) w;
1152         }
1153
1154     cvFree( &buf );
1155     if( ioFlags )
1156         cvFree( &buffer );
1157     return CV_NO_ERR;
1158 }
1159
1160 /*F///////////////////////////////////////////////////////////////////////////////////////
1161 //    Name: cvCalcCovarMatrixEx
1162 //    Purpose: The function calculates a covariance matrix for a group of input objects
1163 //             (images, vectors, etc.).
1164 //    Context:
1165 //    Parameters:  nObjects    - number of source objects
1166 //                 input       - pointer either to array of input objects
1167 //                               or to read callback function (depending on ioFlags)
1168 //                 ioFlags     - input/output flags (see Notes to
1169 //                               cvCalcEigenObjects function)
1170 //                 ioBufSize   - input/output buffer size
1171 //                 userData    - pointer to the structure which contains all necessary
1172 //                               data for the callback functions
1173 //                 avg         - averaged object
1174 //                 covarMatrix - covariance matrix (output parameter; must be allocated
1175 //                               before call)
1176 //
1177 //    Notes:  See Notes to cvCalcEigenObjects function
1178 //F*/
1179
1180 CV_IMPL void
1181 cvCalcCovarMatrixEx( int  nObjects, void*  input, int  ioFlags,
1182                      int  ioBufSize, uchar*  buffer, void*  userData,
1183                      IplImage* avg, float*  covarMatrix )
1184 {
1185     float *avg_data;
1186     int avg_step = 0;
1187     CvSize avg_size;
1188     int i;
1189
1190     CV_FUNCNAME( "cvCalcCovarMatrixEx" );
1191
1192     __BEGIN__;
1193
1194     cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &avg_size );
1195     if( avg->depth != IPL_DEPTH_32F )
1196         CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1197     if( avg->nChannels != 1 )
1198         CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1199
1200     if( ioFlags == CV_EIGOBJ_NO_CALLBACK )
1201     {
1202         IplImage **images = (IplImage **) (((CvInput *) & input)->data);
1203         uchar **objects = (uchar **) cvAlloc( sizeof( uchar * ) * nObjects );
1204         int img_step = 0, old_step = 0;
1205         CvSize img_size = avg_size, old_size = avg_size;
1206
1207         if( objects == NULL )
1208             CV_ERROR( CV_StsBadArg, "Insufficient memory" );
1209
1210         for( i = 0; i < nObjects; i++ )
1211         {
1212             IplImage *img = images[i];
1213             uchar *img_data;
1214
1215             cvGetImageRawData( img, &img_data, &img_step, &img_size );
1216             if( img->depth != IPL_DEPTH_8U )
1217                 CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1218             if( img_size != avg_size || img_size != old_size )
1219                 CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
1220             if( img->nChannels != 1 )
1221                 CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1222             if( i > 0 && img_step != old_step )
1223                 CV_ERROR( CV_StsBadArg, "Different steps of objects" );
1224
1225             old_step = img_step;
1226             old_size = img_size;
1227             objects[i] = img_data;
1228         }
1229
1230         CV_CALL( icvCalcCovarMatrixEx_8u32fR( nObjects,
1231                                               (void*) objects,
1232                                               img_step,
1233                                               CV_EIGOBJ_NO_CALLBACK,
1234                                               0,
1235                                               NULL,
1236                                               NULL,
1237                                               avg_data,
1238                                               avg_step,
1239                                               avg_size,
1240                                               covarMatrix ));
1241         cvFree( &objects );
1242     }
1243
1244     else
1245
1246     {
1247         CV_CALL( icvCalcCovarMatrixEx_8u32fR( nObjects,
1248                                               input,
1249                                               avg_step / 4,
1250                                               ioFlags,
1251                                               ioBufSize,
1252                                               buffer,
1253                                               userData,
1254                                               avg_data,
1255                                               avg_step,
1256                                               avg_size,
1257                                               covarMatrix ));
1258     }
1259
1260     __END__;
1261 }
1262
1263 /*F///////////////////////////////////////////////////////////////////////////////////////
1264 //    Name: cvCalcEigenObjects
1265 //    Purpose: The function calculates an orthonormal eigen basis and a mean (averaged)
1266 //             object for a group of input objects (images, vectors, etc.).
1267 //    Context:
1268 //    Parameters: nObjects  - number of source objects
1269 //                input     - pointer either to array of input objects
1270 //                            or to read callback function (depending on ioFlags)
1271 //                output    - pointer either to output eigen objects
1272 //                            or to write callback function (depending on ioFlags)
1273 //                ioFlags   - input/output flags (see Notes)
1274 //                ioBufSize - input/output buffer size
1275 //                userData  - pointer to the structure which contains all necessary
1276 //                            data for the callback functions
1277 //                calcLimit - determines the calculation finish conditions
1278 //                avg       - averaged object (has the same size as ROI)
1279 //                eigVals   - pointer to corresponding eigen values (array of <nObjects>
1280 //                            elements in descending order)
1281 //
1282 //    Notes: 1. input/output data (that is, input objects and eigen ones) may either
1283 //              be allocated in the RAM or be read from/written to the HDD (or any
1284 //              other device) by read/write callback functions. It depends on the
1285 //              value of ioFlags paramater, which may be the following:
1286 //                  CV_EIGOBJ_NO_CALLBACK, or 0;
1287 //                  CV_EIGOBJ_INPUT_CALLBACK;
1288 //                  CV_EIGOBJ_OUTPUT_CALLBACK;
1289 //                  CV_EIGOBJ_BOTH_CALLBACK, or
1290 //                            CV_EIGOBJ_INPUT_CALLBACK | CV_EIGOBJ_OUTPUT_CALLBACK.
1291 //              The callback functions as well as the user data structure must be
1292 //              developed by the user.
1293 //
1294 //           2. If ioBufSize = 0, or it's too large, the function dermines buffer size
1295 //              itself.
1296 //
1297 //           3. Depending on calcLimit parameter, calculations are finished either if
1298 //              eigenfaces number comes up to certain value or the relation of the
1299 //              current eigenvalue and the largest one comes down to certain value
1300 //              (or any of the above conditions takes place). The calcLimit->type value
1301 //              must be CV_TERMCRIT_NUMB, CV_TERMCRIT_EPS or
1302 //              CV_TERMCRIT_NUMB | CV_TERMCRIT_EPS. The function returns the real
1303 //              values calcLimit->max_iter and calcLimit->epsilon.
1304 //
1305 //           4. eigVals may be equal to NULL (if you don't need eigen values in further).
1306 //
1307 //F*/
1308 CV_IMPL void
1309 cvCalcEigenObjects( int       nObjects,
1310                     void*     input,
1311                     void*     output,
1312                     int       ioFlags,
1313                     int       ioBufSize,
1314                     void*     userData,
1315                     CvTermCriteria* calcLimit,
1316                     IplImage* avg,
1317                     float*    eigVals )
1318 {
1319     float *avg_data;
1320     int avg_step = 0;
1321     CvSize avg_size;
1322     int i;
1323     int nEigens = nObjects - 1;
1324
1325     CV_FUNCNAME( "cvCalcEigenObjects" );
1326
1327     __BEGIN__;
1328
1329     cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &avg_size );
1330     if( avg->depth != IPL_DEPTH_32F )
1331         CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1332     if( avg->nChannels != 1 )
1333         CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1334
1335     if( nEigens > calcLimit->max_iter && calcLimit->type != CV_TERMCRIT_EPS )
1336         nEigens = calcLimit->max_iter;
1337
1338     switch (ioFlags)
1339     {
1340     case CV_EIGOBJ_NO_CALLBACK:
1341         {
1342             IplImage **objects = (IplImage **) (((CvInput *) & input)->data);
1343             IplImage **eigens = (IplImage **) (((CvInput *) & output)->data);
1344             uchar **objs = (uchar **) cvAlloc( sizeof( uchar * ) * nObjects );
1345             float **eigs = (float **) cvAlloc( sizeof( float * ) * nEigens );
1346             int obj_step = 0, old_step = 0;
1347             int eig_step = 0, oldeig_step = 0;
1348             CvSize obj_size = avg_size, old_size = avg_size,
1349
1350                 eig_size = avg_size, oldeig_size = avg_size;
1351
1352             if( objects == NULL || eigens == NULL )
1353                 CV_ERROR( CV_StsBadArg, "Insufficient memory" );
1354
1355             for( i = 0; i < nObjects; i++ )
1356             {
1357                 IplImage *img = objects[i];
1358                 uchar *obj_data;
1359
1360                 cvGetImageRawData( img, &obj_data, &obj_step, &obj_size );
1361                 if( img->depth != IPL_DEPTH_8U )
1362                     CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1363                 if( obj_size != avg_size || obj_size != old_size )
1364                     CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
1365                 if( img->nChannels != 1 )
1366                     CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1367                 if( i > 0 && obj_step != old_step )
1368                     CV_ERROR( CV_StsBadArg, "Different steps of objects" );
1369
1370                 old_step = obj_step;
1371                 old_size = obj_size;
1372                 objs[i] = obj_data;
1373             }
1374             for( i = 0; i < nEigens; i++ )
1375             {
1376                 IplImage *eig = eigens[i];
1377                 float *eig_data;
1378
1379                 cvGetImageRawData( eig, (uchar **) & eig_data, &eig_step, &eig_size );
1380                 if( eig->depth != IPL_DEPTH_32F )
1381                     CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1382                 if( eig_size != avg_size || eig_size != oldeig_size )
1383                     CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
1384                 if( eig->nChannels != 1 )
1385                     CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1386                 if( i > 0 && eig_step != oldeig_step )
1387                     CV_ERROR( CV_StsBadArg, "Different steps of objects" );
1388
1389                 oldeig_step = eig_step;
1390                 oldeig_size = eig_size;
1391                 eigs[i] = eig_data;
1392             }
1393             CV_CALL( icvCalcEigenObjects_8u32fR( nObjects, (void*) objs, obj_step,
1394                                                  (void*) eigs, eig_step, obj_size,
1395                                                  ioFlags, ioBufSize, userData,
1396                                                  calcLimit, avg_data, avg_step, eigVals ));
1397             cvFree( &objs );
1398             cvFree( &eigs );
1399             break;
1400         }
1401
1402     case CV_EIGOBJ_OUTPUT_CALLBACK:
1403         {
1404             IplImage **objects = (IplImage **) (((CvInput *) & input)->data);
1405             uchar **objs = (uchar **) cvAlloc( sizeof( uchar * ) * nObjects );
1406             int obj_step = 0, old_step = 0;
1407             CvSize obj_size = avg_size, old_size = avg_size;
1408
1409             if( objects == NULL )
1410                 CV_ERROR( CV_StsBadArg, "Insufficient memory" );
1411
1412             for( i = 0; i < nObjects; i++ )
1413             {
1414                 IplImage *img = objects[i];
1415                 uchar *obj_data;
1416
1417                 cvGetImageRawData( img, &obj_data, &obj_step, &obj_size );
1418                 if( img->depth != IPL_DEPTH_8U )
1419                     CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1420                 if( obj_size != avg_size || obj_size != old_size )
1421                     CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
1422                 if( img->nChannels != 1 )
1423                     CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1424                 if( i > 0 && obj_step != old_step )
1425                     CV_ERROR( CV_StsBadArg, "Different steps of objects" );
1426
1427                 old_step = obj_step;
1428                 old_size = obj_size;
1429                 objs[i] = obj_data;
1430             }
1431             CV_CALL( icvCalcEigenObjects_8u32fR( nObjects,
1432                                                  (void*) objs,
1433                                                  obj_step,
1434                                                  output,
1435                                                  avg_step,
1436                                                  obj_size,
1437                                                  ioFlags,
1438                                                  ioBufSize,
1439                                                  userData,
1440                                                  calcLimit,
1441                                                  avg_data,
1442                                                  avg_step,
1443                                                  eigVals   ));
1444             cvFree( &objs );
1445             break;
1446         }
1447
1448     case CV_EIGOBJ_INPUT_CALLBACK:
1449         {
1450             IplImage **eigens = (IplImage **) (((CvInput *) & output)->data);
1451             float **eigs = (float**) cvAlloc( sizeof( float* ) * nEigens );
1452             int eig_step = 0, oldeig_step = 0;
1453             CvSize eig_size = avg_size, oldeig_size = avg_size;
1454
1455             if( eigens == NULL )
1456                 CV_ERROR( CV_StsBadArg, "Insufficient memory" );
1457
1458             for( i = 0; i < nEigens; i++ )
1459             {
1460                 IplImage *eig = eigens[i];
1461                 float *eig_data;
1462
1463                 cvGetImageRawData( eig, (uchar **) & eig_data, &eig_step, &eig_size );
1464                 if( eig->depth != IPL_DEPTH_32F )
1465                     CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1466                 if( eig_size != avg_size || eig_size != oldeig_size )
1467                     CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
1468                 if( eig->nChannels != 1 )
1469                     CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1470                 if( i > 0 && eig_step != oldeig_step )
1471                     CV_ERROR( CV_StsBadArg, "Different steps of objects" );
1472
1473                 oldeig_step = eig_step;
1474                 oldeig_size = eig_size;
1475                 eigs[i] = eig_data;
1476             }
1477             CV_CALL( icvCalcEigenObjects_8u32fR( nObjects,
1478                                                  input,
1479                                                  avg_step / 4,
1480                                                  (void*) eigs,
1481                                                  eig_step,
1482                                                  eig_size,
1483                                                  ioFlags,
1484                                                  ioBufSize,
1485                                                  userData,
1486                                                  calcLimit,
1487                                                  avg_data,
1488                                                  avg_step,
1489                                                  eigVals   ));
1490             cvFree( &eigs );
1491             break;
1492         }
1493     case CV_EIGOBJ_INPUT_CALLBACK | CV_EIGOBJ_OUTPUT_CALLBACK:
1494
1495         CV_CALL( icvCalcEigenObjects_8u32fR( nObjects,
1496                                              input,
1497                                              avg_step / 4,
1498                                              output,
1499                                              avg_step,
1500                                              avg_size,
1501                                              ioFlags,
1502                                              ioBufSize,
1503                                              userData,
1504                                              calcLimit,
1505                                              avg_data,
1506                                              avg_step,
1507                                              eigVals   ));
1508         break;
1509
1510     default:
1511         CV_ERROR( CV_StsBadArg, "Unsupported i/o flag" );
1512     }
1513
1514     __END__;
1515 }
1516
1517 /*--------------------------------------------------------------------------------------*/
1518 /*F///////////////////////////////////////////////////////////////////////////////////////
1519 //    Name: cvCalcDecompCoeff
1520 //    Purpose: The function calculates one decomposition coefficient of input object
1521 //             using previously calculated eigen object and the mean (averaged) object
1522 //    Context:
1523 //    Parameters:  obj     - input object
1524 //                 eigObj  - eigen object
1525 //                 avg     - averaged object
1526 //
1527 //    Returns: decomposition coefficient value or large negative value (if error)
1528 //
1529 //    Notes:
1530 //F*/
1531
1532 CV_IMPL double
1533 cvCalcDecompCoeff( IplImage * obj, IplImage * eigObj, IplImage * avg )
1534 {
1535     double coeff = DBL_MAX;
1536
1537     uchar *obj_data;
1538     float *eig_data;
1539     float *avg_data;
1540     int obj_step = 0, eig_step = 0, avg_step = 0;
1541     CvSize obj_size, eig_size, avg_size;
1542
1543     CV_FUNCNAME( "cvCalcDecompCoeff" );
1544
1545     __BEGIN__;
1546
1547     cvGetImageRawData( obj, &obj_data, &obj_step, &obj_size );
1548     if( obj->depth != IPL_DEPTH_8U )
1549         CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1550     if( obj->nChannels != 1 )
1551         CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1552
1553     cvGetImageRawData( eigObj, (uchar **) & eig_data, &eig_step, &eig_size );
1554     if( eigObj->depth != IPL_DEPTH_32F )
1555         CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1556     if( eigObj->nChannels != 1 )
1557         CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1558
1559     cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &avg_size );
1560     if( avg->depth != IPL_DEPTH_32F )
1561         CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1562     if( avg->nChannels != 1 )
1563         CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1564
1565     if( obj_size != eig_size || obj_size != avg_size )
1566         CV_ERROR( CV_StsBadArg, "different sizes of images" );
1567
1568     coeff = icvCalcDecompCoeff_8u32fR( obj_data, obj_step,
1569                                        eig_data, eig_step,
1570                                        avg_data, avg_step, obj_size );
1571
1572     __END__;
1573
1574     return coeff;
1575 }
1576
1577 /*--------------------------------------------------------------------------------------*/
1578 /*F///////////////////////////////////////////////////////////////////////////////////////
1579 //    Names: cvEigenDecomposite
1580 //    Purpose: The function calculates all decomposition coefficients for input object
1581 //             using previously calculated eigen objects basis and the mean (averaged)
1582 //             object
1583 //
1584 //    Parameters:  obj         - input object
1585 //                 nEigObjs    - number of eigen objects
1586 //                 eigInput    - pointer either to array of pointers to eigen objects
1587 //                               or to read callback function (depending on ioFlags)
1588 //                 ioFlags     - input/output flags
1589 //                 userData    - pointer to the structure which contains all necessary
1590 //                               data for the callback function
1591 //                 avg         - averaged object
1592 //                 coeffs      - calculated coefficients (output data)
1593 //
1594 //    Notes:   see notes for cvCalcEigenObjects function
1595 //F*/
1596
1597 CV_IMPL void
1598 cvEigenDecomposite( IplImage* obj,
1599                     int       nEigObjs,
1600                     void*     eigInput,
1601                     int       ioFlags,
1602                     void*     userData,
1603                     IplImage* avg,
1604                     float*    coeffs )
1605 {
1606     float *avg_data;
1607     uchar *obj_data;
1608     int avg_step = 0, obj_step = 0;
1609     CvSize avg_size, obj_size;
1610     int i;
1611
1612     CV_FUNCNAME( "cvEigenDecomposite" );
1613
1614     __BEGIN__;
1615
1616     cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &avg_size );
1617     if( avg->depth != IPL_DEPTH_32F )
1618         CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1619     if( avg->nChannels != 1 )
1620         CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1621
1622     cvGetImageRawData( obj, &obj_data, &obj_step, &obj_size );
1623     if( obj->depth != IPL_DEPTH_8U )
1624         CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1625     if( obj->nChannels != 1 )
1626         CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1627
1628     if( obj_size != avg_size )
1629         CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
1630
1631     if( ioFlags == CV_EIGOBJ_NO_CALLBACK )
1632     {
1633         IplImage **eigens = (IplImage **) (((CvInput *) & eigInput)->data);
1634         float **eigs = (float **) cvAlloc( sizeof( float * ) * nEigObjs );
1635         int eig_step = 0, old_step = 0;
1636         CvSize eig_size = avg_size, old_size = avg_size;
1637
1638         if( eigs == NULL )
1639             CV_ERROR( CV_StsBadArg, "Insufficient memory" );
1640
1641         for( i = 0; i < nEigObjs; i++ )
1642         {
1643             IplImage *eig = eigens[i];
1644             float *eig_data;
1645
1646             cvGetImageRawData( eig, (uchar **) & eig_data, &eig_step, &eig_size );
1647             if( eig->depth != IPL_DEPTH_32F )
1648                 CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1649             if( eig_size != avg_size || eig_size != old_size )
1650                 CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
1651             if( eig->nChannels != 1 )
1652                 CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1653             if( i > 0 && eig_step != old_step )
1654                 CV_ERROR( CV_StsBadArg, "Different steps of objects" );
1655
1656             old_step = eig_step;
1657             old_size = eig_size;
1658             eigs[i] = eig_data;
1659         }
1660
1661         CV_CALL( icvEigenDecomposite_8u32fR( obj_data,
1662                                              obj_step,
1663                                              nEigObjs,
1664                                              (void*) eigs,
1665                                              eig_step,
1666                                              ioFlags,
1667                                              userData,
1668                                              avg_data,
1669                                              avg_step,
1670                                              obj_size,
1671                                              coeffs   ));
1672         cvFree( &eigs );
1673     }
1674
1675     else
1676
1677     {
1678         CV_CALL( icvEigenDecomposite_8u32fR( obj_data,
1679                                              obj_step,
1680                                              nEigObjs,
1681                                              eigInput,
1682                                              avg_step,
1683                                              ioFlags,
1684                                              userData,
1685                                              avg_data,
1686                                              avg_step,
1687                                              obj_size,
1688                                              coeffs   ));
1689     }
1690
1691     __END__;
1692 }
1693
1694 /*--------------------------------------------------------------------------------------*/
1695 /*F///////////////////////////////////////////////////////////////////////////////////////
1696 //    Name: cvEigenProjection
1697 //    Purpose: The function calculates object projection to the eigen sub-space (restores
1698 //             an object) using previously calculated eigen objects basis, mean (averaged)
1699 //             object and decomposition coefficients of the restored object
1700 //    Context:
1701 //    Parameters:  nEigObjs    - number of eigen objects
1702 //                 eigInput    - pointer either to array of pointers to eigen objects
1703 //                               or to read callback function (depending on ioFlags)
1704 //                 ioFlags     - input/output flags
1705 //                 userData    - pointer to the structure which contains all necessary
1706 //                               data for the callback function
1707 //                 coeffs      - array of decomposition coefficients
1708 //                 avg         - averaged object
1709 //                 proj        - object projection (output data)
1710 //
1711 //    Notes:   see notes for cvCalcEigenObjects function
1712 //F*/
1713
1714 CV_IMPL void
1715 cvEigenProjection( void*     eigInput,
1716                    int       nEigObjs,
1717                    int       ioFlags,
1718                    void*     userData,
1719                    float*    coeffs,
1720                    IplImage* avg,
1721                    IplImage* proj )
1722 {
1723     float *avg_data;
1724     uchar *proj_data;
1725     int avg_step = 0, proj_step = 0;
1726     CvSize avg_size, proj_size;
1727     int i;
1728
1729     CV_FUNCNAME( "cvEigenProjection" );
1730
1731     __BEGIN__;
1732
1733     cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &avg_size );
1734     if( avg->depth != IPL_DEPTH_32F )
1735         CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1736     if( avg->nChannels != 1 )
1737         CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1738
1739     cvGetImageRawData( proj, &proj_data, &proj_step, &proj_size );
1740     if( proj->depth != IPL_DEPTH_8U )
1741         CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1742     if( proj->nChannels != 1 )
1743         CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1744
1745     if( proj_size != avg_size )
1746         CV_ERROR( CV_StsBadArg, "Different sizes of projects" );
1747
1748     if( ioFlags == CV_EIGOBJ_NO_CALLBACK )
1749     {
1750         IplImage **eigens = (IplImage**) (((CvInput *) & eigInput)->data);
1751         float **eigs = (float**) cvAlloc( sizeof( float * ) * nEigObjs );
1752         int eig_step = 0, old_step = 0;
1753         CvSize eig_size = avg_size, old_size = avg_size;
1754
1755         if( eigs == NULL )
1756             CV_ERROR( CV_StsBadArg, "Insufficient memory" );
1757
1758         for( i = 0; i < nEigObjs; i++ )
1759         {
1760             IplImage *eig = eigens[i];
1761             float *eig_data;
1762
1763             cvGetImageRawData( eig, (uchar **) & eig_data, &eig_step, &eig_size );
1764             if( eig->depth != IPL_DEPTH_32F )
1765                 CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
1766             if( eig_size != avg_size || eig_size != old_size )
1767                 CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
1768             if( eig->nChannels != 1 )
1769                 CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
1770             if( i > 0 && eig_step != old_step )
1771                 CV_ERROR( CV_StsBadArg, "Different steps of objects" );
1772
1773             old_step = eig_step;
1774             old_size = eig_size;
1775             eigs[i] = eig_data;
1776         }
1777
1778         CV_CALL( icvEigenProjection_8u32fR( nEigObjs,
1779                                             (void*) eigs,
1780                                             eig_step,
1781                                             ioFlags,
1782                                             userData,
1783                                             coeffs,
1784                                             avg_data,
1785                                             avg_step,
1786                                             proj_data,
1787                                             proj_step,
1788                                             avg_size   ));
1789         cvFree( &eigs );
1790     }
1791
1792     else
1793
1794     {
1795         CV_CALL( icvEigenProjection_8u32fR( nEigObjs,
1796                                             eigInput,
1797                                             avg_step,
1798                                             ioFlags,
1799                                             userData,
1800                                             coeffs,
1801                                             avg_data,
1802                                             avg_step,
1803                                             proj_data,
1804                                             proj_step,
1805                                             avg_size   ));
1806     }
1807
1808     __END__;
1809 }
1810
1811 /* End of file. */