modules/legacy/src/eigenobjects.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 "precomp.hpp"
  43
  44 static CvStatus
  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         ioBufSize = MIN( icvDefaultBufferSize(), n );
 543     }
 544
 545     /* memory allocation (if necesseay) */
 546
 547     if( ioFlags & CV_EIGOBJ_INPUT_CALLBACK )
 548     {
 549         buf = (uchar *) cvAlloc( sizeof( uchar ) * objSize.width );
 550         if( buf == NULL )
 551             return CV_OUTOFMEM_ERR;
 552     }
 553
 554     if( ioFlags )
 555     {
 556         buffer = (void *) cvAlloc( ioBufSize );
 557         if( buffer == NULL )
 558         {
 559             if( buf )
 560                 cvFree( &buf );
 561             return CV_OUTOFMEM_ERR;
 562         }
 563     }
 564
 565     /* Calculation of averaged object */
 566     bf = avg;
 567     for( i = 0; i < avgSize.height; i++, bf += avgStep )
 568         for( j = 0; j < avgSize.width; j++ )
 569             bf[j] = 0.f;
 570
 571     for( i = 0; i < nObjects; i++ )
 572     {
 573         int k, l;
 574         uchar *bu = (ioFlags & CV_EIGOBJ_INPUT_CALLBACK) ? buf : ((uchar **) input)[i];
 575
 576         if( ioFlags & CV_EIGOBJ_INPUT_CALLBACK )
 577         {
 578             CvCallback read_callback = ((CvInput *) & input)->callback;
 579
 580             r = (CvStatus)read_callback( i, (void *) buf, userData );
 581             if( r )
 582             {
 583                 if( buffer )
 584                     cvFree( &buffer );
 585                 if( buf )
 586                     cvFree( &buf );
 587                 return r;
 588             }
 589         }
 590
 591         bf = avg;
 592         for( k = 0; k < avgSize.height; k++, bf += avgStep, bu += objStep1 )
 593             for( l = 0; l < avgSize.width; l++ )
 594                 bf[l] += bu[l];
 595     }
 596
 597     bf = avg;
 598     for( i = 0; i < avgSize.height; i++, bf += avgStep )
 599         for( j = 0; j < avgSize.width; j++ )
 600             bf[j] *= m;
 601
 602     /* Calculation of covariance matrix */
 603     c = (float *) cvAlloc( sizeof( float ) * nObjects * nObjects );
 604
 605     if( c == NULL )
 606     {
 607         if( buffer )
 608             cvFree( &buffer );
 609         if( buf )
 610             cvFree( &buf );
 611         return CV_OUTOFMEM_ERR;
 612     }
 613
 614     r = icvCalcCovarMatrixEx_8u32fR( nObjects, input, objStep1, ioFlags, ioBufSize,
 615                                       (uchar *) buffer, userData, avg, 4 * avgStep, size, c );
 616     if( r )
 617     {
 618         cvFree( &c );
 619         if( buffer )
 620             cvFree( &buffer );
 621         if( buf )
 622             cvFree( &buf );
 623         return r;
 624     }
 625
 626     /* Calculation of eigenvalues & eigenvectors */
 627     ev = (float *) cvAlloc( sizeof( float ) * nObjects * nObjects );
 628
 629     if( ev == NULL )
 630     {
 631         cvFree( &c );
 632         if( buffer )
 633             cvFree( &buffer );
 634         if( buf )
 635             cvFree( &buf );
 636         return CV_OUTOFMEM_ERR;
 637     }
 638
 639     if( eigVals == NULL )
 640     {
 641         eigVals = (float *) cvAlloc( sizeof( float ) * nObjects );
 642
 643         if( eigVals == NULL )
 644         {
 645             cvFree( &c );
 646             cvFree( &ev );
 647             if( buffer )
 648                 cvFree( &buffer );
 649             if( buf )
 650                 cvFree( &buf );
 651             return CV_OUTOFMEM_ERR;
 652         }
 653         iev = 1;
 654     }
 655
 656     r = icvJacobiEigens_32f( c, ev, eigVals, nObjects, 0.0f );
 657     cvFree( &c );
 658     if( r )
 659     {
 660         cvFree( &ev );
 661         if( buffer )
 662             cvFree( &buffer );
 663         if( buf )
 664             cvFree( &buf );
 665         if( iev )
 666             cvFree( &eigVals );
 667         return r;
 668     }
 669
 670     /* Eigen objects number determination */
 671     if( calcLimit->type != CV_TERMCRIT_NUMBER )
 672     {
 673         for( i = 0; i < m1; i++ )
 674             if( fabs( eigVals[i] / eigVals[0] ) < calcLimit->epsilon )
 675                 break;
 676         m1 = calcLimit->max_iter = i;
 677     }
 678     else
 679         m1 = calcLimit->max_iter;
 680     calcLimit->epsilon = (float) fabs( eigVals[m1 - 1] / eigVals[0] );
 681
 682     for( i = 0; i < m1; i++ )
 683         eigVals[i] = (float) (1.0 / sqrt( (double)eigVals[i] ));
 684
 685     /* ----------------- Calculation of eigenobjects ----------------------- */
 686     if( ioFlags & CV_EIGOBJ_OUTPUT_CALLBACK )
 687     {
 688         int nio, ngr, igr;
 689
 690         nio = ioBufSize / (4 * eigSize.width);  /* number of eigen objects in buffer */
 691         ngr = m1 / nio;         /* number of io groups */
 692         if( nObjects % nio )
 693             ngr += 1;
 694
 695         for( igr = 0; igr < ngr; igr++ )
 696         {
 697             int io, ie, imin = igr * nio, imax = imin + nio;
 698
 699             if( imax > m1 )
 700                 imax = m1;
 701
 702             for(int k = 0; k < eigSize.width * (imax - imin); k++ )
 703                 ((float *) buffer)[k] = 0.f;
 704
 705             for( io = 0; io < nObjects; io++ )
 706             {
 707                 uchar *bu = ioFlags & CV_EIGOBJ_INPUT_CALLBACK ? buf : ((uchar **) input)[io];
 708
 709                 if( ioFlags & CV_EIGOBJ_INPUT_CALLBACK )
 710                 {
 711                     CvCallback read_callback = ((CvInput *) & input)->callback;
 712
 713                     r = (CvStatus)read_callback( io, (void *) buf, userData );
 714                     if( r )
 715                     {
 716                         cvFree( &ev );
 717                         if( iev )
 718                             cvFree( &eigVals );
 719                         if( buffer )
 720                             cvFree( &buffer );
 721                         if( buf )
 722                             cvFree( &buf );
 723                         return r;
 724                     }
 725                 }
 726
 727                 for( ie = imin; ie < imax; ie++ )
 728                 {
 729                     int k, l;
 730                     uchar *bv = bu;
 731                     float e = ev[ie * nObjects + io] * eigVals[ie];
 732                     float *be = ((float *) buffer) + ((ie - imin) * eigStep);
 733
 734                     bf = avg;
 735                     for( k = 0; k < size.height; k++, bv += objStep1,
 736                          bf += avgStep, be += eigStep1 )
 737                     {
 738                         for( l = 0; l < size.width - 3; l += 4 )
 739                         {
 740                             float f = bf[l];
 741                             uchar v = bv[l];
 742
 743                             be[l] += e * (v - f);
 744                             f = bf[l + 1];
 745                             v = bv[l + 1];
 746                             be[l + 1] += e * (v - f);
 747                             f = bf[l + 2];
 748                             v = bv[l + 2];
 749                             be[l + 2] += e * (v - f);
 750                             f = bf[l + 3];
 751                             v = bv[l + 3];
 752                             be[l + 3] += e * (v - f);
 753                         }
 754                         for( ; l < size.width; l++ )
 755                             be[l] += e * (bv[l] - bf[l]);
 756                     }
 757                 }
 758             }                   /* io */
 759
 760             for( ie = imin; ie < imax; ie++ )   /* calculated eigen objects writting */
 761             {
 762                 CvCallback write_callback = ((CvInput *) & output)->callback;
 763                 float *be = ((float *) buffer) + ((ie - imin) * eigStep);
 764
 765                 r = (CvStatus)write_callback( ie, (void *) be, userData );
 766                 if( r )
 767                 {
 768                     cvFree( &ev );
 769                     if( iev )
 770                         cvFree( &eigVals );
 771                     if( buffer )
 772                         cvFree( &buffer );
 773                     if( buf )
 774                         cvFree( &buf );
 775                     return r;
 776                 }
 777             }
 778         }                       /* igr */
 779     }
 780
 781     else
 782     {
 783         int k, p, l;
 784
 785         for( i = 0; i < m1; i++ )       /* e.o. annulation */
 786         {
 787             float *be = ((float **) output)[i];
 788
 789             for( p = 0; p < eigSize.height; p++, be += eigStep )
 790                 for( l = 0; l < eigSize.width; l++ )
 791                     be[l] = 0.0f;
 792         }
 793
 794         for( k = 0; k < nObjects; k++ )
 795         {
 796             uchar *bv = (ioFlags & CV_EIGOBJ_INPUT_CALLBACK) ? buf : ((uchar **) input)[k];
 797
 798             if( ioFlags & CV_EIGOBJ_INPUT_CALLBACK )
 799             {
 800                 CvCallback read_callback = ((CvInput *) & input)->callback;
 801
 802                 r = (CvStatus)read_callback( k, (void *) buf, userData );
 803                 if( r )
 804                 {
 805                     cvFree( &ev );
 806                     if( iev )
 807                         cvFree( &eigVals );
 808                     if( buffer )
 809                         cvFree( &buffer );
 810                     if( buf )
 811                         cvFree( &buf );
 812                     return r;
 813                 }
 814             }
 815
 816             for( i = 0; i < m1; i++ )
 817             {
 818                 float v = eigVals[i] * ev[i * nObjects + k];
 819                 float *be = ((float **) output)[i];
 820                 uchar *bu = bv;
 821
 822                 bf = avg;
 823
 824                 for( p = 0; p < size.height; p++, bu += objStep1,
 825                      bf += avgStep, be += eigStep1 )
 826                 {
 827                     for( l = 0; l < size.width - 3; l += 4 )
 828                     {
 829                         float f = bf[l];
 830                         uchar u = bu[l];
 831
 832                         be[l] += v * (u - f);
 833                         f = bf[l + 1];
 834                         u = bu[l + 1];
 835                         be[l + 1] += v * (u - f);
 836                         f = bf[l + 2];
 837                         u = bu[l + 2];
 838                         be[l + 2] += v * (u - f);
 839                         f = bf[l + 3];
 840                         u = bu[l + 3];
 841                         be[l + 3] += v * (u - f);
 842                     }
 843                     for( ; l < size.width; l++ )
 844                         be[l] += v * (bu[l] - bf[l]);
 845                 }
 846             }                   /* i */
 847         }                       /* k */
 848     }                           /* else */
 849
 850     cvFree( &ev );
 851     if( iev )
 852         cvFree( &eigVals );
 853     else
 854         for( i = 0; i < m1; i++ )
 855             eigVals[i] = 1.f / (eigVals[i] * eigVals[i]);
 856     if( buffer )
 857         cvFree( &buffer );
 858     if( buf )
 859         cvFree( &buf );
 860     return CV_NO_ERR;
 861 }
 862
 863 /* --- End of icvCalcEigenObjects_8u32fR --- */
 864
 865 /*F///////////////////////////////////////////////////////////////////////////////////////
 866 //    Name: icvCalcDecompCoeff_8u32fR
 867 //    Purpose: The function calculates one decomposition coefficient of input object
 868 //             using previously calculated eigen object and the mean (averaged) object
 869 //    Context:
 870 //    Parameters:  obj     - input object
 871 //                 objStep - its step (in bytes)
 872 //                 eigObj  - pointer to eigen object
 873 //                 eigStep - its step (in bytes)
 874 //                 avg     - pointer to averaged object
 875 //                 avgStep - its step (in bytes)
 876 //                 size    - ROI size of each source object
 877 //
 878 //    Returns: decomposition coefficient value or large negative value (if error)
 879 //
 880 //    Notes:
 881 //F*/
 882 static float CV_STDCALL
 883 icvCalcDecompCoeff_8u32fR( uchar* obj, int objStep,
 884                            float *eigObj, int eigStep,
 885                            float *avg, int avgStep, CvSize size )
 886 {
 887     int i, k;
 888     float w = 0.0f;
 889
 890     if( size.width > objStep || 4 * size.width > eigStep
 891         || 4 * size.width > avgStep || size.height < 1 )
 892         return -1.0e30f;
 893     if( obj == NULL || eigObj == NULL || avg == NULL )
 894         return -1.0e30f;
 895
 896     eigStep /= 4;
 897     avgStep /= 4;
 898
 899     if( size.width == objStep && size.width == eigStep && size.width == avgStep )
 900     {
 901         size.width *= size.height;
 902         size.height = 1;
 903         objStep = eigStep = avgStep = size.width;
 904     }
 905
 906     for( i = 0; i < size.height; i++, obj += objStep, eigObj += eigStep, avg += avgStep )
 907     {
 908         for( k = 0; k < size.width - 4; k += 4 )
 909         {
 910             float o = (float) obj[k];
 911             float e = eigObj[k];
 912             float a = avg[k];
 913
 914             w += e * (o - a);
 915             o = (float) obj[k + 1];
 916             e = eigObj[k + 1];
 917             a = avg[k + 1];
 918             w += e * (o - a);
 919             o = (float) obj[k + 2];
 920             e = eigObj[k + 2];
 921             a = avg[k + 2];
 922             w += e * (o - a);
 923             o = (float) obj[k + 3];
 924             e = eigObj[k + 3];
 925             a = avg[k + 3];
 926             w += e * (o - a);
 927         }
 928         for( ; k < size.width; k++ )
 929             w += eigObj[k] * ((float) obj[k] - avg[k]);
 930     }
 931
 932     return w;
 933 }
 934
 935 /*F///////////////////////////////////////////////////////////////////////////////////////
 936 //    Names: icvEigenDecomposite_8u32fR
 937 //    Purpose: The function calculates all decomposition coefficients for input object
 938 //             using previously calculated eigen objects basis and the mean (averaged)
 939 //             object
 940 //    Context:
 941 //    Parameters:  obj         - input object
 942 //                 objStep     - its step (in bytes)
 943 //                 nEigObjs    - number of eigen objects
 944 //                 eigInput    - pointer either to array of pointers to eigen objects
 945 //                               or to read callback function (depending on ioFlags)
 946 //                 eigStep     - eigen objects step (in bytes)
 947 //                 ioFlags     - input/output flags
 948 //                 iserData    - pointer to the structure which contains all necessary
 949 //                               data for the callback function
 950 //                 avg         - pointer to averaged object
 951 //                 avgStep     - its step (in bytes)
 952 //                 size        - ROI size of each source object
 953 //                 coeffs      - calculated coefficients (output data)
 954 //
 955 //    Returns: icv status
 956 //
 957 //    Notes:   see notes for icvCalcEigenObjects_8u32fR function
 958 //F*/
 959 static CvStatus CV_STDCALL
 960 icvEigenDecomposite_8u32fR( uchar * obj, int objStep, int nEigObjs,
 961                             void *eigInput, int eigStep, int ioFlags,
 962                             void *userData, float *avg, int avgStep,
 963                             CvSize size, float *coeffs )
 964 {
 965     int i;
 966
 967     if( nEigObjs < 2 )
 968         return CV_BADFACTOR_ERR;
 969     if( ioFlags < 0 || ioFlags > 1 )
 970         return CV_BADFACTOR_ERR;
 971     if( size.width > objStep || 4 * size.width > eigStep ||
 972         4 * size.width > avgStep || size.height < 1 )
 973         return CV_BADSIZE_ERR;
 974     if( obj == NULL || eigInput == NULL || coeffs == NULL || avg == NULL )
 975         return CV_NULLPTR_ERR;
 976     if( !ioFlags )
 977         for( i = 0; i < nEigObjs; i++ )
 978             if( ((uchar **) eigInput)[i] == NULL )
 979                 return CV_NULLPTR_ERR;
 980
 981     if( ioFlags )               /* callback */
 982
 983     {
 984         float *buffer;
 985         CvCallback read_callback = ((CvInput *) & eigInput)->callback;
 986
 987         eigStep = 4 * size.width;
 988
 989         /* memory allocation */
 990         buffer = (float *) cvAlloc( sizeof( float ) * size.width * size.height );
 991
 992         if( buffer == NULL )
 993             return CV_OUTOFMEM_ERR;
 994
 995         for( i = 0; i < nEigObjs; i++ )
 996         {
 997             float w;
 998             CvStatus r = (CvStatus)read_callback( i, (void *) buffer, userData );
 999
1000             if( r )
1001             {
1002                 cvFree( &buffer );
1003                 return r;
1004             }
1005             w = icvCalcDecompCoeff_8u32fR( obj, objStep, buffer,
1006                                             eigStep, avg, avgStep, size );
1007             if( w < -1.0e29f )
1008             {
1009                 cvFree( &buffer );
1010                 return CV_NOTDEFINED_ERR;
1011             }
1012             coeffs[i] = w;
1013         }
1014         cvFree( &buffer );
1015     }
1016
1017     else
1018         /* no callback */
1019         for( i = 0; i < nEigObjs; i++ )
1020         {
1021             float w = icvCalcDecompCoeff_8u32fR( obj, objStep, ((float **) eigInput)[i],
1022                                                   eigStep, avg, avgStep, size );
1023
1024             if( w < -1.0e29f )
1025                 return CV_NOTDEFINED_ERR;
1026             coeffs[i] = w;
1027         }
1028
1029     return CV_NO_ERR;
1030 }
1031
1032
1033 /*F///////////////////////////////////////////////////////////////////////////////////////
1034 //    Names: icvEigenProjection_8u32fR
1035 //    Purpose: The function calculates object projection to the eigen sub-space (restores
1036 //             an object) using previously calculated eigen objects basis, mean (averaged)
1037 //             object and decomposition coefficients of the restored object
1038 //    Context:
1039 //    Parameters:  nEigObjs - Number of eigen objects
1040 //                 eigens   - Array of pointers to eigen objects
1041 //                 eigStep  - Eigen objects step (in bytes)
1042 //                 coeffs   - Previously calculated decomposition coefficients
1043 //                 avg      - Pointer to averaged object
1044 //                 avgStep  - Its step (in bytes)
1045 //                 rest     - Pointer to restored object
1046 //                 restStep - Its step (in bytes)
1047 //                 size     - ROI size of each object
1048 //
1049 //    Returns: CV status
1050 //
1051 //    Notes:
1052 //F*/
1053 static CvStatus CV_STDCALL
1054 icvEigenProjection_8u32fR( int nEigObjs, void *eigInput, int eigStep,
1055                            int ioFlags, void *userData, float *coeffs,
1056                            float *avg, int avgStep, uchar * rest,
1057                            int restStep, CvSize size )
1058 {
1059     int i, j, k;
1060     float *buf;
1061     float *buffer = NULL;
1062     float *b;
1063     CvCallback read_callback = ((CvInput *) & eigInput)->callback;
1064
1065     if( size.width > avgStep || 4 * size.width > eigStep || size.height < 1 )
1066         return CV_BADSIZE_ERR;
1067     if( rest == NULL || eigInput == NULL || avg == NULL || coeffs == NULL )
1068         return CV_NULLPTR_ERR;
1069     if( ioFlags < 0 || ioFlags > 1 )
1070         return CV_BADFACTOR_ERR;
1071     if( !ioFlags )
1072         for( i = 0; i < nEigObjs; i++ )
1073             if( ((uchar **) eigInput)[i] == NULL )
1074                 return CV_NULLPTR_ERR;
1075     eigStep /= 4;
1076     avgStep /= 4;
1077
1078     if( size.width == restStep && size.width == eigStep && size.width == avgStep )
1079     {
1080         size.width *= size.height;
1081         size.height = 1;
1082         restStep = eigStep = avgStep = size.width;
1083     }
1084
1085     buf = (float *) cvAlloc( sizeof( float ) * size.width * size.height );
1086
1087     if( buf == NULL )
1088         return CV_OUTOFMEM_ERR;
1089     b = buf;
1090     for( i = 0; i < size.height; i++, avg += avgStep, b += size.width )
1091         for( j = 0; j < size.width; j++ )
1092             b[j] = avg[j];
1093
1094     if( ioFlags )
1095     {
1096         buffer = (float *) cvAlloc( sizeof( float ) * size.width * size.height );
1097
1098         if( buffer == NULL )
1099         {
1100             cvFree( &buf );
1101             return CV_OUTOFMEM_ERR;
1102         }
1103         eigStep = size.width;
1104     }
1105
1106     for( k = 0; k < nEigObjs; k++ )
1107     {
1108         float *e = ioFlags ? buffer : ((float **) eigInput)[k];
1109         float c = coeffs[k];
1110
1111         if( ioFlags )           /* read eigen object */
1112         {
1113             CvStatus r = (CvStatus)read_callback( k, (void *) buffer, userData );
1114
1115             if( r )
1116             {
1117                 cvFree( &buf );
1118                 cvFree( &buffer );
1119                 return r;
1120             }
1121         }
1122
1123         b = buf;
1124         for( i = 0; i < size.height; i++, e += eigStep, b += size.width )
1125         {
1126             for( j = 0; j < size.width - 3; j += 4 )
1127             {
1128                 float b0 = c * e[j];
1129                 float b1 = c * e[j + 1];
1130                 float b2 = c * e[j + 2];
1131                 float b3 = c * e[j + 3];
1132
1133                 b[j] += b0;
1134                 b[j + 1] += b1;
1135                 b[j + 2] += b2;
1136                 b[j + 3] += b3;
1137             }
1138             for( ; j < size.width; j++ )
1139                 b[j] += c * e[j];
1140         }
1141     }
1142
1143     b = buf;
1144     for( i = 0; i < size.height; i++, avg += avgStep, b += size.width, rest += restStep )
1145         for( j = 0; j < size.width; j++ )
1146         {
1147             int w = cvRound( b[j] );
1148
1149             w = !(w & ~255) ? w : w < 0 ? 0 : 255;
1150             rest[j] = (uchar) w;
1151         }
1152
1153     cvFree( &buf );
1154     if( ioFlags )
1155         cvFree( &buffer );
1156     return CV_NO_ERR;
1157 }
1158
1159 /*F///////////////////////////////////////////////////////////////////////////////////////
1160 //    Name: cvCalcCovarMatrixEx
1161 //    Purpose: The function calculates a covariance matrix for a group of input objects
1162 //             (images, vectors, etc.).
1163 //    Context:
1164 //    Parameters:  nObjects    - number of source objects
1165 //                 input       - pointer either to array of input objects
1166 //                               or to read callback function (depending on ioFlags)
1167 //                 ioFlags     - input/output flags (see Notes to
1168 //                               cvCalcEigenObjects function)
1169 //                 ioBufSize   - input/output buffer size
1170 //                 userData    - pointer to the structure which contains all necessary
1171 //                               data for the callback functions
1172 //                 avg         - averaged object
1173 //                 covarMatrix - covariance matrix (output parameter; must be allocated
1174 //                               before call)
1175 //
1176 //    Notes:  See Notes to cvCalcEigenObjects function
1177 //F*/
1178
1179 CV_IMPL void
1180 cvCalcCovarMatrixEx( int  nObjects, void*  input, int  ioFlags,
1181                      int  ioBufSize, uchar*  buffer, void*  userData,
1182                      IplImage* avg, float*  covarMatrix )
1183 {
1184     float *avg_data;
1185     int avg_step = 0;
1186     CvSize avg_size;
1187     int i;
1188
1189     CV_FUNCNAME( "cvCalcCovarMatrixEx" );
1190
1191     __BEGIN__;
1192
1193     cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &avg_size );
1194     if( avg->depth != IPL_DEPTH_32F )
1195         CV_ERROR( CV_BadDepth, "Unsupported format" );
1196     if( avg->nChannels != 1 )
1197         CV_ERROR( CV_BadNumChannels, "Unsupported format" );
1198
1199     if( ioFlags == CV_EIGOBJ_NO_CALLBACK )
1200     {
1201         IplImage **images = (IplImage **) (((CvInput *) & input)->data);
1202         uchar **objects = (uchar **) cvAlloc( sizeof( uchar * ) * nObjects );
1203         int img_step = 0, old_step = 0;
1204         CvSize img_size = avg_size, old_size = avg_size;
1205
1206         if( objects == NULL )
1207             CV_ERROR( CV_StsBadArg, "Insufficient memory" );
1208
1209         for( i = 0; i < nObjects; i++ )
1210         {
1211             IplImage *img = images[i];
1212             uchar *img_data;
1213
1214             cvGetImageRawData( img, &img_data, &img_step, &img_size );
1215             if( img->depth != IPL_DEPTH_8U )
1216                 CV_ERROR( CV_BadDepth, "Unsupported format" );
1217             if( img_size != avg_size || img_size != old_size )
1218                 CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
1219             if( img->nChannels != 1 )
1220                 CV_ERROR( CV_BadNumChannels, "Unsupported format" );
1221             if( i > 0 && img_step != old_step )
1222                 CV_ERROR( CV_StsBadArg, "Different steps of objects" );
1223
1224             old_step = img_step;
1225             old_size = img_size;
1226             objects[i] = img_data;
1227         }
1228
1229         CV_CALL( icvCalcCovarMatrixEx_8u32fR( nObjects,
1230                                               (void*) objects,
1231                                               img_step,
1232                                               CV_EIGOBJ_NO_CALLBACK,
1233                                               0,
1234                                               NULL,
1235                                               NULL,
1236                                               avg_data,
1237                                               avg_step,
1238                                               avg_size,
1239                                               covarMatrix ));
1240         cvFree( &objects );
1241     }
1242
1243     else
1244
1245     {
1246         CV_CALL( icvCalcCovarMatrixEx_8u32fR( nObjects,
1247                                               input,
1248                                               avg_step / 4,
1249                                               ioFlags,
1250                                               ioBufSize,
1251                                               buffer,
1252                                               userData,
1253                                               avg_data,
1254                                               avg_step,
1255                                               avg_size,
1256                                               covarMatrix ));
1257     }
1258
1259     __END__;
1260 }
1261
1262 /*F///////////////////////////////////////////////////////////////////////////////////////
1263 //    Name: cvCalcEigenObjects
1264 //    Purpose: The function calculates an orthonormal eigen basis and a mean (averaged)
1265 //             object for a group of input objects (images, vectors, etc.).
1266 //    Context:
1267 //    Parameters: nObjects  - number of source objects
1268 //                input     - pointer either to array of input objects
1269 //                            or to read callback function (depending on ioFlags)
1270 //                output    - pointer either to output eigen objects
1271 //                            or to write callback function (depending on ioFlags)
1272 //                ioFlags   - input/output flags (see Notes)
1273 //                ioBufSize - input/output buffer size
1274 //                userData  - pointer to the structure which contains all necessary
1275 //                            data for the callback functions
1276 //                calcLimit - determines the calculation finish conditions
1277 //                avg       - averaged object (has the same size as ROI)
1278 //                eigVals   - pointer to corresponding eigen values (array of <nObjects>
1279 //                            elements in descending order)
1280 //
1281 //    Notes: 1. input/output data (that is, input objects and eigen ones) may either
1282 //              be allocated in the RAM or be read from/written to the HDD (or any
1283 //              other device) by read/write callback functions. It depends on the
1284 //              value of ioFlags paramater, which may be the following:
1285 //                  CV_EIGOBJ_NO_CALLBACK, or 0;
1286 //                  CV_EIGOBJ_INPUT_CALLBACK;
1287 //                  CV_EIGOBJ_OUTPUT_CALLBACK;
1288 //                  CV_EIGOBJ_BOTH_CALLBACK, or
1289 //                            CV_EIGOBJ_INPUT_CALLBACK | CV_EIGOBJ_OUTPUT_CALLBACK.
1290 //              The callback functions as well as the user data structure must be
1291 //              developed by the user.
1292 //
1293 //           2. If ioBufSize = 0, or it's too large, the function dermines buffer size
1294 //              itself.
1295 //
1296 //           3. Depending on calcLimit parameter, calculations are finished either if
1297 //              eigenfaces number comes up to certain value or the relation of the
1298 //              current eigenvalue and the largest one comes down to certain value
1299 //              (or any of the above conditions takes place). The calcLimit->type value
1300 //              must be CV_TERMCRIT_NUMB, CV_TERMCRIT_EPS or
1301 //              CV_TERMCRIT_NUMB | CV_TERMCRIT_EPS. The function returns the real
1302 //              values calcLimit->max_iter and calcLimit->epsilon.
1303 //
1304 //           4. eigVals may be equal to NULL (if you don't need eigen values in further).
1305 //
1306 //F*/
1307 CV_IMPL void
1308 cvCalcEigenObjects( int       nObjects,
1309                     void*     input,
1310                     void*     output,
1311                     int       ioFlags,
1312                     int       ioBufSize,
1313                     void*     userData,
1314                     CvTermCriteria* calcLimit,
1315                     IplImage* avg,
1316                     float*    eigVals )
1317 {
1318     float *avg_data;
1319     int avg_step = 0;
1320     CvSize avg_size;
1321     int i;
1322     int nEigens = nObjects - 1;
1323
1324     CV_FUNCNAME( "cvCalcEigenObjects" );
1325
1326     __BEGIN__;
1327
1328     cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &avg_size );
1329     if( avg->depth != IPL_DEPTH_32F )
1330         CV_ERROR( CV_BadDepth, "Unsupported format" );
1331     if( avg->nChannels != 1 )
1332         CV_ERROR( CV_BadNumChannels, "Unsupported format" );
1333
1334     if( nEigens > calcLimit->max_iter && calcLimit->type != CV_TERMCRIT_EPS )
1335         nEigens = calcLimit->max_iter;
1336
1337     switch (ioFlags)
1338     {
1339     case CV_EIGOBJ_NO_CALLBACK:
1340         {
1341             IplImage **objects = (IplImage **) (((CvInput *) & input)->data);
1342             IplImage **eigens = (IplImage **) (((CvInput *) & output)->data);
1343             uchar **objs = (uchar **) cvAlloc( sizeof( uchar * ) * nObjects );
1344             float **eigs = (float **) cvAlloc( sizeof( float * ) * nEigens );
1345             int obj_step = 0, old_step = 0;
1346             int eig_step = 0, oldeig_step = 0;
1347             CvSize obj_size = avg_size, old_size = avg_size,
1348
1349                 eig_size = avg_size, oldeig_size = avg_size;
1350
1351             if( objects == NULL || eigens == NULL )
1352                 CV_ERROR( CV_StsBadArg, "Insufficient memory" );
1353
1354             for( i = 0; i < nObjects; i++ )
1355             {
1356                 IplImage *img = objects[i];
1357                 uchar *obj_data;
1358
1359                 cvGetImageRawData( img, &obj_data, &obj_step, &obj_size );
1360                 if( img->depth != IPL_DEPTH_8U )
1361                     CV_ERROR( CV_BadDepth, "Unsupported format" );
1362                 if( obj_size != avg_size || obj_size != old_size )
1363                     CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
1364                 if( img->nChannels != 1 )
1365                     CV_ERROR( CV_BadNumChannels, "Unsupported format" );
1366                 if( i > 0 && obj_step != old_step )
1367                     CV_ERROR( CV_StsBadArg, "Different steps of objects" );
1368
1369                 old_step = obj_step;
1370                 old_size = obj_size;
1371                 objs[i] = obj_data;
1372             }
1373             for( i = 0; i < nEigens; i++ )
1374             {
1375                 IplImage *eig = eigens[i];
1376                 float *eig_data;
1377
1378                 cvGetImageRawData( eig, (uchar **) & eig_data, &eig_step, &eig_size );
1379                 if( eig->depth != IPL_DEPTH_32F )
1380                     CV_ERROR( CV_BadDepth, "Unsupported format" );
1381                 if( eig_size != avg_size || eig_size != oldeig_size )
1382                     CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
1383                 if( eig->nChannels != 1 )
1384                     CV_ERROR( CV_BadNumChannels, "Unsupported format" );
1385                 if( i > 0 && eig_step != oldeig_step )
1386                     CV_ERROR( CV_StsBadArg, "Different steps of objects" );
1387
1388                 oldeig_step = eig_step;
1389                 oldeig_size = eig_size;
1390                 eigs[i] = eig_data;
1391             }
1392             CV_CALL( icvCalcEigenObjects_8u32fR( nObjects, (void*) objs, obj_step,
1393                                                  (void*) eigs, eig_step, obj_size,
1394                                                  ioFlags, ioBufSize, userData,
1395                                                  calcLimit, avg_data, avg_step, eigVals ));
1396             cvFree( &objs );
1397             cvFree( &eigs );
1398             break;
1399         }
1400
1401     case CV_EIGOBJ_OUTPUT_CALLBACK:
1402         {
1403             IplImage **objects = (IplImage **) (((CvInput *) & input)->data);
1404             uchar **objs = (uchar **) cvAlloc( sizeof( uchar * ) * nObjects );
1405             int obj_step = 0, old_step = 0;
1406             CvSize obj_size = avg_size, old_size = avg_size;
1407
1408             if( objects == NULL )
1409                 CV_ERROR( CV_StsBadArg, "Insufficient memory" );
1410
1411             for( i = 0; i < nObjects; i++ )
1412             {
1413                 IplImage *img = objects[i];
1414                 uchar *obj_data;
1415
1416                 cvGetImageRawData( img, &obj_data, &obj_step, &obj_size );
1417                 if( img->depth != IPL_DEPTH_8U )
1418                     CV_ERROR( CV_BadDepth, "Unsupported format" );
1419                 if( obj_size != avg_size || obj_size != old_size )
1420                     CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
1421                 if( img->nChannels != 1 )
1422                     CV_ERROR( CV_BadNumChannels, "Unsupported format" );
1423                 if( i > 0 && obj_step != old_step )
1424                     CV_ERROR( CV_StsBadArg, "Different steps of objects" );
1425
1426                 old_step = obj_step;
1427                 old_size = obj_size;
1428                 objs[i] = obj_data;
1429             }
1430             CV_CALL( icvCalcEigenObjects_8u32fR( nObjects,
1431                                                  (void*) objs,
1432                                                  obj_step,
1433                                                  output,
1434                                                  avg_step,
1435                                                  obj_size,
1436                                                  ioFlags,
1437                                                  ioBufSize,
1438                                                  userData,
1439                                                  calcLimit,
1440                                                  avg_data,
1441                                                  avg_step,
1442                                                  eigVals   ));
1443             cvFree( &objs );
1444             break;
1445         }
1446
1447     case CV_EIGOBJ_INPUT_CALLBACK:
1448         {
1449             IplImage **eigens = (IplImage **) (((CvInput *) & output)->data);
1450             float **eigs = (float**) cvAlloc( sizeof( float* ) * nEigens );
1451             int eig_step = 0, oldeig_step = 0;
1452             CvSize eig_size = avg_size, oldeig_size = avg_size;
1453
1454             if( eigens == NULL )
1455                 CV_ERROR( CV_StsBadArg, "Insufficient memory" );
1456
1457             for( i = 0; i < nEigens; i++ )
1458             {
1459                 IplImage *eig = eigens[i];
1460                 float *eig_data;
1461
1462                 cvGetImageRawData( eig, (uchar **) & eig_data, &eig_step, &eig_size );
1463                 if( eig->depth != IPL_DEPTH_32F )
1464                     CV_ERROR( CV_BadDepth, "Unsupported format" );
1465                 if( eig_size != avg_size || eig_size != oldeig_size )
1466                     CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
1467                 if( eig->nChannels != 1 )
1468                     CV_ERROR( CV_BadNumChannels, "Unsupported format" );
1469                 if( i > 0 && eig_step != oldeig_step )
1470                     CV_ERROR( CV_StsBadArg, "Different steps of objects" );
1471
1472                 oldeig_step = eig_step;
1473                 oldeig_size = eig_size;
1474                 eigs[i] = eig_data;
1475             }
1476             CV_CALL( icvCalcEigenObjects_8u32fR( nObjects,
1477                                                  input,
1478                                                  avg_step / 4,
1479                                                  (void*) eigs,
1480                                                  eig_step,
1481                                                  eig_size,
1482                                                  ioFlags,
1483                                                  ioBufSize,
1484                                                  userData,
1485                                                  calcLimit,
1486                                                  avg_data,
1487                                                  avg_step,
1488                                                  eigVals   ));
1489             cvFree( &eigs );
1490             break;
1491         }
1492     case CV_EIGOBJ_INPUT_CALLBACK | CV_EIGOBJ_OUTPUT_CALLBACK:
1493
1494         CV_CALL( icvCalcEigenObjects_8u32fR( nObjects,
1495                                              input,
1496                                              avg_step / 4,
1497                                              output,
1498                                              avg_step,
1499                                              avg_size,
1500                                              ioFlags,
1501                                              ioBufSize,
1502                                              userData,
1503                                              calcLimit,
1504                                              avg_data,
1505                                              avg_step,
1506                                              eigVals   ));
1507         break;
1508
1509     default:
1510         CV_ERROR( CV_StsBadArg, "Unsupported i/o flag" );
1511     }
1512
1513     __END__;
1514 }
1515
1516 /*--------------------------------------------------------------------------------------*/
1517 /*F///////////////////////////////////////////////////////////////////////////////////////
1518 //    Name: cvCalcDecompCoeff
1519 //    Purpose: The function calculates one decomposition coefficient of input object
1520 //             using previously calculated eigen object and the mean (averaged) object
1521 //    Context:
1522 //    Parameters:  obj     - input object
1523 //                 eigObj  - eigen object
1524 //                 avg     - averaged object
1525 //
1526 //    Returns: decomposition coefficient value or large negative value (if error)
1527 //
1528 //    Notes:
1529 //F*/
1530
1531 CV_IMPL double
1532 cvCalcDecompCoeff( IplImage * obj, IplImage * eigObj, IplImage * avg )
1533 {
1534     double coeff = DBL_MAX;
1535
1536     uchar *obj_data;
1537     float *eig_data;
1538     float *avg_data;
1539     int obj_step = 0, eig_step = 0, avg_step = 0;
1540     CvSize obj_size, eig_size, avg_size;
1541
1542     CV_FUNCNAME( "cvCalcDecompCoeff" );
1543
1544     __BEGIN__;
1545
1546     cvGetImageRawData( obj, &obj_data, &obj_step, &obj_size );
1547     if( obj->depth != IPL_DEPTH_8U )
1548         CV_ERROR( CV_BadDepth, "Unsupported format" );
1549     if( obj->nChannels != 1 )
1550         CV_ERROR( CV_BadNumChannels, "Unsupported format" );
1551
1552     cvGetImageRawData( eigObj, (uchar **) & eig_data, &eig_step, &eig_size );
1553     if( eigObj->depth != IPL_DEPTH_32F )
1554         CV_ERROR( CV_BadDepth, "Unsupported format" );
1555     if( eigObj->nChannels != 1 )
1556         CV_ERROR( CV_BadNumChannels, "Unsupported format" );
1557
1558     cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &avg_size );
1559     if( avg->depth != IPL_DEPTH_32F )
1560         CV_ERROR( CV_BadDepth, "Unsupported format" );
1561     if( avg->nChannels != 1 )
1562         CV_ERROR( CV_BadNumChannels, "Unsupported format" );
1563
1564     if( obj_size != eig_size || obj_size != avg_size )
1565         CV_ERROR( CV_StsBadArg, "different sizes of images" );
1566
1567     coeff = icvCalcDecompCoeff_8u32fR( obj_data, obj_step,
1568                                        eig_data, eig_step,
1569                                        avg_data, avg_step, obj_size );
1570
1571     __END__;
1572
1573     return coeff;
1574 }
1575
1576 /*--------------------------------------------------------------------------------------*/
1577 /*F///////////////////////////////////////////////////////////////////////////////////////
1578 //    Names: cvEigenDecomposite
1579 //    Purpose: The function calculates all decomposition coefficients for input object
1580 //             using previously calculated eigen objects basis and the mean (averaged)
1581 //             object
1582 //
1583 //    Parameters:  obj         - input object
1584 //                 nEigObjs    - number of eigen objects
1585 //                 eigInput    - pointer either to array of pointers to eigen objects
1586 //                               or to read callback function (depending on ioFlags)
1587 //                 ioFlags     - input/output flags
1588 //                 userData    - pointer to the structure which contains all necessary
1589 //                               data for the callback function
1590 //                 avg         - averaged object
1591 //                 coeffs      - calculated coefficients (output data)
1592 //
1593 //    Notes:   see notes for cvCalcEigenObjects function
1594 //F*/
1595
1596 CV_IMPL void
1597 cvEigenDecomposite( IplImage* obj,
1598                     int       nEigObjs,
1599                     void*     eigInput,
1600                     int       ioFlags,
1601                     void*     userData,
1602                     IplImage* avg,
1603                     float*    coeffs )
1604 {
1605     float *avg_data;
1606     uchar *obj_data;
1607     int avg_step = 0, obj_step = 0;
1608     CvSize avg_size, obj_size;
1609     int i;
1610
1611     CV_FUNCNAME( "cvEigenDecomposite" );
1612
1613     __BEGIN__;
1614
1615     cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &avg_size );
1616     if( avg->depth != IPL_DEPTH_32F )
1617         CV_ERROR( CV_BadDepth, "Unsupported format" );
1618     if( avg->nChannels != 1 )
1619         CV_ERROR( CV_BadNumChannels, "Unsupported format" );
1620
1621     cvGetImageRawData( obj, &obj_data, &obj_step, &obj_size );
1622     if( obj->depth != IPL_DEPTH_8U )
1623         CV_ERROR( CV_BadDepth, "Unsupported format" );
1624     if( obj->nChannels != 1 )
1625         CV_ERROR( CV_BadNumChannels, "Unsupported format" );
1626
1627     if( obj_size != avg_size )
1628         CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
1629
1630     if( ioFlags == CV_EIGOBJ_NO_CALLBACK )
1631     {
1632         IplImage **eigens = (IplImage **) (((CvInput *) & eigInput)->data);
1633         float **eigs = (float **) cvAlloc( sizeof( float * ) * nEigObjs );
1634         int eig_step = 0, old_step = 0;
1635         CvSize eig_size = avg_size, old_size = avg_size;
1636
1637         if( eigs == NULL )
1638             CV_ERROR( CV_StsBadArg, "Insufficient memory" );
1639
1640         for( i = 0; i < nEigObjs; i++ )
1641         {
1642             IplImage *eig = eigens[i];
1643             float *eig_data;
1644
1645             cvGetImageRawData( eig, (uchar **) & eig_data, &eig_step, &eig_size );
1646             if( eig->depth != IPL_DEPTH_32F )
1647                 CV_ERROR( CV_BadDepth, "Unsupported format" );
1648             if( eig_size != avg_size || eig_size != old_size )
1649                 CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
1650             if( eig->nChannels != 1 )
1651                 CV_ERROR( CV_BadNumChannels, "Unsupported format" );
1652             if( i > 0 && eig_step != old_step )
1653                 CV_ERROR( CV_StsBadArg, "Different steps of objects" );
1654
1655             old_step = eig_step;
1656             old_size = eig_size;
1657             eigs[i] = eig_data;
1658         }
1659
1660         CV_CALL( icvEigenDecomposite_8u32fR( obj_data,
1661                                              obj_step,
1662                                              nEigObjs,
1663                                              (void*) eigs,
1664                                              eig_step,
1665                                              ioFlags,
1666                                              userData,
1667                                              avg_data,
1668                                              avg_step,
1669                                              obj_size,
1670                                              coeffs   ));
1671         cvFree( &eigs );
1672     }
1673
1674     else
1675
1676     {
1677         CV_CALL( icvEigenDecomposite_8u32fR( obj_data,
1678                                              obj_step,
1679                                              nEigObjs,
1680                                              eigInput,
1681                                              avg_step,
1682                                              ioFlags,
1683                                              userData,
1684                                              avg_data,
1685                                              avg_step,
1686                                              obj_size,
1687                                              coeffs   ));
1688     }
1689
1690     __END__;
1691 }
1692
1693 /*--------------------------------------------------------------------------------------*/
1694 /*F///////////////////////////////////////////////////////////////////////////////////////
1695 //    Name: cvEigenProjection
1696 //    Purpose: The function calculates object projection to the eigen sub-space (restores
1697 //             an object) using previously calculated eigen objects basis, mean (averaged)
1698 //             object and decomposition coefficients of the restored object
1699 //    Context:
1700 //    Parameters:  nEigObjs    - number of eigen objects
1701 //                 eigInput    - pointer either to array of pointers to eigen objects
1702 //                               or to read callback function (depending on ioFlags)
1703 //                 ioFlags     - input/output flags
1704 //                 userData    - pointer to the structure which contains all necessary
1705 //                               data for the callback function
1706 //                 coeffs      - array of decomposition coefficients
1707 //                 avg         - averaged object
1708 //                 proj        - object projection (output data)
1709 //
1710 //    Notes:   see notes for cvCalcEigenObjects function
1711 //F*/
1712
1713 CV_IMPL void
1714 cvEigenProjection( void*     eigInput,
1715                    int       nEigObjs,
1716                    int       ioFlags,
1717                    void*     userData,
1718                    float*    coeffs,
1719                    IplImage* avg,
1720                    IplImage* proj )
1721 {
1722     float *avg_data;
1723     uchar *proj_data;
1724     int avg_step = 0, proj_step = 0;
1725     CvSize avg_size, proj_size;
1726     int i;
1727
1728     CV_FUNCNAME( "cvEigenProjection" );
1729
1730     __BEGIN__;
1731
1732     cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &avg_size );
1733     if( avg->depth != IPL_DEPTH_32F )
1734         CV_ERROR( CV_BadDepth, "Unsupported format" );
1735     if( avg->nChannels != 1 )
1736         CV_ERROR( CV_BadNumChannels, "Unsupported format" );
1737
1738     cvGetImageRawData( proj, &proj_data, &proj_step, &proj_size );
1739     if( proj->depth != IPL_DEPTH_8U )
1740         CV_ERROR( CV_BadDepth, "Unsupported format" );
1741     if( proj->nChannels != 1 )
1742         CV_ERROR( CV_BadNumChannels, "Unsupported format" );
1743
1744     if( proj_size != avg_size )
1745         CV_ERROR( CV_StsBadArg, "Different sizes of projects" );
1746
1747     if( ioFlags == CV_EIGOBJ_NO_CALLBACK )
1748     {
1749         IplImage **eigens = (IplImage**) (((CvInput *) & eigInput)->data);
1750         float **eigs = (float**) cvAlloc( sizeof( float * ) * nEigObjs );
1751         int eig_step = 0, old_step = 0;
1752         CvSize eig_size = avg_size, old_size = avg_size;
1753
1754         if( eigs == NULL )
1755             CV_ERROR( CV_StsBadArg, "Insufficient memory" );
1756
1757         for( i = 0; i < nEigObjs; i++ )
1758         {
1759             IplImage *eig = eigens[i];
1760             float *eig_data;
1761
1762             cvGetImageRawData( eig, (uchar **) & eig_data, &eig_step, &eig_size );
1763             if( eig->depth != IPL_DEPTH_32F )
1764                 CV_ERROR( CV_BadDepth, "Unsupported format" );
1765             if( eig_size != avg_size || eig_size != old_size )
1766                 CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
1767             if( eig->nChannels != 1 )
1768                 CV_ERROR( CV_BadNumChannels, "Unsupported format" );
1769             if( i > 0 && eig_step != old_step )
1770                 CV_ERROR( CV_StsBadArg, "Different steps of objects" );
1771
1772             old_step = eig_step;
1773             old_size = eig_size;
1774             eigs[i] = eig_data;
1775         }
1776
1777         CV_CALL( icvEigenProjection_8u32fR( nEigObjs,
1778                                             (void*) eigs,
1779                                             eig_step,
1780                                             ioFlags,
1781                                             userData,
1782                                             coeffs,
1783                                             avg_data,
1784                                             avg_step,
1785                                             proj_data,
1786                                             proj_step,
1787                                             avg_size   ));
1788         cvFree( &eigs );
1789     }
1790
1791     else
1792
1793     {
1794         CV_CALL( icvEigenProjection_8u32fR( nEigObjs,
1795                                             eigInput,
1796                                             avg_step,
1797                                             ioFlags,
1798                                             userData,
1799                                             coeffs,
1800                                             avg_data,
1801                                             avg_step,
1802                                             proj_data,
1803                                             proj_step,
1804                                             avg_size   ));
1805     }
1806
1807     __END__;
1808 }
1809
1810 /* End of file. */