sdet_edgel_regions.cxx

Go to the documentation of this file.

// This is brl/bseg/sdet/sdet_edgel_regions.cxx
#include "sdet_edgel_regions.h"
//:
// \file

#include <vcl_iostream.h>
#include <vcl_algorithm.h> // vcl_find()
#include <vbl/vbl_qsort.h> //for sorting labels

#include <vxl_config.h>
#include <vil1/vil1_memory_image_of.h>
#include <vul/vul_timer.h>

#include <vsol/vsol_box_2d.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_edgel_chain.h>
#include <vdgl/vdgl_interpolator.h>
#include <vdgl/vdgl_digital_curve.h>
#include <vtol/vtol_intensity_face.h>

#include <bdgl/bdgl_curve_algs.h>
#include <btol/btol_edge_algs.h>
#include <sdet/sdet_region_edge.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

//: Print the region label array
void sdet_edgel_regions::print_region_array()
{
  vcl_cout << vcl_endl << vcl_endl;
  for (unsigned int y = 0; y<ys_; y++)
  {
    for (unsigned int x = 0; x<xs_; x++)
      if (region_label_array_[y][x]==EDGE
          //&&edge_boundary_array_[Y(y)][X(x)]->IsVertex()
         )
        vcl_cout << '*';
      else
        vcl_cout << region_label_array_[y][x] << ' ';
    vcl_cout << vcl_endl;
  }
  vcl_cout << vcl_endl << vcl_endl;
}

//: Print the contents of the forward equivalence index
void sdet_edgel_regions::print_region_equivalence()
{
  vcl_cout << "\nLabel 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 << " == ";
    vcl_vector<unsigned int>* labels = (*rpf_iterator).second;
    if (labels)
    {
      for (vcl_vector<unsigned int>::iterator lit = labels->begin();
           lit != labels->end(); lit++)
        vcl_cout << *lit << ' ';
      vcl_cout << vcl_endl;
    }
  }
  vcl_cout << vcl_endl;
}

//: Print the contents of the reverse equivalence index
void sdet_edgel_regions::print_reverse_region_equivalence()
{
  vcl_cout << "\nReverse 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 << " == ";
    vcl_vector<unsigned int>* labels = (*rpf_iterator).second;
    if (labels)
    {
      for (vcl_vector<unsigned int>::iterator lit = labels->begin();
           lit != labels->end(); lit++)
        vcl_cout << *lit << ' ';
      vcl_cout << vcl_endl;
    }
  }
  vcl_cout << vcl_endl;
}

//: Print the reduced equivalence relation
void sdet_edgel_regions::print_base_equivalence()
{
  vcl_cout << "\nBase 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 sdet_edgel_regions::print_intensity_data()
{
  for (vcl_vector<vtol_intensity_face_sptr>::iterator fit =faces_->begin();
       fit != faces_->end(); fit++)
  {
    vtol_intensity_face_sptr const & f = *fit;
    vcl_cout << "IntFaceAt(" << f->Xo() << ' ' << f->Yo() << "):\n";
    //      f->PrintFit();
  }
}

//--------------------------------------------------------------
//: The region array can be offset from the corner of the ROI.
//  This method transforms from the ROI coordinate
//  system to the region label array coordinate system.
unsigned int sdet_edgel_regions::X(unsigned int x)
{
  return x-xo_;
}

//: The region array can be offset from the corner of the ROI.
//  This method transforms from the ROI coordinate
//  system to the region label array coordinate system.
unsigned int sdet_edgel_regions::Y(unsigned int y)
{
  return y-yo_;
}

//---------------------------------------------------------
//: Transforms the image x coordinate to the array coordinate with a scale factor
float sdet_edgel_regions::Xf(float x)
{
  float xu = (x-xo_)*s_;
  return xu;
}

//---------------------------------------------------------
//:
float sdet_edgel_regions::Yf(float y)
{
  float yu = (y-yo_)*s_;
  return yu;
}

//----------------------------------------------------------
// Default constructor
sdet_edgel_regions::sdet_edgel_regions(int array_scale,bool verbose,
                                       bool debug)
{
  debug_ = debug;
  verbose_ = verbose;
  s_ = array_scale; //Some boundaries can approach less than 1/2 pixel proximity
#ifdef DEBUG
  debug_data_ = debug_ ? new topo_debug_data : 0;
#endif
  image_source_ = false;
  buf_source_ = false;
  buf_ = NULL;
  xo_=0;
  yo_=0;
  xend_ = 0;
  yend_ = 0;
  xs_ = 0;
  ys_ = 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
sdet_edgel_regions::~sdet_edgel_regions()
{
  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_;
  //fix leaks
  delete faces_;
  delete [] intensity_face_index_;

  for (vcl_map<unsigned int, vcl_vector<vtol_edge_2d_sptr>* >::iterator
       mit = region_edge_adjacency_.begin();
       mit != region_edge_adjacency_.end(); mit++)
  {
    if ((*mit).second)
      (*mit).second->clear();
    delete (*mit).second;
  }

  for (vcl_map<unsigned int, vcl_vector<unsigned int>* >::iterator
       mit = equivalence_set_.begin(); mit != equivalence_set_.end(); mit++)
    delete (*mit).second;

  for (vcl_map<unsigned int, vcl_vector<unsigned int>* >::iterator
       mit = region_pairs_reverse_.begin();
       mit != region_pairs_reverse_.end(); mit++)
    delete (*mit).second;

  for (vcl_map<unsigned int, vcl_vector<unsigned int>* >::iterator
       mit = region_pairs_forward_.begin();
       mit != region_pairs_forward_.end(); mit++)
    delete (*mit).second;
}

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

//-----------------------------------------------------------------
//: The key process loop.
//  Carries out the steps:
//  - Connected components
//  - Edge-label assignment
//  - Collect region boundaries
//  - Construct vtol_intensity_faces
//  - Calculate intensity fit
bool sdet_edgel_regions::
compute_edgel_regions(vil1_image const& image,
                      vcl_vector<vtol_edge_2d_sptr>& sgrp,
                      vcl_vector<vtol_intensity_face_sptr>& faces)
{
#if 0
  image_ = image;
  buf_ = NULL;
  image_source_=true;
  buf_source_ = false;
#endif
  //change to uniform use of bufferxy for ease of interfacing to different
  //image types
  image_ = 0;
  image_source_=false;
  buf_source_ = false;
  buf_ = new gevd_bufferxy(image);
  if(buf_)
    buf_source_ = true;
  else
    return false;
  return compute_edgel_regions(sgrp, faces);
}
bool sdet_edgel_regions::
compute_edgel_regions(vil_image_resource_sptr const& image,
                      vcl_vector<vtol_edge_2d_sptr>& sgrp,
                      vcl_vector<vtol_intensity_face_sptr>& faces)
{
#if 0
  image_ = image;
  buf_ = NULL;
  image_source_=true;
  buf_source_ = false;
#endif
  //change to uniform use of bufferxy for ease of interfacing to different
  //image types
  image_ = 0;
  image_source_=false;
  buf_source_ = false;
  buf_ = new gevd_bufferxy(image);
  if(buf_)
    buf_source_ = true;
  else
    return false;
  return compute_edgel_regions(sgrp, faces);
}

bool sdet_edgel_regions::
compute_edgel_regions(vcl_vector<vtol_edge_2d_sptr>& sgrp,
                      vcl_vector<vtol_intensity_face_sptr>& faces)
{
#ifdef DEBUG
  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 (debug_)
    this->print_region_array();

  unsigned int y, x;
  //Propagate connected components
  // the -1 accounts for the 2x2 label classification neighborhood
  for (y=0; y<ys_-1; y++)
    for (x=0; x<xs_-1; x++)
      this->UpdateConnectedNeighborhood(x,y);

  if (debug_)
    this->print_region_array();
  //Resolve region label equivalence
  if (debug_)
  {
    this->print_region_equivalence();
    this->print_reverse_region_equivalence();
  }
  this->GrowEquivalenceClasses();
  this->PropagateEquivalence();
  if (debug_)
    this->print_base_equivalence();
  this->ApplyRegionEquivalence();
  if (debug_)
    vcl_cout << "After Label Equivalence\n";
  if (debug_)
    this->print_region_array();

  //Associate Edge(s) with region labels (-1?)JLM
  for (y=0; y<ys_-1; y++)
    for (x=0; x<xs_-1; x++)
      this->AssignEdgeLabels(x, y);

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

  //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 (debug_)
    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);

  if (verbose_)
  vcl_cout << "Compute Edgel Regions Total(" << faces_->size() << ") in "
           << t.real() << " msecs.\n"<< vcl_flush;
  return true;
}


//----------------------------------------------------------
//: assign equivalence of region label b to region label a
//
bool sdet_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 sdet_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;
}

//------------------------------------------------------------
//: Return label image (255/0) indicating edgels.
//  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 sdet_edgel_regions::GetEdgeImage(vcl_vector<vtol_edge_2d_sptr>& sg)
{
  if (!this->InitRegionArray(sg)) return NULL;
  unsigned char no_edge = 0, edge = 255;

  vil1_memory_image_of<vxl_byte> image(xs_,ys_);

  for (unsigned int y = 0; y<ys_; y++)
    for (unsigned int x = 0; x<xs_; 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 sdet_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 equivalence_set_.
//  cur_label is the equivalence set being formed. label is the table key of equivalences
//  to be merged. The labels equivalent to label are in the input map, tab.
bool sdet_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;
  //If we can't find the label in tab then there are no equivalences to be merged
  hashi = tab.find(label);
  if (hashi == tab.end()) {
    return false;
  }
  //We did find label, and labels is a set of equivalent labels
  vcl_vector<unsigned int>* labels = hashi->second;

  //If the set is empty then nothing to do
  int len = labels->size();
  if (!len) {
    return false;
  }

  //The set of labels equivalent to cur_label
  vcl_vector<unsigned int>* array;

  hashi = equivalence_set_.find(cur_label);
  if ( hashi == equivalence_set_.end())
  {
    //We didn't find any equivalent labels for cur_label so we initialize a new one
    //and insert it into equivalence_set
    array = new vcl_vector<unsigned int>;
    equivalence_set_[cur_label] = array;
  }
  else
  {
    //otherwise there is already a growing set of equivalent labels
    array = hashi->second;
  }

  //We look through the set of labels found in map for label
  for (int i = 0; i<len; i++)
  {
    unsigned int l = (*labels)[i];
    bool found = false;
    //see if l is in the current set of equivalent labels for cur_label
    for (unsigned int j=0 ; j < array->size() ; j++)
    {
      if ((*array)[j] == l) found = true;
    }
    //If label l is not already there then add it
    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 sdet_edgel_regions::get_next_label(vcl_vector<unsigned int>* labels,
                                        unsigned int& label)
{
  //If the set labels is null then
  //just return the next larger label (if less than max_region_label_)
  unsigned int tmp = label+1;
  if (!labels)
    if (tmp<max_region_label_)
    {
      label = tmp;
      return true;
    }

  //The set is not empty, so search for a label not found in the
  //set.

  for (unsigned int i = tmp; i<max_region_label_; i++)
  {
    //if we don't find i we can use it as the new label
    if (vcl_find(labels->begin(), labels->end(), i) == labels->end())
    {
      label = i;
      return true;
    }
  }
  return false;
}

//----------------------------------------------------------------
//: Form equivalence classes by transitive closure on each label.
//  The resulting label equivalence is stored in the map, equivalence_set_.
//
//  The idea is to add labels to the current equivalence set, cur_set,
//  from either forward or reverse equivalence classes until no
//  new equivalences 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 entries are removed from
//  forward and reverse forward and reverse equivalence maps as they
//  are used.  When the cur_set is completely closed, i.e. no new labels
//  can be added, then a new label not in cur_set is used to
//  seed a new equivalence class.
//
void sdet_edgel_regions::GrowEquivalenceClasses()
{
  if ((max_region_label_-min_region_label_) < 2)
    return;
  //start with the minimum label
  unsigned int cur_label = min_region_label_;
  //the iterator for equivalence_set_
  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;
    int old_len;
    while (merging)
    {
      old_len = len;
      merging = false;

      //find label equivalence in the forward map
      bool find_forward =
        this->merge_equivalence(region_pairs_forward_, cur_label, i);
      if (find_forward)
      {
        delete region_pairs_forward_[i];//memory leak otherwise
        region_pairs_forward_.erase(i);//remove hits from the forward map
      }
      //find label equivalence in the reverse map
      bool find_reverse =
        this->merge_equivalence(region_pairs_reverse_, cur_label, i);
      if (find_reverse)
      {
        delete region_pairs_reverse_[i];//memory leak otherwise
        region_pairs_reverse_.erase(i);//remove hits from the reverse map
      }
      //At this point we may have established or added to the equivalence set
      //for cur_label stored in equivalence_set[cur_label]
      //we check if we have
      mei = equivalence_set_.find(cur_label);
      if (mei == equivalence_set_.end()) continue;
      //If we don't find cur_label, it means that we couldn't find label i
      //in either region_pairs_forward or region_pairs_reverse, so we need
      //to get a new i.

      //We did find a growing equivalence set, cur_set, but it might be null
      //or empty
      cur_set = equivalence_set_[cur_label];
      if (!cur_set) continue;
      len = cur_set->size();
      if (!len) continue;

      //Now we check cur_set has actually been extended.
      //If it has we have to sort the labels so that we can find
      //the next larger label in cur_set easily
      if (len > old_len)
      {
        i = cur_label;
        vbl_qsort_ascending(*cur_set);//vcl_sort didn't seem to work!
      }

      //Get the next larger label from cur_set
      //so that we can insert its equivalent labels
      for (int j = 0; j<len&&!merging; j++)
        if ((*cur_set)[j]>i)
        {
          i = (*cur_set)[j];
          merging = true;
        }
      //If we reach here with merging = false then we have finished the
      //current equivalence class
    }
    //Find the next largest label that isn't in cur_set to seed the
    //next equivalence class
    if (!get_next_label(cur_set, cur_label)) return;
  }
}

//------------------------------------------------------------
//: Check if the vtol_edge list sg (spatial group) contains edge(s) - NYI
bool sdet_edgel_regions::GroupContainsEdges(vcl_vector<vtol_edge_2d_sptr>& /* sg */)
{
#if 0
  CoolString type(sg.GetSpatialGroupName());
  return type == CoolString("EdgelGroup") ||
    type == CoolString("FittedEdgeGroup");
#endif
  return true; // TODO
}


//----------------------------------------------------------
//: 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 sdet_edgel_regions::insert_edgel(float pre_x, float pre_y,
                                      float xedgel, float yedgel,
                                      sdet_region_edge_sptr const& re)
{
  //Transform image pixel coordinates to array coordinates
  float pre_xf = Xf(pre_x), pre_yf = Yf(pre_y),
    xedgelf = Xf(xedgel), yedgelf = Yf(yedgel);

  float xf, yf;
  unsigned int xinterp, yinterp;

  bool edge_insert = false;
  bool init = true, done = false;
  while (bdgl_curve_algs::line_gen(pre_xf, pre_yf, xedgelf, yedgelf,
                                   init, done, xf, yf))
  {
    xinterp = (unsigned int)xf;
    yinterp = (unsigned int)yf;
    if (out_of_bounds(xinterp, yinterp))
    {
      vcl_cout << "In sdet_edgel_regions::insert_edgel - out of bounds "
               << "at array(" << xinterp << ' ' << yinterp << ")\n";
      continue;
    }

    region_label_array_[yinterp][xinterp] = EDGE;

    if (!re) continue;
    sdet_region_edge_sptr const& old_re =
      edge_boundary_array_[yinterp][xinterp];
    if (old_re&&old_re->get_edge())
    {
      edge_boundary_array_[yinterp][xinterp] = re;
      edge_insert = true;
    }
    if (!old_re)
    {
      edge_boundary_array_[yinterp][xinterp] = re;
      edge_insert = true;
    }
  }
  return edge_insert;
}

int sdet_edgel_regions::bytes_per_pix()
{
  int bypp = 1;
  if (image_source_)
  {
    int bytes_per_comp = (image_.bits_per_component()+7)/8;
    bypp = image_.components() * bytes_per_comp;
  }
  if (buf_source_)
    bypp = buf_->GetBytesPixel();
  return bypp;
}

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

bool sdet_edgel_regions::InitRegionArray(vcl_vector< vtol_edge_2d_sptr>& sg)
{
  if (!this->GroupContainsEdges(sg))
    return false;
  vsol_box_2d b = btol_edge_algs::bounding_box(sg);
  //Get the bounds for the arrays from the input edges
  xo_ = (unsigned int)b.get_min_x();
  yo_ = (unsigned int)b.get_min_y();

  xend_ = (unsigned int)b.get_max_x();
  yend_ = (unsigned int)b.get_max_y();

  xs_ = (this->GetXSize()-1)*s_+1;
  ys_ = (this->GetYSize()-1)*s_+1;

  //JLM debug
  vcl_cout << "bounding box(" << xo_ << ' ' << yo_ << ' ' << xs_ << ' ' << ys_
           << ")\n";

  //Construct the edge and label arrays
  edge_boundary_array_ = vbl_array_2d<sdet_region_edge_sptr>(ys_, xs_);

  region_label_array_ = vbl_array_2d<unsigned int>(ys_, xs_);

  //initialize
  region_label_array_.fill(UNLABELED);
  edge_boundary_array_.fill(NULL);

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

  //Initialize the buffers
  if (ubuf_)
    for (int x = 0; x<this->GetXSize(); ++x)
      ubuf_[x] = 0;

  if (sbuf_)
    for (int x = 0; x<this->GetXSize(); ++x)
      sbuf_[x] = 0;

  //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 sdet_region_edge enables the link between a region label and
    //an Edge from the input segmentation
    sdet_region_edge_sptr const & re = new sdet_region_edge(e);
    region_edges_[e->get_id()] = re;

    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.
    sdet_region_edge_sptr const& vre = new sdet_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 sdet_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 sdet_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 sdet_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 neighborhood with the state of the region array for a given location.
//  \verbatim
//    The Neighborhood    The states are:
//         ul ur          UNLABELED, EDGE, LABEL
//         ll lr              0       1      2
//  \endverbatim
//  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
sdet_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 sdet_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:
//  \verbatim
//  aa ->   aa
//  xx ->   aa
//  \endverbatim
//  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 sdet_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;
   case 80:
    //ee ee
    //xx bb
    ll = max_region_label_++;
    lr = ll;
    return;
   case 81:
    //ee ee
    //xe be
    ll = max_region_label_++;
    return;
   case 88:
    //ee ee
    //ax aa
    lr = ll;
    return;
   case 89:
    //ee ee
    //ae ae
    return;
   case 84:
    //ee ee
    //ex eb
    lr = max_region_label_++;
    return;
   case 85:
    //ee ee
    //ee ee
    return;
   default:
    vcl_cout << "In sdet_edgel_regions::UpdateNeigborhood(..)"
             << "impossible pattern(" << x << ' ' << y << ") = " << (int)nhood << '\n'
             << int(label_code(ul)) << ' ' << int(label_code(ur)) << '\n'
             << int(label_code(ll)) << ' ' << int(label_code(lr)) << "\n\n";
  }
}

static bool reg_edges_neq(sdet_region_edge_sptr