vcs.maemo.org Git - opencv/blob - tests/cv/src/aeigenobjects.inc

   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 /*______________________________________________________________________________________*/
  43 /*                                                                                      */
  44 /*              Test functions for the Eigen Objects functions group                    */
  45 /*______________________________________________________________________________________*/
  46
  47 #include "cvtest.h"
  48
  49 static int _cvCalcCovarMatrix_8u32fR_q( int      nObjects,
  50                                      uchar**  objects,
  51                                      int      objStep,
  52                                      float*   avg,
  53                                      int      avgStep,
  54                                      CvSize size,
  55                                      float*   covarMatrix )
  56 {
  57     int i, j;
  58
  59     if ( nObjects < 2 )                                          return CV_BADFACTOR_ERR;
  60     if ( size.width > objStep || 4*size.width > avgStep || size.height < 1)
  61                                                                  return CV_BADSIZE_ERR;
  62     if ( objects == NULL || avg == NULL || covarMatrix == NULL ) return CV_NULLPTR_ERR;
  63
  64     avgStep /= 4;
  65
  66     for(i=0; i<nObjects; i++)
  67     {
  68         uchar* bu = objects[i];
  69         for(j=i; j<nObjects; j++)
  70         {
  71             int     ij = i*nObjects + j, k, l;
  72             float    w = 0.f;
  73             float*   a = avg;
  74             uchar* bu1 = bu;
  75             uchar* bu2 = objects[j];
  76
  77             for(k=0; k<size.height; k++, bu1 += objStep, bu2 += objStep, a += avgStep)
  78             {
  79                 for(l = 0; l < size.width-3; l+=4)
  80                 {
  81                     float  f = a  [l];
  82                     uchar u1 = bu1[l];
  83                     uchar u2 = bu2[l];
  84                     w += ( u1 - f ) * ( u2 - f );
  85                     f  = a  [l+1];
  86                     u1 = bu1[l+1];
  87                     u2 = bu2[l+1];
  88                     w += ( u1 - f ) * ( u2 - f );
  89                     f  = a  [l+2];
  90                     u1 = bu1[l+2];
  91                     u2 = bu2[l+2];
  92                     w += ( u1 - f ) * ( u2 - f );
  93                     f  = a  [l+3];
  94                     u1 = bu1[l+3];
  95                     u2 = bu2[l+3];
  96                     w += ( u1 - f ) * ( u2 - f );
  97                 }
  98                 for(; l < size.width; l++)
  99                 {
 100                     float  f = a  [l];
 101                     uchar u1 = bu1[l];
 102                     uchar u2 = bu2[l];
 103                     w += ( u1 - f ) * ( u2 - f );
 104                 }
 105             }
 106
 107             covarMatrix[ij] = w;
 108             ij = j*nObjects + i;
 109             covarMatrix[ij] = w;
 110         }
 111     }
 112
 113     return CV_NO_ERR;
 114 }   /* end of _cvCalcCovarMatrix_8u32fR */
 115
 116
 117 /* copy of _cvJacobiEigen_32f */
 118 int _cvJacobiEigens_32f ( float* A,
 119                                float* V,
 120                                float* E,
 121                                int    n,
 122                                float  eps )
 123 {
 124     int i, j, k, ind;
 125     float *AA = A, *VV = V;
 126     double Amax, anorm=0, ax;
 127
 128     if ( A == NULL || V == NULL || E == NULL ) return CV_NULLPTR_ERR;
 129     if ( n <= 0 )                              return CV_BADSIZE_ERR;
 130     if (eps < 1.0e-7f )  eps = 1.0e-7f;
 131
 132     /*-------- Prepare --------*/
 133     for(i=0; i<n; i++, VV+=n, AA+=n)
 134     {
 135         for(j=0; j<i; j++)
 136         {
 137             double Am = AA[j];
 138             anorm += Am*Am;
 139         }
 140         for(j=0; j<n; j++) VV[j] = 0.f;
 141         VV[i] = 1.f;
 142     }
 143
 144     anorm = sqrt( anorm + anorm );
 145     ax    = anorm*eps/n;
 146     Amax  = anorm;
 147
 148     while ( Amax > ax )
 149     {
 150         Amax /= n;
 151         do  /* while (ind) */
 152         {
 153             int p, q;
 154             float *V1  = V, *A1  = A;
 155             ind = 0;
 156             for(p=0; p<n-1; p++, A1+=n, V1+=n)
 157             {
 158                 float *A2 = A + n*(p+1), *V2 = V + n*(p+1);
 159                 for(q=p+1; q<n; q++, A2+=n, V2+=n)
 160                 {
 161                     double x, y, c, s, c2, s2, a;
 162                     float *A3, Apq=A1[q], App, Aqq, Aip, Aiq, Vpi, Vqi;
 163                     if( fabs( Apq ) < Amax ) continue;
 164
 165                     ind=1;
 166
 167                     /*---- Calculation of rotation angle's sine & cosine ----*/
 168                     App = A1[p];
 169                     Aqq = A2[q];
 170                     y   = 5.0e-1*(App - Aqq);
 171                     x = -Apq / sqrt(Apq*Apq + y*y);
 172                     if(y<0.0) x = -x;
 173                     s = x / sqrt(2.0*(1.0 + sqrt(1.0 - x*x)));
 174                     s2 = s*s;
 175                     c  = sqrt(1.0 - s2);
 176                     c2 = c*c;
 177                     a  = 2.0*Apq*c*s;
 178
 179                     /*---- Apq annulation ----*/
 180                     A3 = A;
 181                     for(i=0; i<p; i++, A3+=n)
 182                     {
 183                         Aip = A3[p];
 184                         Aiq = A3[q];
 185                         Vpi = V1[i];
 186                         Vqi = V2[i];
 187                         A3[p] = (float)(Aip*c - Aiq*s);
 188                         A3[q] = (float)(Aiq*c + Aip*s);
 189                         V1[i] = (float)(Vpi*c - Vqi*s);
 190                         V2[i] = (float)(Vqi*c + Vpi*s);
 191                     }
 192                     for(; i<q; i++, A3+=n)
 193                     {
 194                         Aip = A1[i];
 195                         Aiq = A3[q];
 196                         Vpi = V1[i];
 197                         Vqi = V2[i];
 198                         A1[i] = (float)(Aip*c - Aiq*s);
 199                         A3[q] = (float)(Aiq*c + Aip*s);
 200                         V1[i] = (float)(Vpi*c - Vqi*s);
 201                         V2[i] = (float)(Vqi*c + Vpi*s);
 202                     }
 203                     for(; i<n; i++)
 204                     {
 205                         Aip = A1[i];
 206                         Aiq = A2[i];
 207                         Vpi = V1[i];
 208                         Vqi = V2[i];
 209                         A1[i] = (float)(Aip*c - Aiq*s);
 210                         A2[i] = (float)(Aiq*c + Aip*s);
 211                         V1[i] = (float)(Vpi*c - Vqi*s);
 212                         V2[i] = (float)(Vqi*c + Vpi*s);
 213                     }
 214                     A1[p] = (float)(App*c2 + Aqq*s2 - a);
 215                     A2[q] = (float)(App*s2 + Aqq*c2 + a);
 216                     A1[q] = A2[p] = 0.0f;
 217                 } /*q*/
 218             }     /*p*/
 219         } while (ind);
 220         Amax /= n;
 221     }   /* while ( Amax > ax ) */
 222
 223     for(i=0, k=0; i<n; i++, k+=n+1) E[i] = A[k];
 224     /*printf(" M = %d\n", M);*/
 225
 226     /* -------- ordering --------*/
 227     for(i=0; i<n; i++)
 228     {
 229         int m = i;
 230         float Em = (float)fabs(E[i]);
 231         for(j=i+1; j<n; j++)
 232         {
 233             float Ej = (float)fabs(E[j]);
 234             m  = ( Em < Ej ) ?  j :  m;
 235             Em = ( Em < Ej ) ? Ej : Em;
 236         }
 237         if( m != i )
 238         {
 239             int l;
 240             float b = E[i];
 241             E[i] = E[m];
 242             E[m] = b;
 243             for(j=0, k=i*n, l=m*n; j<n; j++, k++, l++)
 244             {
 245                 b    = V[k];
 246                 V[k] = V[l];
 247                 V[l] = b;
 248             }
 249         }
 250     }
 251
 252     return CV_NO_ERR;
 253 }
 254
 255
 256
 257 /*______________________________________________________________________________________*/
 258
 259 int  _cvCalcEigenObjects_8u32fR_q( int      nObjects,
 260                                       uchar**  objects,
 261                                       int      objStep,
 262                                       float**  eigObjs,
 263                                       int      eigStep,
 264                                       CvSize size,
 265                                       float*   eigVals,
 266                                       float*   avg,
 267                                       int      avgStep,
 268                                       int*     nEigObjs,
 269                                       double*  eps      )
 270 {
 271     int i, j, k, l;
 272     uchar *bu;
 273     float *c=0, *ev=0, *bf, *bf1, *bf2, m;
 274     int  r;
 275
 276
 277     if ( nObjects < 2 )                                         return CV_BADFACTOR_ERR;
 278     if ( size.width > objStep || 4*size.width > eigStep ||
 279          4*size.width > avgStep || size.height < 1)             return CV_BADSIZE_ERR;
 280     if ( objects == NULL || eigObjs == NULL || eigVals == NULL ||
 281          avg == NULL || nEigObjs == NULL || eps == NULL )       return CV_NULLPTR_ERR;
 282     for( i=0; i<nObjects;  i++ ) if( objects[i] == NULL ) return CV_NULLPTR_ERR;
 283     for( i=0; i<*nEigObjs; i++ ) if( eigObjs[i] == NULL ) return CV_NULLPTR_ERR;
 284
 285     eigStep /= 4;
 286     avgStep /= 4;
 287
 288 /* Calculation of averaged object */
 289     bf = avg;
 290     for(i = 0; i < size.height; i++, bf += avgStep)
 291         for(j = 0; j < size.width; j++)
 292             bf[j] = 0.f;
 293     for(k = 0; k < nObjects; k++)
 294     {
 295         bu = objects[k];
 296         bf = avg;
 297         for(i = 0; i < size.height; i++, bu +=objStep, bf += avgStep)
 298             for(j = 0; j < size.width; j++)
 299                 bf[j] += bu[j];
 300     }
 301     m = 1.0f/(float)nObjects;
 302     bf = avg;
 303     for(i = 0; i < size.height; i++, bf += avgStep)
 304         for(j = 0; j < size.width; j++)
 305             bf[j] *= m;
 306
 307 /* Calculation of covariance matrix */
 308     c  = (float*)cvAlloc ( sizeof(float)*nObjects*nObjects );
 309     if(c==NULL) return CV_OUTOFMEM_ERR;
 310
 311     r = _cvCalcCovarMatrix_8u32fR_q ( nObjects, objects, objStep,
 312                                      avg, 4*avgStep, size, c );
 313     if(r) { cvFree( &c );  return r; }
 314
 315 /* Calculation of eigenvalues & eigenvectors */
 316     ev = (float*)cvAlloc ( sizeof(float)*nObjects*nObjects );
 317     if(ev==NULL) { cvFree( &c );  return CV_OUTOFMEM_ERR; }
 318
 319     _cvJacobiEigens_32f( c, ev, eigVals, nObjects, 0.0f );
 320     cvFree( &c );
 321
 322     for(i=0; i<*nEigObjs; i++) if( fabs(eigVals[i]/eigVals[0]) < *eps ) break;
 323     *nEigObjs = i;
 324     *eps = fabs(eigVals[*nEigObjs-1]/eigVals[0]);
 325
 326 /* Calculation of eigenobjects */
 327     bf2 = ev;
 328     for(i=0; i<*nEigObjs; i++, bf2+=nObjects)
 329     {
 330         float  e = (float)(1.0/sqrt(eigVals[i]));
 331         float* u = eigObjs[i];
 332
 333         bf  = u;
 334         for(l=0; l<size.height; l++, bf+=eigStep)
 335             for(j=0; j<size.width; j++) bf[j] = 0.0f;
 336
 337         for(k=0; k<nObjects; k++)
 338         {
 339             float v = e*bf2[k];
 340             bf  = u;
 341             bu  = objects[k];
 342             bf1 = avg;
 343             for(l=0; l<size.height; l++, bf+=eigStep, bf1+=avgStep, bu+=objStep)
 344                 for(j=0; j<size.width; j++) bf[j] += v * ((float)bu[j] - bf1[j]);
 345         }
 346     }
 347
 348     cvFree( &ev );
 349     return CV_NO_ERR;
 350 } /* --- End of _cvCalcEigenObjects_8u32fR --- */
 351 /*______________________________________________________________________________________*/
 352
 353 float _cvCalcDecompCoeff_8u32fR_q( uchar*  obj,
 354                                    int     objStep,
 355                                    float*  eigObj,
 356                                    int     eigStep,
 357                                    float*  avg,
 358                                    int     avgStep,
 359                                    CvSize size )
 360 {
 361     int i, k;
 362     float w = 0.0f;
 363
 364     if ( size.width > objStep || 4*size.width > eigStep
 365          || 4*size.width > avgStep || size.height < 1)  return -1.0e30f;
 366     if ( obj == NULL || eigObj == NULL || avg == NULL ) return -1.0e30f;
 367
 368     eigStep /= 4;
 369     avgStep /= 4;
 370
 371     for(i = 0; i < size.height; i++, obj += objStep, eigObj += eigStep, avg += avgStep)
 372         for(k = 0; k < size.width; k++)
 373             w += eigObj[k]*( (float)obj[k] - avg[k] );
 374
 375     return w;
 376 }
 377 /*______________________________________________________________________________________*/
 378
 379 int  _cvEigenDecomposite_8u32fR_q( uchar*  obj,
 380                                           int     objStep,
 381                                           int     nEigObjs,
 382                                           float** eigObjs,
 383                                           int     eigStep,
 384                                           float*  avg,
 385                                           int     avgStep,
 386                                           CvSize size,
 387                                           float*  coeffs )
 388 {
 389     int i;
 390
 391     if ( nEigObjs < 2 )                                    return CV_BADFACTOR_ERR;
 392     if ( size.width > objStep || 4*size.width > eigStep ||
 393          4*size.width > avgStep || size.height < 1)        return CV_BADSIZE_ERR;
 394     if ( obj == NULL || eigObjs == NULL || coeffs == NULL || avg == NULL)
 395         return CV_NULLPTR_ERR;
 396
 397     for(i=0; i<nEigObjs; i++)
 398     {
 399         float w = _cvCalcDecompCoeff_8u32fR_q( obj, objStep, eigObjs[i], eigStep,
 400                                               avg, avgStep, size );
 401         if( w < -1.0e29f ) return CV_NOTDEFINED_ERR;
 402         coeffs[i] = w;
 403     }
 404     return CV_NO_ERR;
 405 }
 406 /*______________________________________________________________________________________*/
 407
 408 int  _cvEigenProjection_8u32fR_q( int     nEigObjs,
 409                                          float** eigens,
 410                                          int     eigStep,
 411                                          float*  coeffs,
 412                                          float*  avg,
 413                                          int     avgStep,
 414                                          uchar*  rest,
 415                                          int     restStep,
 416                                          CvSize size )
 417 {
 418     int i, j, k;
 419
 420     if ( size.width > avgStep || 4*size.width > eigStep || size.height < 1)
 421                                                         return CV_BADSIZE_ERR;
 422     if ( rest == NULL || eigens == NULL || avg == NULL || coeffs == NULL )
 423                                                         return CV_NULLPTR_ERR;
 424     eigStep /= 4;
 425     avgStep /= 4;
 426
 427     for(i = 0; i < size.height; i++, rest+=restStep, avg+=avgStep)
 428     {
 429         int ij = i*eigStep;
 430         for(j = 0; j < size.width; j++, ij++)
 431         {
 432             float w = avg[j];
 433             for(k=0; k<nEigObjs-3; k+=4)
 434             {
 435                 float* b = eigens[k];
 436                 w += coeffs[k  ] * b[ij];
 437                 b = eigens [k+1];
 438                 w += coeffs[k+1] * b[ij];
 439                 b = eigens [k+2];
 440                 w += coeffs[k+2] * b[ij];
 441                 b = eigens [k+3];
 442                 w += coeffs[k+3] * b[ij];
 443             }
 444             for(; k<nEigObjs; k++)
 445             {
 446                 float* b = eigens[k];
 447                 w += coeffs[k] * b[ij];
 448             }
 449             w = w<-0.499999f ? -0.499999f : w>255.499f ? 255.499f : w;
 450             rest[j] = (uchar)cvRound( w );
 451         }
 452     }
 453     return CV_NO_ERR;
 454 }
 455 /*______________________________________________________________________________________*/
 456
 457 /*   << End  of  file >>  */