Update to 2.0.0 tree from current Fremantle build
[opencv] / src / cvaux / cveigenobjects.cpp
diff --git a/src/cvaux/cveigenobjects.cpp b/src/cvaux/cveigenobjects.cpp
new file mode 100644 (file)
index 0000000..041aa98
--- /dev/null
@@ -0,0 +1,1811 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "_cvaux.h"
+
+CvStatus CV_STDCALL
+icvJacobiEigens_32f(float *A, float *V, float *E, int n, float eps)
+{
+    int i, j, k, ind;
+    float *AA = A, *VV = V;
+    double Amax, anorm = 0, ax;
+
+    if( A == NULL || V == NULL || E == NULL )
+        return CV_NULLPTR_ERR;
+    if( n <= 0 )
+        return CV_BADSIZE_ERR;
+    if( eps < 1.0e-7f )
+        eps = 1.0e-7f;
+
+    /*-------- Prepare --------*/
+    for( i = 0; i < n; i++, VV += n, AA += n )
+    {
+        for( j = 0; j < i; j++ )
+        {
+            double Am = AA[j];
+
+            anorm += Am * Am;
+        }
+        for( j = 0; j < n; j++ )
+            VV[j] = 0.f;
+        VV[i] = 1.f;
+    }
+
+    anorm = sqrt( anorm + anorm );
+    ax = anorm * eps / n;
+    Amax = anorm;
+
+    while( Amax > ax )
+    {
+        Amax /= n;
+        do                      /* while (ind) */
+        {
+            int p, q;
+            float *V1 = V, *A1 = A;
+
+            ind = 0;
+            for( p = 0; p < n - 1; p++, A1 += n, V1 += n )
+            {
+                float *A2 = A + n * (p + 1), *V2 = V + n * (p + 1);
+
+                for( q = p + 1; q < n; q++, A2 += n, V2 += n )
+                {
+                    double x, y, c, s, c2, s2, a;
+                    float *A3, Apq = A1[q], App, Aqq, Aip, Aiq, Vpi, Vqi;
+
+                    if( fabs( Apq ) < Amax )
+                        continue;
+
+                    ind = 1;
+
+                    /*---- Calculation of rotation angle's sine & cosine ----*/
+                    App = A1[p];
+                    Aqq = A2[q];
+                    y = 5.0e-1 * (App - Aqq);
+                    x = -Apq / sqrt( (double)Apq * Apq + (double)y * y );
+                    if( y < 0.0 )
+                        x = -x;
+                    s = x / sqrt( 2.0 * (1.0 + sqrt( 1.0 - (double)x * x )));
+                    s2 = s * s;
+                    c = sqrt( 1.0 - s2 );
+                    c2 = c * c;
+                    a = 2.0 * Apq * c * s;
+
+                    /*---- Apq annulation ----*/
+                    A3 = A;
+                    for( i = 0; i < p; i++, A3 += n )
+                    {
+                        Aip = A3[p];
+                        Aiq = A3[q];
+                        Vpi = V1[i];
+                        Vqi = V2[i];
+                        A3[p] = (float) (Aip * c - Aiq * s);
+                        A3[q] = (float) (Aiq * c + Aip * s);
+                        V1[i] = (float) (Vpi * c - Vqi * s);
+                        V2[i] = (float) (Vqi * c + Vpi * s);
+                    }
+                    for( ; i < q; i++, A3 += n )
+                    {
+                        Aip = A1[i];
+                        Aiq = A3[q];
+                        Vpi = V1[i];
+                        Vqi = V2[i];
+                        A1[i] = (float) (Aip * c - Aiq * s);
+                        A3[q] = (float) (Aiq * c + Aip * s);
+                        V1[i] = (float) (Vpi * c - Vqi * s);
+                        V2[i] = (float) (Vqi * c + Vpi * s);
+                    }
+                    for( ; i < n; i++ )
+                    {
+                        Aip = A1[i];
+                        Aiq = A2[i];
+                        Vpi = V1[i];
+                        Vqi = V2[i];
+                        A1[i] = (float) (Aip * c - Aiq * s);
+                        A2[i] = (float) (Aiq * c + Aip * s);
+                        V1[i] = (float) (Vpi * c - Vqi * s);
+                        V2[i] = (float) (Vqi * c + Vpi * s);
+                    }
+                    A1[p] = (float) (App * c2 + Aqq * s2 - a);
+                    A2[q] = (float) (App * s2 + Aqq * c2 + a);
+                    A1[q] = A2[p] = 0.0f;
+                }               /*q */
+            }                   /*p */
+        }
+        while( ind );
+        Amax /= n;
+    }                           /* while ( Amax > ax ) */
+
+    for( i = 0, k = 0; i < n; i++, k += n + 1 )
+        E[i] = A[k];
+    /*printf(" M = %d\n", M); */
+
+    /* -------- ordering -------- */
+    for( i = 0; i < n; i++ )
+    {
+        int m = i;
+        float Em = (float) fabs( E[i] );
+
+        for( j = i + 1; j < n; j++ )
+        {
+            float Ej = (float) fabs( E[j] );
+
+            m = (Em < Ej) ? j : m;
+            Em = (Em < Ej) ? Ej : Em;
+        }
+        if( m != i )
+        {
+            int l;
+            float b = E[i];
+
+            E[i] = E[m];
+            E[m] = b;
+            for( j = 0, k = i * n, l = m * n; j < n; j++, k++, l++ )
+            {
+                b = V[k];
+                V[k] = V[l];
+                V[l] = b;
+            }
+        }
+    }
+
+    return CV_NO_ERR;
+}
+
+/*F///////////////////////////////////////////////////////////////////////////////////////
+//    Name: icvCalcCovarMatrixEx_8u32fR
+//    Purpose: The function calculates a covariance matrix for a group of input objects
+//             (images, vectors, etc.). ROI supported.
+//    Context:
+//    Parameters:  nObjects    - number of source objects
+//                 objects     - array of pointers to ROIs of the source objects
+//                 imgStep     - full width of each source object row in bytes
+//                 avg         - pointer to averaged object
+//                 avgStep     - full width of averaged object row in bytes
+//                 size        - ROI size of each source and averaged objects
+//                 covarMatrix - covariance matrix (output parameter; must be allocated
+//                               before call)
+//
+//    Returns: CV_NO_ERR or error code
+//
+//    Notes:   
+//F*/
+static CvStatus  CV_STDCALL
+icvCalcCovarMatrixEx_8u32fR( int nObjects, void *input, int objStep1,
+                             int ioFlags, int ioBufSize, uchar* buffer,
+                             void *userData, float *avg, int avgStep,
+                             CvSize size, float *covarMatrix )
+{
+    int objStep = objStep1;
+
+    /* ---- TEST OF PARAMETERS ---- */
+
+    if( nObjects < 2 )
+        return CV_BADFACTOR_ERR;
+    if( ioFlags < 0 || ioFlags > 3 )
+        return CV_BADFACTOR_ERR;
+    if( ioFlags && ioBufSize < 1024 )
+        return CV_BADFACTOR_ERR;
+    if( ioFlags && buffer == NULL )
+        return CV_NULLPTR_ERR;
+    if( input == NULL || avg == NULL || covarMatrix == NULL )
+        return CV_NULLPTR_ERR;
+    if( size.width > objStep || 4 * size.width > avgStep || size.height < 1 )
+        return CV_BADSIZE_ERR;
+
+    avgStep /= 4;
+
+    if( ioFlags & CV_EIGOBJ_INPUT_CALLBACK )    /* ==== USE INPUT CALLBACK ==== */
+    {
+        int nio, ngr, igr, n = size.width * size.height, mm = 0;
+        CvCallback read_callback = ((CvInput *) & input)->callback;
+        uchar *buffer2;
+
+        objStep = n;
+        nio = ioBufSize / n;    /* number of objects in buffer */
+        ngr = nObjects / nio;   /* number of io groups */
+        if( nObjects % nio )
+            mm = 1;
+        ngr += mm;
+
+        buffer2 = (uchar *)cvAlloc( sizeof( uchar ) * n );
+        if( buffer2 == NULL )
+            return CV_OUTOFMEM_ERR;
+
+        for( igr = 0; igr < ngr; igr++ )
+        {
+            int k, l;
+            int io, jo, imin = igr * nio, imax = imin + nio;
+            uchar *bu1 = buffer, *bu2;
+
+            if( imax > nObjects )
+                imax = nObjects;
+
+            /* read igr group */
+            for( io = imin; io < imax; io++, bu1 += n )
+            {
+                CvStatus r;
+
+                r = (CvStatus)read_callback( io, (void *) bu1, userData );
+                if( r )
+                    return r;
+            }
+
+            /* diagonal square calc */
+            bu1 = buffer;
+            for( io = imin; io < imax; io++, bu1 += n )
+            {
+                bu2 = bu1;
+                for( jo = io; jo < imax; jo++, bu2 += n )
+                {
+                    float w = 0.f;
+                    float *fu = avg;
+                    int ij = 0;
+
+                    for( k = 0; k < size.height; k++, fu += avgStep )
+                        for( l = 0; l < size.width; l++, ij++ )
+                        {
+                            float f = fu[l], u1 = bu1[ij], u2 = bu2[ij];
+
+                            w += (u1 - f) * (u2 - f);
+                        }
+                    covarMatrix[io * nObjects + jo] = covarMatrix[jo * nObjects + io] = w;
+                }
+            }
+
+            /* non-diagonal elements calc */
+            for( jo = imax; jo < nObjects; jo++ )
+            {
+                CvStatus r;
+
+                bu1 = buffer;
+                bu2 = buffer2;
+
+                /* read jo object */
+                r = (CvStatus)read_callback( jo, (void *) bu2, userData );
+                if( r )
+                    return r;
+
+                for( io = imin; io < imax; io++, bu1 += n )
+                {
+                    float w = 0.f;
+                    float *fu = avg;
+                    int ij = 0;
+
+                    for( k = 0; k < size.height; k++, fu += avgStep )
+                    {
+                        for( l = 0; l < size.width - 3; l += 4, ij += 4 )
+                        {
+                            float f = fu[l];
+                            uchar u1 = bu1[ij];
+                            uchar u2 = bu2[ij];
+
+                            w += (u1 - f) * (u2 - f);
+                            f = fu[l + 1];
+                            u1 = bu1[ij + 1];
+                            u2 = bu2[ij + 1];
+                            w += (u1 - f) * (u2 - f);
+                            f = fu[l + 2];
+                            u1 = bu1[ij + 2];
+                            u2 = bu2[ij + 2];
+                            w += (u1 - f) * (u2 - f);
+                            f = fu[l + 3];
+                            u1 = bu1[ij + 3];
+                            u2 = bu2[ij + 3];
+                            w += (u1 - f) * (u2 - f);
+                        }
+                        for( ; l < size.width; l++, ij++ )
+                        {
+                            float f = fu[l], u1 = bu1[ij], u2 = bu2[ij];
+
+                            w += (u1 - f) * (u2 - f);
+                        }
+                    }
+                    covarMatrix[io * nObjects + jo] = covarMatrix[jo * nObjects + io] = w;
+                }
+            }
+        }                       /* igr */
+
+        cvFree( &buffer2 );
+    }                           /* if() */
+
+    else
+        /* ==== NOT USE INPUT CALLBACK ==== */
+    {
+        int i, j;
+        uchar **objects = (uchar **) (((CvInput *) & input)->data);
+
+        for( i = 0; i < nObjects; i++ )
+        {
+            uchar *bu = objects[i];
+
+            for( j = i; j < nObjects; j++ )
+            {
+                int k, l;
+                float w = 0.f;
+                float *a = avg;
+                uchar *bu1 = bu;
+                uchar *bu2 = objects[j];
+
+                for( k = 0; k < size.height;
+                     k++, bu1 += objStep, bu2 += objStep, a += avgStep )
+                {
+                    for( l = 0; l < size.width - 3; l += 4 )
+                    {
+                        float f = a[l];
+                        uchar u1 = bu1[l];
+                        uchar u2 = bu2[l];
+
+                        w += (u1 - f) * (u2 - f);
+                        f = a[l + 1];
+                        u1 = bu1[l + 1];
+                        u2 = bu2[l + 1];
+                        w += (u1 - f) * (u2 - f);
+                        f = a[l + 2];
+                        u1 = bu1[l + 2];
+                        u2 = bu2[l + 2];
+                        w += (u1 - f) * (u2 - f);
+                        f = a[l + 3];
+                        u1 = bu1[l + 3];
+                        u2 = bu2[l + 3];
+                        w += (u1 - f) * (u2 - f);
+                    }
+                    for( ; l < size.width; l++ )
+                    {
+                        float f = a[l];
+                        uchar u1 = bu1[l];
+                        uchar u2 = bu2[l];
+
+                        w += (u1 - f) * (u2 - f);
+                    }
+                }
+
+                covarMatrix[i * nObjects + j] = covarMatrix[j * nObjects + i] = w;
+            }
+        }
+    }                           /* else */
+
+    return CV_NO_ERR;
+}
+
+/*======================== end of icvCalcCovarMatrixEx_8u32fR ===========================*/
+
+
+static int
+icvDefaultBufferSize( void )
+{
+    return 10 * 1024 * 1024;
+}
+
+/*F///////////////////////////////////////////////////////////////////////////////////////
+//    Name: icvCalcEigenObjects_8u32fR
+//    Purpose: The function calculates an orthonormal eigen basis and a mean (averaged)
+//             object for a group of input objects (images, vectors, etc.). ROI supported.
+//    Context:
+//    Parameters: nObjects  - number of source objects
+//                input     - pointer either to array of pointers to input objects
+//                            or to read callback function (depending on ioFlags)
+//                imgStep   - full width of each source object row in bytes
+//                output    - pointer either to array of pointers to output eigen objects
+//                            or to write callback function (depending on ioFlags)
+//                eigStep   - full width of each eigenobject row in bytes
+//                size      - ROI size of each source object
+//                ioFlags   - input/output flags (see Notes)
+//                ioBufSize - input/output buffer size
+//                userData  - pointer to the structure which contains all necessary
+//                            data for the callback functions
+//                calcLimit - determines the calculation finish conditions
+//                avg       - pointer to averaged object (has the same size as ROI)
+//                avgStep   - full width of averaged object row in bytes
+//                eigVals   - pointer to corresponding eigenvalues (array of <nObjects>
+//                            elements in descending order)
+//
+//    Returns: CV_NO_ERR or error code
+//
+//    Notes: 1. input/output data (that is, input objects and eigen ones) may either
+//              be allocated in the RAM or be read from/written to the HDD (or any
+//              other device) by read/write callback functions. It depends on the
+//              value of ioFlags paramater, which may be the following:
+//                  CV_EIGOBJ_NO_CALLBACK, or 0;
+//                  CV_EIGOBJ_INPUT_CALLBACK;
+//                  CV_EIGOBJ_OUTPUT_CALLBACK;
+//                  CV_EIGOBJ_BOTH_CALLBACK, or
+//                            CV_EIGOBJ_INPUT_CALLBACK | CV_EIGOBJ_OUTPUT_CALLBACK.
+//              The callback functions as well as the user data structure must be
+//              developed by the user.
+//
+//           2. If ioBufSize = 0, or it's too large, the function dermines buffer size
+//              itself.
+//
+//           3. Depending on calcLimit parameter, calculations are finished either if
+//              eigenfaces number comes up to certain value or the relation of the
+//              current eigenvalue and the largest one comes down to certain value
+//              (or any of the above conditions takes place). The calcLimit->type value
+//              must be CV_TERMCRIT_NUMB, CV_TERMCRIT_EPS or
+//              CV_TERMCRIT_NUMB | CV_TERMCRIT_EPS. The function returns the real
+//              values calcLimit->max_iter and calcLimit->epsilon.
+//
+//           4. eigVals may be equal to NULL (if you don't need eigen values in further).
+//
+//F*/
+static CvStatus CV_STDCALL
+icvCalcEigenObjects_8u32fR( int nObjects, void* input, int objStep,
+                            void* output, int eigStep, CvSize size,
+                            int  ioFlags, int ioBufSize, void* userData,
+                            CvTermCriteria* calcLimit, float* avg,
+                            int    avgStep, float  *eigVals )
+{
+    int i, j, n, iev = 0, m1 = nObjects - 1, objStep1 = objStep, eigStep1 = eigStep / 4;
+    CvSize objSize, eigSize, avgSize;
+    float *c = 0;
+    float *ev = 0;
+    float *bf = 0;
+    uchar *buf = 0;
+    void *buffer = 0;
+    float m = 1.0f / (float) nObjects;
+    CvStatus r;
+
+    if( m1 > calcLimit->max_iter && calcLimit->type != CV_TERMCRIT_EPS )
+        m1 = calcLimit->max_iter;
+
+    /* ---- TEST OF PARAMETERS ---- */
+
+    if( nObjects < 2 )
+        return CV_BADFACTOR_ERR;
+    if( ioFlags < 0 || ioFlags > 3 )
+        return CV_BADFACTOR_ERR;
+    if( input == NULL || output == NULL || avg == NULL )
+        return CV_NULLPTR_ERR;
+    if( size.width > objStep || 4 * size.width > eigStep ||
+        4 * size.width > avgStep || size.height < 1 )
+        return CV_BADSIZE_ERR;
+    if( !(ioFlags & CV_EIGOBJ_INPUT_CALLBACK) )
+        for( i = 0; i < nObjects; i++ )
+            if( ((uchar **) input)[i] == NULL )
+                return CV_NULLPTR_ERR;
+    if( !(ioFlags & CV_EIGOBJ_OUTPUT_CALLBACK) )
+        for( i = 0; i < m1; i++ )
+            if( ((float **) output)[i] == NULL )
+                return CV_NULLPTR_ERR;
+
+    avgStep /= 4;
+    eigStep /= 4;
+
+    if( objStep == size.width && eigStep == size.width && avgStep == size.width )
+    {
+        size.width *= size.height;
+        size.height = 1;
+        objStep = objStep1 = eigStep = eigStep1 = avgStep = size.width;
+    }
+    objSize = eigSize = avgSize = size;
+
+    if( ioFlags & CV_EIGOBJ_INPUT_CALLBACK )
+    {
+        objSize.width *= objSize.height;
+        objSize.height = 1;
+        objStep = objSize.width;
+        objStep1 = size.width;
+    }
+
+    if( ioFlags & CV_EIGOBJ_OUTPUT_CALLBACK )
+    {
+        eigSize.width *= eigSize.height;
+        eigSize.height = 1;
+        eigStep = eigSize.width;
+        eigStep1 = size.width;
+    }
+
+    n = objSize.height * objSize.width * (ioFlags & CV_EIGOBJ_INPUT_CALLBACK) +
+        2 * eigSize.height * eigSize.width * (ioFlags & CV_EIGOBJ_OUTPUT_CALLBACK);
+
+    /* Buffer size determination */
+    if( ioFlags )
+    {
+        int size = icvDefaultBufferSize();
+        ioBufSize = MIN( size, n );
+    }
+
+    /* memory allocation (if necesseay) */
+
+    if( ioFlags & CV_EIGOBJ_INPUT_CALLBACK )
+    {
+        buf = (uchar *) cvAlloc( sizeof( uchar ) * objSize.width );
+        if( buf == NULL )
+            return CV_OUTOFMEM_ERR;
+    }
+
+    if( ioFlags )
+    {
+        buffer = (void *) cvAlloc( ioBufSize );
+        if( buffer == NULL )
+        {
+            if( buf )
+                cvFree( &buf );
+            return CV_OUTOFMEM_ERR;
+        }
+    }
+
+    /* Calculation of averaged object */
+    bf = avg;
+    for( i = 0; i < avgSize.height; i++, bf += avgStep )
+        for( j = 0; j < avgSize.width; j++ )
+            bf[j] = 0.f;
+
+    for( i = 0; i < nObjects; i++ )
+    {
+        int k, l;
+        uchar *bu = (ioFlags & CV_EIGOBJ_INPUT_CALLBACK) ? buf : ((uchar **) input)[i];
+
+        if( ioFlags & CV_EIGOBJ_INPUT_CALLBACK )
+        {
+            CvCallback read_callback = ((CvInput *) & input)->callback;
+
+            r = (CvStatus)read_callback( i, (void *) buf, userData );
+            if( r )
+            {
+                if( buffer )
+                    cvFree( &buffer );
+                if( buf )
+                    cvFree( &buf );
+                return r;
+            }
+        }
+
+        bf = avg;
+        for( k = 0; k < avgSize.height; k++, bf += avgStep, bu += objStep1 )
+            for( l = 0; l < avgSize.width; l++ )
+                bf[l] += bu[l];
+    }
+
+    bf = avg;
+    for( i = 0; i < avgSize.height; i++, bf += avgStep )
+        for( j = 0; j < avgSize.width; j++ )
+            bf[j] *= m;
+
+    /* Calculation of covariance matrix */
+    c = (float *) cvAlloc( sizeof( float ) * nObjects * nObjects );
+
+    if( c == NULL )
+    {
+        if( buffer )
+            cvFree( &buffer );
+        if( buf )
+            cvFree( &buf );
+        return CV_OUTOFMEM_ERR;
+    }
+
+    r = icvCalcCovarMatrixEx_8u32fR( nObjects, input, objStep1, ioFlags, ioBufSize,
+                                      (uchar *) buffer, userData, avg, 4 * avgStep, size, c );
+    if( r )
+    {
+        cvFree( &c );
+        if( buffer )
+            cvFree( &buffer );
+        if( buf )
+            cvFree( &buf );
+        return r;
+    }
+
+    /* Calculation of eigenvalues & eigenvectors */
+    ev = (float *) cvAlloc( sizeof( float ) * nObjects * nObjects );
+
+    if( ev == NULL )
+    {
+        cvFree( &c );
+        if( buffer )
+            cvFree( &buffer );
+        if( buf )
+            cvFree( &buf );
+        return CV_OUTOFMEM_ERR;
+    }
+
+    if( eigVals == NULL )
+    {
+        eigVals = (float *) cvAlloc( sizeof( float ) * nObjects );
+
+        if( eigVals == NULL )
+        {
+            cvFree( &c );
+            cvFree( &ev );
+            if( buffer )
+                cvFree( &buffer );
+            if( buf )
+                cvFree( &buf );
+            return CV_OUTOFMEM_ERR;
+        }
+        iev = 1;
+    }
+
+    r = icvJacobiEigens_32f( c, ev, eigVals, nObjects, 0.0f );
+    cvFree( &c );
+    if( r )
+    {
+        cvFree( &ev );
+        if( buffer )
+            cvFree( &buffer );
+        if( buf )
+            cvFree( &buf );
+        if( iev )
+            cvFree( &eigVals );
+        return r;
+    }
+
+    /* Eigen objects number determination */
+    if( calcLimit->type != CV_TERMCRIT_NUMBER )
+    {
+        for( i = 0; i < m1; i++ )
+            if( fabs( eigVals[i] / eigVals[0] ) < calcLimit->epsilon )
+                break;
+        m1 = calcLimit->max_iter = i;
+    }
+    else
+        m1 = calcLimit->max_iter;
+    calcLimit->epsilon = (float) fabs( eigVals[m1 - 1] / eigVals[0] );
+
+    for( i = 0; i < m1; i++ )
+        eigVals[i] = (float) (1.0 / sqrt( (double)eigVals[i] ));
+
+    /* ----------------- Calculation of eigenobjects ----------------------- */
+    if( ioFlags & CV_EIGOBJ_OUTPUT_CALLBACK )
+    {
+        int nio, ngr, igr;
+
+        nio = ioBufSize / (4 * eigSize.width);  /* number of eigen objects in buffer */
+        ngr = m1 / nio;         /* number of io groups */
+        if( nObjects % nio )
+            ngr += 1;
+
+        for( igr = 0; igr < ngr; igr++ )
+        {
+            int i, io, ie, imin = igr * nio, imax = imin + nio;
+
+            if( imax > m1 )
+                imax = m1;
+
+            for( i = 0; i < eigSize.width * (imax - imin); i++ )
+                ((float *) buffer)[i] = 0.f;
+
+            for( io = 0; io < nObjects; io++ )
+            {
+                uchar *bu = ioFlags & CV_EIGOBJ_INPUT_CALLBACK ? buf : ((uchar **) input)[io];
+
+                if( ioFlags & CV_EIGOBJ_INPUT_CALLBACK )
+                {
+                    CvCallback read_callback = ((CvInput *) & input)->callback;
+
+                    r = (CvStatus)read_callback( io, (void *) buf, userData );
+                    if( r )
+                    {
+                        cvFree( &ev );
+                        if( iev )
+                            cvFree( &eigVals );
+                        if( buffer )
+                            cvFree( &buffer );
+                        if( buf )
+                            cvFree( &buf );
+                        return r;
+                    }
+                }
+
+                for( ie = imin; ie < imax; ie++ )
+                {
+                    int k, l;
+                    uchar *bv = bu;
+                    float e = ev[ie * nObjects + io] * eigVals[ie];
+                    float *be = ((float *) buffer) + ((ie - imin) * eigStep);
+
+                    bf = avg;
+                    for( k = 0; k < size.height; k++, bv += objStep1,
+                         bf += avgStep, be += eigStep1 )
+                    {
+                        for( l = 0; l < size.width - 3; l += 4 )
+                        {
+                            float f = bf[l];
+                            uchar v = bv[l];
+
+                            be[l] += e * (v - f);
+                            f = bf[l + 1];
+                            v = bv[l + 1];
+                            be[l + 1] += e * (v - f);
+                            f = bf[l + 2];
+                            v = bv[l + 2];
+                            be[l + 2] += e * (v - f);
+                            f = bf[l + 3];
+                            v = bv[l + 3];
+                            be[l + 3] += e * (v - f);
+                        }
+                        for( ; l < size.width; l++ )
+                            be[l] += e * (bv[l] - bf[l]);
+                    }
+                }
+            }                   /* io */
+
+            for( ie = imin; ie < imax; ie++ )   /* calculated eigen objects writting */
+            {
+                CvCallback write_callback = ((CvInput *) & output)->callback;
+                float *be = ((float *) buffer) + ((ie - imin) * eigStep);
+
+                r = (CvStatus)write_callback( ie, (void *) be, userData );
+                if( r )
+                {
+                    cvFree( &ev );
+                    if( iev )
+                        cvFree( &eigVals );
+                    if( buffer )
+                        cvFree( &buffer );
+                    if( buf )
+                        cvFree( &buf );
+                    return r;
+                }
+            }
+        }                       /* igr */
+    }
+
+    else
+    {
+        int k, p, l;
+
+        for( i = 0; i < m1; i++ )       /* e.o. annulation */
+        {
+            float *be = ((float **) output)[i];
+
+            for( p = 0; p < eigSize.height; p++, be += eigStep )
+                for( l = 0; l < eigSize.width; l++ )
+                    be[l] = 0.0f;
+        }
+
+        for( k = 0; k < nObjects; k++ )
+        {
+            uchar *bv = (ioFlags & CV_EIGOBJ_INPUT_CALLBACK) ? buf : ((uchar **) input)[k];
+
+            if( ioFlags & CV_EIGOBJ_INPUT_CALLBACK )
+            {
+                CvCallback read_callback = ((CvInput *) & input)->callback;
+
+                r = (CvStatus)read_callback( k, (void *) buf, userData );
+                if( r )
+                {
+                    cvFree( &ev );
+                    if( iev )
+                        cvFree( &eigVals );
+                    if( buffer )
+                        cvFree( &buffer );
+                    if( buf )
+                        cvFree( &buf );
+                    return r;
+                }
+            }
+
+            for( i = 0; i < m1; i++ )
+            {
+                float v = eigVals[i] * ev[i * nObjects + k];
+                float *be = ((float **) output)[i];
+                uchar *bu = bv;
+
+                bf = avg;
+
+                for( p = 0; p < size.height; p++, bu += objStep1,
+                     bf += avgStep, be += eigStep1 )
+                {
+                    for( l = 0; l < size.width - 3; l += 4 )
+                    {
+                        float f = bf[l];
+                        uchar u = bu[l];
+
+                        be[l] += v * (u - f);
+                        f = bf[l + 1];
+                        u = bu[l + 1];
+                        be[l + 1] += v * (u - f);
+                        f = bf[l + 2];
+                        u = bu[l + 2];
+                        be[l + 2] += v * (u - f);
+                        f = bf[l + 3];
+                        u = bu[l + 3];
+                        be[l + 3] += v * (u - f);
+                    }
+                    for( ; l < size.width; l++ )
+                        be[l] += v * (bu[l] - bf[l]);
+                }
+            }                   /* i */
+        }                       /* k */
+    }                           /* else */
+
+    cvFree( &ev );
+    if( iev )
+        cvFree( &eigVals );
+    else
+        for( i = 0; i < m1; i++ )
+            eigVals[i] = 1.f / (eigVals[i] * eigVals[i]);
+    if( buffer )
+        cvFree( &buffer );
+    if( buf )
+        cvFree( &buf );
+    return CV_NO_ERR;
+}
+
+/* --- End of icvCalcEigenObjects_8u32fR --- */
+
+/*F///////////////////////////////////////////////////////////////////////////////////////
+//    Name: icvCalcDecompCoeff_8u32fR
+//    Purpose: The function calculates one decomposition coefficient of input object
+//             using previously calculated eigen object and the mean (averaged) object
+//    Context:
+//    Parameters:  obj     - input object
+//                 objStep - its step (in bytes)
+//                 eigObj  - pointer to eigen object
+//                 eigStep - its step (in bytes)
+//                 avg     - pointer to averaged object
+//                 avgStep - its step (in bytes)
+//                 size    - ROI size of each source object
+//
+//    Returns: decomposition coefficient value or large negative value (if error)
+//
+//    Notes:
+//F*/
+static float CV_STDCALL
+icvCalcDecompCoeff_8u32fR( uchar* obj, int objStep,
+                           float *eigObj, int eigStep,
+                           float *avg, int avgStep, CvSize size )
+{
+    int i, k;
+    float w = 0.0f;
+
+    if( size.width > objStep || 4 * size.width > eigStep
+        || 4 * size.width > avgStep || size.height < 1 )
+        return -1.0e30f;
+    if( obj == NULL || eigObj == NULL || avg == NULL )
+        return -1.0e30f;
+
+    eigStep /= 4;
+    avgStep /= 4;
+
+    if( size.width == objStep && size.width == eigStep && size.width == avgStep )
+    {
+        size.width *= size.height;
+        size.height = 1;
+        objStep = eigStep = avgStep = size.width;
+    }
+
+    for( i = 0; i < size.height; i++, obj += objStep, eigObj += eigStep, avg += avgStep )
+    {
+        for( k = 0; k < size.width - 4; k += 4 )
+        {
+            float o = (float) obj[k];
+            float e = eigObj[k];
+            float a = avg[k];
+
+            w += e * (o - a);
+            o = (float) obj[k + 1];
+            e = eigObj[k + 1];
+            a = avg[k + 1];
+            w += e * (o - a);
+            o = (float) obj[k + 2];
+            e = eigObj[k + 2];
+            a = avg[k + 2];
+            w += e * (o - a);
+            o = (float) obj[k + 3];
+            e = eigObj[k + 3];
+            a = avg[k + 3];
+            w += e * (o - a);
+        }
+        for( ; k < size.width; k++ )
+            w += eigObj[k] * ((float) obj[k] - avg[k]);
+    }
+
+    return w;
+}
+
+/*F///////////////////////////////////////////////////////////////////////////////////////
+//    Names: icvEigenDecomposite_8u32fR
+//    Purpose: The function calculates all decomposition coefficients for input object
+//             using previously calculated eigen objects basis and the mean (averaged)
+//             object
+//    Context:
+//    Parameters:  obj         - input object
+//                 objStep     - its step (in bytes)
+//                 nEigObjs    - number of eigen objects
+//                 eigInput    - pointer either to array of pointers to eigen objects
+//                               or to read callback function (depending on ioFlags)
+//                 eigStep     - eigen objects step (in bytes)
+//                 ioFlags     - input/output flags
+//                 iserData    - pointer to the structure which contains all necessary
+//                               data for the callback function
+//                 avg         - pointer to averaged object
+//                 avgStep     - its step (in bytes)
+//                 size        - ROI size of each source object
+//                 coeffs      - calculated coefficients (output data)
+//
+//    Returns: icv status
+//
+//    Notes:   see notes for icvCalcEigenObjects_8u32fR function
+//F*/
+static CvStatus CV_STDCALL
+icvEigenDecomposite_8u32fR( uchar * obj, int objStep, int nEigObjs,
+                            void *eigInput, int eigStep, int ioFlags,
+                            void *userData, float *avg, int avgStep,
+                            CvSize size, float *coeffs )
+{
+    int i;
+
+    if( nEigObjs < 2 )
+        return CV_BADFACTOR_ERR;
+    if( ioFlags < 0 || ioFlags > 1 )
+        return CV_BADFACTOR_ERR;
+    if( size.width > objStep || 4 * size.width > eigStep ||
+        4 * size.width > avgStep || size.height < 1 )
+        return CV_BADSIZE_ERR;
+    if( obj == NULL || eigInput == NULL || coeffs == NULL || avg == NULL )
+        return CV_NULLPTR_ERR;
+    if( !ioFlags )
+        for( i = 0; i < nEigObjs; i++ )
+            if( ((uchar **) eigInput)[i] == NULL )
+                return CV_NULLPTR_ERR;
+
+    if( ioFlags )               /* callback */
+
+    {
+        float *buffer;
+        CvCallback read_callback = ((CvInput *) & eigInput)->callback;
+
+        eigStep = 4 * size.width;
+
+        /* memory allocation */
+        buffer = (float *) cvAlloc( sizeof( float ) * size.width * size.height );
+
+        if( buffer == NULL )
+            return CV_OUTOFMEM_ERR;
+
+        for( i = 0; i < nEigObjs; i++ )
+        {
+            float w;
+            CvStatus r = (CvStatus)read_callback( i, (void *) buffer, userData );
+
+            if( r )
+            {
+                cvFree( &buffer );
+                return r;
+            }
+            w = icvCalcDecompCoeff_8u32fR( obj, objStep, buffer,
+                                            eigStep, avg, avgStep, size );
+            if( w < -1.0e29f )
+            {
+                cvFree( &buffer );
+                return CV_NOTDEFINED_ERR;
+            }
+            coeffs[i] = w;
+        }
+        cvFree( &buffer );
+    }
+
+    else
+        /* no callback */
+        for( i = 0; i < nEigObjs; i++ )
+        {
+            float w = icvCalcDecompCoeff_8u32fR( obj, objStep, ((float **) eigInput)[i],
+                                                  eigStep, avg, avgStep, size );
+
+            if( w < -1.0e29f )
+                return CV_NOTDEFINED_ERR;
+            coeffs[i] = w;
+        }
+
+    return CV_NO_ERR;
+}
+
+
+/*F///////////////////////////////////////////////////////////////////////////////////////
+//    Names: icvEigenProjection_8u32fR
+//    Purpose: The function calculates object projection to the eigen sub-space (restores
+//             an object) using previously calculated eigen objects basis, mean (averaged)
+//             object and decomposition coefficients of the restored object
+//    Context:
+//    Parameters:  nEigObjs - Number of eigen objects
+//                 eigens   - Array of pointers to eigen objects
+//                 eigStep  - Eigen objects step (in bytes)
+//                 coeffs   - Previously calculated decomposition coefficients
+//                 avg      - Pointer to averaged object
+//                 avgStep  - Its step (in bytes)
+//                 rest     - Pointer to restored object
+//                 restStep - Its step (in bytes)
+//                 size     - ROI size of each object
+//
+//    Returns: CV status
+//
+//    Notes:
+//F*/
+static CvStatus CV_STDCALL
+icvEigenProjection_8u32fR( int nEigObjs, void *eigInput, int eigStep,
+                           int ioFlags, void *userData, float *coeffs,
+                           float *avg, int avgStep, uchar * rest,
+                           int restStep, CvSize size )
+{
+    int i, j, k;
+    float *buf;
+    float *buffer = NULL;
+    float *b;
+    CvCallback read_callback = ((CvInput *) & eigInput)->callback;
+
+    if( size.width > avgStep || 4 * size.width > eigStep || size.height < 1 )
+        return CV_BADSIZE_ERR;
+    if( rest == NULL || eigInput == NULL || avg == NULL || coeffs == NULL )
+        return CV_NULLPTR_ERR;
+    if( ioFlags < 0 || ioFlags > 1 )
+        return CV_BADFACTOR_ERR;
+    if( !ioFlags )
+        for( i = 0; i < nEigObjs; i++ )
+            if( ((uchar **) eigInput)[i] == NULL )
+                return CV_NULLPTR_ERR;
+    eigStep /= 4;
+    avgStep /= 4;
+
+    if( size.width == restStep && size.width == eigStep && size.width == avgStep )
+    {
+        size.width *= size.height;
+        size.height = 1;
+        restStep = eigStep = avgStep = size.width;
+    }
+
+    buf = (float *) cvAlloc( sizeof( float ) * size.width * size.height );
+
+    if( buf == NULL )
+        return CV_OUTOFMEM_ERR;
+    b = buf;
+    for( i = 0; i < size.height; i++, avg += avgStep, b += size.width )
+        for( j = 0; j < size.width; j++ )
+            b[j] = avg[j];
+
+    if( ioFlags )
+    {
+        buffer = (float *) cvAlloc( sizeof( float ) * size.width * size.height );
+
+        if( buffer == NULL )
+        {
+            cvFree( &buf );
+            return CV_OUTOFMEM_ERR;
+        }
+        eigStep = size.width;
+    }
+
+    for( k = 0; k < nEigObjs; k++ )
+    {
+        float *e = ioFlags ? buffer : ((float **) eigInput)[k];
+        float c = coeffs[k];
+
+        if( ioFlags )           /* read eigen object */
+        {
+            CvStatus r = (CvStatus)read_callback( k, (void *) buffer, userData );
+
+            if( r )
+            {
+                cvFree( &buf );
+                cvFree( &buffer );
+                return r;
+            }
+        }
+
+        b = buf;
+        for( i = 0; i < size.height; i++, e += eigStep, b += size.width )
+        {
+            for( j = 0; j < size.width - 3; j += 4 )
+            {
+                float b0 = c * e[j];
+                float b1 = c * e[j + 1];
+                float b2 = c * e[j + 2];
+                float b3 = c * e[j + 3];
+
+                b[j] += b0;
+                b[j + 1] += b1;
+                b[j + 2] += b2;
+                b[j + 3] += b3;
+            }
+            for( ; j < size.width; j++ )
+                b[j] += c * e[j];
+        }
+    }
+
+    b = buf;
+    for( i = 0; i < size.height; i++, avg += avgStep, b += size.width, rest += restStep )
+        for( j = 0; j < size.width; j++ )
+        {
+            int w = cvRound( b[j] );
+
+            w = !(w & ~255) ? w : w < 0 ? 0 : 255;
+            rest[j] = (uchar) w;
+        }
+
+    cvFree( &buf );
+    if( ioFlags )
+        cvFree( &buffer );
+    return CV_NO_ERR;
+}
+
+/*F///////////////////////////////////////////////////////////////////////////////////////
+//    Name: cvCalcCovarMatrixEx
+//    Purpose: The function calculates a covariance matrix for a group of input objects
+//             (images, vectors, etc.).
+//    Context:
+//    Parameters:  nObjects    - number of source objects
+//                 input       - pointer either to array of input objects
+//                               or to read callback function (depending on ioFlags)
+//                 ioFlags     - input/output flags (see Notes to
+//                               cvCalcEigenObjects function)
+//                 ioBufSize   - input/output buffer size
+//                 userData    - pointer to the structure which contains all necessary
+//                               data for the callback functions
+//                 avg         - averaged object
+//                 covarMatrix - covariance matrix (output parameter; must be allocated
+//                               before call)
+//
+//    Notes:  See Notes to cvCalcEigenObjects function
+//F*/
+
+CV_IMPL void
+cvCalcCovarMatrixEx( int  nObjects, void*  input, int  ioFlags,
+                     int  ioBufSize, uchar*  buffer, void*  userData,
+                     IplImage* avg, float*  covarMatrix )
+{
+    float *avg_data;
+    int avg_step = 0;
+    CvSize avg_size;
+    int i;
+
+    CV_FUNCNAME( "cvCalcCovarMatrixEx" );
+
+    __BEGIN__;
+
+    cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &avg_size );
+    if( avg->depth != IPL_DEPTH_32F )
+        CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
+    if( avg->nChannels != 1 )
+        CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
+
+    if( ioFlags == CV_EIGOBJ_NO_CALLBACK )
+    {
+        IplImage **images = (IplImage **) (((CvInput *) & input)->data);
+        uchar **objects = (uchar **) cvAlloc( sizeof( uchar * ) * nObjects );
+        int img_step = 0, old_step = 0;
+        CvSize img_size = avg_size, old_size = avg_size;
+
+        if( objects == NULL )
+            CV_ERROR( CV_StsBadArg, "Insufficient memory" );
+
+        for( i = 0; i < nObjects; i++ )
+        {
+            IplImage *img = images[i];
+            uchar *img_data;
+
+            cvGetImageRawData( img, &img_data, &img_step, &img_size );
+            if( img->depth != IPL_DEPTH_8U )
+                CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
+            if( img_size != avg_size || img_size != old_size )
+                CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
+            if( img->nChannels != 1 )
+                CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
+            if( i > 0 && img_step != old_step )
+                CV_ERROR( CV_StsBadArg, "Different steps of objects" );
+
+            old_step = img_step;
+            old_size = img_size;
+            objects[i] = img_data;
+        }
+
+        CV_CALL( icvCalcCovarMatrixEx_8u32fR( nObjects,
+                                              (void*) objects,
+                                              img_step,
+                                              CV_EIGOBJ_NO_CALLBACK,
+                                              0,
+                                              NULL,
+                                              NULL,
+                                              avg_data,
+                                              avg_step,
+                                              avg_size,
+                                              covarMatrix ));
+        cvFree( &objects );
+    }
+
+    else
+
+    {
+        CV_CALL( icvCalcCovarMatrixEx_8u32fR( nObjects,
+                                              input,
+                                              avg_step / 4,
+                                              ioFlags,
+                                              ioBufSize,
+                                              buffer,
+                                              userData,
+                                              avg_data,
+                                              avg_step,
+                                              avg_size,
+                                              covarMatrix ));
+    }
+
+    __END__;
+}
+
+/*F///////////////////////////////////////////////////////////////////////////////////////
+//    Name: cvCalcEigenObjects
+//    Purpose: The function calculates an orthonormal eigen basis and a mean (averaged)
+//             object for a group of input objects (images, vectors, etc.).
+//    Context:
+//    Parameters: nObjects  - number of source objects
+//                input     - pointer either to array of input objects
+//                            or to read callback function (depending on ioFlags)
+//                output    - pointer either to output eigen objects
+//                            or to write callback function (depending on ioFlags)
+//                ioFlags   - input/output flags (see Notes)
+//                ioBufSize - input/output buffer size
+//                userData  - pointer to the structure which contains all necessary
+//                            data for the callback functions
+//                calcLimit - determines the calculation finish conditions
+//                avg       - averaged object (has the same size as ROI)
+//                eigVals   - pointer to corresponding eigen values (array of <nObjects>
+//                            elements in descending order)
+//
+//    Notes: 1. input/output data (that is, input objects and eigen ones) may either
+//              be allocated in the RAM or be read from/written to the HDD (or any
+//              other device) by read/write callback functions. It depends on the
+//              value of ioFlags paramater, which may be the following:
+//                  CV_EIGOBJ_NO_CALLBACK, or 0;
+//                  CV_EIGOBJ_INPUT_CALLBACK;
+//                  CV_EIGOBJ_OUTPUT_CALLBACK;
+//                  CV_EIGOBJ_BOTH_CALLBACK, or
+//                            CV_EIGOBJ_INPUT_CALLBACK | CV_EIGOBJ_OUTPUT_CALLBACK.
+//              The callback functions as well as the user data structure must be
+//              developed by the user.
+//
+//           2. If ioBufSize = 0, or it's too large, the function dermines buffer size
+//              itself.
+//
+//           3. Depending on calcLimit parameter, calculations are finished either if
+//              eigenfaces number comes up to certain value or the relation of the
+//              current eigenvalue and the largest one comes down to certain value
+//              (or any of the above conditions takes place). The calcLimit->type value
+//              must be CV_TERMCRIT_NUMB, CV_TERMCRIT_EPS or
+//              CV_TERMCRIT_NUMB | CV_TERMCRIT_EPS. The function returns the real
+//              values calcLimit->max_iter and calcLimit->epsilon.
+//
+//           4. eigVals may be equal to NULL (if you don't need eigen values in further).
+//
+//F*/
+CV_IMPL void
+cvCalcEigenObjects( int       nObjects,
+                    void*     input,
+                    void*     output,
+                    int       ioFlags,
+                    int       ioBufSize,
+                    void*     userData,
+                    CvTermCriteria* calcLimit,
+                    IplImage* avg, 
+                    float*    eigVals )
+{
+    float *avg_data;
+    int avg_step = 0;
+    CvSize avg_size;
+    int i;
+    int nEigens = nObjects - 1;
+
+    CV_FUNCNAME( "cvCalcEigenObjects" );
+
+    __BEGIN__;
+
+    cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &avg_size );
+    if( avg->depth != IPL_DEPTH_32F )
+        CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
+    if( avg->nChannels != 1 )
+        CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
+
+    if( nEigens > calcLimit->max_iter && calcLimit->type != CV_TERMCRIT_EPS )
+        nEigens = calcLimit->max_iter;
+
+    switch (ioFlags)
+    {
+    case CV_EIGOBJ_NO_CALLBACK:
+        {
+            IplImage **objects = (IplImage **) (((CvInput *) & input)->data);
+            IplImage **eigens = (IplImage **) (((CvInput *) & output)->data);
+            uchar **objs = (uchar **) cvAlloc( sizeof( uchar * ) * nObjects );
+            float **eigs = (float **) cvAlloc( sizeof( float * ) * nEigens );
+            int obj_step = 0, old_step = 0;
+            int eig_step = 0, oldeig_step = 0;
+            CvSize obj_size = avg_size, old_size = avg_size,
+
+                eig_size = avg_size, oldeig_size = avg_size;
+
+            if( objects == NULL || eigens == NULL )
+                CV_ERROR( CV_StsBadArg, "Insufficient memory" );
+
+            for( i = 0; i < nObjects; i++ )
+            {
+                IplImage *img = objects[i];
+                uchar *obj_data;
+
+                cvGetImageRawData( img, &obj_data, &obj_step, &obj_size );
+                if( img->depth != IPL_DEPTH_8U )
+                    CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
+                if( obj_size != avg_size || obj_size != old_size )
+                    CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
+                if( img->nChannels != 1 )
+                    CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
+                if( i > 0 && obj_step != old_step )
+                    CV_ERROR( CV_StsBadArg, "Different steps of objects" );
+
+                old_step = obj_step;
+                old_size = obj_size;
+                objs[i] = obj_data;
+            }
+            for( i = 0; i < nEigens; i++ )
+            {
+                IplImage *eig = eigens[i];
+                float *eig_data;
+
+                cvGetImageRawData( eig, (uchar **) & eig_data, &eig_step, &eig_size );
+                if( eig->depth != IPL_DEPTH_32F )
+                    CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
+                if( eig_size != avg_size || eig_size != oldeig_size )
+                    CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
+                if( eig->nChannels != 1 )
+                    CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
+                if( i > 0 && eig_step != oldeig_step )
+                    CV_ERROR( CV_StsBadArg, "Different steps of objects" );
+
+                oldeig_step = eig_step;
+                oldeig_size = eig_size;
+                eigs[i] = eig_data;
+            }
+            CV_CALL( icvCalcEigenObjects_8u32fR( nObjects, (void*) objs, obj_step,
+                                                 (void*) eigs, eig_step, obj_size,
+                                                 ioFlags, ioBufSize, userData,
+                                                 calcLimit, avg_data, avg_step, eigVals ));
+            cvFree( &objs );
+            cvFree( &eigs );
+            break;
+        }
+
+    case CV_EIGOBJ_OUTPUT_CALLBACK:
+        {
+            IplImage **objects = (IplImage **) (((CvInput *) & input)->data);
+            uchar **objs = (uchar **) cvAlloc( sizeof( uchar * ) * nObjects );
+            int obj_step = 0, old_step = 0;
+            CvSize obj_size = avg_size, old_size = avg_size;
+
+            if( objects == NULL )
+                CV_ERROR( CV_StsBadArg, "Insufficient memory" );
+
+            for( i = 0; i < nObjects; i++ )
+            {
+                IplImage *img = objects[i];
+                uchar *obj_data;
+
+                cvGetImageRawData( img, &obj_data, &obj_step, &obj_size );
+                if( img->depth != IPL_DEPTH_8U )
+                    CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
+                if( obj_size != avg_size || obj_size != old_size )
+                    CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
+                if( img->nChannels != 1 )
+                    CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
+                if( i > 0 && obj_step != old_step )
+                    CV_ERROR( CV_StsBadArg, "Different steps of objects" );
+
+                old_step = obj_step;
+                old_size = obj_size;
+                objs[i] = obj_data;
+            }
+            CV_CALL( icvCalcEigenObjects_8u32fR( nObjects,
+                                                 (void*) objs,
+                                                 obj_step,
+                                                 output,
+                                                 avg_step,
+                                                 obj_size,
+                                                 ioFlags,
+                                                 ioBufSize,
+                                                 userData,
+                                                 calcLimit,
+                                                 avg_data,
+                                                 avg_step,
+                                                 eigVals   ));
+            cvFree( &objs );
+            break;
+        }
+
+    case CV_EIGOBJ_INPUT_CALLBACK:
+        {
+            IplImage **eigens = (IplImage **) (((CvInput *) & output)->data);
+            float **eigs = (float**) cvAlloc( sizeof( float* ) * nEigens );
+            int eig_step = 0, oldeig_step = 0;
+            CvSize eig_size = avg_size, oldeig_size = avg_size;
+
+            if( eigens == NULL )
+                CV_ERROR( CV_StsBadArg, "Insufficient memory" );
+
+            for( i = 0; i < nEigens; i++ )
+            {
+                IplImage *eig = eigens[i];
+                float *eig_data;
+
+                cvGetImageRawData( eig, (uchar **) & eig_data, &eig_step, &eig_size );
+                if( eig->depth != IPL_DEPTH_32F )
+                    CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
+                if( eig_size != avg_size || eig_size != oldeig_size )
+                    CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
+                if( eig->nChannels != 1 )
+                    CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
+                if( i > 0 && eig_step != oldeig_step )
+                    CV_ERROR( CV_StsBadArg, "Different steps of objects" );
+
+                oldeig_step = eig_step;
+                oldeig_size = eig_size;
+                eigs[i] = eig_data;
+            }
+            CV_CALL( icvCalcEigenObjects_8u32fR( nObjects,
+                                                 input,
+                                                 avg_step / 4,
+                                                 (void*) eigs,
+                                                 eig_step,
+                                                 eig_size,
+                                                 ioFlags,
+                                                 ioBufSize,
+                                                 userData,
+                                                 calcLimit,
+                                                 avg_data,
+                                                 avg_step,
+                                                 eigVals   ));
+            cvFree( &eigs );
+            break;
+        }
+    case CV_EIGOBJ_INPUT_CALLBACK | CV_EIGOBJ_OUTPUT_CALLBACK:
+
+        CV_CALL( icvCalcEigenObjects_8u32fR( nObjects,
+                                             input,
+                                             avg_step / 4,
+                                             output,
+                                             avg_step,
+                                             avg_size,
+                                             ioFlags,
+                                             ioBufSize,
+                                             userData,
+                                             calcLimit,
+                                             avg_data,
+                                             avg_step,
+                                             eigVals   ));
+        break;
+
+    default:
+        CV_ERROR( CV_StsBadArg, "Unsupported i/o flag" );
+    }
+
+    __END__;
+}
+
+/*--------------------------------------------------------------------------------------*/
+/*F///////////////////////////////////////////////////////////////////////////////////////
+//    Name: cvCalcDecompCoeff
+//    Purpose: The function calculates one decomposition coefficient of input object
+//             using previously calculated eigen object and the mean (averaged) object
+//    Context:
+//    Parameters:  obj     - input object
+//                 eigObj  - eigen object
+//                 avg     - averaged object
+//
+//    Returns: decomposition coefficient value or large negative value (if error)
+//
+//    Notes:
+//F*/
+
+CV_IMPL double
+cvCalcDecompCoeff( IplImage * obj, IplImage * eigObj, IplImage * avg )
+{
+    double coeff = DBL_MAX;
+
+    uchar *obj_data;
+    float *eig_data;
+    float *avg_data;
+    int obj_step = 0, eig_step = 0, avg_step = 0;
+    CvSize obj_size, eig_size, avg_size;
+
+    CV_FUNCNAME( "cvCalcDecompCoeff" );
+
+    __BEGIN__;
+
+    cvGetImageRawData( obj, &obj_data, &obj_step, &obj_size );
+    if( obj->depth != IPL_DEPTH_8U )
+        CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
+    if( obj->nChannels != 1 )
+        CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
+
+    cvGetImageRawData( eigObj, (uchar **) & eig_data, &eig_step, &eig_size );
+    if( eigObj->depth != IPL_DEPTH_32F )
+        CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
+    if( eigObj->nChannels != 1 )
+        CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
+
+    cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &avg_size );
+    if( avg->depth != IPL_DEPTH_32F )
+        CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
+    if( avg->nChannels != 1 )
+        CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
+
+    if( obj_size != eig_size || obj_size != avg_size )
+        CV_ERROR( CV_StsBadArg, "different sizes of images" );
+
+    coeff = icvCalcDecompCoeff_8u32fR( obj_data, obj_step,
+                                       eig_data, eig_step,
+                                       avg_data, avg_step, obj_size );
+
+    __END__;
+    
+    return coeff;
+}
+
+/*--------------------------------------------------------------------------------------*/
+/*F///////////////////////////////////////////////////////////////////////////////////////
+//    Names: cvEigenDecomposite
+//    Purpose: The function calculates all decomposition coefficients for input object
+//             using previously calculated eigen objects basis and the mean (averaged)
+//             object
+//
+//    Parameters:  obj         - input object
+//                 nEigObjs    - number of eigen objects
+//                 eigInput    - pointer either to array of pointers to eigen objects
+//                               or to read callback function (depending on ioFlags)
+//                 ioFlags     - input/output flags
+//                 userData    - pointer to the structure which contains all necessary
+//                               data for the callback function
+//                 avg         - averaged object
+//                 coeffs      - calculated coefficients (output data)
+//
+//    Notes:   see notes for cvCalcEigenObjects function
+//F*/
+
+CV_IMPL void
+cvEigenDecomposite( IplImage* obj,
+                    int       nEigObjs,
+                    void*     eigInput,
+                    int       ioFlags, 
+                    void*     userData, 
+                    IplImage* avg, 
+                    float*    coeffs )
+{
+    float *avg_data;
+    uchar *obj_data;
+    int avg_step = 0, obj_step = 0;
+    CvSize avg_size, obj_size;
+    int i;
+
+    CV_FUNCNAME( "cvEigenDecomposite" );
+
+    __BEGIN__;
+
+    cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &avg_size );
+    if( avg->depth != IPL_DEPTH_32F )
+        CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
+    if( avg->nChannels != 1 )
+        CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
+
+    cvGetImageRawData( obj, &obj_data, &obj_step, &obj_size );
+    if( obj->depth != IPL_DEPTH_8U )
+        CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
+    if( obj->nChannels != 1 )
+        CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
+
+    if( obj_size != avg_size )
+        CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
+
+    if( ioFlags == CV_EIGOBJ_NO_CALLBACK )
+    {
+        IplImage **eigens = (IplImage **) (((CvInput *) & eigInput)->data);
+        float **eigs = (float **) cvAlloc( sizeof( float * ) * nEigObjs );
+        int eig_step = 0, old_step = 0;
+        CvSize eig_size = avg_size, old_size = avg_size;
+
+        if( eigs == NULL )
+            CV_ERROR( CV_StsBadArg, "Insufficient memory" );
+
+        for( i = 0; i < nEigObjs; i++ )
+        {
+            IplImage *eig = eigens[i];
+            float *eig_data;
+
+            cvGetImageRawData( eig, (uchar **) & eig_data, &eig_step, &eig_size );
+            if( eig->depth != IPL_DEPTH_32F )
+                CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
+            if( eig_size != avg_size || eig_size != old_size )
+                CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
+            if( eig->nChannels != 1 )
+                CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
+            if( i > 0 && eig_step != old_step )
+                CV_ERROR( CV_StsBadArg, "Different steps of objects" );
+
+            old_step = eig_step;
+            old_size = eig_size;
+            eigs[i] = eig_data;
+        }
+
+        CV_CALL( icvEigenDecomposite_8u32fR( obj_data,
+                                             obj_step,
+                                             nEigObjs,
+                                             (void*) eigs,
+                                             eig_step,
+                                             ioFlags,
+                                             userData,
+                                             avg_data,
+                                             avg_step,
+                                             obj_size,
+                                             coeffs   ));
+        cvFree( &eigs );
+    }
+
+    else
+
+    {
+        CV_CALL( icvEigenDecomposite_8u32fR( obj_data,
+                                             obj_step,
+                                             nEigObjs,
+                                             eigInput,
+                                             avg_step,
+                                             ioFlags,
+                                             userData,
+                                             avg_data,
+                                             avg_step,
+                                             obj_size,
+                                             coeffs   ));
+    }
+
+    __END__;
+}
+
+/*--------------------------------------------------------------------------------------*/
+/*F///////////////////////////////////////////////////////////////////////////////////////
+//    Name: cvEigenProjection
+//    Purpose: The function calculates object projection to the eigen sub-space (restores
+//             an object) using previously calculated eigen objects basis, mean (averaged)
+//             object and decomposition coefficients of the restored object
+//    Context:
+//    Parameters:  nEigObjs    - number of eigen objects
+//                 eigInput    - pointer either to array of pointers to eigen objects
+//                               or to read callback function (depending on ioFlags)
+//                 ioFlags     - input/output flags
+//                 userData    - pointer to the structure which contains all necessary
+//                               data for the callback function
+//                 coeffs      - array of decomposition coefficients
+//                 avg         - averaged object
+//                 proj        - object projection (output data)
+//
+//    Notes:   see notes for cvCalcEigenObjects function
+//F*/
+
+CV_IMPL void
+cvEigenProjection( void*     eigInput,
+                   int       nEigObjs,
+                   int       ioFlags,
+                   void*     userData,
+                   float*    coeffs, 
+                   IplImage* avg,
+                   IplImage* proj )
+{
+    float *avg_data;
+    uchar *proj_data;
+    int avg_step = 0, proj_step = 0;
+    CvSize avg_size, proj_size;
+    int i;
+
+    CV_FUNCNAME( "cvEigenProjection" );
+
+    __BEGIN__;
+
+    cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &avg_size );
+    if( avg->depth != IPL_DEPTH_32F )
+        CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
+    if( avg->nChannels != 1 )
+        CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
+
+    cvGetImageRawData( proj, &proj_data, &proj_step, &proj_size );
+    if( proj->depth != IPL_DEPTH_8U )
+        CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
+    if( proj->nChannels != 1 )
+        CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
+
+    if( proj_size != avg_size )
+        CV_ERROR( CV_StsBadArg, "Different sizes of projects" );
+
+    if( ioFlags == CV_EIGOBJ_NO_CALLBACK )
+    {
+        IplImage **eigens = (IplImage**) (((CvInput *) & eigInput)->data);
+        float **eigs = (float**) cvAlloc( sizeof( float * ) * nEigObjs );
+        int eig_step = 0, old_step = 0;
+        CvSize eig_size = avg_size, old_size = avg_size;
+
+        if( eigs == NULL )
+            CV_ERROR( CV_StsBadArg, "Insufficient memory" );
+
+        for( i = 0; i < nEigObjs; i++ )
+        {
+            IplImage *eig = eigens[i];
+            float *eig_data;
+
+            cvGetImageRawData( eig, (uchar **) & eig_data, &eig_step, &eig_size );
+            if( eig->depth != IPL_DEPTH_32F )
+                CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
+            if( eig_size != avg_size || eig_size != old_size )
+                CV_ERROR( CV_StsBadArg, "Different sizes of objects" );
+            if( eig->nChannels != 1 )
+                CV_ERROR( CV_BadNumChannels, cvUnsupportedFormat );
+            if( i > 0 && eig_step != old_step )
+                CV_ERROR( CV_StsBadArg, "Different steps of objects" );
+
+            old_step = eig_step;
+            old_size = eig_size;
+            eigs[i] = eig_data;
+        }
+
+        CV_CALL( icvEigenProjection_8u32fR( nEigObjs,
+                                            (void*) eigs,
+                                            eig_step,
+                                            ioFlags,
+                                            userData,
+                                            coeffs,
+                                            avg_data,
+                                            avg_step,
+                                            proj_data,
+                                            proj_step,
+                                            avg_size   ));
+        cvFree( &eigs );
+    }
+
+    else
+
+    {
+        CV_CALL( icvEigenProjection_8u32fR( nEigObjs,
+                                            eigInput,
+                                            avg_step,
+                                            ioFlags,
+                                            userData,
+                                            coeffs,
+                                            avg_data,
+                                            avg_step,
+                                            proj_data,
+                                            proj_step,
+                                            avg_size   ));
+    }
+
+    __END__;
+}
+
+/* End of file. */