Home | History | Annotate | Download | only in src
      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. */
   1812