ContourBoundaryDetector.C

Go to the documentation of this file.
00001 /*!@file Gist/ContourBoundaryDetector.C Detect meaningful contours by
00002    weighting longer countour chains more */
00003 // //////////////////////////////////////////////////////////////////// //
00004 // The iLab Neuromorphic Vision C++ Toolkit - Copyright (C) 2001 by the //
00005 // University of Southern California (USC) and the iLab at USC.         //
00006 // See http://iLab.usc.edu for information about this project.          //
00007 // //////////////////////////////////////////////////////////////////// //
00008 // Major portions of the iLab Neuromorphic Vision Toolkit are protected //
00009 // under the U.S. patent ``Computation of Intrinsic Perceptual Saliency //
00010 // in Visual Environments, and Applications'' by Christof Koch and      //
00011 // Laurent Itti, California Institute of Technology, 2001 (patent       //
00012 // pending; application number 09/912,225 filed July 23, 2001; see      //
00013 // http://pair.uspto.gov/cgi-bin/final/home.pl for current status).     //
00014 // //////////////////////////////////////////////////////////////////// //
00015 // This file is part of the iLab Neuromorphic Vision C++ Toolkit.       //
00016 //                                                                      //
00017 // The iLab Neuromorphic Vision C++ Toolkit is free software; you can   //
00018 // redistribute it and/or modify it under the terms of the GNU General  //
00019 // Public License as published by the Free Software Foundation; either  //
00020 // version 2 of the License, or (at your option) any later version.     //
00021 //                                                                      //
00022 // The iLab Neuromorphic Vision C++ Toolkit is distributed in the hope  //
00023 // that it will be useful, but WITHOUT ANY WARRANTY; without even the   //
00024 // implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR      //
00025 // PURPOSE.  See the GNU General Public License for more details.       //
00026 //                                                                      //
00027 // You should have received a copy of the GNU General Public License    //
00028 // along with the iLab Neuromorphic Vision C++ Toolkit; if not, write   //
00029 // to the Free Software Foundation, Inc., 59 Temple Place, Suite 330,   //
00030 // Boston, MA 02111-1307 USA.                                           //
00031 // //////////////////////////////////////////////////////////////////// //
00032 //
00033 // Primary maintainer for this file: Christian Siagian <siagian@usc.edu>
00034 // $HeadURL: svn://ilab.usc.edu/trunk/saliency/src/Robots/Beobot2/Navigation/FOE_Navigation/MiddleTemporal.C
00035 // $ $Id: $
00036 //
00037 //////////////////////////////////////////////////////////////////////////
00038 
00039 #include "Gist/ContourBoundaryDetector.H"
00040 #include "Image/ColorOps.H"
00041 #include "Image/MathOps.H"
00042 #include "Image/DrawOps.H"
00043 #include  <cstdio>
00044 
00045 // ######################################################################
00046 // ######################################################################
00047 
00048 // ######################################################################
00049 ContourBoundaryDetector::ContourBoundaryDetector()
00050 {
00051   itsWin.reset();
00052 }
00053 
00054 // ######################################################################
00055 ContourBoundaryDetector::~ContourBoundaryDetector()
00056 { }
00057 
00058 // ######################################################################
00059 void ContourBoundaryDetector::computeContourBoundary
00060 (Image<PixRGB<byte> > ima, int r)
00061 {
00062   itsImage = ima;
00063   if(r == -1)
00064     itsRadius = DEFAULT_NEIGHBORHOOD_RAD;
00065 
00066   // execute each step
00067   computeVarianceRidgeBoundaryMap();
00068   computeNmsBoundaryMap();
00069   computeContourBoundaryEdgels();
00070   computeContourBoundaryMap();
00071 }
00072 
00073 // ######################################################################
00074 void ContourBoundaryDetector::computeVarianceRidgeBoundaryMap()
00075 {
00076   // get the CIELab variance of the input image
00077   Image<float> varImg   = getLabStandardDeviation(itsImage, itsRadius);
00078 
00079   // get the gradient of the variance of the image
00080   std::vector<Image<float> > gradImg  = calculateGradient(varImg, itsRadius);
00081 
00082   // get the ridges in the gradient of the variance image
00083   Image<float> ridgeImg = getRidge(gradImg, itsRadius);
00084 
00085   // substract the gradient from the ridge image
00086   itsBoundaryImage  = substractGradImg(ridgeImg, gradImg);
00087 }
00088 
00089 // ######################################################################
00090 void ContourBoundaryDetector::computeNmsBoundaryMap()
00091 {
00092   // run non-max suppression to thin out image
00093   itsBoundaryImageNMS = getNonMaxSuppression(itsBoundaryImage);
00094   //bImgNMS = getNonMaxSuppression(bImgNMS);
00095 }
00096 
00097 // ######################################################################
00098 void ContourBoundaryDetector::computeContourBoundaryEdgels()
00099 {
00100   // compute the contour edgels
00101   itsEdgelBoundaryImage = getContourBoundaryEdgels();
00102 }
00103 
00104 // ######################################################################
00105 void ContourBoundaryDetector::computeContourBoundaryMap()
00106 {
00107   // connect the contour edgels
00108   connectContourEdgels();
00109 
00110   itsContourBoundaryImage = getContourBoundaryImage();
00111 
00112   //displayContourBoundary();
00113 }
00114 
00115 // ######################################################################
00116 Image<float> ContourBoundaryDetector::getVarianceRidgeBoundaryMap
00117 (Image<PixRGB<byte> > ima, int r)
00118 {
00119   itsImage = ima;
00120   if(r == -1)
00121     itsRadius = DEFAULT_NEIGHBORHOOD_RAD;
00122 
00123   computeVarianceRidgeBoundaryMap();
00124 
00125   return itsBoundaryImage;
00126 }
00127 
00128 // ######################################################################
00129 Image<float> ContourBoundaryDetector::getVarianceRidgeBoundaryMap()
00130 {
00131   return itsBoundaryImage;
00132 }
00133 
00134 // ######################################################################
00135 Image<float> ContourBoundaryDetector::getNmsBoundaryMap
00136 (Image<PixRGB<byte> > ima, int r)
00137 {
00138   // assume that we need to recompute everything
00139   itsImage = ima;
00140   if(r == -1)
00141     itsRadius = DEFAULT_NEIGHBORHOOD_RAD;
00142 
00143   computeVarianceRidgeBoundaryMap();
00144   computeNmsBoundaryMap();
00145 
00146   return itsBoundaryImageNMS;
00147 }
00148 
00149 // ######################################################################
00150 Image<float> ContourBoundaryDetector::getNmsBoundaryMap()
00151 {
00152   return itsBoundaryImageNMS;
00153 }
00154 
00155 // ######################################################################
00156 Image<float> ContourBoundaryDetector::getEdgelBoundaryMap
00157 (Image<PixRGB<byte> > ima, int r)
00158 {
00159   // assume that we need to recompute everything
00160   itsImage = ima;
00161   if(r == -1)
00162     itsRadius = DEFAULT_NEIGHBORHOOD_RAD;
00163 
00164   computeVarianceRidgeBoundaryMap();
00165   computeNmsBoundaryMap();
00166   computeContourBoundaryEdgels();
00167 
00168   return itsEdgelBoundaryImage;
00169 }
00170 
00171 // ######################################################################
00172 Image<float> ContourBoundaryDetector::getEdgelBoundaryMap()
00173 {
00174   return itsEdgelBoundaryImage;
00175 }
00176 
00177 // ######################################################################
00178 Image<PixRGB<byte> > ContourBoundaryDetector::getContourBoundaryMap
00179 (Image<PixRGB<byte> > ima, int r)
00180 {
00181   // assume that we need to recompute everything
00182   itsImage = ima;
00183   if(r == -1)
00184     itsRadius = DEFAULT_NEIGHBORHOOD_RAD;
00185 
00186   computeVarianceRidgeBoundaryMap();
00187   computeNmsBoundaryMap();
00188   computeContourBoundaryEdgels();
00189   computeContourBoundaryMap();
00190 
00191   return itsContourBoundaryImage;
00192 }
00193 
00194 // ######################################################################
00195 Image<PixRGB<byte> > ContourBoundaryDetector::getContourBoundaryMap()
00196 {
00197   return itsContourBoundaryImage;
00198 }
00199 
00200 // ######################################################################
00201 Image<float> ContourBoundaryDetector::getLabStandardDeviation
00202 (Image<PixRGB<byte> > ima, int r)
00203 {
00204   // convert to CIElab color space
00205   Image<float> lImg; 
00206   Image<float> aImg;
00207   Image<float> bImg;
00208   getLAB(ima, lImg, aImg, bImg);
00209 
00210   // compute squares of each channel
00211   Image<float> lSqImg = lImg * lImg; 
00212   Image<float> aSqImg = aImg * aImg;
00213   Image<float> bSqImg = bImg * bImg;
00214 
00215   // compute box blur for each channel
00216   Image<float> lBbImg = boxBlur(lImg, r);
00217   Image<float> aBbImg = boxBlur(aImg, r); 
00218   Image<float> bBbImg = boxBlur(bImg, r); 
00219 
00220   Image<float> lSqBbImg = boxBlur(lSqImg, r);
00221   Image<float> aSqBbImg = boxBlur(aSqImg, r); 
00222   Image<float> bSqBbImg = boxBlur(bSqImg, r); 
00223 
00224   //calculate standard deviation of each l a b channel
00225   Image<float> vl = lBbImg*lBbImg - lSqBbImg;
00226   Image<float> va = aBbImg*aBbImg - aSqBbImg;
00227   Image<float> vb = bBbImg*bBbImg - bSqBbImg;
00228   Image<float> varImg = sqCombine(vl, va, vb);
00229 
00230   return varImg;
00231 }
00232 
00233 // ######################################################################
00234 Image<float> ContourBoundaryDetector::boxBlur(Image<float> img, int rad)
00235 {
00236   int w = img.getWidth();
00237   int h = img.getHeight();
00238 
00239   Image<float> ret(w,h,NO_INIT);
00240   
00241   int r = rad;
00242 
00243   // go through each location
00244   for(int i = 0; i < w; i++)
00245     {
00246       for(int j = 0; j < h; j++)
00247         {
00248           float sum   = 0.0;
00249           uint  count = 0; 
00250 
00251           for(int di = -r; di <= r; di++)
00252             {
00253               for(int dj = -r; dj <= r; dj++)
00254                 {
00255                   int ii = i + di;
00256                   int jj = j + dj;
00257 
00258                   if(ii <  0) ii = -ii;
00259                   if(jj <  0) jj = -jj;
00260                   if(ii >= w) ii = 2*w - 2 - ii;
00261                   if(jj >= h) jj = 2*h - 2 - jj;
00262 
00263                   //if(img.coordsOk(ii,jj))
00264                   //  {
00265                       sum += img.getVal(ii,jj);
00266                       count++;
00267                   //  }
00268                 }
00269             }
00270           ret.setVal(i,j, sum/count);
00271         }
00272     }
00273 
00274   return ret;
00275 }
00276 
00277 // ######################################################################
00278 Image<float> ContourBoundaryDetector::sqCombine
00279 (Image<float> a, Image<float> b, Image<float> c)
00280 {
00281   ASSERT(a.isSameSize(b));
00282   ASSERT(b.isSameSize(c));
00283 
00284   int w = a.getWidth();
00285   int h = a.getHeight();
00286   Image<float> ret(w,h,NO_INIT);
00287 
00288   Image<float>::const_iterator aitr = a.begin(), stop = a.end();
00289   Image<float>::const_iterator bitr = b.begin();
00290   Image<float>::const_iterator citr = c.begin();
00291   Image<float>::iterator ritr = ret.beginw();
00292 
00293   while (aitr != stop)
00294     {
00295       float av = *aitr++;
00296       float bv = *bitr++;
00297       float cv = *citr++;
00298 
00299       *ritr++ = pow(av*av + bv*bv + cv*cv, 0.5);        
00300     }
00301  
00302   return ret;
00303 }
00304 
00305 // ######################################################################
00306 std::vector<Image<float> > ContourBoundaryDetector::calculateGradient
00307 (Image<float> varImg, int r)
00308 {  
00309   int w = varImg.getWidth();
00310   int h = varImg.getHeight();
00311 
00312   Image<float> gradImgX(w,h,NO_INIT);
00313   Image<float> gradImgY(w,h,NO_INIT);
00314 
00315   // smooth the variance Image using BoxBlur
00316   Image<float> varBbImg = boxBlur(varImg, r);
00317 
00318   std::vector<float> dx(NUM_GRADIENT_DIRECTIONS);
00319   std::vector<float> dy(NUM_GRADIENT_DIRECTIONS);    
00320   for(uint k = 0; k < NUM_GRADIENT_DIRECTIONS; k++)
00321     {
00322       dx[k] = cos(k*2*M_PI/NUM_GRADIENT_DIRECTIONS); 
00323       dy[k] = sin(k*2*M_PI/NUM_GRADIENT_DIRECTIONS);
00324       LDEBUG("%d %f %f", k, dx[k], dy[k]);
00325     }
00326 
00327   // go through each location
00328   for(int i = 0; i < w; i++)
00329     {
00330       for(int j = 0; j < h; j++)
00331         {
00332           float sumX = 0.0;
00333           float sumY = 0.0;
00334           for(uint k = 0; k < NUM_GRADIENT_DIRECTIONS; k++)
00335             {
00336               int i1 = i + r*dx[k];
00337               int j1 = j + r*dy[k];
00338 
00339               int i2 = i - r*dx[k];             
00340               int j2 = j - r*dy[k];
00341 
00342               if(i1 <  0) i1 = -i1;
00343               if(j1 <  0) j1 = -j1;
00344               if(i1 >= w) i1 = 2*w - 2 - i1;
00345               if(j1 >= h) j1 = 2*h - 2 - j1;
00346               
00347               if(i2 <  0) i2 = -i2;
00348               if(j2 <  0) j2 = -j2;
00349               if(i2 >= w) i2 = 2*w - 2 - i2;
00350               if(j2 >= h) j2 = 2*h - 2 - j2;
00351               
00352               float val = varBbImg.getVal(i1,j1) - varBbImg.getVal(i2,j2);
00353 
00354               sumX +=  val * dx[k];
00355               sumY +=  val * dy[k]; 
00356 
00357               // LINFO("%d %d %d| val: %f, (%f %f) %f %f", 
00358               //            i,j,k, val, val*dx[k], val*dy[k], sumX, sumY);
00359             }
00360 
00361           gradImgX.setVal(i,j, sumX);
00362           gradImgY.setVal(i,j, sumY);
00363         }
00364     }
00365 
00366   std::vector<Image<float> > gradImg(2);
00367   gradImg[0] = gradImgX;
00368   gradImg[1] = gradImgY;
00369 
00370   return gradImg;
00371 }
00372 
00373 // ######################################################################
00374 Image<float> ContourBoundaryDetector::getRidge
00375 (std::vector<Image<float> > gradImg, int r)
00376 {
00377   Image<float> gradImgX = gradImg[0];
00378   Image<float> gradImgY = gradImg[1];
00379 
00380   int w = gradImg[0].getWidth();
00381   int h = gradImg[0].getHeight();
00382   Image<float> ridgeImg(w,h,NO_INIT);
00383 
00384   std::vector<float> dx(NUM_GRADIENT_DIRECTIONS);
00385   std::vector<float> dy(NUM_GRADIENT_DIRECTIONS);    
00386   for(uint k = 0; k < NUM_GRADIENT_DIRECTIONS; k++)
00387     {
00388       dx[k] = cos(k*2*M_PI/NUM_GRADIENT_DIRECTIONS); 
00389       dy[k] = sin(k*2*M_PI/NUM_GRADIENT_DIRECTIONS);
00390     }
00391 
00392   // threshold the gradient image
00393   std::vector<std::vector<Image<float> > > dVin(NUM_GRADIENT_DIRECTIONS); 
00394   for(uint k = 0; k < NUM_GRADIENT_DIRECTIONS; k++)
00395     dVin[k] = std::vector<Image<float> >(2);
00396 
00397   for(uint k = 0; k < NUM_GRADIENT_DIRECTIONS; k++)
00398     {
00399       dVin[k][0] = Image<float>(w,h,NO_INIT);
00400       dVin[k][1] = Image<float>(w,h,NO_INIT);
00401 
00402       for(int i = 0; i < w; i++)
00403         {
00404           for(int j = 0; j < h; j++)
00405             {
00406               float x = 0.0;
00407               float y = 0.0;
00408 
00409               int ii = int(i + r*dx[k]);
00410               int jj = int(j + r*dy[k]); 
00411 
00412               if(ii <  0) ii = -ii;
00413               if(jj <  0) jj = -jj;
00414               if(ii >= w) ii = 2*w - 2 - ii;
00415               if(jj >= h) jj = 2*h - 2 - jj;
00416               
00417               float vX = gradImgX.getVal(ii,jj); 
00418               float vY = gradImgY.getVal(ii,jj);
00419               if((vX*dx[k] + vY*dy[k]) < 0.0)
00420                 {
00421                   x = vX;
00422                   y = vY;
00423                 }
00424               
00425               dVin[k][0].setVal(i,j, x);
00426               dVin[k][1].setVal(i,j, y);
00427               
00428               // if(i == 100 && j == 100)
00429               //        LINFO("[%3d %3d %d %d]: %3d %3d: "
00430               //              "%10.5f %10.5f %10.5f %10.5f %10.5f %10.5f", 
00431               //              i,j,k, r, ii, jj, 
00432               //              vX, vY, dx[k], dy[k], x, y);
00433             }
00434         }
00435     }
00436   itsDVin = dVin;
00437 
00438   std::vector<Image<float> > rDir (NUM_GRADIENT_DIRECTIONS); 
00439   Image<float> rDirIndex(w,h,NO_INIT);
00440 
00441   // calculate the geometric and arithmetic ridge direction
00442   // and sum the two together
00443   for(uint k = 0; k < NUM_GRADIENT_DIRECTIONS; k++)
00444     {
00445       rDir[k] = Image<float>(w,h,NO_INIT);
00446 
00447       for(int i = 0; i < w; i++)
00448         {
00449           for(int j = 0; j < h; j++)
00450             {
00451               float x1 = dVin[k][0].getVal(i,j);
00452               float y1 = dVin[k][1].getVal(i,j);
00453 
00454               uint k2 = k + NUM_RIDGE_DIRECTIONS;
00455               if(k >= NUM_RIDGE_DIRECTIONS) 
00456                 k2 = k - NUM_RIDGE_DIRECTIONS/2;
00457 
00458               float x2 = dVin[k2][0].getVal(i,j);
00459               float y2 = dVin[k2][1].getVal(i,j);
00460 
00461               float gVal = -(x1*x2 + y1*y2);                  
00462               if(gVal < 0.0) 
00463                 gVal = 0.0;
00464               else 
00465                 gVal = pow(gVal, 0.5);
00466 
00467               // float ax = x2 - x1;
00468               // float ay = y2 - y1;
00469               // float aVal =  pow(ax*ax + ay*ay, 0.5)/ 2.0;          
00470               // rDir[k].setVal(i,j, gVal + aVal);
00471 
00472               rDir[k].setVal(i,j, gVal);
00473 
00474               // if(i == 100 && j == 100)
00475               //        LINFO("[%3d %3d %d %d %d]: "
00476               //              "%10.5f %10.5f %10.5f %10.5f || "
00477               //              "%10.5f %10.5f = %10.5f", 
00478               //              i,j,k, k2, r, x1, y1, x2, y2, 
00479               //              gVal, aVal, gVal+aVal);
00480             }
00481         }
00482     }
00483 
00484   itsRidgeDirection = rDir;
00485 
00486   std::vector<Image<float> > rDirMax(NUM_RIDGE_DIRECTIONS); 
00487   for(uint k = 0; k < NUM_RIDGE_DIRECTIONS; k++)    
00488     rDirMax[k] = Image<float>(w,h,ZEROS);
00489 
00490   // get the maximum among the directions
00491   for(int i = 0; i < w; i++)
00492     {
00493       for(int j = 0; j < h; j++)
00494         {
00495           float max = 0.0; int maxk = -1;
00496           for(uint k = 0; k < NUM_RIDGE_DIRECTIONS; k++)
00497             {
00498               float val = rDir[k].getVal(i,j);
00499               if(val > max) { max = val; maxk = k; }
00500             }
00501           ridgeImg.setVal(i,j, max);
00502           rDirIndex.setVal(i,j, maxk);
00503           if(maxk != -1)
00504             rDirMax[maxk].setVal(i,j, max);
00505         }
00506     }
00507 
00508   itsRidgeDirectionIndex = rDirIndex;
00509   itsRidgeDirectionMax   = rDirMax;
00510 
00511   return ridgeImg;
00512 }
00513 
00514 // ######################################################################
00515 Image<float> ContourBoundaryDetector::substractGradImg
00516 (Image<float> ridgeImg, std::vector<Image<float> >  gradImg)
00517 {
00518   Image<float> gradImgX = gradImg[0];
00519   Image<float> gradImgY = gradImg[1];
00520 
00521   Image<float> res = 
00522     //clampedDiff(ridgeImg, sqrt(gradImgX*gradImgX + gradImgY*gradImgY));
00523     ridgeImg - sqrt(gradImgX*gradImgX + gradImgY*gradImgY);
00524 
00525   return res;
00526 }
00527 
00528 // ######################################################################
00529 Image<float> ContourBoundaryDetector::getNonMaxSuppression
00530 (Image<float> bImg)
00531 {
00532   int w = bImg.getWidth();
00533   int h = bImg.getHeight();
00534   Image<float> bImgNMS(w,h,ZEROS);
00535 
00536   // LINFO("ridge image");
00537   // int lb = 152, rb = lb+16, tb = 64, bb = tb+16; 
00538   // for(int j = tb; j < bb; j++)
00539   //   {
00540   //     for(int i = lb; i < rb; i++)
00541   //    {
00542   //      printf("%10.5f ", bImg.getVal(i,j));
00543   //    }
00544   //     printf("\n");
00545   //   }
00546 
00547   // set up kernel for non-max suppression
00548   int wSize = BOUNDARY_STEP_SIZE+1;
00549   std::vector<std::vector<Point2D<int> > > sCoordsL(NUM_RIDGE_DIRECTIONS);
00550   std::vector<std::vector<Point2D<int> > > sCoordsR(NUM_RIDGE_DIRECTIONS);
00551   std::vector<std::vector<Point2D<int> > > cCoords (NUM_RIDGE_DIRECTIONS);
00552 
00553 
00554   for(uint k = 0; k < NUM_RIDGE_DIRECTIONS; k++)
00555     {
00556       cCoords [k] = std::vector<Point2D<int> >();  
00557       sCoordsL[k] = std::vector<Point2D<int> >();  
00558       sCoordsR[k] = std::vector<Point2D<int> >();  
00559     }
00560 
00561   // Non-max suppressed values for each direction
00562   itsRidgeDirectionNMS.clear();
00563   for(uint k = 0; k < NUM_RIDGE_DIRECTIONS; k++)
00564     itsRidgeDirectionNMS.push_back( Image<float>(w,h,ZEROS));  
00565 
00566   for(int i = -wSize/2; i <= wSize/2; i++)
00567     {
00568       for(int j = -wSize/2; j <= wSize/2; j++)
00569         {
00570           if( i == 0) cCoords [0].push_back(Point2D<int>(i,j));
00571           if( i <  0) sCoordsL[0].push_back(Point2D<int>(i,j));
00572           if( i >  0) sCoordsR[0].push_back(Point2D<int>(i,j));
00573 
00574           if(-j == i) cCoords [1].push_back(Point2D<int>(i,j));
00575           if(-j >  i) sCoordsL[1].push_back(Point2D<int>(i,j));
00576           if( i > -j) sCoordsR[1].push_back(Point2D<int>(i,j));   
00577 
00578           if( j == 0) cCoords [2].push_back(Point2D<int>(i,j));
00579           if( j <  0) sCoordsL[2].push_back(Point2D<int>(i,j));
00580           if( j >  0) sCoordsR[2].push_back(Point2D<int>(i,j));
00581           
00582           if( i == j) cCoords [3].push_back(Point2D<int>(i,j));
00583           if( j >  i) sCoordsL[3].push_back(Point2D<int>(i,j));
00584           if( i >  j) sCoordsR[3].push_back(Point2D<int>(i,j));
00585         }
00586     }
00587 
00588   // go through each point
00589   for(int i = 0; i < w; i++)
00590     {
00591       for(int j = 0; j < h; j++)
00592         {
00593           // get the value
00594           float val = bImg.getVal(i,j);
00595           Point2D<int> cpt(i,j);
00596           
00597           // if(i == 160 && j == 70)
00598           //   LINFO("%3d %3d: %f", i, j, val);
00599 
00600           for(uint k = 0; k < NUM_RIDGE_DIRECTIONS; k++)
00601             {
00602               float totalC = 0.0; uint ctC = 0;
00603               for(uint cc = 0; cc < cCoords[k].size(); cc++)
00604                 {
00605                   Point2D<int> pt = cCoords[k][cc] + cpt;  
00606                   if(bImg.coordsOk(pt)) 
00607                     {
00608                       totalC += bImg.getVal(pt);
00609                       ctC++;
00610                     }
00611                   // if(i == 160 && j == 70 && k == 1)
00612                   //   {
00613                   //     LINFO("center: %3d %3d: %f", pt.i, pt.j, bImg.getVal(pt));
00614                   //   }
00615                 }
00616               
00617               float totalL = 0.0; uint ctL = 0;
00618               for(uint cl = 0; cl < sCoordsL[k].size(); cl++)
00619                 {
00620                   Point2D<int> pt = sCoordsL[k][cl] + cpt;  
00621                   if(bImg.coordsOk(pt)) 
00622                     {
00623                       totalL += bImg.getVal(pt); ctL++;
00624                     }
00625                   // if(i == 160 && j == 70 && k == 1)
00626                   //   {
00627                   //     LINFO("left: %3d %3d: %f", pt.i, pt.j, bImg.getVal(pt));
00628                   //   }
00629                 }
00630               
00631               float totalR = 0.0; uint ctR = 0;
00632               for(uint cr = 0; cr < sCoordsR[k].size(); cr++)
00633                 {
00634                   Point2D<int> pt = sCoordsR[k][cr] + cpt;  
00635                   if(bImg.coordsOk(pt)) 
00636                     {
00637                       totalR += bImg.getVal(pt); ctR++;
00638                     }
00639                   // if(i == 160 && j == 70 && k == 1)
00640                   //   {
00641                   //     LINFO("right: %3d %3d: %f", pt.i, pt.j, bImg.getVal(pt));
00642                   //   }
00643                 }
00644               
00645               if(totalC/ctC > totalR/ctR && totalC/ctC > totalL/ctL 
00646                  && val > 0.0)  
00647                 {
00648                   bImgNMS.setVal(i,j, val);
00649 
00650                   // if(i == 160 && j == 70)
00651                   //   {
00652                   //     LINFO("K: %d in: %f %f %f | %d %d %d| %f %f %f --> %f", 
00653                   //        k, totalC, totalR, totalL,
00654                   //        ctC, ctR, ctL,
00655                   //        totalC/ctC, totalR/ctR, totalL/ctL, 
00656                   //        totalC/ctC*2 - totalR/ctR - totalL/ctL);
00657                   //   }
00658                   itsRidgeDirectionNMS[k].setVal
00659                     (i,j, totalC/ctC*2 - totalR/ctR - totalL/ctL);
00660                 }
00661             }
00662         }
00663     }
00664 
00665   return bImgNMS;
00666 }
00667 
00668 // ######################################################################
00669 Image<float> ContourBoundaryDetector::getContourBoundaryEdgels()
00670 {
00671   // NOTE: FIXXXX: TENSOR-VOTING IS PROBABLY A BETTER IDEA
00672 
00673   int w = itsImage.getWidth();
00674   int h = itsImage.getHeight();
00675   Image<float> edgelBoundaryImage(w,h,ZEROS);
00676 
00677   // LINFO("ridge NMS image");
00678   // int lb = 152, rb = lb+16, tb = 64, bb = tb+16; 
00679   // for(int j = tb; j < bb; j++)
00680   //   {
00681   //     for(int i = lb; i < rb; i++)
00682   //    {
00683   //      printf("%10.5f ", itsBoundaryImageNMS.getVal(i,j));
00684   //    }
00685   //     printf("\n");
00686   //   }
00687 
00688   // LINFO("dir ridge NMS ");
00689   
00690   // for(uint k = 0; k < NUM_RIDGE_DIRECTIONS; k++)
00691   //   {
00692   //     LINFO("dir: %d", k);
00693   //     for(int j = tb; j < bb; j++)
00694   //    {
00695   //      for(int i = lb; i < rb; i++)
00696   //        {
00697   //          printf("%10.5f ", itsRidgeDirectionNMS[k].getVal(i,j));
00698   //        }
00699   //      printf("\n");
00700   //    }
00701   //   }
00702 
00703   int step  = BOUNDARY_STEP_SIZE;
00704   //int wSize = BOUNDARY_STEP_SIZE+1;
00705 
00706   // set up the center and surround opponency locations
00707   std::vector<std::vector<Point2D<int> > > cCoords(NUM_RIDGE_DIRECTIONS);
00708   std::vector<std::vector<Point2D<int> > > sCoords(NUM_RIDGE_DIRECTIONS);
00709   std::vector<float> angles;
00710 
00711   for(uint k = 0; k < NUM_RIDGE_DIRECTIONS; k++)
00712     {
00713       cCoords[k] = std::vector<Point2D<int> >();  
00714       sCoords[k] = std::vector<Point2D<int> >();  
00715 
00716       angles.push_back(k*2.0*M_PI/float(NUM_RIDGE_DIRECTIONS));
00717     }
00718 
00719   // fill the center coordinates
00720   for(int i = -step/2; i < step/2; i++)
00721     {
00722       cCoords[0].push_back(Point2D<int>( i, 0));
00723       cCoords[1].push_back(Point2D<int>( i, i));
00724       cCoords[2].push_back(Point2D<int>( 0, i));
00725       cCoords[3].push_back(Point2D<int>( i,-i));
00726     }
00727 
00728   int hstep = step/2;
00729   // fill the surround coordinates (bottom or right)
00730   for(int i = 0; i < hstep; i++)
00731     {
00732       sCoords[0].push_back(Point2D<int>( i+hstep,       0));
00733       sCoords[0].push_back(Point2D<int>( i+hstep,   hstep));
00734       sCoords[0].push_back(Point2D<int>( i+hstep,  -hstep));
00735 
00736       sCoords[1].push_back(Point2D<int>( i+hstep, i+hstep));
00737       sCoords[1].push_back(Point2D<int>( i+step , i      ));
00738       sCoords[1].push_back(Point2D<int>( i      ,-i-step ));
00739 
00740       sCoords[2].push_back(Point2D<int>(       0, i+hstep));
00741       sCoords[2].push_back(Point2D<int>(   hstep, i+hstep));
00742       sCoords[2].push_back(Point2D<int>(  -hstep, i+hstep));
00743 
00744       sCoords[3].push_back(Point2D<int>( i+hstep,-i-hstep));
00745       sCoords[3].push_back(Point2D<int>( i      ,-i-step ));
00746       sCoords[3].push_back(Point2D<int>( i+step ,-i      ));
00747     }
00748 
00749   // fill the surround coordinates (top or right)
00750   for(int i = -hstep; i < 0; i++)
00751     {
00752       sCoords[0].push_back(Point2D<int>( i-hstep,       0));
00753       sCoords[0].push_back(Point2D<int>( i-hstep,   hstep));
00754       sCoords[0].push_back(Point2D<int>( i-hstep,  -hstep));
00755 
00756       sCoords[1].push_back(Point2D<int>( i-hstep, i-hstep));
00757       sCoords[1].push_back(Point2D<int>( i-step , i      ));
00758       sCoords[1].push_back(Point2D<int>( i      , i-step ));
00759 
00760       sCoords[2].push_back(Point2D<int>(       0, i-hstep));
00761       sCoords[2].push_back(Point2D<int>(   hstep, i-hstep));
00762       sCoords[2].push_back(Point2D<int>(  -hstep, i-hstep));
00763 
00764       sCoords[3].push_back(Point2D<int>( i-hstep,-i+hstep));
00765       sCoords[3].push_back(Point2D<int>( i      ,-i+step ));
00766       sCoords[3].push_back(Point2D<int>( i-step ,-i      ));
00767     }
00768 
00769   // for(uint k = 0; k < NUM_RIDGE_DIRECTIONS; k++)
00770   //   {
00771   //     LINFO("c Dir: %d", k);
00772   //     for(uint i = 0; i <  cCoords[k].size(); i++)
00773   //    {
00774   //      LINFO("%3d: %5d %5d", i, cCoords[k][i].i, cCoords[k][i].j);
00775   //    }
00776   //     LINFO("s ");
00777   //     for(uint i = 0; i <  sCoords[k].size(); i++)
00778   //    {
00779   //      LINFO("%3d: %5d %5d", i, sCoords[k][i].i, sCoords[k][i].j);
00780   //      if(i%3 == 2) LINFO(" ");
00781   //    }
00782   //     Raster::waitForKey();
00783   //   }
00784   // Raster::waitForKey();
00785 
00786   // reset the edgel storage
00787   // NOTE: we will keep edgel at index 0 empty
00788   int wEdgel = w/step; 
00789   int hEdgel = h/step;
00790   itsEdgels = 
00791     Image<std::vector<rutz::shared_ptr<Edgel> > >(wEdgel, hEdgel, ZEROS);
00792   
00793   // go through each point
00794   // with the specified step size
00795   int wLimit = (w/step)*step;
00796   int hLimit = (h/step)*step;
00797   for(int j = step; j < hLimit; j+= step)
00798     {
00799       for(int i = step; i < wLimit; i+= step)
00800         {
00801           Point2D<int> cpt(i,j);
00802           
00803           int maxk = -1;
00804           Point2D<int> maxPt(-1,-1);
00805           float maxVal = -1.0F;
00806 
00807           // for each direction
00808           for(uint k = 0; k < NUM_RIDGE_DIRECTIONS; k++)
00809             {
00810               Point2D<int> maxCKpt(-1,-1);
00811               float maxCKval = 0.0;
00812 
00813               // get maximum contour value for the center size
00814               // to make the contour detector phase invariant
00815               for(uint ci = 0; ci < cCoords[k].size(); ci++)
00816                 {
00817                   Point2D<int> pt = cCoords[k][ci] + cpt;  
00818                   if(edgelBoundaryImage.coordsOk(pt)) 
00819                     {
00820                       float val = itsRidgeDirectionNMS[k].getVal(pt); 
00821                       if(maxCKval < val) 
00822                         {
00823                           maxCKval = val; maxCKpt = pt;
00824                         }
00825                     }
00826                 }
00827 
00828               float maxSKval = 0.0;
00829 
00830               // get the maximum value for the surround 
00831               for(uint si = 0; si < sCoords[k].size(); si++)
00832                 {
00833                   Point2D<int> pt = sCoords[k][si] + cpt;  
00834                   if(edgelBoundaryImage.coordsOk(pt)) 
00835                     {
00836                       float val = itsRidgeDirectionNMS[k].getVal(pt); 
00837                       if(maxSKval < val) maxSKval = val;
00838                     }
00839                 }
00840               
00841               // if center > 0 and wins
00842               if(maxCKval > 0.0F && maxCKval > maxSKval)
00843                 {
00844                   if(maxCKval > maxVal) 
00845                     {
00846                       maxPt  = maxCKpt;
00847                       maxVal = maxCKval;
00848                       maxk   = k;
00849                     }
00850                   
00851                   uint iEdgel = i/step;
00852                   uint jEdgel = j/step;
00853 
00854                   // put the new edgel in the right position
00855                   rutz::shared_ptr<Edgel> 
00856                     edgel(new Edgel(maxCKpt, angles[k], k, maxCKval)); 
00857                   
00858                   std::vector<rutz::shared_ptr<Edgel> >
00859                     cEdgelList = itsEdgels.getVal(iEdgel,jEdgel);
00860 
00861                   uint eNum = cEdgelList.size();
00862                   uint place = 0;
00863                   for(uint ce = 0; ce < eNum; ce++)
00864                     {
00865                       if(cEdgelList[ce]->val < maxCKval)
00866                         {
00867                           place = ce; ce = eNum;
00868                         }
00869                       else place = ce+1;
00870                     }
00871                   //LINFO("place: %d  | eNum: %d", place, eNum);
00872                   
00873                   cEdgelList.push_back(edgel);
00874                   
00875                   // for(uint ce = 0; ce < cEdgelList.size(); ce++)
00876                   //   LINFO("%12.5f %3d", cEdgelList[ce]->val,  
00877                   //         cEdgelList[ce]->angleIndex);
00878 
00879                   if(place != eNum)
00880                     {
00881                       // LINFO("move one");
00882                   
00883                       for(int ce = int(eNum-1); ce >= int(place); ce--)
00884                         {
00885                           cEdgelList[ce+1] = cEdgelList[ce];
00886                         }
00887                       
00888                       // for(uint ce = 0; ce < cEdgelList.size(); ce++)
00889                       //   LINFO("%12.5f %3d", cEdgelList[ce]->val,  
00890                       //         cEdgelList[ce]->angleIndex);
00891                       
00892                       cEdgelList[place] = edgel;                  
00893                       
00894                       // LINFO("place the new one properly");
00895 
00896                       // for(uint ce = 0; ce < cEdgelList.size(); ce++)
00897                       //   LINFO("%12.5f %3d", cEdgelList[ce]->val,  
00898                       //         cEdgelList[ce]->angleIndex);
00899                     }
00900                   // else LINFO("last place");
00901 
00902                   itsEdgels.setVal(iEdgel,jEdgel, cEdgelList);
00903 
00904 
00905                   // LINFO("%d %d: size: %d ", i,j, 
00906                   //       int(itsEdgels.getVal(iEdgel,jEdgel).size()));
00907 
00908                   // for(uint ce = 0; ce < cEdgelList.size(); ce++)
00909                   //   LINFO("%12.5f %3d", cEdgelList[ce]->val,  
00910                   //         cEdgelList[ce]->angleIndex);
00911                 }
00912             }
00913 
00914           // if there is a winner
00915           if(maxk != -1)
00916             {
00917               float borderK = 
00918                 fmod((maxk+(NUM_RIDGE_DIRECTIONS/2)),NUM_RIDGE_DIRECTIONS);
00919 
00920               float dx = cos(borderK * M_PI/4.0) * hstep; 
00921               float dy = sin(borderK * M_PI/4.0) * hstep;
00922 
00923               Point2D<int> pt = maxPt;
00924               Point2D<int> p1 = pt + Point2D<int>( dx+.5,  dy+.5); 
00925               Point2D<int> p2 = pt + Point2D<int>(-dx-.5, -dy-.5);
00926 
00927               //uint iEdgel = i/step;
00928               //uint jEdgel = j/step;
00929               //if(iEdgel >= 10 && iEdgel <= 25 && jEdgel >= 1 && jEdgel <= 14)
00930               //   {
00931               //     LINFO("maxk: %d -> %d  -> %f  %f (%f %f %f %f)  |%f %f", 
00932               //           maxk, (maxk+(NUM_RIDGE_DIRECTIONS/2)), 
00933               //           (borderK*M_PI)/float(NUM_RIDGE_DIRECTIONS), borderK, 
00934               //           cos(0), cos(M_PI/4.0), cos(M_PI/2.0), cos(M_PI*.75),
00935               //           dx, dy);
00936               
00937 
00938 
00939               //LINFO("%d %d | %d %d | %d %d::::: %d %d  %d", 
00940               //         pt.i, pt.j, p1.i, p1.j, p2.i, p2.j, iEdgel, jEdgel, maxk);           
00941                   
00942               // draw the straightline contour in the image 
00943               // for visualization
00944               drawLine(edgelBoundaryImage, p1, p2, 255.0F);
00945               //drawDisk(edgelBoundaryImage, pt, 2,  255.0F);
00946               // }      
00947               
00948             }
00949         }
00950     }
00951 
00952   return edgelBoundaryImage;
00953 }
00954 
00955 // ######################################################################
00956 void ContourBoundaryDetector::connectContourEdgels()
00957 {
00958   int w = itsImage.getWidth();
00959   int h = itsImage.getHeight();
00960 
00961   // if contour map don't have a contour on the location
00962   // that means the edgel on that location 
00963   // do not belong to any contours yet
00964 
00965   // reset the contours storage
00966   // NOTE: we will keep edgel at index 0 empty
00967   int step   = BOUNDARY_STEP_SIZE;
00968   int wEdgel = w/step; 
00969   int hEdgel = h/step;
00970   itsContourMap = 
00971     Image<rutz::shared_ptr<Contour> >(wEdgel, hEdgel, ZEROS);
00972 
00973   // for each edgel RF location: we move horizontally first
00974   for(int j = 0; j < hEdgel; j++)
00975     {
00976       for(int i = 0; i < wEdgel; i++)
00977         {
00978           std::vector<rutz::shared_ptr<Edgel> >
00979             cEdgelList = itsEdgels.getVal(i,j);
00980           uint eNum = cEdgelList.size();
00981 
00982           // skip if there is no edgel in the location
00983           if(eNum == 0) continue;
00984 
00985           // skip if the edgel already belong to a contours
00986           // FIXXX: ADD THE JUNCTION RECOGNITION AS WELL HERE
00987           if(itsContourMap.getVal(i,j).is_valid()) continue;
00988 
00989           // start growing:
00990           rutz::shared_ptr<Edgel> edgel = cEdgelList[0];          
00991           rutz::shared_ptr<Contour> contour(new Contour());
00992           contour->edgels.push_back(edgel);
00993           itsContourBoundaries.push_back(contour);
00994           itsContourMap.setVal(i, j, contour);
00995 
00996           //LINFO("edgel: %3d %3d: ang: %d", i,j, edgel->angleIndex);
00997           
00998           // while can add an edgel
00999           // to the right (and/or bottom) of the contour
01000           bool isAddingRight = true;
01001           int ci = i, cj = j;
01002           //LINFO("start: ci: %3d cj: %3d: cd: %3d sign: %d", 
01003           //      ci, cj, edgel->angleIndex, 1);
01004           while(isAddingRight)
01005             { 
01006               isAddingRight = addEdgelToContour(ci, cj, edgel, contour, 1);
01007               if(isAddingRight) 
01008                 itsContourMap.setVal(ci, cj, contour);
01009               //LINFO("[%3d %3d]: %3d %3d %3d", ci,cj, 
01010               //      edgel->pt.i, edgel->pt.j, edgel->angleIndex);
01011             }
01012 
01013           //LINFO("Final Size: %d", 
01014           //      int(itsContourMap.getVal(i,j)->edgels.size()));
01015           //Raster::waitForKey();
01016 
01017           // // while can add an edgel
01018           // // to the left (and/or top) 
01019           // bool isAddingLeft = true;
01020           // while(isAddingLeft)
01021           //   {
01022           //     isAddingLeft = addEdgelToContour(i,j, edgel, contour, -1);  
01023           //   }
01024 
01025           
01026 
01027         }
01028     }
01029 
01030   // go through the contours
01031 
01032   // only link with negative edgels if the two contours are 
01033 
01034   // report the contours from the most salient first
01035   //  using equation that includes both length of contours and magnitude
01036 
01037   // better graphical drawing 
01038   //  --> probably need the stdev gradient info
01039 
01040 
01041 }
01042 
01043 // ######################################################################
01044 bool ContourBoundaryDetector::addEdgelToContour
01045 (int &ci, int &cj, 
01046  rutz::shared_ptr<Edgel> &edgel, 
01047  rutz::shared_ptr<Contour> &contour,
01048  int sign)
01049 {
01050   bool isAdding = false;
01051 
01052   // direction of the current edgel
01053   int cd = edgel->angleIndex;
01054 
01055   // get the neighborhood of edgel 
01056   // depending on the direction
01057   int ni1, nj1, ni2, nj2, ni3, nj3;
01058   getEdgelDirectionNeighborhood
01059     (ci, cj, cd, sign, ni1, nj1, ni2, nj2, ni3, nj3);
01060 
01061   // LINFO("ni1: %3d nj1: %3d", ni1, nj1);
01062   // LINFO("ni2: %3d nj2: %3d", ni2, nj2);
01063   // LINFO("ni3: %3d nj3: %3d", ni3, nj3);
01064 
01065   // only examine if coordinate is within image
01066   // and  there is an edgel in the coordinate 
01067   uint bEdge1 = 0;  uint bEdge2 = 0;  uint bEdge3 = 0;
01068 
01069   if(itsEdgels.coordsOk(ni1, nj1) &&      
01070      itsEdgels.getVal(ni1, nj1).size() != 0) bEdge1 = 1;
01071 
01072   if(itsEdgels.coordsOk(ni2, nj2) &&      
01073      itsEdgels.getVal(ni2, nj2).size() != 0) bEdge2 = 1;
01074 
01075   if(itsEdgels.coordsOk(ni3, nj3) &&      
01076      itsEdgels.getVal(ni3, nj3).size() != 0) bEdge3 = 1;
01077 
01078   uint totalEdgels = bEdge1+bEdge2+bEdge3;
01079 
01080   //LINFO("  have %d neighbors: "
01081   //      "(%d %d %d) (%d %d %d) (%d %d %d)", totalEdgels,
01082   //      ni1, nj1, bEdge1, ni2, nj2, bEdge2, ni3, nj3, bEdge3);
01083 
01084   // if no edgels to continue, we're done 
01085   if(totalEdgels == 0) 
01086     return false;
01087 
01088   // if there is only 1 edgel to continue to 
01089   else if(totalEdgels == 1) 
01090     {
01091       // get the next edgel
01092       int ni = -1; int nj = -1; int nd = -1;
01093       if(bEdge1 == 1)      { ni = ni1; nj = nj1; }
01094       else if(bEdge2 == 1) { ni = ni2; nj = nj2; }
01095       else if(bEdge3 == 1) { ni = ni3; nj = nj3; }      
01096       
01097       std::vector<rutz::shared_ptr<Edgel> >
01098         nEdgelList = itsEdgels.getVal(ni,nj);          
01099       rutz::shared_ptr<Edgel> nEdgel = nEdgelList[0]; 
01100 
01101       // get the next direction 
01102       nd = nEdgel->angleIndex;
01103 
01104       // if unclaimed and legal direction to add
01105       bool legalDirection = 
01106         isLegalDirection(ci, cj, cd, sign, ni, nj, nd);
01107       if(itsContourMap.getVal(ni,nj).is_invalid() && 
01108          legalDirection)
01109         {      
01110           // add and advance the end of contour
01111           isAdding = true;
01112           contour->edgels.push_back(nEdgel);
01113           edgel = nEdgel;
01114           ci = ni; cj = nj; cd = edgel->angleIndex;
01115           
01116           //LINFO("  --> [1]  adding ci: %3d cj: %3d: cd: %3d sign: %d", 
01117           //      ci, cj, cd, sign);
01118 
01119           // FIXXX: maybe there is a second go round 
01120           //        to connect even more
01121           //        so may not want to return right away
01122           return isAdding;
01123         }
01124       else return false;
01125     }
01126   else if(totalEdgels == 2)
01127     {
01128       // FIXXX: have to check for junction example:
01129       //        where there are more than 1 valid direction 
01130 
01131       // priority on what??
01132       int n1i = -1, n1j = -1, n2i = -1, n2j = -1;
01133       if(bEdge1 == 0)
01134         {
01135           n1i = ni2; n1j = nj2; n2i = ni3; n2j = nj3;       
01136         }
01137       else if(bEdge2 == 0)
01138         {
01139           n1i = ni1; n1j = nj1; n2i = ni3; n2j = nj3;       
01140       
01141         }
01142       else if(bEdge3 == 0)
01143         {
01144           n1i = ni1; n1j = nj1; n2i = ni2; n2j = nj2;       
01145         }
01146       
01147       // get the first next edgel
01148       std::vector<rutz::shared_ptr<Edgel> >
01149         nEdgelList1 = itsEdgels.getVal(n1i,n1j);          
01150       rutz::shared_ptr<Edgel> nEdgel1 = nEdgelList1[0]; 
01151       int n1d = nEdgel1->angleIndex;
01152 
01153       // get the second next edgel
01154       std::vector<rutz::shared_ptr<Edgel> >
01155         nEdgelList2 = itsEdgels.getVal(n2i,n2j);          
01156       rutz::shared_ptr<Edgel> nEdgel2 = nEdgelList2[0]; 
01157       int n2d = nEdgel2->angleIndex;
01158 
01159       // just deal with the diagonal first
01160       if(cd == 1)
01161         {
01162           if(n1d == 1 && n2d == 1)
01163             {
01164               // if edgel unclaimed
01165               if(itsContourMap.getVal(n2i,n2j).is_invalid())
01166                 {      
01167                   // add and advance the end of contour
01168                   isAdding = true;
01169                   contour->edgels.push_back(nEdgel2);
01170                   edgel = nEdgel2;
01171                   ci = n2i; cj = n2j;
01172 
01173                   //LINFO("  --> [2a] adding ci: %3d cj: %3d: cd: %3d sign: %d", 
01174                   //      ci, cj, cd, sign);
01175                   // FIXXX: maybe there is a second go round 
01176                   //        to connect even more
01177                   //        so may not want to return right away
01178                   return isAdding;
01179                 }
01180               else
01181                 { 
01182                   //LINFO("  CLAIMED"); 
01183                   return false; 
01184                 }
01185             }
01186           else 
01187             {
01188               //LINFO("  --> [21a]NOT 1(%3d %3d %3d) 2(%3d %3d %3d)", 
01189               //      n1i, n1j, n1d, n2i, n2j, n2d);
01190               return false;
01191             }
01192         }
01193       if(cd == 3)
01194         {
01195           if(n1d == 3 && n2d == 3)
01196             {
01197               // if edgel unclaimed
01198               if(itsContourMap.getVal(n2i,n2j).is_invalid())
01199                 {      
01200                   // add and advance the end of contour
01201                   isAdding = true;
01202                   contour->edgels.push_back(nEdgel2);
01203                   edgel = nEdgel2;
01204                   ci = n2i; cj = n2j; cd = edgel->angleIndex;
01205 
01206                   //LINFO("  --> [23b] adding ci: %3d cj: %3d: cd: %3d sign: %d", 
01207                   //      ci, cj, cd, sign);
01208 
01209                   // FIXXX: maybe there is a second go round 
01210                   //        to connect even more
01211                   //        so may not want to return right away
01212                   return isAdding;
01213                 }
01214               else
01215                 { 
01216                   //LINFO("  CLAIMED"); 
01217                   return false; 
01218                 }
01219             }
01220           else 
01221             {
01222               //LINFO("  --> [2b]NOT 1(%3d %3d %3d) 2(%3d %3d %3d)", 
01223               //      n1i, n1j, n1d, n2i, n2j, n2d);
01224               return false;
01225             }
01226         }
01227 
01228       if(cd == 0)
01229         {
01230           if(n1d == 3 && n2d == 3)
01231             {
01232               // if edgel unclaimed
01233               if(itsContourMap.getVal(n1i,n1j).is_invalid())
01234                 {      
01235                   // add and advance the end of contour
01236                   isAdding = true;
01237                   contour->edgels.push_back(nEdgel1);
01238                   edgel = nEdgel1;
01239                   ci = n1i; cj = n1j; cd = edgel->angleIndex;
01240 
01241                   //LINFO("  --> [20a] adding ci: %3d cj: %3d: cd: %3d sign: %d", 
01242                   //      ci, cj, cd, sign);
01243 
01244                   // FIXXX: maybe there is a second go round 
01245                   //        to connect even more
01246                   //        so may not want to return right away
01247                   return isAdding;
01248                 }
01249             }
01250           if(n1d == 1 && n2d == 1)
01251             {
01252               // if edgel unclaimed
01253               if(itsContourMap.getVal(n1i,n1j).is_invalid())
01254                 {      
01255                   // add and advance the end of contour
01256                   isAdding = true;
01257                   contour->edgels.push_back(nEdgel1);
01258                   edgel = nEdgel1;
01259                   ci = n1i; cj = n1j; cd = edgel->angleIndex;
01260 
01261                   //LINFO("  --> [20b] adding ci: %3d cj: %3d: cd: %3d sign: %d", 
01262                   //      ci, cj, cd, sign);
01263                   // FIXXX: maybe there is a second go round 
01264                   //        to connect even more
01265                   //        so may not want to return right away
01266                   return isAdding;
01267                 }
01268             }
01269 
01270           // FIXXX: pick the better of the two
01271           //LINFO("  Branch:0: junction");
01272           //LINFO("  --> [2_0]NOT 1(%3d %3d %3d) 2(%3d %3d %3d)", 
01273           //      n1i, n1j, n1d, n2i, n2j, n2d);
01274 
01275           return false;
01276         }
01277 
01278       if(cd == 2)
01279         {
01280           if(n1d == 3 && n2d == 3)
01281             {
01282               // if edgel unclaimed
01283               if(itsContourMap.getVal(n1i,n1j).is_invalid())
01284                 {      
01285                   // add and advance the end of contour
01286                   isAdding = true;
01287                   contour->edgels.push_back(nEdgel1);
01288                   edgel = nEdgel1;
01289                   ci = n1i; cj = n1j; cd = edgel->angleIndex;
01290                   
01291                   //LINFO("  --> [20a] adding ci: %3d cj: %3d: cd: %3d sign: %d", 
01292                   //      ci, cj, cd, sign);
01293                   
01294                   // FIXXX: maybe there is a second go round 
01295                   //        to connect even more
01296                   //        so may not want to return right away
01297                   return isAdding;
01298                 }
01299             }          
01300           if(n1d == 1 && n2d == 1)
01301             {
01302               // if edgel unclaimed
01303               if(itsContourMap.getVal(n1i,n1j).is_invalid())
01304                 {      
01305                   // add and advance the end of contour
01306                   isAdding = true;
01307                   contour->edgels.push_back(nEdgel1);
01308                   edgel = nEdgel1;
01309                   ci = n1i; cj = n1j; cd = edgel->angleIndex;
01310 
01311                   //LINFO("  --> [20b] adding ci: %3d cj: %3d: cd: %3d sign: %d", 
01312                   //      ci, cj, cd, sign);
01313                   // FIXXX: maybe there is a second go round 
01314                   //        to connect even more
01315                   //        so may not want to return right away
01316                   return isAdding;
01317                 }
01318             }
01319 
01320           // FIXXX: pick the better of the two
01321           //LINFO("  Branch:2: junction");
01322           //LINFO("  --> [2_2]NOT 1(%3d %3d %3d) 2(%3d %3d %3d)", 
01323           //      n1i, n1j, n1d, n2i, n2j, n2d);
01324           return false;
01325         }
01326 
01327 
01328   // if we can still include with the rules 
01329   // (in this order):
01330   //   1 space lookahead
01331   //   2 space lookahead
01332   //  [3 space lookahead]???
01333   //    note 3 space lookahead means 
01334   //    a gap of 16 pixels to interpolate    
01335   //    --> probably an absolute max
01336   
01337   // in any of these lookahead 
01338   // only consider positive NMS difference
01339   
01340   
01341   // if there is branching: also split (it's a junction)
01342   
01343   // if it's a junction
01344   // stop a contour on the end
01345   // split if it hit the middle
01346   
01347   // split if the NMS difference is very not similar 
01348   // remember merging is easy - but spliting is hard
01349   
01350   // add the edgel to the contour
01351  
01352   //LINFO("ci: %3d cj: %3d", ci, cj);
01353 
01354       
01355       // for now return false
01356       return false;
01357     }
01358   else if(totalEdgels == 3)
01359     {
01360       // more than likely it's too ambiguous to decide what to do
01361       LINFO("3 edgels. May want to take a look at it");
01362       //Raster::waitForKey();
01363       return false;
01364     }
01365 
01366 
01367   // if(itsEdgels.coordsOk(ni1, nj1) &&      
01368   //    itsEdgels.getVal(ni1, nj1).size() != 0 &&
01369   //    itsContourMap.getVal(ni1,nj1).is_invalid())
01370   //   {
01371   //     std::vector<rutz::shared_ptr<Edgel> >
01372   //       nEdgelList1 = itsEdgels.getVal(ni1,nj1);
01373       
01374   //     rutz::shared_ptr<Edgel> nEdgel1 = nEdgelList1[0]; 
01375       
01376   //     //LINFO("e1: %d %d %d %f", nEdgel1->pt.i, nEdgel1->pt.j, nEdgel1->angleIndex, nEdgel1->val);
01377 
01378   //     // if the next edgel can be added to the contour
01379   //     if(nEdgel1->angleIndex == aind) 
01380   //       {
01381   //         // add and advance the end of contour
01382   //         isAdding = true;
01383   //         contour->edgels.push_back(nEdgel1);
01384   //         edgel = nEdgel1;
01385 
01386   //         //LINFO("in1");
01387 
01388   //         ci = ni1; cj = nj1;
01389   //       }
01390   //   }
01391   // else if(itsEdgels.coordsOk(ni2, nj2) &&
01392   //    itsEdgels.getVal(ni2, nj2).size() != 0 &&
01393   //    itsContourMap.getVal(ni2,nj2).is_invalid())
01394   //   {
01395   //     std::vector<rutz::shared_ptr<Edgel> >
01396   //       nEdgelList2 = itsEdgels.getVal(ni2,nj2);
01397       
01398   //     rutz::shared_ptr<Edgel> nEdgel2 = nEdgelList2[0]; 
01399       
01400   //     //LINFO("e2: %d %d %d %f", nEdgel2->pt.i, nEdgel2->pt.j, nEdgel2->angleIndex, nEdgel2->val);
01401 
01402   //     // if the next edgel can be added to the contour
01403   //     if(nEdgel2->angleIndex == aind) 
01404   //       {
01405   //         // add and advance the end of contour
01406   //         isAdding = true;
01407   //         contour->edgels.push_back(nEdgel2);
01408   //         edgel = nEdgel2;
01409 
01410   //         //LINFO("in2");
01411 
01412   //         ci = ni2; cj = nj2;
01413   //       }
01414   //   }
01415   // else if(itsEdgels.coordsOk(ni3, nj3) &&
01416   //    itsEdgels.getVal(ni3, nj3).size() != 0 &&
01417   //    itsContourMap.getVal(ni3,nj3).is_invalid())
01418   //   {
01419   //     std::vector<rutz::shared_ptr<Edgel> >
01420   //       nEdgelList3 = itsEdgels.getVal(ni3,nj3);
01421       
01422   //     rutz::shared_ptr<Edgel> nEdgel3 = nEdgelList3[0]; 
01423       
01424   //     //LINFO("e3: %d %d %d %f", nEdgel3->pt.i, nEdgel3->pt.j, nEdgel3->angleIndex, nEdgel3->val);
01425 
01426   //     // if the next edgel can be added to the contour
01427   //     if(nEdgel3->angleIndex == aind) 
01428   //       {
01429   //         // add and advance the end of contour
01430   //         isAdding = true;
01431   //         contour->edgels.push_back(nEdgel3);
01432   //         edgel = nEdgel3;
01433 
01434   //         //LINFO("in3");
01435 
01436   //         ci = ni3; cj = nj3;
01437   //       }
01438   //   }
01439 
01440 
01441 
01442 
01443   return isAdding;
01444 }
01445 
01446 // ######################################################################
01447 void ContourBoundaryDetector::getEdgelDirectionNeighborhood
01448 (int ci, int cj, int dir, int sign, 
01449  int &ni1, int &nj1, int &ni2, int &nj2, int &ni3, int &nj3)
01450 {
01451   // 0 degrees: | boundary 
01452   if(dir == 0)
01453     {
01454       ni1 = ci;      nj1 = cj + sign;
01455       ni2 = ci-sign; nj2 = cj + sign;
01456       ni3 = ci+sign; nj3 = cj + sign;
01457     }      
01458   // 45 degrees: / boundary
01459   else if(dir == 1)
01460     {
01461       ni1 = ci-sign; nj1 = cj + sign;
01462       ni2 = ci-sign; nj2 = cj;
01463       ni3 = ci;      nj3 = cj + sign;
01464     }
01465 
01466   // 90 degrees: -- boundary
01467   else if(dir == 2)
01468     {
01469       ni1 = ci+sign; nj1 = cj;
01470       ni2 = ci+sign; nj2 = cj - sign;
01471       ni3 = ci+sign; nj3 = cj + sign;
01472     }
01473 
01474   // 135 degrees: \ boundary
01475   else if(dir == 3)
01476     {
01477       ni1 = ci+sign; nj1 = cj + sign;
01478       ni2 = ci+sign; nj2 = cj;
01479       ni3 = ci;      nj3 = cj + sign;
01480     }
01481 }
01482 
01483 // ######################################################################
01484 bool ContourBoundaryDetector::isLegalDirection
01485 (int ci, int cj, int cd, int sign, int ni, int nj, int nd)
01486 {
01487   // important that this operation symmetric!!!
01488   // that is: (a isLegalDirection b) == (b isLegalDirection a)
01489 
01490   int di = sign*(ni - ci);
01491   int dj = sign*(nj - cj);
01492 
01493   // LINFO("%d %d %d - %d - %d %d %d ", ci, cj, cd, sign, ni, nj, nd);
01494   // LINFO("di dj: %d %d", di, dj);
01495 
01496   // 0 degrees: | boundary:  
01497   if(cd == 0)
01498     {
01499       // loc n1
01500       if     (di ==  0 && dj ==  1)
01501         {
01502           return (nd != 2);
01503         }
01504       // loc n2
01505       else if(di == -1 && dj ==  1)
01506         {
01507           return (nd != 3);
01508         }
01509       // loc n3
01510       else if(di ==  1 && dj ==  1)
01511         {
01512           return (nd != 1);
01513         }
01514       else 
01515         {
01516           LFATAL("invalid dir. loc:"
01517                  "%d %d %d - %d - %d %d %d "
01518                  "di dj: %d %d",
01519                  ci, cj, cd, sign, ni, nj, nd, di, dj);
01520           return false;
01521         }
01522     }      
01523   // 45 degrees: / boundary
01524   else if(cd == 1)
01525     {
01526       // loc n1
01527       if     (di == -1 && dj ==  1)
01528         {
01529           return (nd != 3);
01530         }
01531       // loc n2
01532       else if(di == -1 && dj ==  0)
01533         {
01534           return (nd != 0);
01535         }
01536       // loc n3
01537       else if(di ==  0 && dj ==  1)
01538         {
01539           return (nd != 2);
01540         }
01541       else 
01542         {
01543           LFATAL("invalid dir. loc:"
01544                  "%d %d %d - %d - %d %d %d "
01545                  "di dj: %d %d",
01546                  ci, cj, cd, sign, ni, nj, nd, di, dj);
01547           return false;
01548         }
01549     }
01550 
01551   // 90 degrees: -- boundary
01552   else if(cd == 2)
01553     {
01554       // loc n1
01555       if     (di ==  1 && dj ==  0)
01556         {
01557           return (nd != 0);
01558         }
01559       // loc n2
01560       else if(di ==  1 && dj == -1)
01561         {
01562           return (nd != 3);
01563         }
01564       // loc n3
01565       else if(di ==  1 && dj ==  1)
01566         {
01567           return (nd != 1);
01568         }
01569       else 
01570         {
01571           LFATAL("invalid dir. loc:"
01572                  "%d %d %d - %d - %d %d %d "
01573                  "di dj: %d %d",
01574                  ci, cj, cd, sign, ni, nj, nd, di, dj);
01575           return false;
01576         }
01577     }
01578 
01579   // 135 degrees: \ boundary
01580   else if(cd == 3)
01581     {
01582       // loc n1
01583       if     (di ==  1 && dj ==  1)
01584         {
01585           return (nd != 1);
01586         }
01587       // loc n2
01588       else if(di ==  1 && dj ==  0)
01589         {
01590           return (nd != 0);
01591         }
01592       // loc n3
01593       else if(di ==  0 && dj ==  1)
01594         {
01595           return (nd != 2);
01596         }
01597       else 
01598         {
01599           LFATAL("invalid dir. loc:"
01600                  "%d %d %d - %d - %d %d %d "
01601                  "di dj: %d %d",
01602                  ci, cj, cd, sign, ni, nj, nd, di, dj);
01603           return false;
01604         }
01605     }
01606   else LFATAL("invalid direction");
01607   return false;
01608 }
01609 
01610 // ######################################################################
01611 Image<PixRGB<byte> > ContourBoundaryDetector::getContourBoundaryImage()
01612 {
01613   int w = itsImage.getWidth();
01614   int h = itsImage.getHeight();
01615 
01616   Image<PixRGB<byte> > image = itsImage;
01617   float mVal = 32;
01618   float bVal = 255 - mVal;
01619   Image<byte> dImaR, dImaG, dImaB;
01620   getComponents(image, dImaR, dImaG, dImaB);
01621   inplaceNormalize(dImaR, byte(0), byte(mVal));
01622   inplaceNormalize(dImaG, byte(0), byte(mVal));
01623   inplaceNormalize(dImaB, byte(0), byte(mVal));
01624   Image<PixRGB<byte> > dIma  = makeRGB(dImaR,dImaG,dImaB);
01625 
01626   int hstep = BOUNDARY_STEP_SIZE/2;
01627 
01628   // clean image
01629   Image<byte> dCBmapR(w,h,ZEROS);
01630   Image<byte> dCBmapG(w,h,ZEROS);
01631   Image<byte> dCBmapB(w,h,ZEROS);
01632 
01633   for(uint i = 0; i < itsContourBoundaries.size(); i++)
01634     {
01635       rutz::shared_ptr<Contour> contour = itsContourBoundaries[i];
01636 
01637       byte rcVal = bVal*rand()/(RAND_MAX+1.0);
01638       byte gcVal = bVal*rand()/(RAND_MAX+1.0);
01639       byte bcVal = bVal*rand()/(RAND_MAX+1.0);
01640 
01641       if(contour->edgels.size() < 5) continue;
01642 
01643       for(uint j = 0; j < contour->edgels.size(); j++)
01644         {
01645           rutz::shared_ptr<Edgel> edgel = contour->edgels[j];
01646           
01647           uint aind = edgel->angleIndex; 
01648           float baind = 
01649             fmod((aind+(NUM_RIDGE_DIRECTIONS/2)),NUM_RIDGE_DIRECTIONS);
01650           
01651           float dx = cos(baind * M_PI/4.0) * hstep; 
01652           float dy = sin(baind * M_PI/4.0) * hstep;
01653           
01654           Point2D<int> pt = edgel->pt;
01655           Point2D<int> p1 = pt + Point2D<int>( dx+.5,  dy+.5); 
01656           Point2D<int> p2 = pt + Point2D<int>(-dx-.5, -dy-.5);
01657           
01658           //uint iEdgel = i/step;
01659           //uint jEdgel = j/step;
01660           //if(iEdgel >= 10 && iEdgel <= 25 && jEdgel >= 1 && jEdgel <= 14)
01661           //   {
01662           //     LINFO("maxk: %d -> %d  -> %f  %f (%f %f %f %f)  |%f %f", 
01663           //           maxk, (maxk+(NUM_RIDGE_DIRECTIONS/2)), 
01664           //           (borderK*M_PI)/float(NUM_RIDGE_DIRECTIONS), borderK, 
01665           //           cos(0), cos(M_PI/4.0), cos(M_PI/2.0), cos(M_PI*.75),
01666           //           dx, dy);
01667           
01668           //LINFO("%d %d | %d %d | %d %d::::: %d %d  %d", 
01669           //         pt.i, pt.j, p1.i, p1.j, p2.i, p2.j, iEdgel, jEdgel, maxk);       
01670 
01671           // draw the edgel
01672 
01673           // draw the straightline contour in the image 
01674           // for visualization
01675           drawLine(dCBmapR, p1, p2, byte(rcVal));
01676           drawLine(dCBmapG, p1, p2, byte(gcVal));
01677           drawLine(dCBmapB, p1, p2, byte(bcVal));
01678           //drawDisk(itsContourBoundaryImage, pt, 2,  255.0F);
01679         }
01680     }
01681 
01682   Image<PixRGB<byte> > dCBmap = makeRGB(dCBmapR,dCBmapG,dCBmapB);
01683   return Image<PixRGB<byte> >(dIma+dCBmap);
01684 }
01685 
01686 // ######################################################################
01687 void ContourBoundaryDetector::displayContourBoundary()
01688 {
01689   int w = itsImage.getWidth();
01690   int h = itsImage.getHeight();
01691 
01692   if(itsWin.is_invalid())
01693     itsWin.reset(new XWinManaged(Dims(w,h), w+10, 0, "Contours"));
01694   else itsWin->setDims(Dims(w,h));
01695 
01696   // display the contour
01697   itsWin->drawImage(getContourBoundaryImage(),0,0);
01698   //Raster::waitForKey();
01699 }
01700 
01701 // ######################################################################
01702 void ContourBoundaryDetector::displayGradImage
01703 (std::vector<Image<float> > gradImg)
01704 {
01705   Image<float> gradImgX = gradImg[0];
01706   Image<float> gradImgY = gradImg[1];
01707 
01708   uint w = gradImgX.getWidth();
01709   uint h = gradImgX.getHeight();
01710 
01711   Image<PixRGB<byte> > dGradImg(w, h, NO_INIT);
01712 
01713   float mn, mx;
01714   Image<float> temp = gradImgX*gradImgX + gradImgY*gradImgY;
01715   getMinMax(temp, mn,mx);
01716   float max = pow(mx,.5);
01717 
01718   for(uint i = 0; i < w; i++)
01719     {
01720       for(uint j = 0; j < h; j++)
01721         {
01722           float x = gradImgX.getVal(i,j);
01723           float y = gradImgY.getVal(i,j);
01724                     
01725           float dir = atan2(y,x) * 180.0 / M_PI;
01726           if(dir < 0.0) dir = 360.0 + dir;
01727           float mag = pow(x*x + y*y, 0.5);
01728 
01729           // PixHSV<byte> hsvp(byte(dir/360.0 * 255.0),
01730           //                100,
01731           //                byte(max/max* 255.0));                
01732 
01733           PixHSV<byte> hsvp(180,100,50);
01734 
01735           //PixRGB<byte> rgbp(hsvp);
01736           byte m = byte(mag/max * 255.0);
01737           //byte d = byte(dir/360.0 * 255.0);
01738           PixRGB<byte> rgbp(m, m, m);
01739           dGradImg.setVal(i,j, rgbp);
01740 
01741           // PixRGB<byte> a(0,255,255);
01742           // PixHSV<byte> b(a);
01743           // PixRGB<byte> c(b);
01744           // LINFO("%d %d %d ==> %d %d %d ==> %d %d %d",                
01745           //    a.red(), a.green(), a.blue(),
01746           //    b.H(), b.S(), b.V(),
01747           //    c.red(), c.green(), c.blue());
01748           // Raster::waitForKey();
01749 
01750           // LINFO("%f %f || %f %f (%f): %d %d %d ==> %d %d %d", 
01751           //    x,y, dir, mag, mag/max* 255.0,
01752           //    hsvp.H(), hsvp.S(), hsvp.V(),
01753           //    rgbp.red(), rgbp.green(), rgbp.blue());
01754         }
01755     }
01756 
01757   itsWin->drawImage(dGradImg, 0,0);
01758   Raster::waitForKey(); 
01759 }
01760 
01761 
01762 // ######################################################################
01763 /* So things look consistent in everyone's emacs... */
01764 /* Local Variables: */
01765 /* indent-tabs-mode: nil */
01766 /* End: */
Generated on Sun May 8 08:40:39 2011 for iLab Neuromorphic Vision Toolkit by  doxygen 1.6.3