BeoSubPipe.C

Go to the documentation of this file.
00001 /*!@file BeoSub/BeoSubPipe.C find pipe     */
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: Michael Montalbo <montalbo@usc.edu>
00033 // $HeadURL: svn://isvn.usc.edu/software/invt/trunk/saliency/src/BeoSub/BeoSubPipe.C $
00034 // $Id: BeoSubPipe.C 12782 2010-02-05 22:14:30Z irock $
00035 
00036 #include "Image/OpenCVUtil.H"  // must be first to avoid conflicting defs of int64, uint64
00037 
00038 #include "GUI/XWinManaged.H"
00039 
00040 #include "Image/Pixels.H"
00041 #include "Raster/Raster.H"
00042 
00043 #include "Util/Timer.H"
00044 #include "Util/Types.H"
00045 #include "Util/log.H"
00046 
00047 #include "Image/ColorOps.H"
00048 #include "Image/CutPaste.H"
00049 #include "Image/MathOps.H"
00050 #include "Image/DrawOps.H"
00051 #include "Image/FilterOps.H"
00052 #include "Image/Transforms.H"
00053 
00054 #include "BeoSub/hysteresis.H"
00055 #include "VFAT/segmentImageTrackMC.H"
00056 #include "BeoSub/HoughTransform.H"
00057 #include "BeoSub/ColorTracker.H"
00058 //#include "BeoSub/getPipeDir.H"
00059 #include "BeoSub/CannyEdge.H"
00060 #include "BeoSub/IsolateColor.H"
00061 #include "rutz/compat_cmath.h" // for isnan()
00062 
00063 #include "BeoSub/BeoSubPipe.H"
00064 
00065 #include <cstdio>
00066 #include <cstdlib>
00067 #include <cstring>
00068 #include <iostream> //needed for segmentImageTrackMC!
00069 #include <math.h>
00070 #include <vector>
00071 #include <cmath>
00072 
00073 // ######################################################################
00074 BeoSubPipe::BeoSubPipe()
00075 {
00076   houghThreshold = 30;
00077   minThreshold = 10;
00078   maxThreshold = 50;
00079   sigma = .7;
00080   tlow = 0.2;
00081   thigh = .97;
00082   linescale =80;
00083 
00084   // number of consecutive times that a line was found (going up to 8)
00085   foundCount = 0;
00086 
00087   itsWinInitialized = false;
00088   fNum = 0;
00089 }
00090 
00091 // ######################################################################
00092 BeoSubPipe::~BeoSubPipe()
00093 { }
00094 
00095 // ######################################################################
00096 std::vector<LineSegment2D> BeoSubPipe::getHoughLines
00097 (Image<byte> &cameraImage,
00098  Image< PixRGB <byte> > &outputImage)
00099 {
00100   std::vector <LineSegment2D> lines;
00101 
00102 #ifndef HAVE_OPENCV
00103   LFATAL("OpenCV must be installed in order to use this function");
00104 #else
00105 
00106   IplImage *edge = cvCreateImage( cvGetSize(img2ipl(cameraImage)), 8, 1 );
00107   cvCanny( img2ipl(luminance(cameraImage)), edge, 100, 150, 3 );//150,200,3
00108   //  inplacePaste(dispImage, toRGB(edgeImage), Point2D<int>(0,h));
00109 
00110   Image< PixRGB<byte> >houghImage = ipl2gray(edge);
00111   //  tim.reset();
00112   //  rutz::shared_ptr<CvMemStorage*> storage;
00113   //  *storage = cvCreateMemStorage(0);
00114   CvMemStorage* storage = cvCreateMemStorage(0);
00115 
00116   outputImage.clear();
00117 
00118   outputImage = toRGB(houghImage);
00119 
00120   //  CvSeq* cvlines = cvHoughLines2(edge, storage, CV_HOUGH_PROBABILISTIC, 1, CV_PI/180
00121   CvSeq* cvlines = cvHoughLines2(edge, storage, CV_HOUGH_PROBABILISTIC, 1, CV_PI/180, 80 ,30, 10);
00122   //  t = tim.get();
00123   //  LDEBUG("houghTransform takes: %f ms", (float)t/1000.0);
00124 
00125   for(int i = 0; i < cvlines->total; i++ )
00126     {
00127       CvPoint* line = (CvPoint*)cvGetSeqElem(cvlines,i);
00128       Point2D<int> pt1 = Point2D<int>(line[0].x,line[0].y);
00129       Point2D<int> pt2 = Point2D<int>(line[1].x,line[1].y);
00130       //      cvLine( color_dst, line[0], line[1], CV_RGB(255,0,0), 3, 8 );
00131 
00132       lines.push_back(LineSegment2D(pt1,pt2));
00133       drawLine(outputImage, pt1, pt2, PixRGB<byte>(255,0,0));
00134     }
00135 
00136 //     for(int i = 0; i < MIN(cvlines->total,10); i++ )
00137 //     {
00138 //       float* line = (float*)cvGetSeqElem(cvlines,i);
00139 //       float rho = line[0];
00140 //       float theta = line[1];
00141 //       CvPoint pt1, pt2;
00142 //       double a = cos(theta), b = sin(theta);
00143 //       double x0 = a*rho, y0 = b*rho;
00144 //       pt1.x = cvRound(x0 + 1000*(-b));
00145 //       pt1.y = cvRound(y0 + 1000*(a));
00146 //       pt2.x = cvRound(x0 - 1000*(-b));
00147 //       pt2.y = cvRound(y0 - 1000*(a));
00148 //       lines.push_back(LineSegment2D(Point2D<int>(pt1.x,pt1.y),Point2D<int>(pt2.x,pt2.y)));
00149 //       drawLine(outputImage, Point2D<int>(pt1.x,pt1.y), Point2D<int>(pt2.x, pt2.y), PixRGB<byte>(255,0,0));
00150 //     }
00151 
00152     //    delete edge;
00153     //    delete cvlines;
00154 
00155 #endif // HAVE_OPENCV
00156 
00157     return lines;
00158 }
00159 
00160 // ######################################################################
00161 
00162 float BeoSubPipe::pipeOrientation
00163 (Image< PixRGB <byte> > &cameraImage, Image<PixRGB <byte> > &outputImage)
00164 {
00165 
00166   // ------------ STEP 1: SEGMENT ORANGE  --------------------------
00167 
00168   Timer tim(1000000);
00169   uint w = cameraImage.getWidth();
00170   uint h = cameraImage.getHeight();
00171   Image<PixRGB<byte> > dispImage(w*2,h*2,NO_INIT);
00172   // Image<PixRGB<byte> > houghImage(w,h,NO_INIT);
00173 
00174 
00175   if(!itsWinInitialized)
00176     {
00177       itsWin.reset
00178         (new XWinManaged(Dims(w*2,h*2),
00179                          w+10, 0, "Pipe Direction Output"));
00180       itsWinInitialized = true;
00181     }
00182 
00183   Image<byte> orangeIsoImage(w,h, ZEROS);
00184 
00185   inplacePaste(dispImage, cameraImage, Point2D<int>(0,0));
00186 
00187   //Isolate the orange pixels in the image
00188   tim.reset();
00189   float orange = isolateOrange(cameraImage, orangeIsoImage); //, fnb=0;
00190   uint64 t = tim.get();
00191   LINFO("isoOrange[%f] takes: %f ms", orange, (float)t/1000.0);
00192   inplacePaste(dispImage, toRGB(orangeIsoImage), Point2D<int>(w,0));
00193 
00194   // ------------ STEP 2: GET EDGE PIXELS  --------------------------
00195 
00196 
00197   std::vector <LineSegment2D> lines;
00198   unsigned char *edge;
00199   char *dirfilename = NULL;
00200   canny(orangeIsoImage.getArrayPtr(), h, w,
00201         sigma, tlow, thigh, &edge, dirfilename);
00202 
00203   Image<byte> edgeImage(edge, w, h);
00204   inplacePaste(dispImage, toRGB(edgeImage), Point2D<int>(0,h));
00205 
00206   Image< PixRGB<byte> > houghImage = edgeImage;  //houghImage = orangeIsoImage;
00207   tim.reset();
00208   lines = houghTransform(edgeImage, M_PI/180, 1, houghThreshold, houghImage);
00209   t = tim.get();
00210   LINFO("houghTransform takes: %f ms", (float)t/1000.0);
00211 
00212   // ******* ---- ********
00213 
00214 //   std::vector <LineSegment2D> lines;
00215 
00216 //   tim.reset();
00217 //   lines = getHoughLines(orangeIsoImage, outputImage);
00218 //   //lines = houghTransform(edgeImage, M_PI/180, 1, houghThreshold, houghImage);
00219 //   t = tim.get();
00220 //   LINFO("houghTransform takes: %f ms", (float)t/1000.0);
00221 
00222   inplacePaste(dispImage, outputImage, Point2D<int>(0,h));
00223 
00224   // ------------ STEP 3: GET PIPE  --------------------------
00225 
00226     uint totlines=0;
00227   //std::vector<LineSegment2D> centerPipeLines = getPipes(lines);
00228   float pipeDir = getPipeDir(lines, 20, totlines);
00229 
00230   Point2D<int> p(cameraImage.getWidth()/2, cameraImage.getHeight()/2);
00231   Point2D<int> p2((int)(cameraImage.getWidth()/2+cos(pipeDir)*linescale),
00232              (int)(cameraImage.getHeight()/2+sin(pipeDir)*linescale));
00233   outputImage = houghImage;
00234 
00235 //  for(uint i = 0; i < centerPipeLines.size(); i++)
00236 //    {
00237 //       drawLine(outputImage, centerPipeLines[i].point1(), centerPipeLines[i].point2(), PixRGB <byte> (255, 255,0), 3);
00238   drawLine(outputImage, p, p2, PixRGB <byte> (255, 255,0), 3);
00239       //itsWin->setTitle(sformat("p_dir: %10.6f", pipeDir).c_str());
00240       //    }
00241 
00242   // ------------ STEP 4: REPORT (DISPLAY RESULTS)  ---------------
00243 
00244   inplacePaste(dispImage, outputImage, Point2D<int>(w, h));
00245   drawLine(dispImage, Point2D<int>(0,h),Point2D<int>(w*2-1,h),
00246            PixRGB<byte>(255,255,255),1);
00247   drawLine(dispImage, Point2D<int>(w,0),Point2D<int>(w,h*2-1),
00248            PixRGB<byte>(255,255,255),1);
00249   writeText( dispImage, Point2D<int>(0,0), sformat("Frame: %6d", fNum).c_str(),
00250              PixRGB<byte>(255,0,0));
00251   writeText( dispImage, Point2D<int>(w,0), sformat("Segment Color").c_str(),
00252              PixRGB<byte>(255,0,0));
00253   writeText( dispImage, Point2D<int>(0,h), sformat("Detect Edge").c_str(),
00254              PixRGB<byte>(255,0,0));
00255   writeText( dispImage, Point2D<int>(w,h), sformat("Identify Pipe").c_str(),
00256              PixRGB<byte>(255,0,0));
00257 
00258 //   std::string saveFName =  sformat("data/pipeDetect_%07d.ppm", fNum);
00259 //   LINFO("saving: %s",saveFName.c_str());
00260 //   Raster::WriteRGB(dispImage,saveFName);
00261   itsWin->drawImage(dispImage, 0, 0);
00262 
00263   fNum++;
00264 
00265   // decrease when:
00266   // it didn't find many lines, but there is some orange in the image
00267 
00268   // increase when:
00269   // it finds too many lines
00270 
00271   // increase and decrease the hough threshold
00272   // in an attempt to get some lines, but not too many.
00273 //   if (orange > 0 &&
00274 //       totlines < 20 &&
00275 //       houghThreshold > minThreshold) {
00276 //     houghThreshold--;  //make it more relaxed
00277 //     LINFO("Decreasing Hough Threshold");
00278 //   }
00279 //   else if (totlines > 20 && houghThreshold < maxThreshold) {
00280 //     houghThreshold++;  //make it stricter
00281 //     LINFO("Increasing Hough Threshold");
00282 //   }
00283 
00284 //   LINFO("Line Count: %"ZU, lines.size());
00285 //   LINFO("Adj Line Count: %d", totlines);
00286 //   LINFO("Orange: %d", orange);
00287 //   LINFO("Hough Threshold %d", houghThreshold);
00288 //   LINFO("Found Streak: %d", foundCount);
00289 
00290   return 0.0;//pipeDir;
00291 }
00292 
00293 // ######################################################################
00294 std::vector<LineSegment2D> BeoSubPipe::getPipes
00295 (const std::vector<LineSegment2D> lines)
00296 {
00297   uint nLines = lines.size();
00298   if(nLines == 0) { LDEBUG("*** NO LINES ***"); }
00299 
00300   std::vector< std::vector<LineSegment2D> > pipeLines;
00301 
00302   //Go through all the lines
00303   for(uint r = 0; r < nLines; r++)
00304     {
00305       int lnIndex = -1;
00306 
00307       //check to see if the current lines fits into a bucket
00308       for(uint c = 0; c < pipeLines.size(); c++)
00309         {
00310           if(lines[r].angleBetween(pipeLines[c][0]) < 5*(M_PI/180))//convert to radians
00311           {
00312             lnIndex = c;
00313             break;
00314           }
00315 
00316         }
00317 
00318       //if the line fits into a pre-existing bucket, add it to the bucket
00319       if(lnIndex > 0)
00320         {
00321           pipeLines[lnIndex].push_back(lines[r]);
00322           //average the old bucket's value with the new line added
00323           //so as to create a moving bucket
00324           Point2D<int> newPt1 = Point2D<int>(((lines[r].point1().i + pipeLines[lnIndex][0].point1().i)/2),
00325                             ((lines[r].point1().j + pipeLines[lnIndex][0].point1().j)/2));
00326 
00327           Point2D<int> newPt2 = Point2D<int>(((lines[r].point2().i + pipeLines[lnIndex][0].point2().i)/2),
00328                             ((lines[r].point2().j + pipeLines[lnIndex][0].point2().j)/2));
00329 
00330           pipeLines[lnIndex][0] = LineSegment2D(newPt1,newPt2);
00331 
00332         }
00333       //otherwise, create a new bucket
00334       else
00335         {
00336           std::vector<LineSegment2D> newCntrLines;
00337           newCntrLines.push_back(lines[r]);
00338           pipeLines.push_back(newCntrLines);
00339         }
00340     }
00341 
00342   std::vector<LineSegment2D> centerPipeLines;
00343 
00344   Point2D<int> two = Point2D<int>(2,2);
00345 
00346   for(uint c = 0; c < pipeLines.size(); c++)
00347     {
00348       if(pipeLines[c].size() == 2)
00349         {
00350           Point2D<int> endPoint1 = Point2D<int>((pipeLines[c][0].point1()+pipeLines[c][1].point1())/two);
00351           Point2D<int> endPoint2 = Point2D<int>((pipeLines[c][0].point2()+pipeLines[c][1].point2())/two);
00352 
00353           centerPipeLines.push_back(LineSegment2D(endPoint1,endPoint2));
00354         }
00355     }
00356 
00357   return centerPipeLines;
00358 
00359 //   // get the mean
00360 //   float sumLength = 0.0;  float meanLength = 0.0;
00361 //   for (uint i = 0; i < nLines; i++) sumLength += lines[i].length();
00362 //   meanLength = sumLength / nLines;
00363 
00364 //   // get the standard deviation
00365 //   float sumsqr = 0.0;
00366 //   for (uint i = 0; i < nLines; i++)
00367 //     sumsqr += pow((float)(lines[i].length()) - meanLength, 2.0);
00368 //   float stdevLength = sqrt(sumsqr / (nLines-1));
00369 
00370 //   // kick out lines that aren't within the stddev of the mean length
00371 //   LINFO ("Mean: %f   StdDev: %f", meanLength, stdevLength);
00372 //   std::vector<LineSegment2D> nlines;
00373 //   for(uint i = 0; i < nLines; i++)
00374 //     {
00375 //       if (lines[i].length() > (meanLength - stdevLength) &&
00376 //           lines[i].length() < (meanLength + stdevLength)   )
00377 //         nlines.push_back(lines[i]);
00378 //     }
00379 
00380 //   if (nlines.size() == 0)
00381 //     { LINFO("*** NO LINE LENGTH WITHIN STD DEV ***"); return avgPipeAngle; }
00382 
00383 //   std::vector<LineSegment2D> flines = nlines;
00384 //   LINFO("\n\n\n\nSum the Average...");
00385 //   double sumAngle = 0;
00386 //   std::vector<LineSegment2D>::iterator itr  = flines.begin();
00387 //   while(itr < flines.end())
00388 //     {
00389 //       float angle = (*itr).angle();
00390 //       LINFO("Line Angle[%4d,%4d][%4d,%4d]: %f",
00391 //             (*itr).point1().i, (*itr).point1().j,
00392 //             (*itr).point2().i, (*itr).point2().j,
00393 //             angle * 180/M_PI);
00394 
00395 //       // eliminate lines that are close to vertical
00396 //       if(!isnan(angle))
00397 //         {
00398 //           sumAngle += angle; itr++;
00399 //         }
00400 //       else
00401 //         { sumAngle += 90.0;
00402 //           //itr = flines.erase(itr); LEBUG("Drop a line");
00403 //         }
00404 //     }
00405 //   if (flines.size() == 0) return avgPipeAngle;
00406 
00407 //   // also drop off angles that do not fall
00408 //   // within a standard deviation of the average angle
00409 //   float meanAngle = sumAngle / flines.size();
00410 //   float stdAngleSum = 0.0;
00411 //   float stdAngleSqr = 0.0;
00412 //   for(uint i = 0; i < flines.size(); i++)
00413 //     {
00414 //       float angle = flines[i].angle();
00415 //       stdAngleSum = angle - meanAngle;
00416 //       stdAngleSqr += (stdAngleSum * stdAngleSum);
00417 //     }
00418 //   float stdevAngle = stdAngleSqr / flines.size();
00419 //   double stdtemp = sqrt((double)stdevAngle);
00420 //   stdevAngle = (float)stdtemp;
00421 
00422 //   // throw out angles that do not fit within stdev of the angle
00423 //   sumAngle = 0.0;  itr  = flines.begin();
00424 //   while(itr < flines.end())
00425 //     {
00426 //       float angle = (*itr).angle();
00427 //       if(angle >= (meanAngle - stdevAngle) &&
00428 //          angle <= (meanAngle + stdevAngle))
00429 //         { sumAngle += angle; itr++; }
00430 //       else itr = flines.erase(itr);
00431 //     }
00432 
00433 //   if (flines.size() != 0) sumAngle /= flines.size();  else return avgPipeAngle;
00434 //   if (sumAngle > 0) sumAngle -= M_PI;
00435 
00436 //   float avgAngle = sumAngle;
00437 
00438 //   // Smooth out the angles being displayed by averaging the position of the last 30 angles
00439 //   uint angleBuffSize = angleBuff.size();
00440 
00441 //   bool isWithinStdDev = false;
00442 //   // Check if angle is within std dev
00443 //   if(angleBuff.size() < 2
00444 //      || (fabs(avgAngle - avgPipeAngle) <= stdDevPipeAngle ))
00445 //     {
00446 //       stdDevAngleCount = 0;
00447 //       isWithinStdDev = true;
00448 //     }
00449 //   else
00450 //     {
00451 //       stdDevAngleCount++;
00452 //     }
00453 
00454 //   if(stdDevAngleCount > 5)
00455 //     {
00456 //       angleBuff.clear();
00457 //       isWithinStdDev = true;
00458 //     }
00459 
00460 //   if(isWithinStdDev)
00461 //     {
00462 //       if(angleBuffSize >= 30)
00463 //         {
00464 //           angleBuff.pop_front();
00465 //         }
00466 
00467 //       angleBuff.push_back(avgAngle);
00468 
00469 //       angleBuffSize = angleBuff.size();
00470 //     }
00471 
00472 //   if(angleBuffSize > 0)
00473 //     {
00474 //       float sumBuffAngle = 0.0;
00475 //       std::list<float>::iterator it;
00476 
00477 //       for(it = angleBuff.begin(); it != angleBuff.end(); ++it)
00478 //         {
00479 //           sumBuffAngle += (*it);
00480 //         }
00481 
00482 //       avgPipeAngle = sumBuffAngle / angleBuffSize;
00483 
00484 //       //      LINFO("AVG Pipe Angle: %f",avgPipeAngle);
00485 
00486 
00487 //       float sqrStdDevAngle = 0.0;
00488 //       for(it = angleBuff.begin(); it != angleBuff.end(); ++it)
00489 //       {
00490 //         sqrStdDevAngle = ((*it - avgPipeAngle) * (*it - avgPipeAngle));
00491 //       }
00492 
00493 //       float stdevAngle = sqrStdDevAngle / angleBuffSize;
00494 //       double stdTempAngle = sqrt((double)stdevAngle);
00495 //       stdDevPipeAngle = (float)stdTempAngle;
00496 //     }
00497 //   else
00498 //     {
00499 //       avgPipeAngle = avgAngle;
00500 //     }
00501 
00502 // //   LDEBUG("StdDev: %f", stdevAngle * 180/M_PI);
00503 // //   LDEBUG("Mean: %f", meanAngle * 180/M_PI);
00504 // //   LDEBUG("Pipe Direction: %f", ((float)sumAngle) * 180.0/M_PI);
00505 // //   LDEBUG("Line Count: %"ZU, lines.size());
00506 // //   LDEBUG("Adjusted Line Count: %"ZU, flines.size());
00507 //   return avgPipeAngle;
00508 }
00509 
00510 // ######################################################################
00511 float BeoSubPipe::getPipeDir
00512 (const std::vector<LineSegment2D> lines, const float thresh, uint &tot)
00513 {
00514   uint nLines = lines.size();
00515   if(nLines == 0) { LINFO("*** NO LINES ***"); return avgPipeAngle; }
00516 
00517   // get the mean
00518   float sumLength = 0.0;  float meanLength = 0.0;
00519   for (uint i = 0; i < nLines; i++) sumLength += lines[i].length();
00520   meanLength = sumLength / nLines;
00521 
00522   // get the standard deviation
00523   float sumsqr = 0.0;
00524   for (uint i = 0; i < nLines; i++)
00525     sumsqr += pow((float)(lines[i].length()) - meanLength, 2.0);
00526   float stdevLength = sqrt(sumsqr / (nLines-1));
00527 
00528   // kick out lines that aren't within the stddev of the mean length
00529   LINFO ("Mean: %f   StdDev: %f", meanLength, stdevLength);
00530   std::vector<LineSegment2D> nlines;
00531   for(uint i = 0; i < nLines; i++)
00532     {
00533       if (lines[i].length() > (meanLength - stdevLength) &&
00534           lines[i].length() < (meanLength + stdevLength)   )
00535         nlines.push_back(lines[i]);
00536     }
00537 
00538   if (nlines.size() == 0)
00539     { LINFO("*** NO LINE LENGTH WITHIN STD DEV ***"); return avgPipeAngle; }
00540 
00541   std::vector<LineSegment2D> flines = nlines;
00542   LINFO("\n\n\n\nSum the Average...");
00543   double sumAngle = 0;
00544   std::vector<LineSegment2D>::iterator itr  = flines.begin();
00545   while(itr < flines.end())
00546     {
00547       float angle = (*itr).angle();
00548       LINFO("Line Angle[%4d,%4d][%4d,%4d]: %f",
00549             (*itr).point1().i, (*itr).point1().j,
00550             (*itr).point2().i, (*itr).point2().j,
00551             angle * 180/M_PI);
00552 
00553       // eliminate lines that are close to vertical
00554       if(!isnan(angle))
00555         {
00556           sumAngle += angle; itr++;
00557         }
00558       else
00559         { sumAngle += 90.0;
00560           //itr = flines.erase(itr); LEBUG("Drop a line");
00561         }
00562     }
00563   if (flines.size() == 0) return avgPipeAngle;
00564 
00565   // also drop off angles that do not fall
00566   // within a standard deviation of the average angle
00567   float meanAngle = sumAngle / flines.size();
00568   float stdAngleSum = 0.0;
00569   float stdAngleSqr = 0.0;
00570   for(uint i = 0; i < flines.size(); i++)
00571     {
00572      float angle = flines[i].angle();
00573       stdAngleSum = angle - meanAngle;
00574       stdAngleSqr += (stdAngleSum * stdAngleSum);
00575     }
00576   float stdevAngle = stdAngleSqr / flines.size();
00577   double stdtemp = sqrt((double)stdevAngle);
00578   stdevAngle = (float)stdtemp;
00579 
00580   // throw out angles that do not fit within stdev of the angle
00581   sumAngle = 0.0;  itr  = flines.begin();
00582   while(itr < flines.end())
00583     {
00584       float angle = (*itr).angle();
00585       if(angle >= (meanAngle - stdevAngle) &&
00586          angle <= (meanAngle + stdevAngle))
00587         { sumAngle += angle; itr++; }
00588       else itr = flines.erase(itr);
00589     }
00590 
00591   if (flines.size() != 0) sumAngle /= flines.size();  else return avgPipeAngle;
00592   if (sumAngle > 0) sumAngle -= M_PI;
00593 
00594   float avgAngle = sumAngle;
00595 
00596   // Smooth out the angles being displayed by averaging the position of the last 30 angles
00597   uint angleBuffSize = angleBuff.size();
00598 
00599   bool isWithinStdDev = false;
00600   // Check if angle is within std dev
00601   if(angleBuff.size() < 2
00602      || (fabs(avgAngle - avgPipeAngle) <= stdDevPipeAngle ))
00603     {
00604       stdDevAngleCount = 0;
00605       isWithinStdDev = true;
00606     }
00607   else
00608     {
00609       stdDevAngleCount++;
00610     }
00611 
00612   if(stdDevAngleCount > 5)
00613     {
00614       angleBuff.clear();
00615       isWithinStdDev = true;
00616     }
00617 
00618   if(isWithinStdDev)
00619     {
00620       if(angleBuffSize >= 30)
00621         {
00622           angleBuff.pop_front();
00623         }
00624 
00625       angleBuff.push_back(avgAngle);
00626 
00627       angleBuffSize = angleBuff.size();
00628     }
00629 
00630   if(angleBuffSize > 0)
00631     {
00632       float sumBuffAngle = 0.0;
00633       std::list<float>::iterator it;
00634 
00635       for(it = angleBuff.begin(); it != angleBuff.end(); ++it)
00636         {
00637           sumBuffAngle += (*it);
00638         }
00639 
00640       avgPipeAngle = sumBuffAngle / angleBuffSize;
00641 
00642       //      LINFO("AVG Pipe Angle: %f",avgPipeAngle);
00643 
00644 
00645       float sqrStdDevAngle = 0.0;
00646       for(it = angleBuff.begin(); it != angleBuff.end(); ++it)
00647       {
00648         sqrStdDevAngle = ((*it - avgPipeAngle) * (*it - avgPipeAngle));
00649       }
00650 
00651       float stdevAngle = sqrStdDevAngle / angleBuffSize;
00652       double stdTempAngle = sqrt((double)stdevAngle);
00653       stdDevPipeAngle = (float)stdTempAngle;
00654     }
00655   else
00656     {
00657       avgPipeAngle = avgAngle;
00658     }
00659 
00660 //   LDEBUG("StdDev: %f", stdevAngle * 180/M_PI);
00661 //   LDEBUG("Mean: %f", meanAngle * 180/M_PI);
00662 //   LDEBUG("Pipe Direction: %f", ((float)sumAngle) * 180.0/M_PI);
00663 //   LDEBUG("Line Count: %"ZU, lines.size());
00664 //   LDEBUG("Adjusted Line Count: %"ZU, flines.size());
00665   return avgPipeAngle;
00666 }
00667 
00668 // ######################################################################
00669 /* So things look consistent in everyone's emacs... */
00670 /* Local Variables: */
00671 /* indent-tabs-mode: nil */
00672 /* End: */
Generated on Sun May 8 08:40:19 2011 for iLab Neuromorphic Vision Toolkit by  doxygen 1.6.3