gevd_edgel_regions.cxx

Go to the documentation of this file.

// This is gel/gevd/gevd_edgel_regions.cxx
#include "gevd_edgel_regions.h"
//:
// \file

#include <vcl_iostream.h>
#include <vcl_cstdlib.h> // for vcl_abs(int)
#include <vcl_cassert.h>
#include <vcl_algorithm.h> // for vcl_sort() and vcl_find()

#include <vtol/vtol_intensity_face.h>

#include <vxl_config.h>
#include <vil1/vil1_memory_image_of.h>
#include <vcl_cmath.h>     // for sqrt()
#include <vul/vul_timer.h>

#include <vsol/vsol_box_2d.h>
#include <vsol/vsol_box_2d_sptr.h>
#include <vtol/vtol_edge.h>
#include <vtol/vtol_edge_2d.h>
#include <vtol/vtol_one_chain.h>
#include <vtol/vtol_cycle_processor.h>

#include <vdgl/vdgl_digital_curve.h>
#include <vdgl/vdgl_edgel_chain.h>
#include <vdgl/vdgl_interpolator.h>


#define bytePixel(buf,x,y)   (*((unsigned char*)buf->GetElementAddr(x,y)))
#define shortPixel(buf,x,y)  (*((short*)buf->GetElementAddr(x,y)))

#if 0 // not ported to vxl
//----------------------------------------------------------------
//: cache debug info
void edgel_regions::cache_bad_edges(CoolArrayP<Edge*>& bad_edges)
{
  vcl_vector<Edge*> corrupt_edges;
  for (CoolArrayP<Edge*>::iterator eit = bad_edges.begin();
       eit != bad_edges.end(); eit++)
    corrupt_edges.push_back(*eit);
  corrupt_edges_.push_back(corrupt_edges);
}

void edgel_regions::cache_bad_edges(vcl_vector<Edge*>& bad_edges)
{
  vcl_vector<Edge*> corrupt_edges;
  for (vcl_vector<Edge*>::iterator eit = bad_edges.begin();
       eit != bad_edges.end(); eit++)
    corrupt_edges.push_back(*eit);
  corrupt_edges_.push_back(corrupt_edges);
}

void edgel_regions::cache_bad_verts(CoolArrayP<Vertex*>& bad_verts)
{
  vcl_vector<Vertex*> corrupt_verts;
  for (CoolArrayP<Vertex*>::iterator vit = bad_verts.begin();
       vit != bad_verts.end(); vit++)
    corrupt_verts.push_back(*vit);
  corrupt_vertices_.push_back(corrupt_verts);
}
#endif

//--------------------------------------------------------
//: compare function to be used with CoolList<float>::sort() [increasing]
static int increasing_compare(unsigned int const&a, unsigned int const&b)
{
  if (a > b) return 1;
  if (a < b) return -1;
  return 0;
}

//Print the region label array
void gevd_edgel_regions::print_region_array()
{
  vcl_cout << vcl_endl;
  for (unsigned int y = yo_; y<=yend_; y++)
  {
    vcl_cout << vcl_endl; // << setw(2);
    for (unsigned int x = xo_; x<=xend_; x++)
      if (region_label_array_[Y(y)][X(x)]==EDGE)
          //&&edge_boundary_array_[Y(y)][X(x)]->IsVertex()
        vcl_cout << "* ";
      else
        vcl_cout << region_label_array_[Y(y)][X(x)] << ' ';
  }
  vcl_cout << "\n\n\n";
}

//Print the contents of the forward equivalence index
void gevd_edgel_regions::print_region_equivalence()
{
  vcl_cout << vcl_endl << "Label Equivalence:\n"
           << "----------------\n";
  vcl_map<unsigned int, vcl_vector<unsigned int>*>::iterator rpf_iterator;
  for (rpf_iterator= region_pairs_forward_.begin();
       rpf_iterator!=region_pairs_forward_.end(); rpf_iterator++)
  {
    vcl_cout << (*rpf_iterator).first << " == "
             << (*rpf_iterator).second << vcl_endl;
  }
  vcl_cout << vcl_endl;
}

//Print the contents of the reverse equivalence index
void gevd_edgel_regions::print_reverse_region_equivalence()
{
  vcl_cout << vcl_endl << "Reverse Label Equivalence:\n"
           << "----------------\n";
  vcl_map<unsigned int, vcl_vector<unsigned int>*>::iterator rpf_iterator;
  for (rpf_iterator= region_pairs_reverse_.begin();
       rpf_iterator!=region_pairs_reverse_.end(); rpf_iterator++)
    {
      vcl_cout << (*rpf_iterator).first << " == "
               << (*rpf_iterator).second << vcl_endl;
    }
  vcl_cout << vcl_endl;
}

//Print the reduced equivalence relation
void gevd_edgel_regions::print_base_equivalence()
{
  vcl_cout << vcl_endl << "Base Label Equivalence:\n"
           << "----------------\n";

  for (unsigned int i = min_region_label_; i<max_region_label_; i++)
    vcl_cout << i << " == " << this->BaseLabel(i) << vcl_endl;
}

//Print the fitted intensity data for all faces
void gevd_edgel_regions::print_intensity_data()
{
  for (vcl_vector<vtol_intensity_face_sptr>::iterator fit =faces_->begin();
       fit != faces_->end(); fit++)
    {
#if 0
    vtol_intensity_face* f = *fit;
      vcl_cout << "IntFaceAt(" << f->Xo() << ' ' << f->Yo() << "):\n";
      f->PrintFit();
#endif
    }
}

//--------------------------------------------------------------
//: The region array can be offset from the corner of the ROI.
//    The following two methods transform from the ROI coordinate
//    system to the region label array coordinate system
unsigned int gevd_edgel_regions::X(unsigned int x)
{
  return x-xo_;
}

unsigned int gevd_edgel_regions::Y(unsigned int y)
{
  return y-yo_;
}

//---------------------------------------------------------
//: Casts the float x location of a point to an unsigned-int location
unsigned int gevd_edgel_regions::Xf(float x)
{
  unsigned int xu = ((unsigned int)(x));
  return xu;
}

//---------------------------------------------------------
//: Casts the float x location of a point to an unsigned-int location
unsigned int gevd_edgel_regions::Yf(float y)
{
  unsigned int yu = ((unsigned int)(y));
  return yu;
}

//----------------------------------------------------------
//: Default constructor
gevd_edgel_regions::gevd_edgel_regions(bool debug)
{
  verbose_ = false;
  debug_ = debug;
#if 0
  if (debug_)
    debug_data_ = new topo_debug_data;
  else
    debug_data_ = 0;
#endif
    image_ = NULL;
  image_source_ = false;
  buf_source_ = false;
  buf_ = NULL;
  edge_boundary_array_ = NULL;
  region_label_array_ = NULL;
  xo_=0;
  yo_=0;
  xend_ = 0;
  yend_ = 0;
  min_region_label_ = LABEL;
  max_region_label_ = LABEL;
  faces_ = new vcl_vector<vtol_intensity_face_sptr>;
  face_edge_index_ = NULL;
  intensity_face_index_ = NULL;
  failed_insertions_ = new vcl_vector<vtol_edge_2d_sptr>;
  ubuf_ = NULL;
  sbuf_ = NULL;
}

//----------------------------------------------------------
//: Default destructor
gevd_edgel_regions::~gevd_edgel_regions()
{
  if (image_) {
    image_ = NULL;
  }

  unsigned int y;
  if (region_label_array_)
    for (y=yo_; y<=yend_; y++)
      delete [] region_label_array_[Y(y)];
  delete [] region_label_array_;

#if 0
  for (region_pairs_reverse_.reset(); region_pairs_reverse_.next();)
    delete region_pairs_reverse_.value();

  for (region_pairs_forward_.reset(); region_pairs_forward_.next();)
    delete region_pairs_forward_.value();
  for (equivalence_set_.reset(); equivalence_set_.next();)
    delete equivalence_set_.value();

  for (region_edges_.reset(); region_edges_.next();)
  {
    region_edges_.value()->UnProtect();
  }

   if (edge_boundary_array_)
     for (y = yo_; y<=yend_; y++)
     {
       for (x = xo_; x<=xend_; x++)
         if (edge_boundary_array_[Y(y)][X(x)])
           edge_boundary_array_[Y(y)][X(x)]->UnProtect();
       delete [] edge_boundary_array_[Y(y)];
     }
   delete [] edge_boundary_array_;


  for (region_edge_adjacency_.reset(); region_edge_adjacency_.next();)
    delete region_edge_adjacency_.value();

  for (faces_->reset(); faces_->next();)
    faces_->value()->UnProtect();
  delete faces_;

  if (intensity_face_index_)
  {
    for (unsigned int i = 0; i<max_region_label_; i++)
      if (intensity_face_index_[i])
        intensity_face_index_[i]= NULL;
    delete [] intensity_face_index_;
  }
#endif


  if (face_edge_index_)
    {
      for (unsigned int i = 0; i<max_region_label_; i++)
        delete face_edge_index_[i];
      delete [] face_edge_index_;
    }
  delete failed_insertions_;
  delete [] ubuf_;
  delete [] sbuf_;
}

bool gevd_edgel_regions::compute_edgel_regions(gevd_bufferxy* buf,
                                               vcl_vector<vtol_edge_2d_sptr>& sgrp,
                                               vcl_vector<vtol_intensity_face_sptr>& faces)
{
  buf_ = buf;
  image_= NULL;
  buf_source_=true;
  image_source_=false;
  return compute_edgel_regions(sgrp, faces);
}

//-----------------------------------------------------------------
//: The key process loop.
//  Carries out the steps:
//   1) Connected components 2)Edge-label assignment 3)Collect region
//   boundaries 4) Construct vtol_intensity_faces 5)Calculate intensity fit
bool
gevd_edgel_regions::compute_edgel_regions(vil1_image* image,
                                          vcl_vector<vtol_edge_2d_sptr>& sgrp,
                                          vcl_vector<vtol_intensity_face_sptr>& faces)
{
  image_ = image;
  buf_ = NULL;
  image_source_=true;
  buf_source_ = false;
  return compute_edgel_regions(sgrp, faces);
}

bool gevd_edgel_regions::compute_edgel_regions(vcl_vector<vtol_edge_2d_sptr>& sgrp,
                                               vcl_vector<vtol_intensity_face_sptr>& faces)
{
//    corrupt_edges_.clear();
//    corrupt_vertices_.clear();
#if 0
  if (debug_data_)
  {
    debug_data_->clear();
    if (buf_source_)
      debug_data_->set_buf(buf_);
  }
#endif
  vul_timer t;
  if (!sgrp.size())
    return false;

  //Set up the region label array with edge boundaries
  this->InitRegionArray(sgrp);
  if (verbose_)
    this->print_region_array();
  //Propagate connected components
  unsigned int y, x;
  for (y=yo_; y<yend_; y++)
    for (x=xo_; x<xend_; x++)
      this->UpdateConnectedNeighborhood(X(x), Y(y));
  if (verbose_)
    this->print_region_array();
  //Resolve region label equivalence
  this->GrowEquivalenceClasses();
  this->PropagateEquivalence();
  if (verbose_)
    this->print_base_equivalence();
  this->ApplyRegionEquivalence();
  if (verbose_)
    this->print_region_array();

  //Associate Edge(s) with region labels
  for (y=yo_; y<yend_; y++)
    for (x=xo_; x<xend_; x++)
      this->AssignEdgeLabels(X(x), Y(y));

  //Collect Edge(s) bounding each region
  this->CollectEdges();
  vcl_cout << "Propagate Regions and Collect Edges("
           << max_region_label_-min_region_label_ << ") in "
           << t.real() << " msecs.\n";

  //Collect Face-Edge associations
  this->CollectFaceEdges();

  //Construct vtol_intensity_faces
  this->ConstructFaces();
  if (!faces_||!faces_->size())
    return false;

  //Collect intensity data for each region
  this->InsertFaceData();

  if (verbose_)
    this->print_intensity_data();
  //Output the result
  faces.clear();
  for (vcl_vector<vtol_intensity_face_sptr>::iterator fit = faces_->begin();
       fit != faces_->end(); fit++)
  {
    faces.push_back(*fit);
    //        faces_->value()->Protect(); //This caused a big leak
  }
  vcl_cout << "Compute Edgel Regions Total(" << sgrp.size() << ") in "
           << t.real() << " msecs.\n";
  return true;
}


//----------------------------------------------------------
//: assign equivalence of region label b to region label a
//
bool gevd_edgel_regions::InsertRegionEquivalence(unsigned int label_b,
                                                 unsigned int label_a)
{
  bool result = true;
  result = result && this->add_to_forward(label_b, label_a);
  result = result && this->add_to_reverse(label_a, label_b);
  return result;
}

//--------------------------------------------------------------
//: Get the most basic label equivalent to a given label.
//    If the label_map_ is not defined, this function uses the
//    forward label hash table.  Otherwise 0 is returned.
unsigned int gevd_edgel_regions::BaseLabel(unsigned int label)
{
  vcl_map<unsigned int, unsigned int >::iterator lmtest;
  lmtest = label_map_.find(label);
  if ( lmtest != label_map_.end() )
    return lmtest->second;
  else
    return 0;
}

//------------------------------------------------------------
//:
//  Paint the edgels into the region label array and then
//    output an image where the value is 255 if the pixel is an
//    edge, 0 otherwise
vil1_image* gevd_edgel_regions::GetEdgeImage(vcl_vector<vtol_edge_2d_sptr>& sg)
{
  if (!this->InitRegionArray(sg)) return NULL;
  int sizex = this->GetXSize(), sizey = this->GetYSize();
  unsigned char no_edge = 0, edge = 255;
#if 0
  ImageTemplate temp(sizex, sizey, 8, 32, 32);
  MemoryImage* image = new MemoryImage(&temp);
  for (int y = 0; y<sizey; y++)
    for (int x = 0; x<sizex; x++)
      if (region_label_array_[y][x] == EDGE)
        image->PutPixel(&edge, x, y);
      else
        image->PutPixel(&no_edge, x, y);
#endif
  vil1_memory_image_of<vxl_byte> * image = new vil1_memory_image_of<vxl_byte>(sizex,sizey);

  for (int y = 0; y<sizey; y++)
    for (int x = 0; x<sizex; x++)
      if (region_label_array_[y][x] == EDGE)
        (*image)[y][x] = edge;
      else
        (*image)[y][x] = no_edge;
  return image;
}

//-----------------------------------------------------------
//: Populate the label_map_ to reflect the equivalences between labels.
//
void gevd_edgel_regions::PropagateEquivalence()
{
  vcl_map<unsigned int, vcl_vector<unsigned int>* >::iterator  esi;

  for (esi = equivalence_set_.begin(); esi != equivalence_set_.end(); esi++)
  {
    unsigned int base = esi->first;
    vcl_vector<unsigned int>* labels = esi->second;
    if (!labels) continue;
    int len = labels->size();
    if (!len) continue;
    for (int i = 0; i<len; i++)
      label_map_[(*labels)[i]] = base;
  }
  for (unsigned int i = min_region_label_; i<max_region_label_; i++)
    if (label_map_.find(i) == label_map_.end())
      label_map_[i] = i;
}

//---------------------------------------------------------------------
//:
//  Find the set of labels equivalent to label from a given hash table
//    and merge into the appropriate equivalence set. cur_label is the
//    equivalence set being formed. label is the table key of equivalences
//    to be merged.
bool gevd_edgel_regions::
merge_equivalence(vcl_map<unsigned int, vcl_vector<unsigned int>* >& tab,
                  unsigned int cur_label,
                  unsigned int label)
{
  vcl_map<unsigned int, vcl_vector<unsigned int>* >::iterator hashi;

  hashi = tab.find(label);
  if (hashi == tab.end()) {
  return false;
  }
  vcl_vector<unsigned int>* labels = hashi->second;

  int len = labels->size();
  if (!len) {
    return false;
  }
  vcl_vector<unsigned int>* array;

  hashi = equivalence_set_.find(cur_label);
  if ( hashi == equivalence_set_.end())
  {
    array = new vcl_vector<unsigned int>;
    equivalence_set_[cur_label] = array;
  }
  else
  {
    array = hashi->second;
  }
  for (int i = 0; i<len; i++)
  {
    unsigned int l = (*labels)[i];
    bool found = false;
    for (unsigned int j=0 ; j < array->size() ; j++)
    {
      if ((*array)[j] == l) found = true;
    }
    if (!found)
    {
      array->push_back(l);
    }
  }

  return true;
}

//----------------------------------------------------------------
//: Find the next label not accounted for in the current equivalence set.
//    The set of labels is searched to find a label larger than label, but
//    not in the set, labels.
bool gevd_edgel_regions::get_next_label(vcl_vector<unsigned int>* labels,
                                        unsigned int& label)
{
  unsigned int tmp = label+1;
  if (!labels)
    if (tmp<max_region_label_)
    {
      label = tmp;
      return true;
    }
  for (unsigned int i = tmp; i<max_region_label_; i++)
  {
    bool found = false;
    for ( unsigned int j = 0 ; j < labels->size() ; j++ )
    {
      if ((*labels)[j] == i)
      {
        found = true;
      }
    }
    if (!found)
    {
      label = i;
      return true;
    }
  }
  return false;
}

//----------------------------------------------------------------
//: Form equivalence classes by transitive closure on each label.
//    The idea is to add labels to the set, cur_set, by equivalence until no
//    new equivalence groups can be found.  The current equivalence class,
//    cur_set is repeatedly scanned when new labels are added to pick up
//    new equivalence groups.  Old equivalence table entries are removed
//    as they are used.  When the equivalence class corresponding to cur_lable
//    is completely closed, then a new label not in a previously formed
//    equivalence class is used to seed a new equivalence class.
void gevd_edgel_regions::GrowEquivalenceClasses()
{
  if ((max_region_label_-min_region_label_) < 2)
    return;
  unsigned int cur_label = min_region_label_;
  vcl_map<unsigned int, vcl_vector<unsigned int>* >::iterator mei;
  while (true)
    {
      bool merging = true;
      vcl_vector<unsigned int>* cur_set = NULL;
      unsigned int i = cur_label;
      int len = 0;
      while (merging)
      {
        int old_len = len;
        merging = false;
        bool find_forward =
          this->merge_equivalence(region_pairs_forward_, cur_label, i);
        if (find_forward)
          region_pairs_forward_.erase(i);
        bool find_reverse =
          this->merge_equivalence(region_pairs_reverse_, cur_label, i);
        if (find_reverse)
          region_pairs_reverse_.erase(i);
        mei = equivalence_set_.find(cur_label);
        if (mei == equivalence_set_.end()) continue;
        cur_set = equivalence_set_[cur_label];
        if (!cur_set) continue;
        len = cur_set->size();
        if (!len) continue;
        if (len > old_len)
        {
          i = cur_label;
          vcl_sort(cur_set->begin(), cur_set->end(), &increasing_compare);
#if 0
          cur_set->sort(increasing_compare);
          //-tpk- need to convert
#endif
        }

        for (int j = 0; j<len&&!merging; j++)
          if ((*cur_set)[j]>i)
          {
            i = (*cur_set)[j];
            merging = true;
          }
      }
      if (!get_next_label(cur_set, cur_label)) return;
    }
}

//------------------------------------------------------------
//: Check if the SpatialGroup contains Edge(s)
//
bool gevd_edgel_regions::GroupContainsEdges(vcl_vector<vtol_edge_2d_sptr>& /*sg*/)
{
  vcl_cerr << "gevd_edgel_regions::GroupContainsEdges() not yet implemented\n";
#if 0
  vcl_string type(sg.GetSpatialGroupName());
  return type == vcl_string("EdgelGroup") ||
         type == vcl_string("FittedEdgeGroup");
#endif
  return true; // TODO
}

//: Advance along a line and generate contiguous pixels on the line.
//  This function assures that the spatial digitization of the Edge(s)
//  produces completely connected boundaries.
static bool line_gen(float xs, float ys, float xe, float ye,
                     unsigned int& x, unsigned int& y)
{
  assert(xs > 0.0f); assert(ys > 0.0f);
  const float pix_edge = 1.0f; //We are working at scale = 1.0
  static float xi, yi;
  static bool init = true;
  static bool done = false;
  if (init)
    {
      xi = xs;
      yi = ys;
      x = (unsigned int)(xi/pix_edge);
      y = (unsigned int)(yi/pix_edge);
      init = false;
      return true;
    }
  if (done) return false;
  float dx = xe-xs;
  float dy = ye-ys;
  float mag = vcl_sqrt(dx*dx + dy*dy);
  if (mag<pix_edge)//Can't reach the next pixel under any circumstances
    {             //so just output the target, xe, ye.
      x = (unsigned int)xe; y = (unsigned int)ye;
      done = true;
      return true;
    }
  float delta = (0.5f*pix_edge)/mag; //move in 1/2 pixel increments
  //Previous pixel location
  int xp = int(xi/pix_edge);
  int yp = int(yi/pix_edge);
  //Increment along the line until the motion is greater than one pixel
  for (int i = 0; i<3; i++)
    {
      xi += dx*delta;
      yi += dy*delta;
      //Check for end of segment, make sure we emit the end of the segment
      if ((xe>=xs&&xi>xe)||(xe<=xs&&xi<xe)||(ye>=ys&&yi>ye)||(ye<=ys&&yi<ye))
        {
          x = (unsigned int)xe; y = (unsigned int)ye;
          done = true;
          return true;
        }
      //Check if we have advanced by more than .5 pixels
      x = (unsigned int)(xi/pix_edge);
      y = (unsigned int)(yi/pix_edge);
      if (vcl_abs(int(x)-xp)>(.5*pix_edge)||vcl_abs(int(y)-yp)>(.5*pix_edge))
        return true;
    }
  vcl_cout << "in gevd_edgel_regions line_gen: - shouldn't happen\n";
  return false;
}

//----------------------------------------------------------
//: A utility for inserting an edgel into the region_label_array_.
//    An edgel and a previous edgel in the chain are used to interpolate
//    intermediate edgels to take account of pixel quantization
bool gevd_edgel_regions::insert_edgel(float pre_x, float pre_y,
                                      float xedgel, float yedgel, gevd_region_edge* re)
{
  unsigned int xinterp, yinterp;
  bool edge_insert = false;
  while (line_gen(pre_x, pre_y, xedgel, yedgel, xinterp, yinterp))
    {
      if (out_of_bounds(X(xinterp), Y(yinterp)))
        {
          vcl_cout << "In gevd_edgel_regions::insert_edgel - out of bounds "
                   << "at (" << xinterp << ' ' << yinterp << ")\n";
          continue;
        }

      region_label_array_[Y(yinterp)][X(xinterp)] = EDGE;

      if (!re) continue;
      gevd_region_edge* old_re = edge_boundary_array_[Y(yinterp)][X(xinterp)];
      if (old_re&&old_re->get_edge())
        {
        //        old_re->UnProtect();
          edge_boundary_array_[Y(yinterp)][X(xinterp)] = re;
          //      re->Protect();
          edge_insert = true;
        }
      if (!old_re)
        {
          edge_boundary_array_[Y(yinterp)][X(xinterp)] = re;
          //      re->Protect();
          edge_insert = true;
        }
    }
  return edge_insert;
}

int gevd_edgel_regions::bytes_per_pix()
{
  int bypp = 1;
  if (image_source_)
    bypp = (image_->components() / image_->bits_per_component());

  if (buf_source_)
    bypp = buf_->GetBytesPixel();
  return bypp;
}

//-----------------------------------------------------------
//: Initialize the region label array.
//  There are three types of region label symbols:
//  1) UNLABELED - no label has yet been assigned,
//  2) EDGE, the existence of an edgel boundary pixel.
//  3) An unsigned integer which represents an existing region label.

bool gevd_edgel_regions::InitRegionArray(vcl_vector< vtol_edge_2d_sptr>& sg)
{
  if (!this->GroupContainsEdges(sg))
    return false;

  double xmin;
  double ymin;
  double xmax;
  double ymax;
  vsol_box_2d_sptr b;

  vcl_vector<vtol_edge_2d_sptr>::iterator i = sg.begin();
  assert( i != sg.end() );
  b=(*i)->get_bounding_box();
  xmin=b->get_min_x();
  ymin=b->get_min_y();
  xmax=b->get_max_x();
  ymax=b->get_max_y();
  for (++i; i!=sg.end(); ++i)
    {
      b=(*i)->get_bounding_box();
      if (b->get_min_x()<xmin)
        xmin=b->get_min_x();
      if (b->get_min_y()<ymin)
        ymin=b->get_min_y();
      if (b->get_max_x()>xmax)
        xmax=b->get_max_x();
      if (b->get_max_y()>ymax)
        ymax=b->get_max_y();
    }

  //Get the size of the arrays
  xo_ = (unsigned int)xmin;
  yo_ = (unsigned int)ymin;

  xend_ = (unsigned int)xmax;
  yend_ = (unsigned int)ymax;

  unsigned int sizex = xend_ - xo_ + 1;
  unsigned int sizey = yend_ - yo_ + 1;
  //Construct the arrays
  edge_boundary_array_ = new gevd_region_edge**[sizey];
  unsigned int y;
  for (y = yo_; y<=yend_; y++)
    edge_boundary_array_[Y(y)] = new gevd_region_edge*[sizex];

  region_label_array_ = new unsigned int*[sizey];
  for (y = yo_; y<=yend_; y++)
    region_label_array_[Y(y)] = new unsigned int[sizex];

  switch (this->bytes_per_pix())
  {
   case 1:     // Grey scale image with one byte per pixel
    ubuf_ = new unsigned char[sizex];
    break;
   case 2:     // Grey scale image with an unsigned short per pixel
     sbuf_ = new unsigned short[sizex];
     break;
   default:
    vcl_cout<<"In vtol_intensity_face::get_intensity(): bytes/pixel not 1 or 2\n";
  }

  //Initialize the arrays
  unsigned int x;
  if (ubuf_)
    for (x = xo_; x<=xend_; x++)
      ubuf_[X(x)] = 0;

  if (sbuf_)
    for (x = xo_; x<=xend_; x++)
      sbuf_[X(x)] = 0;

  for (y = yo_; y<=yend_; y++)
    for (x = xo_; x<=xend_; x++)
      {
        region_label_array_[Y(y)][X(x)] = UNLABELED;
        edge_boundary_array_[Y(y)][X(x)] = NULL;
      }
  //Insert edgels into arrays.

  int counter=0;
  for (vcl_vector<vtol_edge_2d_sptr >::iterator sgit = sg.begin();
       sgit != sg.end(); sgit++)
    {
      vtol_edge_2d_sptr e = (*sgit);
      if (!e)
        continue;
      e->set_id(counter++);
      vdgl_digital_curve* dc = e->curve()->cast_to_vdgl_digital_curve();
      if (!dc)
        continue;
      //The gevd_region_edge enables the link between a region label and
      //an Edge from the input segmentation
      gevd_region_edge* re = new gevd_region_edge(e);
      region_edges_[e->get_id()] = re;
      //      e->Protect(); re->Protect();
      //      float * ex = dc->GetX();
      //      float * ey = dc->GetY();
      vdgl_edgel_chain_sptr xy= dc->get_interpolator()->get_edgel_chain();
      int n_edgels = xy->size();

      //Insert the vertices at the ends of the edge
      //If there is a gap between the vertex locations and the
      //nearest edgels, the insert_edgel function will generate
      //edge insertions to prevent any gaps in the region boundary
      //There shouldn't be DigitalCurve(s) with this defect but it
      //does seem to occur.
      gevd_region_edge* vre = new gevd_region_edge(NULL);
      float pex1, pey1, pex2, pey2;
      if (n_edgels>0)
        {
          pex1 = float((*xy)[0].x()); pex2 = float((*xy)[n_edgels-1].x());
          pey1 = float((*xy)[0].y()); pey2 = float((*xy)[n_edgels-1].y());
          this->insert_edgel(pex1, pey1,
                             float(e->v1()->cast_to_vertex_2d()->x()),
                             float(e->v1()->cast_to_vertex_2d()->y()),
                             vre);
          this->insert_edgel(pex2, pey2,
                             float(e->v2()->cast_to_vertex_2d()->x()),
                             float(e->v2()->cast_to_vertex_2d()->y()),
                             vre);
        }
      else
        {
          pex1 = float(e->v1()->cast_to_vertex_2d()->x());
          pey1 = float(e->v1()->cast_to_vertex_2d()->y());
          pex2 = float(e->v2()->cast_to_vertex_2d()->x());
          pey2 = float(e->v2()->cast_to_vertex_2d()->y());
          this->insert_edgel(pex1, pey1, pex1, pey1, vre);
          this->insert_edgel(pex1, pey1, pex2, pey2, vre);
        }
      //Insert the edgels between the vertices
      //An edgel cannot be inserted on top of an existing vertex
      //insertion.  A vertex is distinguished by having a NULL Edge.
      //That is, for vertices it is ambiguous what Edge instance
      //exists at that location. Edge insertion fails if no edgels
      //of the Edge can be sucessfully inserted.
      bool edge_insert = false;
      for (int k =0; k < n_edgels; k++)
        {
          bool inserted = this->insert_edgel(pex1, pey1, float((*xy)[k].x()), float((*xy)[k].y()), re);
          edge_insert = edge_insert || inserted;
          pex1 = float((*xy)[k].x());
          pey1 = float((*xy)[k].y());
        }
      if (!edge_insert)
        {
          vcl_cout << "Edge Insert Failed for (" << e->v1() << ' '
                   << e->v2() << ")N: "<< n_edgels << vcl_endl;
          failed_insertions_->push_back(e);
        }
    }
  return true;
}

//---------------------------------------------------------------
//: Get code for a given label.
//
unsigned char gevd_edgel_regions::label_code(unsigned int label)
{
  unsigned char result;
  if (label<min_region_label_)
    result = label;
  else
    result = LABEL;
  return result;
}

//----------------------------------------------------------------
//: Add a new pair to the forward equivalence list.
//    That is, a ==> b. Note that there can be multiple equivalences which are
//    stored in a vcl_vector
bool gevd_edgel_regions::add_to_forward(unsigned int key, unsigned int value)
{
  bool result = true;
  vcl_map<unsigned int, vcl_vector<unsigned int>* >::iterator rpfi;
  rpfi = region_pairs_forward_.find(key);

  if (rpfi !=region_pairs_forward_.end())
    {
    vcl_vector<unsigned int> * vec = region_pairs_forward_[key];
    bool found = false;
    for (unsigned int i =0 ; i < vec->size() ; i++)
      {
      if ((*vec)[i] == value)
        found = true;
      }
    if (!found)
      {
      vec->push_back(value);
      }
    }
  else
    {
      vcl_vector<unsigned int>* larray = new vcl_vector<unsigned int>;
      larray->push_back(value);
       region_pairs_forward_[key]=larray;
    }
  return result;
}

//----------------------------------------------------------------
//: Add a new pair to the reverse equivalence list.
//    That is, b==>a. Note that there can be multiple equivalences which are
//    stored in a vcl_vector
//
bool gevd_edgel_regions::add_to_reverse(unsigned int key, unsigned int value)
{
  bool result = true;
  vcl_map<unsigned int, vcl_vector<unsigned int>* >::iterator rpfi;
  rpfi = region_pairs_reverse_.find(key);

  if (rpfi !=region_pairs_reverse_.end())
    {
    vcl_vector<unsigned int> * vec = region_pairs_reverse_[key];
    bool found = false;
    for (unsigned int i =0 ; i < vec->size() ; i++)
      {
      if ((*vec)[i] == value)
        found = true;
      }
    if (!found)
      {
      vec->push_back(value);
      }
    }
  else
    {
      vcl_vector<unsigned int>* larray = new vcl_vector<unsigned int>;
      larray->push_back(value);
       region_pairs_reverse_[key]=larray;
    }
  return result;
}

//-------------------------------------------------------------------
//: Encode a 2x2 neighbor hood with the state of the region array for a given location.
//    The Neighborhood    The states are:
//         ul ur          UNLABELED, EDGE, LABEL
//         ll lr              0       1      2
//    The encoding maps the 2x2 pattern to an unsigned char.  The layout
//    of the uchar is:  [ul|ur|ll|lr] with 2 bits for the state of each
//    position.
unsigned char
gevd_edgel_regions::EncodeNeighborhood(unsigned int ul, unsigned int ur,
                                       unsigned int ll, unsigned int lr)
{
  unsigned char nhood = 0;

  unsigned char nul = label_code(ul);
  nul = nul<<6;
  nhood |= nul;

  unsigned char nur = label_code(ur);
  nur = nur<<4;
  nhood |= nur;

  unsigned char nll = label_code(ll);
  nll = nll << 2;
  nhood |= nll;

  unsigned char nlr = label_code(lr);
  nhood |= nlr;

  return nhood;
}

//----------------------------------------------------------------------
// --This is the fundamental assignment of label equivalence
//
void gevd_edgel_regions::insert_equivalence(unsigned int ll, unsigned int ur, unsigned int& lr)
{
  if (ll!=ur)//Is a=b?
    {//No. We want the smallest label to be the equivalence base
      //Also we don't want to loose earlier equivalences
      unsigned int smaller_label, larger_label;
      //Get the ordering right
      if (ll>ur)
        {smaller_label = ur; larger_label = ll;}
      else
        {smaller_label = ll; larger_label = ur;}
      this->add_to_forward(larger_label, smaller_label);
      this->add_to_reverse(smaller_label, larger_label);
      lr = smaller_label;
    }
  else //yes, a=b the simple case
    lr = ur;
}

//-------------------------------------------------------------------
//: Propagate connected components.
//    Uses an unsigned char encoding
//    a 2x2 neighborhood to propagate region labels. For example:
//    aa ->   aa
//    xx ->   aa
//    Here the lower two labels are set to the upper label.
//    This method operates directly on the region_label_array_.
//
//    Note that the 2x2 neighborhood is encoded as uchar = [ul|ur|ll|lr]
//                                                          2  2  2  2  bits
void gevd_edgel_regions::UpdateConnectedNeighborhood(unsigned int x, unsigned int y)
{
  unsigned int ul = region_label_array_[y][x];
  unsigned int ur = region_label_array_[y][x+1];
  unsigned int& ll = region_label_array_[y+1][x];
  unsigned int& lr = region_label_array_[y+1][x+1];

  unsigned char nhood = this->EncodeNeighborhood(ul, ur, ll, lr);
  switch(int(nhood))
  {
    case 0:
      //xx xx
      //xx xx
      return;
    case 1:
      //xx xx
      //xe xe
      return;
    case 5:
      //xx xx
      //ee ee
      return;

    case 17:
      //xe xe
      //xe xe
      return;

    case 20:
      //xe xe
      //ex ea
      lr = max_region_label_++;
      return;

    case 21:
      //xe xe
      //ee ee
      return;

    case 22:
      //xe xe
      //ea ea
      return;

    case 68:
      //ex ex
      //ex ex
      return;

    case 72:
      //ex ex
      //ax ax
      return;

    case 133:
      //ax aa
      //ee ee
      return;

    case 136:
      //ax ax
      //ax ax
      return;

    case 160:
      //aa aa
      //xx aa
      ll = ul;
      lr = ul;
      return;

    case 161:
      //aa aa
      //xe ae
      ll = ul;
      return;

    case 168:
      //aa aa
      //ax aa
      insert_equivalence(ll, ur, lr);
      return;

    case 169:
      //aa aa
      //ae ae
      return;

    case 164:
      //aa aa
      //ex ea
      lr = ul;
      return;

    case 165:
      //aa aa
      //ee ee
      return;

    case 144:
      //ae ae
      //xx aa
      ll = ul;
      lr = ul;
      return;

    case 145:
      //ae ae
      //xe ae
      ll = ul;
      return;

    case 152:
      //ae ae
      //ax aa
      lr = ul;
      return;

    case 153:
      //ae ae
      //ae ae
      return;

    case 148:
      //ae ae
      //ex eb
      lr = max_region_label_++;
      return;

    case 149:
      //ae ae
      //ee ee
      return;

    case 96:
      //ea ea
      //xx aa
      ll = ur;
      lr = ur;
      return;

    case 97:
      //ea ea
      //xe be
      ll = max_region_label_++;
      return;

    case 104:
      //ea ea
      //bx aa
      insert_equivalence(ll, ur, lr);
      return;

    case 105:
      //ea ea
      //ae ae
      return;

    case 100:
      //ea ea
      //ex ea
      lr = ur;
      return;

    case 101:
      //ea ea
      //ee ee
      return<