MotionEnergy.C

Go to the documentation of this file.
00001 /*!@file RCBot/Motion/MotionEnergy.C detect motion in an image stream   */
00002 // //////////////////////////////////////////////////////////////////// //
00003 // The iLab Neuromorphic Vision C++ Toolkit - Copyright (C) 2001 by the //
00004 // University of Southern California (USC) and the iLab at USC.         //
00005 // See http://iLab.usc.edu for information about this project.          //
00006 // //////////////////////////////////////////////////////////////////// //
00007 // Major portions of the iLab Neuromorphic Vision Toolkit are protected //
00008 // under the U.S. patent ``Computation of Intrinsic Perceptual Saliency //
00009 // in Visual Environments, and Applications'' by Christof Koch and      //
00010 // Laurent Itti, California Institute of Technology, 2001 (patent       //
00011 // pending; application number 09/912,225 filed July 23, 2001; see      //
00012 // http://pair.uspto.gov/cgi-bin/final/home.pl for current status).     //
00013 // //////////////////////////////////////////////////////////////////// //
00014 // This file is part of the iLab Neuromorphic Vision C++ Toolkit.       //
00015 //                                                                      //
00016 // The iLab Neuromorphic Vision C++ Toolkit is free software; you can   //
00017 // redistribute it and/or modify it under the terms of the GNU General  //
00018 // Public License as published by the Free Software Foundation; either  //
00019 // version 2 of the License, or (at your option) any later version.     //
00020 //                                                                      //
00021 // The iLab Neuromorphic Vision C++ Toolkit is distributed in the hope  //
00022 // that it will be useful, but WITHOUT ANY WARRANTY; without even the   //
00023 // implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR      //
00024 // PURPOSE.  See the GNU General Public License for more details.       //
00025 //                                                                      //
00026 // You should have received a copy of the GNU General Public License    //
00027 // along with the iLab Neuromorphic Vision C++ Toolkit; if not, write   //
00028 // to the Free Software Foundation, Inc., 59 Temple Place, Suite 330,   //
00029 // Boston, MA 02111-1307 USA.                                           //
00030 // //////////////////////////////////////////////////////////////////// //
00031 //
00032 // Primary maintainer for this file: Lior Elazary <lelazary@yahoo.com>
00033 // $HeadURL: svn://isvn.usc.edu/software/invt/trunk/saliency/src/RCBot/Motion/MotionEnergy.C $
00034 // $Id: MotionEnergy.C 13756 2010-08-04 21:57:32Z siagian $
00035 //
00036 
00037 #include "MotionEnergy.H"
00038 
00039 // ######################################################################
00040 // ##### MotionEnergyPyrBuilder Functions:
00041 // ######################################################################
00042 
00043 //debuging
00044 //#define DEBUG_MotionEnergy
00045 
00046 #ifdef DEBUG_MotionEnergy
00047 #include "GUI/DebugWin.H"
00048 #endif
00049 
00050 // ######################################################################
00051 template <class T>
00052 MotionEnergyPyrBuilder<T>::MotionEnergyPyrBuilder(const PyramidType typ,
00053                                                   const float gabor_theta,
00054                                                   const float intens,
00055                                                   const int timeDomainSize,
00056                                                   const float magThreshold ) :
00057   PyrBuilder<T>(),
00058   itsPtype(typ),
00059   itsGaborAngle(gabor_theta),
00060   itsGaborIntens(intens),
00061   itsTimeDomainSize(timeDomainSize),
00062   itsMagThreshold(magThreshold)
00063 {
00064   ASSERT(itsTimeDomainSize == 3 || itsTimeDomainSize == 5);
00065 }
00066 
00067 // ######################################################################
00068 template <class T>
00069 ImageSet<T> MotionEnergyPyrBuilder<T>::build(const Image<T>& image,
00070                                              const int firstlevel,
00071                                              const int depth,
00072                                              PyramidCache<T>* cache)
00073 {
00074   ImageSet<T> result(depth);
00075   return result;
00076 }
00077 
00078 // ######################################################################
00079 template <class T>
00080 void MotionEnergyPyrBuilder<T>::updateMotion(const Image<T>& image,
00081                                              const int depth)
00082 {
00083 
00084   // create a pyramid with the input image
00085   ImageSet<T> pyr = buildPyrGeneric(image, 0, depth, itsPtype,
00086                                     itsGaborAngle, itsGaborIntens);
00087 
00088   imgPyrQ.push_back(pyr); //build the time domain
00089 
00090   //the first frame (push the same frame itsTimeDomainSize-1 times)
00091   if (imgPyrQ.size() == 1){
00092     for(unsigned int i=0; i<itsTimeDomainSize-1; i++)
00093       imgPyrQ.push_back(pyr); //build the time domain
00094   }
00095 
00096   if (imgPyrQ.size() > itsTimeDomainSize)        //limit the image time queue
00097     imgPyrQ.pop_front();
00098 }
00099 
00100 // ######################################################################
00101 // get the motion from the image by using sobel grad detector
00102 // to detect edge orientations in the space and time
00103 template <class T>
00104 Image<float> MotionEnergyPyrBuilder<T>::buildHorizontalMotionLevel(int scale)
00105 {
00106   // if we have too few imgs in the Q, just return an empty image:
00107   if (imgPyrQ.size() < itsTimeDomainSize)
00108     return Image<float>();
00109 
00110   int width  = imgPyrQ[0][scale].getWidth();
00111   int height = imgPyrQ[0][scale].getHeight();
00112 
00113   Image<float> motion(width, height, ZEROS); // motion place holder
00114 
00115   // to make the operation look nice and readable
00116   // the parenthesis around the x are needed
00117   // otherwise a bug happens when using Pix(0,i-1)
00118 #define Pix(t,x) (*(imgT[t] + (y_pos*width)+(x)))
00119 
00120   // build the time domain pointers
00121   typename Image<T>::iterator imgT[itsTimeDomainSize];
00122   for (unsigned int i=0; i < itsTimeDomainSize; i++){
00123     imgT[i] = imgPyrQ[i][scale].beginw();
00124   }
00125 
00126   for (int y_pos = 0; y_pos < height; y_pos++)
00127     {
00128       // start and end before the template runs out
00129       for (unsigned int i = (itsTimeDomainSize/2); i < width-(itsTimeDomainSize/2); i++){
00130         // calculate the sobel Gx and Gy
00131 
00132         float Gx = 0, Gy = 0;
00133         // use a sobel 3x3 or 5x5. based on timedomainsize
00134         switch(itsTimeDomainSize){
00135           // 3x3 sobel operation
00136         case 3:
00137           Gx = -1*Pix(0,i-1) + 0 + 1*Pix(0,i+1) +
00138             -2*Pix(1,i-1) + 0 + 2*Pix(1,i+1) +
00139             -1*Pix(2,i-1) + 0 + 1*Pix(2,i+1);
00140 
00141           Gy = 1*Pix(0,i-1) +  2*Pix(0,i) +  1*Pix(0,i+1) +
00142             0*Pix(1,i-1) +  0*Pix(1,i) +  0*Pix(1,i+1) +
00143             -1*Pix(2,i-1) + -2*Pix(2,i) + -1*Pix(2,i+1);
00144           break;
00145 
00146           // 5x5 sobel operation
00147         case 5:
00148           Gx = -1*Pix(0,i-2) +  -2*Pix(0,i-1) +   0*Pix(0,i) +  2*Pix(0,i+1) +  1*Pix(0,i+2) +
00149             -4*Pix(1,i-2) +  -8*Pix(1,i-1) +   0*Pix(1,i) +  8*Pix(1,i+1) +  4*Pix(1,i+2) +
00150             -6*Pix(2,i-2) + -12*Pix(2,i-1) +   0*Pix(2,i) + 12*Pix(2,i+1) +  6*Pix(2,i+2) +
00151             -4*Pix(3,i-2) +  -8*Pix(3,i-1) +   0*Pix(3,i) +  8*Pix(3,i+1) +  4*Pix(3,i+2) +
00152             -1*Pix(4,i-2) +  -2*Pix(4,i-1) +   0*Pix(4,i) +  2*Pix(4,i+1) +  1*Pix(4,i+2);
00153 
00154           Gy =  1*Pix(0,i-2) +   4*Pix(0,i-1) +   6*Pix(0,i) +  4*Pix(0,i+1) +  1*Pix(0,i+2) +
00155             2*Pix(1,i-2) +   8*Pix(1,i-1) +  12*Pix(1,i) +  8*Pix(1,i+1) +  2*Pix(1,i+2) +
00156             0*Pix(2,i-2) +   0*Pix(2,i-1) +   0*Pix(2,i) +  0*Pix(2,i+1) +  0*Pix(2,i+2) +
00157             -2*Pix(3,i-2) +  -8*Pix(3,i-1) + -12*Pix(3,i) + -8*Pix(3,i+1) + -2*Pix(3,i+2) +
00158             -1*Pix(4,i-2) +  -4*Pix(4,i-1) +  -6*Pix(4,i) + -4*Pix(4,i+1) + -1*Pix(4,i+2);
00159         default:
00160           break;
00161         }
00162 
00163         double mag =  sqrt((Gx*Gx) + (Gy*Gy));
00164         double angle = 0; // this represents the vel of the motion. 90 no motion.
00165 
00166         if (mag > itsMagThreshold){        // try to not process motion as a result of noise
00167           // divide x by y and not y by x as in the orig sobel alg
00168           // this is to get the vector oriented perallel to the edge
00169 
00170           if (Gy == 0){ // watch for divide by 0
00171             angle = atan(0); // + M_PI/2;
00172           } else {
00173             angle = atan((double)(Gx/Gy)); // + (M_PI/2);
00174           }
00175 
00176           // flip angles to be at the top.
00177           // Could also be done by changing the convX and convY
00178           angle *= -1;
00179 
00180           // Have all the angles pointing in the positive direction by
00181           // flipping the angle if its less then 0
00182           if (angle < 0) angle += M_PI;
00183 
00184           // dismiss any angles that violates the motion
00185           if (angle > 0.15 && angle < M_PI-0.15) {
00186             motion.setVal(i, y_pos, (float)(angle) - (M_PI/2));
00187           } else {
00188             motion.setVal(i, y_pos,0);
00189           }
00190         } else {
00191           motion.setVal(i, y_pos, 0);
00192         }
00193       }
00194     }
00195 
00196   return motion;
00197 }
00198 
00199 // ######################################################################
00200 // get the motion from the image by using sobel grad detector
00201 // to detect edge orientations in the space and time
00202 template <class T>
00203 Image<float> MotionEnergyPyrBuilder<T>::buildVerticalMotionLevel(int scale)
00204 {
00205   // if we have too few imgs in the Q, just return an empty image:
00206   if (imgPyrQ.size() < itsTimeDomainSize)
00207     return Image<float>();
00208 
00209   int width  = imgPyrQ[0][scale].getWidth();
00210   int height = imgPyrQ[0][scale].getHeight();
00211 
00212   Image<float> motion(width, height, ZEROS); // motion place holder
00213 
00214   // build the time domain pointers
00215   typename Image<T>::iterator imgT[itsTimeDomainSize];
00216   for (unsigned int i=0; i<itsTimeDomainSize; i++)
00217     {
00218       imgT[i] = imgPyrQ[i][scale].beginw(); 
00219     }
00220 
00221   for (int x_pos = 0; x_pos<width; x_pos++){
00222 
00223 #ifdef DEBUG_MotionEnergy
00224     // build an image of one slice throught the time and y domain
00225     Image<byte> debugImg(itsTimeDomainSize, height, ZEROS);
00226 #endif
00227 
00228     // start and end before the template runs out
00229     for (unsigned int i = (itsTimeDomainSize/2); i < height-(itsTimeDomainSize/2); i++)
00230       {
00231         // calculate the sobel Gx and Gy
00232 
00233         float Gx=0, Gy=0;
00234 
00235         // to make the operation look nice and readable
00236         // the parenthesis around the y are needed because when
00237         // otherwise a bug happends when using Pix2(0,i-1)
00238 #define Pix2(t,y) (*(imgT[t] + ((y)*width)+x_pos))
00239 
00240 #ifdef DEBUG_MotionEnergy
00241         // build the image
00242         for (uint j=0; j<itsTimeDomainSize; j++)
00243           {
00244             debugImg.setVal(j,i,Pix2(j,i));
00245           }
00246 #endif
00247         // use a sobel 3x3 or 5x5. based on timedomainsize
00248         switch(itsTimeDomainSize){
00249           // 3x3 sobel operation
00250         case 3:
00251           Gx = -1*Pix2(0,i+1) + 0 + 1*Pix2(2,i+1)  +
00252                -2*Pix2(0,i)   + 0 + 2*Pix2(2,i)    +
00253                -1*Pix2(0,i-1) + 0 + 1*Pix2(2,i-1);
00254           
00255           Gy =  1*Pix2(0,i+1) +  2*Pix2(1,i+1) +  1*Pix2(2,i+1) +
00256                 0*Pix2(0,i)   +  0*Pix2(1,i)   +  0*Pix2(2,i)   +
00257                -1*Pix2(0,i-1) + -2*Pix2(1,i-1) + -1*Pix2(2,i-1);
00258 
00259           break;
00260 
00261           // 5x5 sobel operation
00262         case 5:
00263           Gx = -1*Pix2(0,i+2) + -2*Pix2(1,i+2) +  0*Pix2(2,i+2) +  2*Pix2(3,i+2) +  1*Pix2(4,i+2) +
00264                -4*Pix2(0,i+1) + -8*Pix2(1,i+1) +  0*Pix2(2,i+1) +  8*Pix2(3,i+1) +  4*Pix2(4,i+1) +
00265                -6*Pix2(0,i)   + -12*Pix2(1,i)  +  0*Pix2(2,i)   + 12*Pix2(3,i)   +  6*Pix2(4,i)   +
00266                -4*Pix2(0,i-1) + -8*Pix2(1,i-1) +  0*Pix2(2,i-1) +  8*Pix2(3,i-1) +  4*Pix2(4,i-1) +
00267                -1*Pix2(0,i-2) + -2*Pix2(1,i-2) +  0*Pix2(2,i-2) +  2*Pix2(3,i-2) +  1*Pix2(4,i-2);
00268 
00269           Gy =  1*Pix2(0,i+2) +  4*Pix2(1,i+2) +   6*Pix2(2,i+2) +  4*Pix2(3,i+2) +   1*Pix2(4,i+2) +
00270                 2*Pix2(0,i+1) +  8*Pix2(1,i+1) +  12*Pix2(2,i+1) +  8*Pix2(3,i+1) +   2*Pix2(4,i+1) +
00271                 0*Pix2(0,i)   +  0*Pix2(1,i)   +   0*Pix2(2,i)   +  0*Pix2(3,i)   +   0*Pix2(4,i)   +
00272                -2*Pix2(0,i-1) + -8*Pix2(1,i-1) + -12*Pix2(2,i-1) + -8*Pix2(3,i-1) +  -2*Pix2(4,i-1) +
00273                -1*Pix2(0,i-2) + -4*Pix2(1,i-2) +  -6*Pix2(2,i-2) + -4*Pix2(3,i-2) +  -1*Pix2(4,i-2);
00274         default:
00275           break;
00276         }
00277 
00278       double mag =  sqrt((Gx*Gx) + (Gy*Gy));
00279       double angle = 0; // this represents the vel of the motion. 90 no motion.
00280 
00281       if (mag > itsMagThreshold){        // try to not process motion as a result of noise
00282         // divide x by y and not y by x as in the orig sobel alg
00283         // this is to get the vector oriented perallel to the edge
00284 
00285         if (Gy == 0){ // watch for divide by 0
00286           angle = atan(0); // + M_PI/2;
00287         } else {
00288           angle = atan((double)(Gx/Gy)); // + (M_PI/2);
00289         }
00290 
00291         angle *= -1; // flip angles to be at the top. Could also be done by changing the
00292         // convX and convY
00293 
00294         // Have all the angles pointing in the positive direction by
00295         // fliping the angle if its less then 0
00296         if (angle < 0) angle += M_PI;
00297 
00298         // dismiss any angles that violate the motion
00299         if (angle > 0.15 && angle < M_PI-0.15) {
00300           motion.setVal(x_pos, i, (float)(angle) - (M_PI/2));
00301         } else {
00302           motion.setVal(x_pos, i, 0);
00303         }
00304       } else {
00305         motion.setVal(x_pos, i, 0);
00306       }
00307 
00308     }
00309 #ifdef DEBUG_MotionEnergy
00310     Image<PixRGB<byte> > img = toRGB(rescale(debugImg, 256, 256));
00311     SHOWIMG(img);
00312     //SHOWIMG(rescale(debugImg,Dims(256, 256)));
00313 #endif
00314   }
00315 
00316   return motion;
00317 }
00318 
00319 // ######################################################################
00320 template <class T>
00321 ImageSet<float> MotionEnergyPyrBuilder<T>::buildHorizontalMotion()
00322 {
00323   const int depth = imgPyrQ[0].size();
00324 
00325   ImageSet<float> result(depth);
00326 
00327   for (int scale = 0; scale < depth; ++scale)
00328     result[scale] = buildHorizontalMotionLevel(scale);
00329 
00330   return result;
00331 }
00332 
00333 // ######################################################################
00334 // get the motion from the image by using sobel grad detector to detect edge orientations
00335 // in the space and time
00336 template <class T>
00337 ImageSet<float> MotionEnergyPyrBuilder<T>::buildVerticalMotion()
00338 {
00339   const int depth = imgPyrQ[0].size();
00340 
00341   ImageSet<float> result(depth);
00342 
00343   for (int scale = 0; scale < depth; scale++)
00344     result[scale] = buildVerticalMotionLevel(scale);
00345 
00346   return result;
00347 }
00348 
00349 // ######################################################################
00350 template <class T>
00351 float MotionEnergyPyrBuilder<T>::DrawVectors(Image<T> &img, Image<float> &motion)
00352 {
00353   //TODO: should check the mag is the same size as dir
00354 
00355   Image<float>::const_iterator mag_ptr = motion.begin();
00356   Image<float>::const_iterator mag_stop = motion.end();
00357 
00358   int inx=0;
00359   int avg_i = 0;
00360   double avg_angle = 0;
00361 
00362   while (mag_ptr != mag_stop)
00363     {
00364       int y = inx/motion.getWidth();
00365       int x = inx - (y*motion.getWidth());
00366 
00367       if (*mag_ptr != 0) {
00368         avg_i++;
00369         //avg_angle += (*mag_ptr+M_PI/2);
00370         avg_angle += (*mag_ptr);
00371 
00372         int scale_x = x * (img.getWidth()/motion.getWidth());
00373         int scale_y = y * (img.getHeight()/motion.getHeight());
00374         drawLine(img, Point2D<int>(scale_x,scale_y),
00375                  Point2D<int>((int)(scale_x+25*cos((*mag_ptr))),
00376                          (int)(scale_y-25*sin((*mag_ptr)))),
00377                  (byte)0);
00378       }
00379       mag_ptr++;
00380       inx++;
00381     }
00382 
00383   if (avg_i > 0){
00384     int xi = img.getWidth()/2;
00385     int yi = img.getHeight()/2;
00386 
00387     drawLine(img,Point2D<int>(xi, yi),
00388              Point2D<int>((int)(xi+75*cos(avg_angle/avg_i)),
00389                      (int)(yi-75*sin(avg_angle/avg_i))),
00390              (byte)0, 3);
00391     return avg_angle/avg_i;
00392   }
00393 
00394   return -999;
00395 }
00396 
00397 // ######################################################################
00398 template <class T>
00399 MotionEnergyPyrBuilder<T>* MotionEnergyPyrBuilder<T>::clone() const
00400 { return new MotionEnergyPyrBuilder<T>(*this); }
00401 
00402 // ######################################################################
00403 template <class T>
00404 void MotionEnergyPyrBuilder<T>::reset()
00405 {
00406   // imgPyrQ.clear();
00407 }
00408 
00409 template class MotionEnergyPyrBuilder<byte>;
00410 
00411 // ######################################################################
00412 /* So things look consistent in everyone's emacs... */
00413 /* Local Variables: */
00414 /* indent-tabs-mode: nil */
00415 /* End: */
Generated on Sun May 8 08:41:17 2011 for iLab Neuromorphic Vision Toolkit by  doxygen 1.6.3