git.maemo.org Git - opencv/blob - tests/cv/src/aoptflowhs.cpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 //
  10 //                        Intel License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
  14 // Third party copyrights are property of their respective owners.
  15 //
  16 // Redistribution and use in source and binary forms, with or without modification,
  17 // are permitted provided that the following conditions are met:
  18 //
  19 //   * Redistribution's of source code must retain the above copyright notice,
  20 //     this list of conditions and the following disclaimer.
  21 //
  22 //   * Redistribution's in binary form must reproduce the above copyright notice,
  23 //     this list of conditions and the following disclaimer in the documentation
  24 //     and/or other materials provided with the distribution.
  25 //
  26 //   * The name of Intel Corporation may not be used to endorse or promote products
  27 //     derived from this software without specific prior written permission.
  28 //
  29 // This software is provided by the copyright holders and contributors "as is" and
  30 // any express or implied warranties, including, but not limited to, the implied
  31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  32 // In no event shall the Intel Corporation or contributors be liable for any direct,
  33 // indirect, incidental, special, exemplary, or consequential damages
  34 // (including, but not limited to, procurement of substitute goods or services;
  35 // loss of use, data, or profits; or business interruption) however caused
  36 // and on any theory of liability, whether in contract, strict liability,
  37 // or tort (including negligence or otherwise) arising in any way out of
  38 // the use of this software, even if advised of the possibility of such damage.
  39 //
  40 //M*/
  41
  42 #include "cvtest.h"
  43
  44 #if 0
  45
  46 /* Testing parameters */
  47 static char FuncName[] = "cvCalcOpticalFlowHS";
  48 static char TestName[] = "Optical flow (Horn & Schunck)";
  49 static char TestClass[] = "Algorithm";
  50
  51 static long lImageWidth;
  52 static long lImageHeight;
  53 static float lambda;
  54
  55 #define EPSILON 0.0001f
  56
  57 static int fmaCalcOpticalFlowHS( void )
  58 {
  59     /* Some Variables */
  60     int            i,j,k;
  61
  62     uchar*         roiA;
  63     uchar*         roiB;
  64
  65     float*         VelocityX;
  66     float*         VelocityY;
  67
  68     float*         auxVelocityX;
  69     float*         auxVelocityY;
  70
  71     float*         DerX;
  72     float*         DerY;
  73     float*         DerT;
  74
  75     long            lErrors = 0;
  76
  77     CvTermCriteria criteria;
  78
  79     int usePrevious;
  80
  81     int Stop = 0;
  82     int iteration = 0;
  83     float epsilon = 0;
  84
  85     static int  read_param = 0;
  86
  87     /* Initialization global parameters */
  88     if( !read_param )
  89     {
  90         read_param = 1;
  91         /* Reading test-parameters */
  92         trslRead( &lImageHeight, "300", "Image height" );
  93         trslRead( &lImageWidth, "300", "Image width" );
  94         trssRead( &lambda, "20", "lambda" );
  95     }
  96
  97     /* initialization - for warning disable */
  98     criteria.epsilon = 0;
  99     criteria.max_iter = 0;
 100     criteria.type  = 1;
 101
 102     /* Allocating memory for all frames */
 103     IplImage* imgA = cvCreateImage( cvSize(lImageWidth,lImageHeight), IPL_DEPTH_8U, 1 );
 104     IplImage* imgB = cvCreateImage( cvSize(lImageWidth,lImageHeight), IPL_DEPTH_8U, 1 );
 105
 106     IplImage* testVelocityX = cvCreateImage( cvSize(lImageWidth,lImageHeight), IPL_DEPTH_32F, 1 );
 107     IplImage* testVelocityY = cvCreateImage( cvSize(lImageWidth,lImageHeight), IPL_DEPTH_32F, 1 );
 108
 109     VelocityX = (float*)cvAlloc(  lImageWidth*lImageHeight * sizeof(float) );
 110     VelocityY = (float*)cvAlloc(  lImageWidth*lImageHeight * sizeof(float) );
 111
 112     auxVelocityX = (float*)cvAlloc(  lImageWidth*lImageHeight * sizeof(float) );
 113     auxVelocityY = (float*)cvAlloc(  lImageWidth*lImageHeight * sizeof(float) );
 114
 115     DerX = (float*)cvAlloc( lImageWidth*lImageHeight * sizeof(float) );
 116     DerY = (float*)cvAlloc( lImageWidth*lImageHeight * sizeof(float) );
 117     DerT = (float*)cvAlloc( lImageWidth*lImageHeight * sizeof(float) );
 118
 119     /* Filling images */
 120     ats1bInitRandom( 0, 255, (uchar*)imgA->imageData, lImageWidth * lImageHeight );
 121     ats1bInitRandom( 0, 255, (uchar*)imgB->imageData, lImageWidth * lImageHeight );
 122
 123     /* set ROI of images */
 124     roiA = (uchar*)imgA->imageData;
 125     roiB = (uchar*)imgB->imageData;
 126
 127     /* example of 3*3 ROI*/
 128     /*roiA[0] = 0;
 129     roiA[1] = 1;
 130     roiA[2] = 2;
 131     roiA[lImageWidth] = 0;
 132     roiA[lImageWidth+1] = 1;
 133     roiA[lImageWidth+2] = 2;
 134     roiA[2*lImageWidth] = 0;
 135     roiA[2*lImageWidth+1] = 1;
 136     roiA[2*lImageWidth+2] = 2;
 137
 138     roiB[0] = 1;
 139     roiB[1] = 2;
 140     roiB[2] = 3;
 141     roiB[lImageWidth] = 1;
 142     roiB[lImageWidth+1] = 2;
 143     roiB[lImageWidth+2] = 3;
 144     roiB[2*lImageWidth] = 1;
 145     roiB[2*lImageWidth+1] = 2;
 146     roiB[2*lImageWidth+2] = 3;*/
 147 /****************************************************************************************\
 148 *                  Calculate derivatives                                                 *
 149 \****************************************************************************************/
 150     for (i=0; i<lImageHeight; i++)
 151     {
 152         for(j=0; j<lImageWidth; j++)
 153         {
 154             int jr,jl,it,ib;
 155
 156             if ( j==lImageWidth-1 )
 157                 jr = lImageWidth-1;
 158             else jr = j + 1;
 159
 160             if ( j==0 )
 161                 jl = 0;
 162             else jl = j - 1;
 163
 164             if ( i==(lImageHeight - 1) )
 165                 ib = lImageHeight - 1;
 166             else ib = i + 1;
 167
 168             if ( i==0 )
 169                 it = 0;
 170             else it = i - 1;
 171
 172             DerX[ i*lImageWidth + j ] = (float)
 173                 (roiA[ (it)*imgA->widthStep + jr ]
 174                 - roiA[ (it)*imgA->widthStep + jl ]
 175                 + 2*roiA[ (i)*imgA->widthStep + jr ]
 176                 - 2*roiA[ (i)*imgA->widthStep + jl ]
 177                 + roiA[ (ib)*imgA->widthStep + jr ]
 178                 - roiA[ (ib)*imgA->widthStep + jl ])/8 ;
 179
 180             DerY[ i*lImageWidth + j ] = (float)
 181                 ( roiA[ (ib)*imgA->widthStep + jl ]
 182                 + 2*roiA[ (ib)*imgA->widthStep + j  ]
 183                 + roiA[ (ib)*imgA->widthStep + jr ]
 184                 - roiA[ (it)*imgA->widthStep + jl ]
 185                 - 2*roiA[ (it)*imgA->widthStep + j  ]
 186                 - roiA[ (it)*imgA->widthStep + jr ])/8  ;
 187
 188             DerT[ i*lImageWidth + j ] = (float)
 189                 (roiB[i*imgB->widthStep + j] - roiA[i*imgA->widthStep + j]);
 190         }
 191     }
 192 for( usePrevious = 0; usePrevious < 2; usePrevious++ )
 193 {
 194 /****************************************************************************************\
 195 *                    Cases                                                               *
 196 \****************************************************************************************/
 197     for ( k = 0; k < 4; k++ )
 198     {
 199         switch (k)
 200         {
 201         case 0:
 202             {
 203                 criteria.type = CV_TERMCRIT_ITER;
 204                 criteria.max_iter = 3;
 205
 206                 trsWrite( ATS_LST|ATS_CON,
 207                          "usePrevious = %d, criteria = ITER, max_iter = %d\n",
 208                          usePrevious, criteria.max_iter);
 209
 210                 break;
 211             }
 212         case 1:
 213             {
 214                 criteria.type = CV_TERMCRIT_EPS;
 215                 criteria.epsilon = 0.001f;
 216                 trsWrite( ATS_LST|ATS_CON,
 217                          "usePrevious = %d, criteria = EPS, epsilon = %f\n",
 218                          usePrevious, criteria.epsilon);
 219
 220                 break;
 221             }
 222         case 2:
 223             {
 224                 criteria.type = CV_TERMCRIT_EPS | CV_TERMCRIT_ITER;
 225                 criteria.epsilon = 0.0001f;
 226                 criteria.max_iter = 3;
 227                 trsWrite( ATS_LST|ATS_CON,
 228                          "usePrevious = %d,"
 229                          "criteria = EPS|ITER,"
 230                          "epsilon = %f, max_iter = %d\n",
 231                          usePrevious, criteria.epsilon, criteria.max_iter);
 232
 233                 break;
 234             }
 235         case 3:
 236             {
 237                 criteria.type = CV_TERMCRIT_EPS | CV_TERMCRIT_ITER;
 238                 criteria.epsilon = 0.00001f;
 239                 criteria.max_iter = 100;
 240                 trsWrite( ATS_LST|ATS_CON,
 241                          "usePrevious = %d,"
 242                          "criteria = EPS|ITER,"
 243                          "epsilon = %f, max_iter = %d\n",
 244                          usePrevious, criteria.epsilon, criteria.max_iter);
 245
 246                 break;
 247             }
 248         }
 249         Stop = 0;
 250
 251         /* Run CVL function */
 252         cvCalcOpticalFlowHS( imgA , imgB, usePrevious,
 253                              testVelocityX, testVelocityY,
 254                              lambda, criteria );
 255
 256         /* Calc by other way */
 257         if (!usePrevious)
 258         {
 259             /* Filling initial velocity with zero */
 260             for (i = 0; i < lImageWidth * lImageHeight; i++ )
 261             {
 262                 VelocityX[i] = 0 ;
 263                 VelocityY[i] = 0 ;
 264             }
 265         }
 266         iteration = 0;
 267         while ( !Stop )
 268         {
 269             float* oldX;
 270             float* oldY;
 271             float* newX;
 272             float* newY;
 273
 274             iteration++;
 275
 276             if ( iteration & 1 )
 277             {
 278                 oldX = VelocityX;
 279                 oldY = VelocityY;
 280                 newX = auxVelocityX;
 281                 newY = auxVelocityY;
 282             }
 283             else
 284             {
 285                 oldX = auxVelocityX;
 286                 oldY = auxVelocityY;
 287                 newX = VelocityX;
 288                 newY = VelocityY;
 289             }
 290
 291             for( i = 0; i < lImageHeight; i++)
 292             {
 293                 for(j = 0; j< lImageWidth; j++)
 294                 {
 295                     float aveX = 0;
 296                     float aveY = 0;
 297                     float dx,dy,dt;
 298
 299                     aveX +=(j==0) ? oldX[ i*lImageWidth + j ] : oldX[ i*lImageWidth + j-1 ];
 300                     aveX +=(j==lImageWidth-1) ? oldX[ i*lImageWidth + j ] :
 301                                               oldX[ i*lImageWidth + j+1 ];
 302                     aveX +=(i==0) ? oldX[ i*lImageWidth + j ] : oldX[ (i-1)*lImageWidth + j ];
 303                     aveX +=(i==lImageHeight-1) ? oldX[ i*lImageWidth + j ] :
 304                                                oldX[ (i+1)*lImageWidth + j ];
 305                     aveX /=4;
 306
 307                     aveY +=(j==0) ? oldY[ i*lImageWidth + j ] : oldY[ i*lImageWidth + j-1 ];
 308                     aveY +=(j==lImageWidth-1) ? oldY[ i*lImageWidth + j ] :
 309                                               oldY[ i*lImageWidth + j+1 ];
 310                     aveY +=(i==0) ? oldY[ i*lImageWidth + j ] : oldY[ (i-1)*lImageWidth + j ];
 311                     aveY +=(i==lImageHeight-1) ? oldY[ i*lImageWidth + j ] :
 312                                                oldY[ (i+1)*lImageWidth + j ];
 313                     aveY /=4;
 314
 315                     dx = DerX[ i*lImageWidth + j ];
 316                     dy = DerY[ i*lImageWidth + j ];
 317                     dt = DerT[ i*lImageWidth + j ];
 318
 319                     /* Horn & Schunck pure formulas */
 320                     newX[ i*lImageWidth + j ] = aveX - ( dx * aveX +
 321                                                        dy * aveY + dt ) * lambda * dx /
 322                                                        (1 + lambda * ( dx*dx + dy*dy ));
 323
 324                     newY[ i*lImageWidth + j ] = aveY - ( dx * aveX +
 325                                                        dy * aveY + dt ) * lambda * dy /
 326                                                        (1 + lambda * ( dx*dx + dy*dy ));
 327                 }
 328             }
 329             /* evaluate epsilon */
 330             epsilon = 0;
 331             for ( i = 0; i < lImageHeight; i++)
 332             {
 333                 for ( j = 0; j < lImageWidth; j++)
 334                 {
 335                     epsilon = MAX((float)fabs(newX[i*lImageWidth + j]
 336                                               - oldX[i*lImageWidth + j]), epsilon );
 337                     epsilon = MAX((float)fabs(newY[i*lImageWidth + j]
 338                                               - oldY[i*lImageWidth + j]), epsilon );
 339                 }
 340             }
 341
 342             switch (criteria.type)
 343             {
 344             case CV_TERMCRIT_ITER:
 345                 Stop = (criteria.max_iter == iteration );break;
 346             case CV_TERMCRIT_EPS:
 347                 Stop = (criteria.epsilon > epsilon );break;
 348             case CV_TERMCRIT_ITER|CV_TERMCRIT_EPS:
 349                 Stop = ( ( criteria.epsilon > epsilon    ) ||
 350                          ( criteria.max_iter == iteration ));
 351                 break;
 352             }
 353             if (Stop)
 354             {
 355                 if ( (newX != VelocityX) && (newY != VelocityY) )
 356                 {
 357                     memcpy( VelocityX, newX, lImageWidth * lImageHeight * sizeof(float) );
 358                     memcpy( VelocityY, newY, lImageWidth * lImageHeight * sizeof(float) );
 359                 }
 360             }
 361         }
 362         trsWrite( ATS_LST|ATS_CON,
 363                          "%d iterations are made\n", iteration );
 364
 365         for( i = 0; i < lImageHeight; i++)
 366         {
 367             for(j = 0; j< lImageWidth; j++)
 368             {
 369                 float tvx = ((float*)(testVelocityX->imageData + i*testVelocityX->widthStep))[j];
 370                 float tvy = ((float*)(testVelocityY->imageData + i*testVelocityY->widthStep))[j];
 371
 372                 if (( fabs( tvx - VelocityX[i*lImageWidth + j])>EPSILON )||
 373                     ( fabs( tvy - VelocityY[i*lImageWidth + j])>EPSILON ) )
 374                 {
 375                     //trsWrite( ATS_LST | ATS_CON, " ValueX %f \n",
 376                     //          testVelocityX[i*lROIWidth + j] );
 377                     //trsWrite( ATS_LST | ATS_CON, " mustX  %f \n",
 378                     //          VelocityX[i*lROIWidth + j] );
 379
 380                     //trsWrite( ATS_LST | ATS_CON, " ValueY %f \n",
 381                     //          testVelocityY[i*lROIWidth + j] );
 382                     //trsWrite( ATS_LST | ATS_CON, " mustY  %f \n",
 383                     //          VelocityY[i*lROIWidth + j] );
 384
 385                     //trsWrite( ATS_LST | ATS_CON, " Coordinates %d %d\n", i, j );
 386
 387                     lErrors++;
 388                 }
 389             }
 390         }
 391     }/* for */
 392     /* Filling initial velocity with zero */
 393     cvZero( testVelocityX );
 394     cvZero( testVelocityY );
 395     for (i = 0; i < lImageWidth * lImageHeight; i++ )
 396     {
 397         VelocityX[i] = 0 ;
 398         VelocityY[i] = 0 ;
 399     }
 400 }
 401
 402     /* Free memory */
 403     cvFree( &VelocityX );
 404     cvFree( &VelocityY );
 405     cvFree( &auxVelocityX );
 406     cvFree( &auxVelocityY );
 407
 408
 409     cvFree( &DerX );
 410     cvFree( &DerY );
 411     cvFree( &DerT );
 412
 413     cvReleaseImage( &imgA );
 414     cvReleaseImage( &imgB );
 415     cvReleaseImage( &testVelocityX );
 416     cvReleaseImage( &testVelocityY );
 417
 418
 419     if( lErrors == 0 ) return trsResult( TRS_OK, "No errors fixed for this text" );
 420     else return trsResult( TRS_FAIL, "Total fixed %d errors", lErrors );
 421 } /*fmaCalcOpticalFlowHS*/
 422
 423 void InitACalcOpticalFlowHS( void )
 424 {
 425     /* Registering test function */
 426     trsReg( FuncName, TestName, TestClass, fmaCalcOpticalFlowHS );
 427 } /* InitACalcOpticalFlowHS */
 428
 429 #endif
 430
 431 /* End of file. */