SalientRegionSegmenter.C

Go to the documentation of this file.

/*!@file Gist/SalientRegionSegmenter.C segment out object depicted by the
salient point. Here we use both region growing and boundary detection
 */
// //////////////////////////////////////////////////////////////////// //
// The iLab Neuromorphic Vision C++ Toolkit - Copyright (C) 2001 by the //
// University of Southern California (USC) and the iLab at USC.         //
// See http://iLab.usc.edu for information about this project.          //
// //////////////////////////////////////////////////////////////////// //
// Major portions of the iLab Neuromorphic Vision Toolkit are protected //
// under the U.S. patent ``Computation of Intrinsic Perceptual Saliency //
// in Visual Environments, and Applications'' by Christof Koch and      //
// Laurent Itti, California Institute of Technology, 2001 (patent       //
// pending; application number 09/912,225 filed July 23, 2001; see      //
// http://pair.uspto.gov/cgi-bin/final/home.pl for current status).     //
// //////////////////////////////////////////////////////////////////// //
// This file is part of the iLab Neuromorphic Vision C++ Toolkit.       //
//                                                                      //
// The iLab Neuromorphic Vision C++ Toolkit is free software; you can   //
// redistribute it and/or modify it under the terms of the GNU General  //
// Public License as published by the Free Software Foundation; either  //
// version 2 of the License, or (at your option) any later version.     //
//                                                                      //
// The iLab Neuromorphic Vision C++ Toolkit is distributed in the hope  //
// that it will be useful, but WITHOUT ANY WARRANTY; without even the   //
// implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR      //
// PURPOSE.  See the GNU General Public License for more details.       //
//                                                                      //
// You should have received a copy of the GNU General Public License    //
// along with the iLab Neuromorphic Vision C++ Toolkit; if not, write   //
// to the Free Software Foundation, Inc., 59 Temple Place, Suite 330,   //
// Boston, MA 02111-1307 USA.                                           //
// //////////////////////////////////////////////////////////////////// //
//
// Primary maintainer for this file: Christian Siagian <siagian@usc.edu>
// $HeadURL: svn://isvn.usc.edu/software/invt/trunk/saliency/src/Gist/SalientRegionSegmenter.C $
// $Id: $
//
//////////////////////////////////////////////////////////////////////////

#include "Gist/SalientRegionSegmenter.H"

#include "Image/CutPaste.H"     // for inplacePaste()
#include "Image/DrawOps.H"
#include "Image/ColorOps.H"
#include "Image/ShapeOps.H"
#include "Image/Kernels.H"
#include  <cstdio>

#define MINIMUM_CONTOUR_LENGTH     5

#define LEN_WEIGHT               .20 // contour boundary length weight
#define SIM_WEIGHT               .20 // similarity with closer side weight 
#define STR_WEIGHT               .20 // strength of boundary weight
#define RAT_WEIGHT               .20 // ratio of in:out weight
#define DIS_WEIGHT               .20 // distance between sal pt and segment

#define NEIGHBORHOOD_RADIUS       12
//#define NEIGHBORHOOD_RADIUS        8
#define MAX_FAR_TO_NEAR_RATIO    5.0
#define MAX_DISTANCE             .25 // quarter of the diagonal

#define NUM_CHANNELS               3 // current number of channels   
#define NUM_L_BINS                25 // # bins for L component of CIELab
#define NUM_A_BINS                25 // # bins for a component of CIELab
#define NUM_B_BINS                25 // # bins for b component of CIELab
#define NUM_HISTOGRAM_DIMS        75 // current total number of bins

#define UNVISITED_REGION           0 // unvisited mean, just looked at
                                     // for the first time 
#define CENTER_REGION              1
#define CENTER_CORE_REGION         2
#define SURROUND_REGION            3
#define UNASSIGNED_REGION          4
#define OTHER_OBJECT_REGION        5 

#define CLOSER_SIDE                1
#define FAR_SIDE                   2

// ######################################################################
// ######################################################################

// ######################################################################
SalientRegionSegmenter::SalientRegionSegmenter()
:
itsContourBoundaryDetector(new ContourBoundaryDetector())
{
  itsWin.reset();

  // These color features are inspired by Martin, PAMI 2004 pb paper
  float bg_smooth_sigma    = 0.1;       // bg histogram smoothing sigma
  float cg_smooth_sigma    = 0.05;      // cg histogram smoothing sigma

  itsLKernel = 
    gaussian<float>(1.0F, NUM_L_BINS*bg_smooth_sigma, 
                    int(3.0F*NUM_L_BINS*bg_smooth_sigma + 0.5F));
  itsAKernel = 
    gaussian<float>(1.0F, NUM_A_BINS*cg_smooth_sigma, 
                    int(3.0F*NUM_A_BINS*cg_smooth_sigma + 0.5F));
  itsBKernel = 
    gaussian<float>(1.0F, NUM_B_BINS*cg_smooth_sigma, 
                    int(3.0F*NUM_B_BINS*cg_smooth_sigma + 0.5F));

  float suml = sum(itsLKernel); itsLKernel = itsLKernel/suml;
  float suma = sum(itsAKernel); itsAKernel = itsAKernel/suma;
  float sumb = sum(itsBKernel); itsBKernel = itsBKernel/sumb;

  // LINFO("check kernels");

  // for(int i = 0; i < itsLKernel.getWidth(); i++)
  //   {
  //     for(int j = 0; j < itsLKernel.getHeight(); j++)
  //       printf("%10.3f ", itsLKernel.getVal(i,j));
  //     printf("\n");
  //   }
  // Raster::waitForKey();

  // for(int i = 0; i < aKernel.getWidth(); i++)
  //   {
  //     for(int j = 0; j < aKernel.getHeight(); j++)
  //       printf("%10.3f ", aKernel.getVal(i,j));
  //     printf("\n");
  //   }
  // Raster::waitForKey();

  // for(int i = 0; i < bKernel.getWidth(); i++)
  //   {
  //     for(int j = 0; j < bKernel.getHeight(); j++)
  //       printf("%10.3f ", bKernel.getVal(i,j));
  //     printf("\n");
  //   }
  // Raster::waitForKey();
}

// ######################################################################
SalientRegionSegmenter::~SalientRegionSegmenter()
{ }

// ######################################################################
void SalientRegionSegmenter::setImage(Image<PixRGB<byte> > image)
{ 
  // reset everything
  itsSalientPoints.clear();

  itsImage = image;

  // this is the feature histogram in Martin's pb PAMI 2004
  setImageFeatureHistogramValues();

  uint w = itsImage.getWidth();
  uint h = itsImage.getHeight();

  if(itsWin.is_invalid())
    itsWin.reset(new XWinManaged(Dims(4*w,2*h), 0, 0, 
                                 "SalReg Segmenter"));
  else itsWin->setDims(Dims(4*w,2*h));

  // calculate the contour boundary and region map
  computeBoundaryAndRegionInformation();

  // keep track of what histograms has been computed
  uint hstep = BOUNDARY_STEP_SIZE/2;
  uint wG = (w-hstep)/BOUNDARY_STEP_SIZE;
  uint hG = (h-hstep)/BOUNDARY_STEP_SIZE;
  itsComputedHistograms = 
    Image<rutz::shared_ptr<Histogram> >(wG, hG, ZEROS);
}

// ######################################################################
void SalientRegionSegmenter::setImageFeatureHistogramValues()
{
  // convert to CIElab color space
  Image<float> lImg; 
  Image<float> aImg;
  Image<float> bImg;
  getNormalizedLAB(itsImage, lImg, aImg, bImg);
  
  // convert to bin numbers
  itsImageHistogramEntries.clear();
  itsImageHistogramEntries.push_back(quantize_values(lImg, NUM_L_BINS));
  itsImageHistogramEntries.push_back(quantize_values(aImg, NUM_A_BINS));
  itsImageHistogramEntries.push_back(quantize_values(bImg, NUM_B_BINS));

  // uint w = itsImage.getWidth();
  // uint h = itsImage.getHeight();
  // for(uint i = 0; i < w; i++)
  //   for(uint j = 0; j < h; j++)
  //     {
  //       if(lImg.getVal(i,j) > 1.0) LINFO("limg: %3d %3d", i,j);
  //       if(aImg.getVal(i,j) > 1.0) LINFO("aimg: %3d %3d", i,j);
  //       if(bImg.getVal(i,j) > 1.0) LINFO("bimg: %3d %3d", i,j);
  //     }
}

// ######################################################################
Image<int> SalientRegionSegmenter::quantize_values
(Image<float> image, int num_bins, bool normalize)
{
  // FIXXX: May want to move to MathOps later
  // if(normalize)
  //   inplaceNormalize(image, 0.0F, 1.0F);

  uint w = itsImage.getWidth();
  uint h = itsImage.getHeight();
  Image<int> result(w,h,NO_INIT);

  Image<float>::iterator iptr = image.beginw();  

  Image<int>::iterator aptr = result.beginw(), stop = result.endw();  
  float nbins = num_bins;
  while(aptr != stop)
    {
      float in = *iptr++;

      // bins will be 0 to num_bins-1
      int bin = int(in*nbins);
      if(bin == num_bins) bin = num_bins - 1;
      *aptr++ = bin;

      //LINFO("val: %f --> bin: %d", in, bin);
    }

  return result;
}

// ######################################################################
void SalientRegionSegmenter::computeBoundaryAndRegionInformation()
{
  computeBoundaryInformation(); 
  computeRegionInformation(); 

  displayBoundaryAndRegionInformation();
}

// ######################################################################
void SalientRegionSegmenter::computeBoundaryInformation()
{
  // compute contour boundary map
  Timer timer(1000000); timer.reset();
  itsContourBoundaryDetector->computeContourBoundary(itsImage);
  LINFO("time: %f ms", timer.get()/1000.0);

  itsContourBoundaries = 
    itsContourBoundaryDetector->getContourBoundaries();

  float maxVal = 0.0; 
  for(uint i = 0; i < itsContourBoundaries.size(); i++)
    {
      rutz::shared_ptr<Contour> contour = 
        itsContourBoundaries[i];
      
      // go through each contour
      uint cLen = contour->edgels.size();
      for(uint j = 0; j < cLen; j++)
        {
          float val = contour->edgels[j]->val;
          if(maxVal < val) maxVal = val;
        }
    }

  itsMaximumEdgelStrength = maxVal;
}

// ######################################################################
void SalientRegionSegmenter::computeRegionInformation()
{
  // default parameters for the graph-based segmentation
  // NOTE: in original Felzenswalb-Huttenlocher paper: 
  // k = 500 --> here k = 25 create over-segmented object
  float sigma = .5; uint k = 25; uint min_size = 20; int num_ccs;

  itsInitialRegionImage = 
    SuperPixelSegment
    (itsImage, sigma, k, min_size, num_ccs, &itsInitialRegions);

  itsInitialRegionColorImage = 
    SuperPixelDebugImage(itsInitialRegions, itsImage);
 
  itsRegionAdjecencyList = getRegionAdjecencyList(); 
  //  displayRegionAdjecencyList();
}

// ######################################################################
std::vector<std::vector<uint> > 
SalientRegionSegmenter::getRegionAdjecencyList()
{
  uint nRegions = itsInitialRegions.size();
  std::vector<std::vector<uint> >  adjList(nRegions);
  for(uint i = 0; i < nRegions; i++)
    adjList[i] = std::vector<uint>();

  // find regions:
  // FIXXX: there may be a better way

  // check out all the 
  for(uint i = 0; i < nRegions; i++)
    {
      for(uint j = 0; j < itsInitialRegions[i].size(); j++)
        {          
          Point2D<int> pt = itsInitialRegions[i][j];          
          uint label = i;

          for(int di = -1; di <= 1; di++)
            {
              for(int dj = -1; dj <= -1; dj++)
                {
                  Point2D<int> dpt(di,dj);
                  Point2D<int> pt2 = pt + dpt;
                  if(!(di == 0 && dj == 0) &&  
                     itsInitialRegionImage.coordsOk(pt2))
                    {
                      //if the neighbor is different
                      uint label2 = itsInitialRegionImage.getVal(pt2);

                      // check if the label is already 
                      // in the adjecency list
                      if(label != label2)
                        {
                          bool in = false;
                          for(uint l = 0; l < adjList[i].size(); l++)
                            {
                              if(adjList[i][l] == label2) 
                                { in = true; l = adjList[i].size(); }
                            }
                          if(!in) adjList[i].push_back(label2);
                        }
                    }
                }
            }
        }
    }  
  return adjList;
}

// ######################################################################
void SalientRegionSegmenter::displayRegionAdjecencyList()
{
  // display stuff
  Image<PixRGB<byte> > ima = itsImage;

  Image<float> fIma(luminance(ima));
  Image<float> tempFIma = fIma;
  inplaceNormalize(tempFIma, 0.0f,255.0f);  
  Image<byte>  bIma(tempFIma); 

  // get non-max suppressed boundary map image
  float mVal = 32;
  float bVal = 255 - mVal;

  Image<byte> dImaR, dImaG, dImaB;
  getComponents(ima, dImaR, dImaG, dImaB);
  inplaceNormalize(dImaR, byte(0), byte(mVal));
  inplaceNormalize(dImaG, byte(0), byte(mVal));
  inplaceNormalize(dImaB, byte(0), byte(mVal));
  Image<PixRGB<byte> > dIma  = makeRGB(dImaR,dImaG,dImaB);

  uint w = ima.getWidth();
  uint h = ima.getHeight();
  Image<PixRGB<byte> > dispIma(4*w,2*h,ZEROS);

  inplacePaste (dispIma, ima, Point2D<int>(0,0));

  uint nRegions = itsInitialRegions.size();
  for(uint i = 0; i < nRegions; i++)
    {     
      LINFO("region: %d (%d pix): has %d neighbors", 
            i, int(itsInitialRegions[i].size()), 
            int(itsRegionAdjecencyList[i].size()));
      Image<float> regionMap(w,h, ZEROS);
      
      for(uint j = 0; j < itsInitialRegions[i].size(); j++)
        {
          Point2D<int> pt = itsInitialRegions[i][j]; 
          regionMap.setVal(pt, 1.0);

          //LINFO("  pt[%3d]: %3d %3d ", j, pt.i, pt.j);
        }

      for(uint j = 0; j < itsRegionAdjecencyList[i].size(); j++)
        {
          uint label = itsRegionAdjecencyList[i][j];
          LINFO("  [%3d] label: %d has: %d pixels", 
                j, label,  int(itsInitialRegions[label].size()));
          float randVal = rand()/(RAND_MAX+1.0);
          for(uint k = 0; k < itsInitialRegions[label].size(); k++)
            {
              Point2D<int> ptNeighbor = itsInitialRegions[label][k]; 
              regionMap.setVal(ptNeighbor, randVal);
              //LINFO("  pt[%3d]: %3d %3d ", j, ptNeighbor.i, ptNeighbor.j);
            }
        }
      
      inplaceNormalize(regionMap, 0.0f,bVal);
      Image<byte> dBmapc(regionMap);
      Image<PixRGB<byte> > dBmap = toRGB(dBmapc);
  
      inplacePaste
        (dispIma, Image<PixRGB<byte> >(dIma+dBmap), Point2D<int>(w,0));

      itsWin->drawImage(dispIma, 0,0);

      Raster::waitForKey();
    }
}


// ######################################################################
void SalientRegionSegmenter::displayBoundaryAndRegionInformation()
{
  // display stuff
  Image<PixRGB<byte> > ima = itsImage;
  uint w = ima.getWidth();
  uint h = ima.getHeight();

  Image<float> fIma(luminance(ima));
  Image<float> tempFIma = fIma;
  inplaceNormalize(tempFIma, 0.0f,255.0f);  
  Image<byte>  bIma(tempFIma); 

  // get non-max suppressed boundary map image
  float mVal = 32;
  float bVal = 255 - mVal;
  Image<byte> dImaR, dImaG, dImaB;
  getComponents(ima, dImaR, dImaG, dImaB);
  inplaceNormalize(dImaR, byte(0), byte(mVal));
  inplaceNormalize(dImaG, byte(0), byte(mVal));
  inplaceNormalize(dImaB, byte(0), byte(mVal));
  Image<PixRGB<byte> > dIma  = makeRGB(dImaR,dImaG,dImaB);

  Image<float> tboundaryMap = 
    itsContourBoundaryDetector->getVarianceRidgeBoundaryMap();
  Image<float> boundaryMap = 
    clampedDiff(tboundaryMap, Image<float>(w,h,ZEROS));
  inplaceNormalize(boundaryMap, 0.0f,bVal);
  Image<byte> dBmapc(boundaryMap);
  Image<PixRGB<byte> > tdBmap = toRGB(dBmapc);
  Image<PixRGB<byte> > dBmap(dIma+tdBmap);
  //inplacePaste(dispIma, Image<PixRGB<byte> >(dIma+dBmap), Point2D<int>(w,0));

  // the non-max suppressed boundary map image
  Image<float> dBmapNMStemp = 
    itsContourBoundaryDetector->getNmsBoundaryMap();
  inplaceNormalize(dBmapNMStemp, 0.0F, 255.0F);
  Image<PixRGB<byte> > dBmapNMS = 
    toRGB(Image<byte>(dBmapNMStemp));

  // the contour boundary edgel map image
  Image<float> dCBEmapTemp = 
    itsContourBoundaryDetector->getEdgelBoundaryMap();
  inplaceNormalize(dCBEmapTemp, 0.0F, bVal);
  Image<PixRGB<byte> > tdCBEmap = 
    toRGB(Image<byte>(dCBEmapTemp));
  Image<PixRGB<byte> > dCBEmap(dIma+tdCBEmap);

  // get the contour boundary map
  Image<PixRGB<byte> > dCBmapTemp =
    itsContourBoundaryDetector->getContourBoundaryMap();
  Image<PixRGB<byte> > dCBmap(dIma+dCBmapTemp);

  // add the receptive field
  // 300, 140
  //uint cl = 156; uint ct = 68; uint cr = 8;
  //uint sl = cl - 4; uint st = ct - 4; uint sr = 16;

  uint cl = 90; uint ct = 14; uint cr = 8;
  uint sl = cl - 4; uint st = ct - 4; uint sr = 16;

  // to the image
  drawRect(ima,
           Rectangle(Point2D<int>(cl, ct), Dims(cr, cr)), 
           PixRGB<byte>(255,0,0));
  drawRect(ima,
           Rectangle(Point2D<int>(sl, st), Dims(sr, sr)), 
           PixRGB<byte>(255,255,0));
  //Raster::WriteRGB(ima, std::string("starfish_image.png"));

  // to the VRD
  drawRect(dBmap,
           Rectangle(Point2D<int>(cl, ct), Dims(cr, cr)), 
           PixRGB<byte>(255,0,0));
  drawRect(dBmap,
           Rectangle(Point2D<int>(sl, st), Dims(sr, sr)), 
           PixRGB<byte>(255,255,0));
  //Raster::WriteRGB(dBmap, std::string("starfish_vrd_image.png"));

  // to the NMS VRD
  drawRect(dBmapNMS, 
           Rectangle(Point2D<int>(cl, ct), Dims(cr, cr)), 
           PixRGB<byte>(255,0,0));
  drawRect(dBmapNMS, 
           Rectangle(Point2D<int>(sl, st), Dims(sr, sr)), 
           PixRGB<byte>(255,255,0));
  //Raster::WriteRGB(dBmapNMS, std::string("starfish_nms_vrd_image.png"));

  // to the edgel map
  drawRect(dCBEmap,
           Rectangle(Point2D<int>(cl, ct), Dims(cr, cr)), 
           PixRGB<byte>(255,0,0));
  drawRect(dCBEmap,
           Rectangle(Point2D<int>(sl, st), Dims(sr, sr)), 
           PixRGB<byte>(255,255,0));
  //Raster::WriteRGB(dCBEmap, std::string("starfish_cbe_image.png"));
  
  // to the contour map
  drawRect(dCBmap, 
           Rectangle(Point2D<int>(cl, ct), Dims(cr, cr)), 
           PixRGB<byte>(255,0,0));
  drawRect(dCBmap, 
           Rectangle(Point2D<int>(sl, st), Dims(sr, sr)), 
           PixRGB<byte>(255,255,0));
  //Raster::WriteRGB(dCBmap, std::string("starfish_cb_image.png"));

  // setup the display map
  Image<PixRGB<byte> > dispIma(4*w,2*h,ZEROS);
  inplacePaste(dispIma, ima,      Point2D<int>(  0, 0));
  inplacePaste(dispIma, dBmap,    Point2D<int>(  w, 0));
  inplacePaste(dispIma, dBmapNMS, Point2D<int>(  0, h));
  inplacePaste(dispIma, dCBEmap,  Point2D<int>(  w, h));
  inplacePaste(dispIma, dCBmap,   Point2D<int>(2*w, 0));
  inplacePaste(dispIma, itsInitialRegionColorImage, Point2D<int>(2*w, h));
  //Raster::WriteRGB(itsInitialRegionColorImage, 
  //                 std::string("starfish_FelzHutt_500.png"));

  itsWin->drawImage(dispIma, 0,0);
  Raster::waitForKey();
}

// ######################################################################
Image<float> SalientRegionSegmenter::getSalientRegion(Point2D<int> pt)
{ 
  // FIXXX: we may want to move the salient points 
  //        away from the contours
  //        This may be a problem with the orientation channel  
  //        May be scale dependent
  // from salient point get most appropriate salient region 
  // FIXXX: for orientation saliency the border is salient, 
  //        we move the point perpendicular to the area 
  //        (pick it randomly or the one with
  //        higher average salient energy) We might hit object 
  //        or background (it's ok, either way we are going to 
  //        reduce the number of unassign region
  pt = Point2D<int>(155,37); 
  LINFO("\n\n\n\n\n hard coding salient location: %d %d "
        "TAKE IT OUT LATER!\n\n\n", pt.i, pt.j);

  // Round it to the closest 8x8 grid location
  uint hstep = BOUNDARY_STEP_SIZE/2;
  uint i = ((pt.i+hstep)/BOUNDARY_STEP_SIZE) * BOUNDARY_STEP_SIZE; 
  uint j = ((pt.j+hstep)/BOUNDARY_STEP_SIZE) * BOUNDARY_STEP_SIZE; 
  itsSalientPoints.push_back(Point2D<int>(i,j));

  LINFO("adjusted point: %d %d", i,j);

  // calculate the region first
  computeSalientRegion();

  return itsSalientRegionImage;
}

// ######################################################################
void SalientRegionSegmenter::computeSalientRegion()
{ 
  // we have original resolution
  // and an 8x8 grid overlay resolution
  // --> this can be refined later or make MULTISCALE
  uint w = itsImage.getWidth();
  uint h = itsImage.getHeight();

  uint hstep = BOUNDARY_STEP_SIZE/2;
  uint wG = (w-hstep)/BOUNDARY_STEP_SIZE;
  uint hG = (h-hstep)/BOUNDARY_STEP_SIZE;

  Point2D<int> currentSalPoint = 
    itsSalientPoints[itsSalientPoints.size()-1];
  LINFO("START at salient point %d %d", 
        currentSalPoint.i, currentSalPoint.j);

  Point2D<int> currentSalPointGrid
    (currentSalPoint.i/BOUNDARY_STEP_SIZE, 
     currentSalPoint.j/BOUNDARY_STEP_SIZE);
 
  // assignment map on the grid overlay resolution
  Image<int> assignmentMap(wG, hG, ZEROS);

  // the image borders are out of bounds
  for(uint i = 0; i < wG; i++) 
    {
      assignmentMap.setVal(i,0, OTHER_OBJECT_REGION);
      assignmentMap.setVal(i,1, OTHER_OBJECT_REGION);
    }
  for(uint j = 0; j < hG; j++) 
    {
      assignmentMap.setVal(0,j, OTHER_OBJECT_REGION);
      assignmentMap.setVal(1,j, OTHER_OBJECT_REGION);      
    }
  
  for(uint i = 0; i < wG; i++) 
    {
      if(hG > 0)
        assignmentMap.setVal(i,hG-1, OTHER_OBJECT_REGION);
      if(hG > 1)
        assignmentMap.setVal(i,hG-2, OTHER_OBJECT_REGION);
    }
  for(uint j = 0; j < hG; j++) 
    {
      if(wG > 0)
        assignmentMap.setVal(0,j, OTHER_OBJECT_REGION);
      if(wG > 1)
        assignmentMap.setVal(1,j, OTHER_OBJECT_REGION);      
    }

  // list of grid nodes that are not yet assigned 
  // to either object, background, or unassigned list
  std::vector<Point2D<int> > currentObjectLocations;
  std::vector<Point2D<int> > currentObjectCoreLocations;
  std::vector<Point2D<int> > currentBackgroundLocations;
  std::vector<Point2D<int> > currentUnassignedLocations;

  // assignment map on the grid overlay resolution
  Image<int> regionAssignmentMap(w, h, ZEROS);

  // list of initial oversegmented regions that are  
  // assigned to either object or background
  std::vector<int> currentObjectRegionList;
  //std::vector<int> currentObjectCoreRegionList;
  std::vector<int> currentBackgroundRegionList;
  //std::vector<int> currentUnassignedRegionList;

  // add the start region to the core list
  currentObjectCoreLocations.push_back(currentSalPointGrid);  
  //assignmentMap.setVal(currentSalPointGrid, CENTER_CORE_REGION);

  // create boundary fragments (also call segments): 
  //   divide contour to segments of 5 edgels: about 30 - 40 pix long
  //   --> only consider contours of 5 edgels or more.
  // sort on longest contour, most similar chi-sq, closest, most salient, 
  std::list<Segment> sortedSegments = fillSegmentQueue(currentSalPoint);

  // current starting point to grow from
  Point2D<int> currPoint     = currentSalPoint;
  Point2D<int> currPointGrid = currentSalPointGrid;

  LINFO("START PT: %3d %3d --> %3d %3d", 
        currPoint.i, currPoint.j, currPointGrid.i, currPointGrid.j);

  std::vector<std::pair<Segment,Point2D<int> > > segmentPtCombo;

  // add the starting salient point 
  // to unassigned list just to start things up
  currentUnassignedLocations.push_back(currPointGrid); 
  assignmentMap.setVal(currPointGrid, UNASSIGNED_REGION);

  // while there are unassigned locations
  // around the growing center group
  while(currentUnassignedLocations.size() != 0)
    {      
      // if no segments qualified -- we break out 
      filterSegments(sortedSegments);
      if(sortedSegments.size() == 0) break;

      // pick best segment (on metric above)
      Segment segment = sortedSegments.front();
      sortedSegments.pop_front();

      // make sure we do not repeat the same combo twice
      segmentPtCombo.push_back
        (std::pair<Segment,Point2D<int> >(segment, currPoint));

      // compute histogram of distribution of points within 12pix
      Histogram hist1 = getHistogramDistribution(currPoint);      
      
      // compute distribution of both side of the boundary 
      // closer is pts2, farther is pts3
      std::vector<Point2D<int> > pts2;
      std::vector<Point2D<int> > pts3;
      getSegmentPointDistributions(segment, currPoint, pts2, pts3);
      Histogram hist2 = getHistogramDistribution(pts2);
      Histogram hist3 = getHistogramDistribution(pts3);
      
      // compute the Chi-sq difference: this is the threshold for growing
      float threshDiff = hist1.getChiSqDiff(hist2);
      LINFO("THRESH Diff: %f", threshDiff);

      LINFO("{%3d} %10.3f, %10.3f, %10.3f, %10.3f, %10.3f) = %10.3f", 
            segment.orgIndex,
            segment.lenVal, segment.simVal, segment.strVal, 
            segment.ratVal, segment.disVal, segment.pVal  );

      // grow CENTER in the direction of the boundary 
      // to connect the two 
      int index = segment.segStartIndex + MINIMUM_CONTOUR_LENGTH/2; 
      rutz::shared_ptr<Contour> contour = 
        itsContourBoundaries[segment.contourIndex];  
      rutz::shared_ptr<Edgel> edgel = contour->edgels[index];
      //Point2D<int> segCenterPoint = edgel->pt; 
      Point2D<int> segCenterPointGrid
        ((edgel->pt.i+hstep)/BOUNDARY_STEP_SIZE,
         (edgel->pt.j+hstep)/BOUNDARY_STEP_SIZE );

      std::vector<Point2D<int> > points = 
        getLine(currPointGrid, segCenterPointGrid);

      LINFO("segment points");
      for(uint i = segment.segStartIndex; i <= segment.segEndIndex; i++) 
        {
          rutz::shared_ptr<Edgel> tedgel = contour->edgels[i];
          Point2D<int> tptGrid
            ((tedgel->pt.i+hstep)/BOUNDARY_STEP_SIZE,
             (tedgel->pt.j+hstep)/BOUNDARY_STEP_SIZE );
           
          LINFO("[<<%3d>>] %d %d --> %d %d", 
                i, tptGrid.i, tptGrid.j,
                tedgel->pt.i,tedgel->pt.j);
        }

      LINFO("LINE: %3d %3d to %3d %3d", 
            currPointGrid.i, currPointGrid.j, 
            segCenterPointGrid.i, segCenterPointGrid.j);
      for(uint i = 0; i < points.size(); i++)
        {
          LINFO("[%3d]: %3d %3d", i, points[i].i, points[i].j);          
        }

      // go through each starting grow location
      Image<bool> isComparedWith(wG, hG, ZEROS);
      for(uint i = 0; i < points.size(); i++)
        {          
          Point2D<int> startPt = points[i];          
          LINFO("current starting point: %3d %3d", startPt.i, startPt.j);

          std::vector<Point2D<int> >::iterator culItr = 
            currentUnassignedLocations.begin();

          if(assignmentMap.getVal(startPt) == UNVISITED_REGION) 
            {
              // // find the current starting growing point
              // uint pos = 0;
              // for(uint i = 0; i < currentUnassignedLocations.size(); i++)
              //   if(startPt == currentUnassignedLocations[i]) 
              //     { pos = i; i = currentUnassignedLocations.size(); }
              // culItr = currentUnassignedLocations.begin();
              // currentUnassignedLocations.erase(culItr+pos); 
              
              // and move it to the front of the unassigned list
              // even if that point is not there previously
              culItr = currentUnassignedLocations.begin();
              currentUnassignedLocations.insert(culItr, startPt);
              
              assignmentMap.setVal(startPt, UNASSIGNED_REGION);
            }
          else LINFO("skip adding starting point -- already visited");
          Raster::waitForKey();

          // start from the front of the queue 
          uint currIndex = 0;

          // while there is still neighbor that we have not visited
          // as a result of the current salreg-salcontour combo
          while(currIndex < currentUnassignedLocations.size())
            {
              // get location
              Point2D<int> ptGrid = currentUnassignedLocations[currIndex];
              Point2D<int> pt(ptGrid.i*BOUNDARY_STEP_SIZE, 
                              ptGrid.j*BOUNDARY_STEP_SIZE);

              LINFO("curr pt to check: %3d %3d --> %3d %3d", 
                    ptGrid.i, ptGrid.j, pt.i, pt.j);
             
              // if not yet checked
              if(!isComparedWith.getVal(ptGrid)) 
                {
                  LINFO("was NOT compared with yet");

                  // check if the distribution in the location 
                  // goes past the threshold
                  Histogram thist = getHistogramDistribution(pt);
                  float diff = hist1.getChiSqDiff(thist);
                  LINFO("difference:(%3d %3d) & (%3d %3d): %f", 
                        currPoint.i, currPoint.j, pt.i, pt.j, diff);
                  
                  // if it is below threshold:
                  if(diff < threshDiff)
                    {
                      LINFO("pass %d", currIndex);

                      // assign current location to object list
                      currentObjectLocations.push_back(ptGrid);
                      assignmentMap.setVal(ptGrid, CENTER_REGION);
                      
                      // add its neighbors to unassigned locations list
                      uint bSize = currentUnassignedLocations.size();
                      addUnvisitedNeighbors
                        (currentUnassignedLocations, ptGrid, 
                         assignmentMap);
                      uint aSize = currentUnassignedLocations.size();
                      uint addN = aSize - bSize;

                      // take out the location from unassigned list
                      culItr = currentUnassignedLocations.begin();
                      currentUnassignedLocations.erase(culItr+currIndex+addN);

                      // LINFO("after deletion");
                      // for(uint i = 0; i < currentUnassignedLocations.size(); i++)
                      //   {
                      //     LINFO("[%3d]: %3d %3d", i, 
                      //           currentUnassignedLocations[i].i,
                      //           currentUnassignedLocations[i].j);          
                      //   }
                      // Raster::waitForKey();

                      currIndex = 0;

                      // display current situation
                      displayGrowMap(assignmentMap, currPoint,
                                     currentObjectLocations,
                                     currentObjectCoreLocations,
                                     currentBackgroundLocations,
                                     currentUnassignedLocations,
                                     isComparedWith, segment);
                    }
                  else
                    {
                      LINFO("didn't pass: %d", currIndex);

                      // otherwise just advance the index
                      currIndex++;
                    }

                  // make note that this location is already compared
                  isComparedWith.setVal(ptGrid, true);
                }
              // otherwise just advance the index
              else 
                {
                  LINFO("already checked: %d", currIndex);

                  currIndex++;
                }              
           }
        }

      LINFO("Done with segment {%3d}", segment.orgIndex);
      Raster::waitForKey();

      // assign regions next to the segment
      assignRegionsAroundSegment
        (segment, pts2, pts3, 
         assignmentMap,
         currentObjectLocations,
         currentBackgroundLocations,
         currentUnassignedLocations,
         regionAssignmentMap,
         currentObjectRegionList, 
         currentBackgroundRegionList);

      
      // GROW SURROUND also check if any of the unassigned neighbors 
      // can be added to the BG list
      // --> if closer to the other side of contour 

      // ==> note that Chi-sq thresholding need 3 regions
      //    ---> grow surround wrt/covering center 
      //         (this is Song Chun Zhu's negative cures)
      //         ---> that is, grow center and grow surround

      
      // display current situation
      LINFO("after the boundary assignment");
      displayGrowMap(assignmentMap, currPoint,
                     currentObjectLocations,
                     currentObjectCoreLocations,
                     currentBackgroundLocations,
                     currentUnassignedLocations,
                     isComparedWith, segment);

      // get the next region
      // FIXXX ADD THIS STEP LATER              




      // stop if all surroundings of the center list are assigned
      // bool allAssigned = false;
      // for(uint i = 0; i < ...; i++)
      //   {
      //     // check its neighbor if neither are
      //     // center nor surround assigned 
      //     allAssigned = true;
      //   }
      //
      // if(allAssigned) 
      //   {      
      //      ==> there is no need to go through all the boundaries 
      //      because it's already sorted
      //   }
      // else
      //   {
      //      add the neighbor that passes the threshold
      //      of the original seed to the core    
      //
      //      recalculate closest most similar boundary from regions
      //      somehow promote the ones that go through unassigned region
      //   }

      // currPoint     =
      // currPointGrid =

      // dead zone (no boundary and very non-similar colors) 
      // will be visited once core grows to it
    }

  Raster::waitForKey();
}

// ######################################################################
std::list<Segment> SalientRegionSegmenter::fillSegmentQueue
(Point2D<int> pt)
{
  std::list<Segment> sortedSegments;

  // for each contour
  int index = 0;
  for(uint i = 0; i < itsContourBoundaries.size(); i++)
    {
      rutz::shared_ptr<Contour> contour = 
        itsContourBoundaries[i];
      
      uint cLen = contour->edgels.size();

      //LINFO("Contour[%3d]: %3d", i, cLen);

      // do we take out the short contours???
      if(cLen < MINIMUM_CONTOUR_LENGTH) continue;

      // go through each contour
      float step = float(MINIMUM_CONTOUR_LENGTH)/2.0F;
      uint numSegments = uint(floor(float(cLen)/step));
      if((cLen%MINIMUM_CONTOUR_LENGTH) == 0) numSegments--;
      LINFO("Contour[%3d]: %3d: numSegments: %d", 
            i, cLen, numSegments);

      for(uint j = 0; j < numSegments; j++)
        {
          uint start = uint(j*step);
          uint end   = start + MINIMUM_CONTOUR_LENGTH - 1;

          if(end >= cLen)
            {
              end   = cLen - 1;
              start = end + 1 - MINIMUM_CONTOUR_LENGTH;
            }
         
          Segment segment(index, i, start, end);
          
          LINFO("[i: %3d c: %3d s: %3d e: %3d]", index, i, start, end);
          
          setPriority(segment, pt);
          sortedSegments.push_back(segment);

          index++;
        }
    }

  Raster::waitForKey();

  // sort the segments
  LINFO("Before Sorting");
  std::list<Segment>::iterator 
    itr  = sortedSegments.begin(),
    stop = sortedSegments.end();
  // while (itr != stop)
  //   {
  //     LINFO("%6.3f*%10.3f + %6.3f*%10.3f + %6.3f*%10.3f + "
  //           "%6.3f*%10.3f + %6.3f*%10.3f = %10.3f", 
  //           LEN_WEIGHT, (*itr).lenVal,
  //           SIM_WEIGHT, (*itr).simVal,
  //           STR_WEIGHT, (*itr).strVal,
  //           RAT_WEIGHT, (*itr).ratVal,
  //           DIS_WEIGHT, (*itr).disVal, (*itr).pVal);
  //     itr++;
  //   }
  // Raster::waitForKey();

  sortedSegments.sort();

  LINFO("After  Sorting");
  
  itr  = sortedSegments.begin();
  stop = sortedSegments.end();
  while (itr != stop)
    {
      LINFO("[%3d] %6.3f*%10.3f + %6.3f*%10.3f + %6.3f*%10.3f + "
            "%6.3f*%10.3f + %6.3f*%10.3f = %10.3f", (*itr).orgIndex,
            LEN_WEIGHT, (*itr).lenVal,
            SIM_WEIGHT, (*itr).simVal,
            STR_WEIGHT, (*itr).strVal,
            RAT_WEIGHT, (*itr).ratVal,
            DIS_WEIGHT, (*itr).disVal, (*itr).pVal);
      itr++;
    }
  Raster::waitForKey();

  return sortedSegments;
}

// ######################################################################
// FIXXX: what happen in the second reprioritization
void SalientRegionSegmenter::setPriority
(Segment &segment, Point2D<int> pt)
{
  uint w = itsImage.getWidth();
  uint h = itsImage.getHeight();
  float diag = pow(w*w + h*h, 0.5);

  // selection priority is:
  // LEN_WEIGHT * length of the contour boundary +
  // SIM_WEIGHT * similarity with the closer side +
  // STR_WEIGHT * strength of boundary +
  // RAT_WEIGHT * ratio of in:out +
  // DIS_WEIGHT * proximity to the point

  // each factor has a max of 1.0

  rutz::shared_ptr<Contour> contour = 
    itsContourBoundaries[segment.contourIndex];  
  uint cLen = contour->edgels.size();

  // length of the contour boundary 
  float tempCLen = float(cLen);
  float maxLen = float(MINIMUM_CONTOUR_LENGTH * 3); 
  if(tempCLen > maxLen) tempCLen = maxLen;
  float lenVal = tempCLen/maxLen;
  segment.lenVal = lenVal; 

  LINFO("1. clen: %d --> lenVal: %f", cLen, lenVal);

  // similarity with the closer side 
  Histogram hist1 = getHistogramDistribution(pt);
  std::vector<Point2D<int> > pts2;
  std::vector<Point2D<int> > pts3;
  getSegmentPointDistributions(segment, pt, pts2, pts3);

  Histogram hist2 = getHistogramDistribution(pts2);
  Histogram hist3 = getHistogramDistribution(pts3);
  float diff12 = hist1.getChiSqDiff(hist2);

  float tempDiff =  diff12;
  if(tempDiff > 1.0) tempDiff = 1.0; 
  float simVal = 1.0F - tempDiff;
  segment.simVal = simVal;

  LINFO("2. diff12: %f --> simVal: %f", diff12, simVal);

  // strength of boundary
  float val = 0.0; 
  uint start = segment.segStartIndex;
  uint end   = segment.segEndIndex;
  for(uint i = start; i <= end; i++)
    {
      val += contour->edgels[i]->val;
    }
  val /= float(end - start + 1);
  float maxStr = .5 * itsMaximumEdgelStrength; 
  float strVal = val/maxStr; if(strVal > 1.0F) strVal = 1.0;
  segment.strVal = strVal;

  LINFO("3. strength of boundary: %f / %f = strVal: %f", 
        val, itsMaximumEdgelStrength, strVal);
    
  // ratio of near:far side
  float diff13  =  hist1.getChiSqDiff(hist3);
  float ratio  =  MAX_FAR_TO_NEAR_RATIO; 
  if(diff12 > 0.0) ratio = diff13/diff12;
  float nratio = MAX_FAR_TO_NEAR_RATIO;
  if(ratio < MAX_FAR_TO_NEAR_RATIO) nratio = ratio; 
  float ratVal = nratio/MAX_FAR_TO_NEAR_RATIO;
  segment.ratVal = ratVal;

  LINFO("4. ratio in&out: %f/%f = %f --> ratVal: %f", 
        diff13, diff12, ratio, ratVal);

  // distance between segment and salient point
  // while NOT going through a surround area
  int index   = segment.segStartIndex + MINIMUM_CONTOUR_LENGTH/2; 
  float dist  = pt.distance(contour->edgels[index]->pt);
  float mdist = MAX_DISTANCE * diag;
  float disVal  = 1.0; 
  if(dist > mdist) disVal = 1.0F - (dist - mdist)/(diag - mdist);
  segment.disVal = disVal;

  LINFO("5. distance: dist: %f/mdist: %f --> disVal: %f", 
        dist, mdist, disVal);

  segment.pVal =  
    LEN_WEIGHT * lenVal +
    SIM_WEIGHT * simVal +
    STR_WEIGHT * strVal +
    RAT_WEIGHT * ratVal +
    DIS_WEIGHT * disVal;

  LINFO("Total: %f", segment.pVal);
  //displayDistributionSetup(segment, pt, pts2, pts3);
}

// ######################################################################
Histogram SalientRegionSegmenter::getHistogramDistribution
(Point2D<int> pt)
{
  // FIXXXX NEED TO CHECK IF ONE IS AVAILABLE (already computed)  

  Histogram lHist(NUM_L_BINS);
  Histogram aHist(NUM_A_BINS);
  Histogram bHist(NUM_L_BINS);

  int rad = NEIGHBORHOOD_RADIUS;
  int numPoints = 0;
  for(int i = -rad; i <= rad; i++) 
    for(int j = -rad; j <= rad; j++) 
      {
        bool isWithinCircle = ((i*i+j*j) <= (rad*rad));

        Point2D<int> pt2 = pt + Point2D<int>(i,j);
        if(itsImage.coordsOk(pt2) && isWithinCircle)
          {
            // L component
            int l_val = itsImageHistogramEntries[0].getVal(pt2);
            lHist.addValue(l_val, 1.0F);

            // a component
            int a_val = itsImageHistogramEntries[1].getVal(pt2);
            aHist.addValue(a_val, 1.0F);

            // b component
            int b_val = itsImageHistogramEntries[2].getVal(pt2);
            bHist.addValue(b_val, 1.0F);

            numPoints++;
          }
      }
  //LINFO("num points: %d", numPoints);

  // maybe later want to add the weight for each channel
  //(Point2D<int> pt, std::vector<float> &dist, float weight)

  // smooth the histogram individually
  lHist.smooth(itsLKernel);
  aHist.smooth(itsAKernel);
  bHist.smooth(itsBKernel);

  // combine the histograms
  std::vector<Histogram> histograms;
  histograms.push_back(lHist);
  histograms.push_back(aHist);
  histograms.push_back(bHist);
  Histogram histogram(histograms);
  
  // normalize with the number of points
  histogram.divide(numPoints);

  return histogram;
}

// ######################################################################
Histogram SalientRegionSegmenter::getHistogramDistribution
(std::vector<Point2D<int> > pts)
{
  Histogram lHist(NUM_L_BINS);
  Histogram aHist(NUM_A_BINS);
  Histogram bHist(NUM_L_BINS);

  int numPoints = 0;
  for(uint i = 0; i < pts.size(); i++) 
      {
        Point2D<int> pt = pts[i];
        if(itsImage.coordsOk(pt))
          {
            // L component
            int l_val = itsImageHistogramEntries[0].getVal(pt);
            lHist.addValue(l_val, 1.0F);

            // a component
            int a_val = itsImageHistogramEntries[1].getVal(pt);
            aHist.addValue(a_val, 1.0F);

            // b component
            int b_val = itsImageHistogramEntries[2].getVal(pt);
            bHist.addValue(b_val, 1.0F);

            numPoints++;
          }
      }

  LINFO("num points: %d", numPoints);
  
  // maybe later want to add the weight for each channel
  //(Point2D<int> pt, std::vector<float> &dist, float weight)

  // smooth the histogram individually
  lHist.smooth(itsLKernel);
  aHist.smooth(itsAKernel);
  bHist.smooth(itsBKernel);

  // combine the histograms
  std::vector<Histogram> histograms;
  histograms.push_back(lHist);
  histograms.push_back(aHist);
  histograms.push_back(bHist);
  Histogram histogram(histograms);
  
  // normalize with the number of points
  histogram.divide(numPoints);

  return histogram;
}

// ######################################################################
void SalientRegionSegmenter::getSegmentPointDistributions
(Segment segment, Point2D<int> pt, 
 std::vector<Point2D<int> > &pts1, std::vector<Point2D<int> > &pts2)
{
  // get the contour divider end points 
  uint start = segment.segStartIndex;
  uint end   = segment.segEndIndex;

  rutz::shared_ptr<Contour> contour = 
    itsContourBoundaries[segment.contourIndex];  

  rutz::shared_ptr<Edgel> edgels = contour->edgels[start];
  rutz::shared_ptr<Edgel> edgele = contour->edgels[end  ];
           
  uint  ainds  = edgels->angleIndex; 
  float bainds = 
    fmod((ainds+(NUM_RIDGE_DIRECTIONS/2)),NUM_RIDGE_DIRECTIONS);
  uint  ainde  = edgele->angleIndex; 
  float bainde = 
    fmod((ainde+(NUM_RIDGE_DIRECTIONS/2)),NUM_RIDGE_DIRECTIONS);

  int step  = BOUNDARY_STEP_SIZE;
  int hstep = step/2;
  float dxs = cos(bainds * M_PI/4.0) * hstep; 
  float dys = sin(bainds * M_PI/4.0) * hstep;

  float dxe = cos(bainde * M_PI/4.0) * hstep; 
  float dye = sin(bainde * M_PI/4.0) * hstep;

  Point2D<float> p1 = edgels->pt + Point2D<float>(-dxs-.5, -dys-.5);
  Point2D<float> p2 = edgele->pt + Point2D<float>( dxe+.5,  dye+.5); 
  
  // estimate the contour average angle
  float ang  = atan2(p2.j-p1.j, p2.i-p1.i);

  // the corrected angle would be to move it to PI/2
  // FIXXX: this could be incorrect 
  //        if we branch the contour to both left and right
  float cAng = M_PI/2.0F - ang;
  if(ang < 0.0F && ang > M_PI) 
    { LFATAL("out of Bounds Angle -- Fix it"); }

  if(cAng == 0.0F)
    {
      LINFO("Deal with ZERO corrected angle");
      Raster::waitForKey();
    }

  float cCAng = cos(cAng);
  float sCAng = sin(cAng);

  // get the midpoint of the angle line estimator
  Point2D<float> mid((p1.i+p2.i)/2.0, (p1.j+p2.j)/2.0);  

  // LINFO("c[%3d]:(%3d:a: %f)-(%3d:a: %3f) ", 
  //       segment.contourIndex, start, bainds, end, bainde);
  // LINFO("hstep: %d [%f %f] [%f %f]", 
  //       hstep, dxs, dys, dxe, dye);
  
  // LINFO("p1:(%f %f) p2:(%f %f) ", 
  //       p1.i, p1.j, p2.i, p2.j);
  // LINFO("ang: %f --> cAng: %f (c: %f s: %f) mid[%f %f] ", 
  //       ang, cAng, cCAng, sCAng, mid.i, mid.j);


  // rotate about that point to vertical
  float tbR = 0.0; float bbR = 0.0; 
  float lbR = 0.0; float rbR = 0.0;
  std::vector<Point2D<float> > sPt;
  std::vector<Point2D<float> > ePt;
  std::vector<Point2D<float> > srPt;
  std::vector<Point2D<float> > erPt;
  std::vector<Point2D<float> > mrPt;
  for(uint k = start; k <= end; k++)
    {
      Point2D<int> pt = contour->edgels[k]->pt;
      float i = pt.i - mid.i;
      float j = pt.j - mid.j; 
      float iR =  i*cCAng - j*sCAng;
      float jR =  i*sCAng + j*cCAng;

      uint aind   = contour->edgels[k]->angleIndex; 
      float baind = 
        fmod((aind+(NUM_RIDGE_DIRECTIONS/2)),NUM_RIDGE_DIRECTIONS);
      float dx = cos(baind * M_PI/4.0) * hstep; 
      float dy = sin(baind * M_PI/4.0) * hstep;

      Point2D<float> ps = pt + Point2D<float>(-dx-.5, -dy-.5);
      Point2D<float> pe = pt + Point2D<float>( dx+.5,  dy+.5); 
      
      float si = ps.i - mid.i;
      float sj = ps.j - mid.j; 
      float siR =  si*cCAng - sj*sCAng;
      float sjR =  si*sCAng + sj*cCAng;
      
      float ei = pe.i - mid.i;
      float ej = pe.j - mid.j; 
      float eiR =  ei*cCAng - ej*sCAng;
      float ejR =  ei*sCAng + ej*cCAng;

      sPt.push_back(ps);
      ePt.push_back(pe);
      srPt.push_back(Point2D<float>(siR,sjR));
      erPt.push_back(Point2D<float>(eiR,ejR));
      mrPt.push_back(Point2D<float>( iR, jR));

      // initialize everything on the first iteration
      if(k == start)
        {
          tbR = sjR;
          bbR = sjR;
          lbR = siR;
          rbR = siR;
        }

      if(tbR > sjR) tbR = sjR;
      if(tbR > ejR) tbR = ejR;

      if(bbR < sjR) bbR = sjR;
      if(bbR < ejR) bbR = ejR;

      if(lbR > siR) lbR = siR;
      if(lbR > eiR) lbR = eiR;

      if(rbR < siR) rbR = siR;
      if(rbR < eiR) rbR = eiR;

      // LINFO("[[%d]]", k);
      // LINFO("M[%3d %3d] --> offset[%f %f] --> [%f %f]", 
      //       pt.i, pt.j, i ,j, iR, jR);

      // LINFO("S[%f %f] --> offset[%f %f] --> [%f %f]", 
      //       ps.i, ps.j, si ,sj, siR, sjR);
      // LINFO("E[%f %f] --> offset[%f %f] --> [%f %f] \n", 
      //       pe.i, pe.j, ei ,ej, eiR, ejR);
    }

  //LINFO("[t: %f b: %f l: %f r: %f] ", tbR, bbR, lbR, rbR);

  // count the pixels in the left and right within the bounding box  
  float areaL = 0; float areaR = 0;
  float w = rbR - lbR;
  float h = bbR - tbR;
  if(w > 0.0F && h > 0.0F)
    {
      for(uint k = start; k <= end; k++)
        {
          int i = k-start;
          // we will use a trapezoid area 
          // to approximate number of pixels
          
          float tl = rbR - srPt[i].i;
          float bl = rbR - erPt[i].i;
          
          float h = erPt[i].j - srPt[i].j;
          
          float trapArea = (tl+bl)*h/2.0F;
          float boxArea  = w*h;
          areaL += (boxArea-trapArea);
          areaR += trapArea;
          
          // LINFO("tl = %f = %f - %f||| bl = %f = %f - %f ",
          //       tl, rbR, srPt[i].i, bl, rbR, erPt[i].i);
          // LINFO("h = %f = %f - %f | w = %f = %f - %f",
          //       h, erPt[i].j, srPt[i].j, w, rbR, lbR);

          // LINFO("[%d] b: %f b-t: %f, t: %f", 
          //       k, boxArea, boxArea-trapArea, trapArea);
        }
    }
  // LINFO("areaL: %f, areaR: %f ==> h: %f", areaL, areaR, h);

  // extend the box so that the number of pixels on both sides 
  // are close to the area circle with radius = NEIGHBORHOOD_RADIUS
  uint nrad = NEIGHBORHOOD_RADIUS;
  float minArea = nrad * nrad * M_PI; 
  float extL = 0.0; float extR = 0.0;

  //LINFO("lbR: %f areaL: %f", lbR, areaL);
  //LINFO("rbR: %f areaR: %f", rbR, areaR);
  if(areaL < minArea)
    {
      extL = lbR - (minArea - areaL)/h;
    }

  if(areaR < minArea)
    {
      extR = rbR + (minArea - areaR)/h;
    }

  //LINFO("add L:%f ->  %f add R: %f -> %f", 
  //      (minArea - areaL), (minArea - areaL)/h,
  //      (minArea - areaR), (minArea - areaR)/h );

  //LINFO("Sanity Check: min Area: %f, Area: ((%f) - (%f))*%f = %f", 
  //      minArea, extR, extL, h, (extR - extL)*h);

  // figure out which side the pt is
  // FIXXX: may want to be more careful here
  int pti = pt.i - mid.i;
  int ptj = pt.j - mid.j; 
  float ptiR =  pti*cCAng - ptj*sCAng;
  //float ptjR =  pti*sCAng + ptj*cCAng;
  float isLeft = true; if(ptiR > 0.0) isLeft = false;
  //LINFO("isLeft: %d", int(isLeft));

  // create a box
  uint finalW = extR-extL;
  uint finalH = h;
  uint maxDim = finalW;
  if(finalH > finalW) maxDim = finalH;
  Image<byte> temp1(maxDim*2, maxDim*2, ZEROS);
  Image<byte> temp2(maxDim*2, maxDim*2, ZEROS);

  // LINFO("finalW: %d finalH: %d ==> maxDim: %d", finalW, finalH, maxDim);
  //uint ttt1 = 0; uint ttt2 = 0;
  for(uint k = start; k <= end; k++)
    {
      uint ind = k - start;

      float t  = srPt[ind].j; float b  =  erPt[ind].j;      
      float h  = b - t;

      float tH = srPt[ind].i; float bH =  erPt[ind].i;
      float w  = bH - tH;

      //LINFO("ind: %d h: %f = %f - %f w: %f = %f - %f", 
      //      ind, h, b, t, w, bH, tH);

      // FIXXX: there may be a vertical gap between indexes

      // get the pixels on the left and right side
      // LINFO("VERT range: %d to %d", int(t), int(b));
      //uint tt1 = 0; uint tt2 = 0;
      for(int j = int(t); j <= int(b + .5F); j++)
        {
          // compute the dividing edge location
          float div = (w/h*(j - b)) + tH;
          if(w < 0.0F) div = (w/h*(j - t)) + tH;

          float leftE  = extL+maxDim;
          float rightE = extR+maxDim;

          // LINFO("horizontal divider: %f, range: %f to %f", 
          //      div+maxDim, leftE, rightE);

          //uint t1 = 0; uint t2 = 0;
          uint jShift = j + maxDim;
          for(int i = int(leftE); i <= int(rightE+.5F); i++) 
            {              
              if(i < (div+maxDim)) 
                {
                  if(isLeft) temp1.setVal(i,jShift, 128); //t1++; }
                  else       temp2.setVal(i,jShift, 128); //t2++; }
                }
              else
                {
                  if(isLeft) temp2.setVal(i,jShift, 128); //t2++; }
                  else       temp1.setVal(i,jShift, 128); //t1++; }
                }

              //LINFO("i:%d %d: div: %f --> %d == %3d %3d", 
              //      i, jShift, div+maxDim, int(i<(div+maxDim)), t1, t2);
            }
          //LINFO("[%3d %3d] t1: %d t2: %d", k, j, t1, t2);
          //tt1 += t1;
          //tt2 += t2;
        }
      //LINFO("[%3d] t1: %d t2: %d", k, tt1, tt2);
      //ttt1 += tt1;
      //ttt2 += tt2;
    }
  //LINFO("ttt1: %d ttt2: %d", ttt1, ttt2);

  //LINFO("sum1: %f sum2: %f", 
  //      float(sum(temp1))/128.0F, float(sum(temp2))/128.0F);

  // re-rotate the bounding box
  Image<byte> temp1R = rotate(temp1, maxDim, maxDim, -cAng);
  Image<byte> temp2R = rotate(temp2, maxDim, maxDim, -cAng);

  //LINFO("md[%d] mid: %f %f", maxDim, mid.i, mid.j);

  // figure out which side is next to the point
  int tw = temp1.getWidth();
  int th = temp1.getHeight();
  for(int i = 0; i < tw; i++)
    for(int j = 0; j < th; j++)
      {
        if(temp1R.getVal(i,j) > 64)
          { 
            Point2D<int> pt(i+int(mid.i)-maxDim,
                            j+int(mid.j)-maxDim );
            pts1.push_back(pt); 
            //LINFO("1 %d %d --> %d %d", i,j,pt.i, pt.j);
            //Raster::waitForKey();
          }

        if(temp2R.getVal(i,j) > 64) 
          {
            Point2D<int> pt(i+int(mid.i)-maxDim,
                            j+int(mid.j)-maxDim );
            pts2.push_back(pt); 
            //LINFO("2 %d %d --> %d %d", i,j,pt.i, pt.j);
            //Raster::waitForKey();
          }
      }

  //LINFO("Final Size: %d %d", int(pts1.size()), int(pts2.size()));
}

// ######################################################################
void SalientRegionSegmenter::filterSegments(std::list<Segment> &segments)
{
  // FIXXX ==> what happens if  best segment is not good enough???
  //  --> need a threshold if the best one is worth pursuing.
  //  --> if not then stop
  //      --> the object is segmented perfectly by Graph-based segmented
}

// ######################################################################
void SalientRegionSegmenter::addUnvisitedNeighbors
(std::vector<Point2D<int> > &currentUnassignedLocations, 
 Point2D<int> currPointGrid, Image<int> &assignmentMap)
{
  //--> make sure no duplicates 
  //    the grid area is unassigned to the list,
  //    as well as surround or previous object
  //--> this will keep growing maximally 
  //    on the current salreg-salcontour combo              

  // LINFO("BEFORE unassigned");
  // for(uint i = 0; i < currentUnassignedLocations.size(); i++)
  //   {
  //     LINFO("[%3d]: %3d %3d", i, 
  //           currentUnassignedLocations[i].i,
  //           currentUnassignedLocations[i].j);          
  //   }
  // Raster::waitForKey();


  // add to the front of the list

  for(int j = 1; j >= -1; j--)
    for(int i = 1; i >= -1; i--)
      if(!(i == 0 && j== 0))
        {
          Point2D<int> pt(currPointGrid.i+i, currPointGrid.j+j);
          
          // if the location is not yet visited
          if(assignmentMap.coordsOk(pt) && 
             assignmentMap.getVal(pt) == UNVISITED_REGION) 
            {
              LINFO("pt(%3d %3d): %d", pt.i,pt.j, 
                    assignmentMap.getVal(pt));

              std::vector<Point2D<int> >::iterator itr = 
                currentUnassignedLocations.begin();
              currentUnassignedLocations.insert(itr, pt);

              assignmentMap.setVal(pt, UNASSIGNED_REGION);
            }
        }

  // LINFO("AFTER  unassigned");
  // for(uint i = 0; i < currentUnassignedLocations.size(); i++)
  //   {
  //     LINFO("[%3d]: %3d %3d", i, 
  //           currentUnassignedLocations[i].i,
  //           currentUnassignedLocations[i].j);          
  //   }
  // Raster::waitForKey();

}

// ######################################################################
std::vector<Point2D<int> > SalientRegionSegmenter::getLine
(Point2D<int> p1, Point2D<int> p2)
{
  std::vector<Point2D<int> > points;

  // ray tracing algorithm
  // from Graphics Gems / Paul Heckbert
  int dx = p2.i - p1.i, ax = abs(dx) << 1, sx = signOf(dx);
  int dy = p2.j - p1.j, ay = abs(dy) << 1, sy = signOf(dy);
  int x = p1.i, y = p1.j;

  const int w = itsImage.getWidth();
  const int h = itsImage.getHeight();

  if (ax > ay)
    {
      int d = ay - (ax >> 1);
      for (;;)
        {
          // assignment
          if (x >= 0 && x < w && y >= 0 && y < h)
            points.push_back(Point2D<int>(x,y));

          if (x == p2.i) return points;
          if (d >= 0) { y += sy; d -= ax; }
          x += sx; d += ay;
        }
    }
  else
    {
      int d = ax - (ay >> 1);
      for (;;)
        {
          // assignment
          if (x >= 0 && x < w && y >= 0 && y < h)
            points.push_back(Point2D<int>(x,y));

          if (y == p2.j) return points;
          if (d >= 0) { x += sx; d -= ay; }
          y += sy; d += ax;
        }
    }
  
  return points;
}

// ######################################################################
void SalientRegionSegmenter::assignRegionsAroundSegment
(Segment segment, 
 std::vector<Point2D<int> > pts2, 
 std::vector<Point2D<int> > pts3, 
 Image<int>                 &assignmentMap,
 std::vector<Point2D<int> > &currentObjectLocations,
 std::vector<Point2D<int> > &currentBackgroundLocations,
 std::vector<Point2D<int> > &currentUnassignedLocations,
 Image<int>       &regionAssignmentMap,
 std::vector<int> &currentObjectRegionList, 
 std::vector<int> &currentBackgroundRegionList)
{
  uint w = itsImage.getWidth();
  uint h = itsImage.getHeight();

  uint hstep = BOUNDARY_STEP_SIZE/2;
  //uint wG = (w-hstep)/BOUNDARY_STEP_SIZE;
  //uint hG = (h-hstep)/BOUNDARY_STEP_SIZE;

  Image<int> image(w,h, ZEROS);
  for(uint i = 0; i < pts2.size(); i++) 
    {
      if(image.coordsOk(pts2[i]))
        image.setVal(pts2[i], CLOSER_SIDE);
    }
  for(uint i = 0; i < pts3.size(); i++) 
    {
      if(image.coordsOk(pts3[i]))
        image.setVal(pts3[i], FAR_SIDE);
    }
  std::vector<std::vector<Point2D<int> > > 
    neighbor(NUM_RIDGE_DIRECTIONS);

  for(uint k = 0; k < NUM_RIDGE_DIRECTIONS; k++)
    {      
      neighbor[k] = std::vector<Point2D<int> >();

      if(k == 0)
        {
          neighbor[k].push_back(Point2D<int>( 1,  0));
          neighbor[k].push_back(Point2D<int>(-1,  0));
        }
      else if(k == 1)
        {
          neighbor[k].push_back(Point2D<int>(-1, -1));
          neighbor[k].push_back(Point2D<int>( 0, -1));
          neighbor[k].push_back(Point2D<int>(-1,  0));
          neighbor[k].push_back(Point2D<int>( 1,  0));
          neighbor[k].push_back(Point2D<int>( 0,  1));
          neighbor[k].push_back(Point2D<int>( 1,  1));
        }
      else if(k == 2)
        {
          neighbor[k].push_back(Point2D<int>( 0,  1));
          neighbor[k].push_back(Point2D<int>( 0, -1));
        }
      else if(k == 3)
        {
          neighbor[k].push_back(Point2D<int>( 0, -1));
          neighbor[k].push_back(Point2D<int>( 1, -1));
          neighbor[k].push_back(Point2D<int>(-1,  0));
          neighbor[k].push_back(Point2D<int>( 1,  0));
          neighbor[k].push_back(Point2D<int>(-1,  1));
          neighbor[k].push_back(Point2D<int>( 0,  1));
        }
    }
  
  // go through all the edgels' grid location first
  for(uint e = segment.segStartIndex; e <= segment.segEndIndex; e++) 
    {
      rutz::shared_ptr<Contour> contour = 
        itsContourBoundaries[segment.contourIndex];  
      rutz::shared_ptr<Edgel> edgel = contour->edgels[e];           
      Point2D<int> pt = edgel->pt;
      Point2D<int> ptGrid
        ((pt.i+hstep)/BOUNDARY_STEP_SIZE,
         (pt.j+hstep)/BOUNDARY_STEP_SIZE );      
      LINFO("[[[%3d]]] %d %d --> %d %d", e, pt.i, pt.j, ptGrid.i, ptGrid.j);
      
      uint label = assignmentMap.getVal(ptGrid);

      // do nothing on center, surround, other objects 
      if(label == CENTER_REGION)       LINFO(" C_REGION : do nothing");
      if(label == SURROUND_REGION)     LINFO(" S_REGION : do nothing");
      if(label == OTHER_OBJECT_REGION) LINFO("OO_REGION : do nothing");
      
      if(label == UNVISITED_REGION)
        {
          // add to background
          currentBackgroundLocations.push_back(ptGrid);
          assignmentMap.setVal(ptGrid, SURROUND_REGION);         

          LINFO("UNVISITED adding %3d %3d to the background", 
                ptGrid.i, ptGrid.j);
        }

      if(label == UNASSIGNED_REGION)
        {

          // LINFO("before changing on edgel UNASSIGNED_REGION");
          // for(uint i = 0; i < currentUnassignedLocations.size(); i++)
          //   {
          //     LINFO("[%3d]: %3d %3d", i, 
          //           currentUnassignedLocations[i].i,
          //           currentUnassignedLocations[i].j);          
          //   }

          // LINFO("Background");
          // for(uint i = 0; i < currentBackgroundLocations.size(); i++)
          //   {
          //     LINFO("[%3d]: %3d %3d", i, 
          //           currentBackgroundLocations[i].i,
          //           currentBackgroundLocations[i].j);          
          //   }

          // Raster::waitForKey();
          

          // find it is in the unassigned list
          uint pos = 0;
          for(uint i = 0; i < currentUnassignedLocations.size(); i++)
            if(ptGrid == currentUnassignedLocations[i]) 
              { pos = i; i = currentUnassignedLocations.size(); }

          // delete it
          std::vector<Point2D<int> >::iterator itr = 
            currentUnassignedLocations.begin();
          currentUnassignedLocations.erase(itr+pos); 
          
          // add to background list
          currentBackgroundLocations.push_back(ptGrid);
          assignmentMap.setVal(ptGrid, SURROUND_REGION);    

          // LINFO("after changing on edgel");
          // for(uint i = 0; i < currentUnassignedLocations.size(); i++)
          //   {
          //     LINFO("[%3d]: %3d %3d", i, 
          //           currentUnassignedLocations[i].i,
          //           currentUnassignedLocations[i].j);          
          //   }


          // LINFO("Background");
          // for(uint i = 0; i < currentBackgroundLocations.size(); i++)
          //   {
          //     LINFO("[%3d]: %3d %3d", i, 
          //           currentBackgroundLocations[i].i,
          //           currentBackgroundLocations[i].j);          
          //   }

          // Raster::waitForKey();
          
      
        }
    }

  // LINFO("before checking edgel neighbors");
  // for(uint i = 0; i < currentUnassignedLocations.size(); i++)
  //   {
  //     LINFO("[%3d]: %3d %3d", i, 
  //           currentUnassignedLocations[i].i,
  //           currentUnassignedLocations[i].j);          
  //   }
  
  
  // LINFO("Background");
  // for(uint i = 0; i < currentBackgroundLocations.size(); i++)
  //   {
  //     LINFO("[%3d]: %3d %3d", i, 
  //           currentBackgroundLocations[i].i,
  //           currentBackgroundLocations[i].j);          
  //   }
  
  // Raster::waitForKey();
 
  // then go through its neighbors
  for(uint e = segment.segStartIndex; e <= segment.segEndIndex; e++) 
    {
      rutz::shared_ptr<Contour> contour = 
        itsContourBoundaries[segment.contourIndex];  
      rutz::shared_ptr<Edgel> edgel = contour->edgels[e];           
      Point2D<int> pt = edgel->pt;
      Point2D<int> ptGrid
        ((pt.i+hstep)/BOUNDARY_STEP_SIZE,
         (pt.j+hstep)/BOUNDARY_STEP_SIZE );
      
      uint dir = edgel->angleIndex; 

      LINFO("[<<%3d>>] %d %d --> %d %d (%d)", 
            e,  pt.i, pt.j, ptGrid.i, ptGrid.j, dir);

      // classify current segment grid location as       
      for(uint i = 0; i < neighbor[dir].size(); i++)
        {
          Point2D<int> ptGrid2 = ptGrid + neighbor[dir][i];
          Point2D<int> pt2(ptGrid2.i*BOUNDARY_STEP_SIZE,
                           ptGrid2.j*BOUNDARY_STEP_SIZE );

          LINFO("    N[%3d]: %d %d --> %d %d", 
                i, ptGrid2.i, ptGrid2.j, pt2.i, pt2.j);

          if(!assignmentMap.coordsOk(ptGrid2))
            {
              LINFO("out of bounds"); continue;
            }

          uint label2 = assignmentMap.getVal(ptGrid2);
          LINFO("LABEL2: %d", label2);

          // do nothing on center, surround, other objects 
          if(label2 == CENTER_REGION)       LINFO(" C_REGION : do nothing");
          if(label2 == SURROUND_REGION)     LINFO(" S_REGION : do nothing");
          if(label2 == OTHER_OBJECT_REGION) LINFO("OO_REGION : do nothing");
          
          if(label2 == UNVISITED_REGION)
            {
              // assign depending on which side of boundary
              int side = image.getVal(pt2);

              // LINFO("UNVISITED Before Background");
              // for(uint i = 0; i < currentBackgroundLocations.size(); i++)
              //   {
              //     LINFO("[%3d]: %3d %3d", i, 
              //           currentBackgroundLocations[i].i,
              //           currentBackgroundLocations[i].j);          
              //   }
              // for(uint i = 0; i < currentObjectLocations.size(); i++)
              //   {
              //     LINFO("[%3d]: %3d %3d", i, 
              //           currentObjectLocations[i].i,
              //           currentObjectLocations[i].j);          
              //   }


              if(side == CLOSER_SIDE)
                {
                  // add to background
                  currentObjectLocations.push_back(ptGrid2);
                  assignmentMap.setVal(ptGrid2, CENTER_REGION);

                  LINFO("  adding %3d %3d to the center", 
                        ptGrid2.i, ptGrid2.j);
                }
              else if(side == FAR_SIDE)
                {
                  // add to background
                  currentBackgroundLocations.push_back(ptGrid2);
                  assignmentMap.setVal(ptGrid2, SURROUND_REGION);

                  LINFO("  adding %3d %3d to the background", 
                        ptGrid2.i, ptGrid2.j);
                }

              // LINFO("UNVISITED after  Background");
              // for(uint i = 0; i < currentBackgroundLocations.size(); i++)
              //   {
              //     LINFO("[%3d]: %3d %3d", i, 
              //           currentBackgroundLocations[i].i,
              //           currentBackgroundLocations[i].j);          
              //   }
              // for(uint i = 0; i < currentObjectLocations.size(); i++)
              //   {
              //     LINFO("[%3d]: %3d %3d", i, 
              //           currentObjectLocations[i].i,
              //           currentObjectLocations[i].j);          
              //   }

            }
          
          if(label2 == UNASSIGNED_REGION)
            {
              // check if point is one either side of boundary          
              int side = image.getVal(pt2);
              LINFO("UN Side: %d", side);

              // LINFO("UNASSIGNED neighbors before");
              // for(uint i = 0; i < currentUnassignedLocations.size(); i++)
              //   {
              //     LINFO("[%3d]: %3d %3d", i, 
              //           currentUnassignedLocations[i].i,
              //           currentUnassignedLocations[i].j);          
              //   }

              if(side == CLOSER_SIDE || side == FAR_SIDE)
                {                
                  // find the location in the unassigned list
                  uint pos = 0;
                  for(uint i = 0; i < currentUnassignedLocations.size(); i++)
                    if(ptGrid2 == currentUnassignedLocations[i]) 
                      { pos = i; i = currentUnassignedLocations.size(); }
                  
                  // delete it
                  std::vector<Point2D<int> >::iterator itr = 
                    currentUnassignedLocations.begin();
                  currentUnassignedLocations.erase(itr+pos); 
                }
              else LINFO("not on the close or far side");

              // LINFO("UNASSIGNED neighbors after ");
              // for(uint i = 0; i < currentUnassignedLocations.size(); i++)
              //   {
              //     LINFO("[%3d]: %3d %3d", i, 
              //           currentUnassignedLocations[i].i,
              //           currentUnassignedLocations[i].j);          
              //   }

              // LINFO("UNASSIGNED before Background");
              // for(uint i = 0; i < currentBackgroundLocations.size(); i++)
              //   {
              //     LINFO("[%3d]: %3d %3d", i, 
              //           currentBackgroundLocations[i].i,
              //           currentBackgroundLocations[i].j);          
              //   }
              // for(uint i = 0; i < currentObjectLocations.size(); i++)
              //   {
              //     LINFO("[%3d]: %3d %3d", i, 
              //           currentObjectLocations[i].i,
              //           currentObjectLocations[i].j);          
              //   }


              if(side == CLOSER_SIDE)
                {
                  // add to background
                  currentObjectLocations.push_back(ptGrid2);
                  assignmentMap.setVal(ptGrid2, CENTER_REGION);

                  LINFO("  adding %3d %3d to the center", 
                        ptGrid2.i, ptGrid2.j);
                }
              else if(side == FAR_SIDE)
                {
                  // add to background
                  currentBackgroundLocations.push_back(ptGrid2);
                  assignmentMap.setVal(ptGrid2, SURROUND_REGION);

                  LINFO("  adding %3d %3d to the background", 
                        ptGrid2.i, ptGrid2.j);
                }


              // LINFO("UNASSIGNED after  Background");
              // for(uint i = 0; i < currentBackgroundLocations.size(); i++)
              //   {
              //     LINFO("[%3d]: %3d %3d", i, 
              //           currentBackgroundLocations[i].i,
              //           currentBackgroundLocations[i].j);          
              //   }
              // for(uint i = 0; i < currentObjectLocations.size(); i++)
              //   {
              //     LINFO("[%3d]: %3d %3d", i, 
              //           currentObjectLocations[i].i,
              //           currentObjectLocations[i].j);          
              //   }

            }
        }
    }

  // get the oversegmented regions
  
  // assign regions next to the segment      
  


  // // region index from the underlying Felz-Hutt region
  // int regIndex = itsInitialRegionImage.getVal(currentSalPoint); 

  // FIXXX: LATER include superpixel underneath the region
  // calculate the overlap on the region: 
  // anything over 50% of total bounding box is good
   // --> ok because the region is similar appearance-wise 


  // =============================================
      
  //currentObjectRegList.push_back(regIndex);
  //for(uint i = 0; i < itsRegionAdjecencyList[regIndex].size(); i++)
  //  currentUnassignedRegionsRegList.push_back
  //    (itsRegionAdjecencyList[regIndex][i]);




      // // find the current starting growing point
      // uint pos = 0;
      // for(uint i = 0; i < currentUnassignedLocations.size(); i++)
      //   if(startPt == currentUnassignedLocations[i]) 
      //     { pos = i; i = currentUnassignedLocations.size(); }
      // culItr = currentUnassignedLocations.begin();
      // currentUnassignedLocations.erase(culItr+pos); 
      
      // and move it to the front of the unassigned list
      // even if that point is not there previously

              // culItr = currentUnassignedLocations.begin();
              // currentUnassignedLocations.insert(culItr, startPt);
 

              //         // assign current location to object list
              //         currentObjectLocations.push_back(ptGrid);
              //         assignmentMap.setVal(ptGrid, CENTER_REGION);
                      
              //         // add its neighbors to unassigned locations list
              //         uint bSize = currentUnassignedLocations.size();
              //         addUnvisitedNeighbors
              //           (currentUnassignedLocations, ptGrid, 
              //            assignmentMap);
              //         uint aSize = currentUnassignedLocations.size();
              //         uint addN = aSize - bSize;

              //         // take out the location from unassigned list
              //         culItr = currentUnassignedLocations.begin();
              //         currentUnassignedLocations.erase(culItr+currIndex+addN);




}

// ######################################################################
void SalientRegionSegmenter::displayDistributionSetup
(Segment segment, Point2D<int> pt, 
 std::vector<Point2D<int> > pts2, 
 std::vector<Point2D<int> > pts3)
{
  // display stuff
  Image<PixRGB<byte> > ima = itsImage;
  uint w = ima.getWidth();
  uint h = ima.getHeight();

  Image<float> fIma(luminance(ima));
  Image<float> tempFIma = fIma;
  inplaceNormalize(tempFIma, 0.0f,255.0f);  
  Image<byte>  bIma(tempFIma); 

  // get the image
  float mVal = 32;
  float bVal = 255 - mVal;
  Image<byte> dImaR, dImaG, dImaB;
  getComponents(ima, dImaR, dImaG, dImaB);
  inplaceNormalize(dImaR, byte(0), byte(mVal));
  inplaceNormalize(dImaG, byte(0), byte(mVal));
  inplaceNormalize(dImaB, byte(0), byte(mVal));
  Image<PixRGB<byte> > dIma  = makeRGB(dImaR,dImaG,dImaB);

  Image<PixRGB<byte> > dispIma(4*w,2*h,ZEROS);

  inplacePaste (dispIma, ima, Point2D<int>(0,0));

  Image<float> stateMapR(w,h, ZEROS);
  Image<float> stateMapG(w,h, ZEROS);
  Image<float> stateMapB(w,h, ZEROS);

  // draw core 
  std::vector<Point2D<int> > cores;
  std::vector<Point2D<int> > center;
  std::vector<Point2D<int> > surround;
  std::vector<Point2D<int> > unassigned;

 cores.push_back(Point2D<int>(152, 40));

 // // unassigned.push_back(Point2D<int>(144, 32));
 // // unassigned.push_back(Point2D<int>(152, 32));
 // // unassigned.push_back(Point2D<int>(160, 32));

 // // unassigned.push_back(Point2D<int>(144, 40));
 // // unassigned.push_back(Point2D<int>(160, 40));

 // // unassigned.push_back(Point2D<int>(144, 48));
 // // unassigned.push_back(Point2D<int>(152, 48));
 // // unassigned.push_back(Point2D<int>(160, 48));

 // cores.push_back(Point2D<int>(144, 32));
 // cores.push_back(Point2D<int>(152, 32));
 // cores.push_back(Point2D<int>(160, 32));

 // cores.push_back(Point2D<int>(144, 40));
 // cores.push_back(Point2D<int>(160, 40));

 // cores.push_back(Point2D<int>(144, 48));
 // cores.push_back(Point2D<int>(152, 48));
 // cores.push_back(Point2D<int>(160, 48));

 // center.push_back(Point2D<int>(136, 24));
 // center.push_back(Point2D<int>(144, 24));
 // center.push_back(Point2D<int>(152, 24));
 // center.push_back(Point2D<int>(136, 32));
 // center.push_back(Point2D<int>(136, 40));

 // center.push_back(Point2D<int>(128, 16));
 // center.push_back(Point2D<int>(136, 16));
 // center.push_back(Point2D<int>(128, 24));
 // center.push_back(Point2D<int>(128, 32));

 // center.push_back(Point2D<int>(120,  8));
 // center.push_back(Point2D<int>(120, 16));
 // center.push_back(Point2D<int>(120, 24));

 // center.push_back(Point2D<int>( 96,  8));
 // center.push_back(Point2D<int>(104,  8));
 // center.push_back(Point2D<int>(112,  8));
 // center.push_back(Point2D<int>(104, 16));
 // center.push_back(Point2D<int>(112, 16));

 // unassigned.push_back(Point2D<int>( 88,  8));
 // surround  .push_back(Point2D<int>( 88, 16));
 // unassigned.push_back(Point2D<int>( 96, 16));
 // surround  .push_back(Point2D<int>( 96, 24));
 // unassigned.push_back(Point2D<int>(104, 24));

 // unassigned.push_back(Point2D<int>(112, 24));
 // surround  .push_back(Point2D<int>(112, 32));
 // unassigned.push_back(Point2D<int>(120, 32));
 // surround  .push_back(Point2D<int>(120, 40));
 // unassigned.push_back(Point2D<int>(128, 40));

 // surround  .push_back(Point2D<int>(128, 48));
 // unassigned.push_back(Point2D<int>(136, 48));
 // surround  .push_back(Point2D<int>(136, 56));
 // unassigned.push_back(Point2D<int>(144, 56));
 // unassigned.push_back(Point2D<int>(152, 56));
 // unassigned.push_back(Point2D<int>(160, 56));
 // unassigned.push_back(Point2D<int>(168, 56));

 // unassigned.push_back(Point2D<int>(128,  8));
 // unassigned.push_back(Point2D<int>(136,  8));

 // unassigned.push_back(Point2D<int>(144, 16));
 // unassigned.push_back(Point2D<int>(152, 16));
 // unassigned.push_back(Point2D<int>(160, 16));

 // unassigned.push_back(Point2D<int>(160, 24));
 // unassigned.push_back(Point2D<int>(168, 24));

 // unassigned.push_back(Point2D<int>(168, 32));
 // unassigned.push_back(Point2D<int>(168, 40));
 // unassigned.push_back(Point2D<int>(168, 48));

  for(uint i = 0; i < cores.size(); i++)
    {
      Rectangle r = Rectangle::tlbrO
        (cores[i].j-4, cores[i].i-4, 
         cores[i].j+4, cores[i].i+4);
      drawFilledRect(stateMapB, r, 0.5F);
    }

  for(uint i = 0; i < center.size(); i++)
    {
      Rectangle r = Rectangle::tlbrO
        (center[i].j-4, center[i].i-4, 
         center[i].j+4, center[i].i+4);
      drawFilledRect(stateMapR, r, 0.5F);
    }

  for(uint i = 0; i < surround.size(); i++)
    {
      Rectangle r = Rectangle::tlbrO
        (surround[i].j-4, surround[i].i-4, 
         surround[i].j+4, surround[i].i+4);
      drawFilledRect(stateMapG, r, 0.5F);
    }

  for(uint i = 0; i < unassigned.size(); i++)
    {
      Rectangle r = Rectangle::tlbrO
        (unassigned[i].j-4, unassigned[i].i-4, 
         unassigned[i].j+4, unassigned[i].i+4);
      drawFilledRect(stateMapR, r, 0.5F);
      drawFilledRect(stateMapG, r, 0.5F);
    }


  // draw grid
  int grad = 8;
  for(uint i = grad; i < w; i++)
    for(uint j = grad; j < h; j++)
      {
        if(i%grad == 0 && j%grad == 0)
          {
            drawCross(stateMapR, Point2D<int>(i,j), 0.5F, 1);
            drawCross(stateMapG, Point2D<int>(i,j), 0.5F, 1);
            drawCross(stateMapB, Point2D<int>(i,j), 0.5F, 1);
          }
      }

  // draw closer population
  //LINFO("pts2 size: %d", int(pts2.size()));
  for(uint i = 0; i < pts2.size(); i++)
    {
      //LINFO("pts2[%3d]: %3d %3d", i, pts2[i].i, pts2[i].j);
      if(stateMapR.coordsOk(pts2[i]))
        {
          stateMapR.setVal(pts2[i], 0.25F);
          //stateMapG.setVal(pts2[i], 0.25F);
          //stateMapB.setVal(pts2[i], 0.25F);
        }
    }

  // draw farther population
  //LINFO("pts3 size: %d", int(pts3.size()));
  for(uint i = 0; i < pts3.size(); i++)
    {
      //LINFO("pts3[%3d]: %3d %3d", i, pts3[i].i, pts3[i].j);
      if(stateMapR.coordsOk(pts3[i]))
        {
          //stateMapR.setVal(pts3[i], 0.25F);
          stateMapG.setVal(pts3[i], 0.25F);
          //stateMapB.setVal(pts3[i], 0.25F);
        }
    } 
  
  // draw the point 
  int rad = NEIGHBORHOOD_RADIUS;
  for(int i = -rad; i <= rad; i++) 
    for(int j = -rad; j <= rad; j++) 
      {  
        bool isWithinCircle = ((i*i+j*j) <= (rad*rad));

        Point2D<int> pt2 = pt + Point2D<int>(i,j);
        if(itsImage.coordsOk(pt2) && isWithinCircle)
          {
            stateMapR.setVal(pt2, 0.5F);
            //stateMapG.setVal(pt2, 0.25F);
            //stateMapB.setVal(pt2, 0.25F);
          }
      }

  drawDisk(stateMapR, pt, 2,  1.0F);
  drawDisk(stateMapG, pt, 2,  1.0F);
  drawDisk(stateMapB, pt, 2,  1.0F);

  // drawLine(stateMapR, pt, Point2D<int>(90,14),  0.5F);
  // drawLine(stateMapG, pt, Point2D<int>(90,14),  0.5F);
  // drawLine(stateMapB, pt, Point2D<int>(90,14),  0.5F);
  
  // draw the contour 
  rutz::shared_ptr<Contour> contour = 
    itsContourBoundaries[segment.contourIndex];  

  int hstep = BOUNDARY_STEP_SIZE/2;
  uint start = segment.segStartIndex;
  uint end   = segment.segEndIndex;
  //for(uint i = start; i <= end; i++)
  //  {
  for(uint i = 0; i < contour->edgels.size(); i++)
    {
      rutz::shared_ptr<Edgel> edgel = contour->edgels[i];
          
      uint aind = edgel->angleIndex; 
      float baind = 
        fmod((aind+(NUM_RIDGE_DIRECTIONS/2)),NUM_RIDGE_DIRECTIONS);
      
      float dx = cos(baind * M_PI/4.0) * hstep; 
      float dy = sin(baind * M_PI/4.0) * hstep;
      
      Point2D<int> pt = edgel->pt;
      Point2D<int> p1 = pt + Point2D<int>( dx+.5,  dy+.5); 
      Point2D<int> p2 = pt + Point2D<int>(-dx-.5, -dy-.5);
      
      if(i >= start && i <= end) 
        {
          drawLine(stateMapR, p1, p2, 1.0F);
          drawLine(stateMapG, p1, p2, 1.0F);
          drawLine(stateMapB, p1, p2, 1.0F);
        }
      else
        {
          drawLine(stateMapR, p1, p2, 0.5F);
          drawLine(stateMapG, p1, p2, 0.5F);
          drawLine(stateMapB, p1, p2, 0.5F);
        }
      //drawDisk(stateMap, pt, 2,  1.0F);
    }

  inplaceNormalize(stateMapR, 0.0f,bVal);
  inplaceNormalize(stateMapG, 0.0f,bVal);
  inplaceNormalize(stateMapB, 0.0f,bVal);
  Image<byte> dSMapR(stateMapR);
  Image<byte> dSMapG(stateMapG);
  Image<byte> dSMapB(stateMapB);
  Image<PixRGB<byte> > dSmap = makeRGB(dSMapR, dSMapG, dSMapB);
  inplacePaste
    (dispIma, Image<PixRGB<byte> >(dIma+dSmap), Point2D<int>(w,0));

  // Raster::WriteRGB(crop(Image<PixRGB<byte> >(dIma+dSmap), 
  //                       Point2D<int>(50,0), Dims(150,100)),
  //                  std::string("starfish_algorithm_4.pnm"));

  itsWin->drawImage(dispIma, 0,0);
  Raster::waitForKey();
}

// ######################################################################
void SalientRegionSegmenter::displayGrowMap
  (Image<int> assignmentMap, Point2D<int> currPoint,
   std::vector<Point2D<int> > currentObjectLocations,
   std::vector<Point2D<int> > currentObjectCoreLocations,
   std::vector<Point2D<int> > currentBackgroundLocations,
   std::vector<Point2D<int> > currentUnassignedLocations,
   Image<bool> isComparedWith, Segment segment)
{
  // display stuff
  Image<PixRGB<byte> > ima = itsImage;
  uint w = ima.getWidth();
  uint h = ima.getHeight();

  Image<float> fIma(luminance(ima));
  Image<float> tempFIma = fIma;
  inplaceNormalize(tempFIma, 0.0f,255.0f);  
  Image<byte>  bIma(tempFIma); 

  // get the image
  float mVal = 32;
  float bVal = 255 - mVal;
  Image<byte> dImaR, dImaG, dImaB;
  getComponents(ima, dImaR, dImaG, dImaB);
  inplaceNormalize(dImaR, byte(0), byte(mVal));
  inplaceNormalize(dImaG, byte(0), byte(mVal));
  inplaceNormalize(dImaB, byte(0), byte(mVal));
  Image<PixRGB<byte> > dIma  = makeRGB(dImaR,dImaG,dImaB);

  Image<PixRGB<byte> > dispIma(4*w,2*h,ZEROS);

  inplacePaste (dispIma, ima, Point2D<int>(0,0));

  Image<float> stateMapR(w,h, ZEROS);
  Image<float> stateMapG(w,h, ZEROS);
  Image<float> stateMapB(w,h, ZEROS);

  std::vector<Point2D<int> >::iterator oitr = 
    currentObjectLocations.begin();
  std::vector<Point2D<int> >::iterator citr = 
    currentObjectCoreLocations.begin();
  std::vector<Point2D<int> >::iterator sitr = 
    currentBackgroundLocations.begin();
  std::vector<Point2D<int> >::iterator uitr = 
    currentUnassignedLocations.begin();
  
  for(uint i = 0; i < currentObjectLocations.size(); i++)
    {
      Point2D<int> pt((*oitr).i*BOUNDARY_STEP_SIZE,
                      (*oitr).j*BOUNDARY_STEP_SIZE );
      Rectangle r = Rectangle::tlbrO(pt.j-4, pt.i-4, pt.j+4, pt.i+4);
      drawFilledRect(stateMapR, r, 0.5F);
      oitr++;
    }

  for(uint i = 0; i < currentObjectCoreLocations.size(); i++)
    {
      Point2D<int> pt((*citr).i*BOUNDARY_STEP_SIZE,
                      (*citr).j*BOUNDARY_STEP_SIZE );
      Rectangle r = Rectangle::tlbrO(pt.j-4, pt.i-4, pt.j+4, pt.i+4);
      drawFilledRect(stateMapB, r, 0.5F);
      citr++;
    }

  for(uint i = 0; i < currentBackgroundLocations.size(); i++)
    {
      Point2D<int> pt((*sitr).i*BOUNDARY_STEP_SIZE,
                      (*sitr).j*BOUNDARY_STEP_SIZE );
      Rectangle r = Rectangle::tlbrO(pt.j-4, pt.i-4, pt.j+4, pt.i+4);
      drawFilledRect(stateMapG, r, 0.5F);
      sitr++;
    }

  for(uint i = 0; i < currentUnassignedLocations.size(); i++)
    {
      Point2D<int> pt((*uitr).i*BOUNDARY_STEP_SIZE,
                      (*uitr).j*BOUNDARY_STEP_SIZE );
      Rectangle r = Rectangle::tlbrO(pt.j-4, pt.i-4, pt.j+4, pt.i+4);
      drawFilledRect(stateMapR, r, 0.5F);
      drawFilledRect(stateMapG, r, 0.5F);
      uitr++;
    }

  // draw grid
  int grad = 8;
  for(uint i = grad; i < w; i++)
    for(uint j = grad; j < h; j++)
      {
        if(i%grad == 0 && j%grad == 0)
          {
            drawCross(stateMapR, Point2D<int>(i,j), 0.5F, 1);
            drawCross(stateMapG, Point2D<int>(i,j), 0.5F, 1);
            drawCross(stateMapB, Point2D<int>(i,j), 0.5F, 1);
          }
      }
  
  // draw the contour 
  rutz::shared_ptr<Contour> contour = 
    itsContourBoundaries[segment.contourIndex];  

  int hstep = BOUNDARY_STEP_SIZE/2;
  uint start = segment.segStartIndex;
  uint end   = segment.segEndIndex;
  //for(uint i = start; i <= end; i++)
  //  {
  for(uint i = 0; i < contour->edgels.size(); i++)
    {
      rutz::shared_ptr<Edgel> edgel = contour->edgels[i];
          
      uint aind = edgel->angleIndex; 
      float baind = 
        fmod((aind+(NUM_RIDGE_DIRECTIONS/2)),NUM_RIDGE_DIRECTIONS);
      
      float dx = cos(baind * M_PI/4.0) * hstep; 
      float dy = sin(baind * M_PI/4.0) * hstep;
      
      Point2D<int> pt = edgel->pt;
      Point2D<int> p1 = pt + Point2D<int>( dx+.5,  dy+.5); 
      Point2D<int> p2 = pt + Point2D<int>(-dx-.5, -dy-.5);
      
      if(i >= start && i <= end) 
        {
          drawLine(stateMapR, p1, p2, 1.0F);
          drawLine(stateMapG, p1, p2, 1.0F);
          drawLine(stateMapB, p1, p2, 1.0F);
        }
      else
        {
          drawLine(stateMapR, p1, p2, 0.5F);
          drawLine(stateMapG, p1, p2, 0.5F);
          drawLine(stateMapB, p1, p2, 0.5F);
        }
      //drawDisk(stateMap, pt, 2,  1.0F);
    }

  inplaceNormalize(stateMapR, 0.0f,bVal);
  inplaceNormalize(stateMapG, 0.0f,bVal);
  inplaceNormalize(stateMapB, 0.0f,bVal);
  Image<byte> dSMapR(stateMapR);
  Image<byte> dSMapG(stateMapG);
  Image<byte> dSMapB(stateMapB);
  Image<PixRGB<byte> > dSmap = makeRGB(dSMapR, dSMapG, dSMapB);
  inplacePaste
    (dispIma, Image<PixRGB<byte> >(dIma+dSmap), Point2D<int>(w,0));

  // Raster::WriteRGB(crop(Image<PixRGB<byte> >(dIma+dSmap), 
  //                       Point2D<int>(50,0), Dims(150,100)),
  //                  std::string("starfish_algorithm_4.pnm"));

  itsWin->drawImage(dispIma, 0,0);
  Raster::waitForKey();
}

// ######################################################################
/* So things look consistent in everyone's emacs... */
/* Local Variables: */
/* indent-tabs-mode: nil */
/* End: */