gevd_contour.cxx

Go to the documentation of this file.

// This is gel/gevd/gevd_contour.cxx
#include "gevd_contour.h"
//:
// \file
#include <vcl_iostream.h>
#include <vcl_cstdlib.h>   // for vcl_abs(int)
#include <vcl_vector.h>
#include <vcl_algorithm.h> // for vcl_max()
#include <vxl_config.h>
#include <vnl/vnl_math.h> // for sqrt(2)
#include <vdgl/vdgl_digital_curve.h>
#include <vdgl/vdgl_edgel_chain.h>
#include <vdgl/vdgl_interpolator.h>
#include <vdgl/vdgl_interpolator_linear.h>
#include <vtol/vtol_vertex_2d.h>
#include <vtol/vtol_edge_2d.h>
#include <gevd/gevd_pixel.h>
#ifdef DEBUG
# include <vul/vul_timer.h>
#endif

const int INVALID = -1;

// Use 8 directions, with 45 degree angle in between them.

const vxl_byte TWOPI = 8, /* FULLPI = 4, */ HALFPI = 2 /* , QUARTERPI = 1 */;
//const vxl_byte DIR0 = 8, DIR1 = 9, DIR2 = 10, DIR3 = 11;
const int DIS[] = { 1, 1, 0,-1,-1,-1, 0, 1, // 8-connected neighbors
                    1, 1, 0,-1,-1,-1, 0, 1, // wrapped by 2PI to
                    1, 1, 0,-1,-1,-1, 0, 1};// avoid modulo operations.
const int DJS[] = { 0, 1, 1, 1, 0,-1,-1,-1,
                    0, 1, 1, 1, 0,-1,-1,-1,
                    0, 1, 1, 1, 0,-1,-1,-1};

//const int RDS[] = {0,-1, 1,-2, 2,-3, 3,-4, 4,-5, 5}; // radial search
const int RIS[] = { 1, 0,-1, 0, // spiral search for 4/8-connected
                    1,-1,-1, 1, // neighbors
                    2, 0,-2, 0,
                    2, 1,-1,-2,-2,-1, 1, 2,
                    2,-2,-2, 2,
                    3, 0,-3, 0,
                    3, 1,-1,-3,-3,-1, 1, 3,
                    3, 2,-2,-3,-3,-2, 2, 3,
                    4, 0,-4, 0};
const int RJS[] = { 0, 1, 0,-1, // rotate CW, increasing radius
                    1, 1,-1,-1,
                    0, 2, 0,-2,
                    1, 2, 2, 1,-1,-2,-2,-1,
                    2, 2,-2,-2,
                    0, 3, 0,-3,
                    1, 3, 3, 1,-1,-3,-3,-1,
                    2, 3, 3, 2,-2,-3,-3,-2,
                    0, 4, 0,-4};
const int RNS[] = { 4, 8, 12, 20, 24, 28, 36, 44, 48}; // at distinct r
const float RGS[] = { 1.f, 1.414213f, 2.f, 2.236067f, 2.828427f, // values of gap
                      3.f, 3.162277f, 3.605551f, 4.f};

// - win32 - moved to here for MSVC++
const int MINLENGTH = 3;        // minimum number of pixels for a chain
const int FRAME = 4;            // border of image

bool gevd_contour::talkative = true;    // By default contour is not silent.


//: Save parameters and create workspace for detecting contours.
// Each contour must have at least 1 pixel above min_strength,
// and its number of internal pixels must be above min_length.
// This is a heuristic hysteresis scheme that prunes weak or short
// isolated chains.
// To join a weaker contour to a stronger contour, a junction must
// have a change in response above min_jump on the stronger contour.
// This way, only strong junctions are detected.
gevd_contour::gevd_contour(float min_strength, int min_length,
                           float min_jump, float max_gap)
  : minStrength(min_strength), minLength(min_length),
    minJump(min_jump), maxSpiral(0),
    edgeMap(), vertexMap()
{
  if (minStrength < 0) {
    vcl_cerr << "gevd_contour::gevd_contour -- negative min_strength: "
             << minStrength << ". Reset to 0.\n";
    minStrength = 0;
  }
  if (minLength < MINLENGTH) {
    vcl_cerr << "gevd_contour::gevd_contour -- too small min_length: "
             << minLength << ". Reset to " << MINLENGTH << ".\n";
    minLength = MINLENGTH;
  }
  if (minJump < 0) {
    vcl_cerr << "gevd_contour::gevd_contour -- negative min_jump: "
             << minJump << ". Reset to 0.\n";
    minJump = 0;
  }
  if (minJump > minStrength) {
    vcl_cerr << "gevd_contour::gevd_contour -- too large min_jump: "
             << min_jump << ". Reset to " << minStrength << ".\n";
    minJump = minStrength;
  }
  if (max_gap < 1) {
    vcl_cerr << "gevd_contour::gevd_contour -- too small max_gap: "
             << max_gap << ". Reset to 1.\n";
    max_gap = 1;
  }
  if (max_gap > FRAME) {
    vcl_cerr << "gevd_contour::gevd_contour -- too large max_gap: "
             << max_gap << ". Reset to " << FRAME << vcl_endl;
    max_gap = FRAME;
  }
  for (int i = 0; i < 9; i++)   // find number of neighbors to search
    if (max_gap <= RGS[i])      // for given gap radius
      maxSpiral= i+1;
}


//: Free space allocated for detecting contours.
gevd_contour::~gevd_contour()
{
  delete edgeMap;               // space shared by LinkJunction/Chain
  delete vertexMap;
}


//: Find network of linked edges and vertices, from 8-connected edge elements.
// The contours must be less than 2 pixel wide,
// for example found from non maximum suppression.
// Isolated edgels and short segments are erased.
bool
gevd_contour::FindNetwork(gevd_bufferxy& edgels,
                          const int njunction,
                          const int* junctionx, const int* junctiony,
                          vcl_vector<vtol_edge_2d_sptr>*& edges,
                          vcl_vector<vtol_vertex_2d_sptr >*& vertices)
{
  // make sure that if no edges are found that edges and vertices
  // get values, to avoid seg faults, WAH
  if (!edges)
    edges = new vcl_vector<vtol_edge_2d_sptr>;
  else
    edges->clear();
  if (!vertices)
    vertices = new vcl_vector<vtol_vertex_2d_sptr >;
  else
    vertices->clear();


  if (talkative)
    vcl_cout << "*** Link edge elements into connected edges/vertices.\n";

  // 1. Setup lookup maps based on (x,y) integer location.
  vertexMap = new vbl_array_2d<vtol_vertex_2d_sptr>(edgels.GetSizeX(), edgels.GetSizeY());
  vertexMap->fill(NULL);
  edgeMap = new vbl_array_2d<vtol_edge_2d_sptr>(edgels.GetSizeX(), edgels.GetSizeY());
  edgeMap->fill(NULL);

  // 2. Collect 4/8-connected pixels into chains
  int n; // = vcl_max(10*njunction, // preallocated size from junctions or
  //       edgels.GetSizeX()*edgels.GetSizeY()/100); // image size
  vcl_vector<vtol_edge_2d_sptr> *edges2 = new vcl_vector<vtol_edge_2d_sptr>;
  n = this->FindChains(edgels, // link pixels into chains
                       njunction, // also use junction pixels
                       junctionx, junctiony,
                       *edges2);
  if (!n)
    return false;               // empty network

  // 3. Sort chains longest first.

  if (edges2->size() < 1000)     // avoid O(nlogn) for very large n
    vcl_sort (edges2->begin(), edges2->end());
  //    edges2.sort(&gevd_contour::LengthCmp); // sort longest/strongest first

  // renumber with order in array
  for (unsigned int i= 0; i< edges2->size(); i++)
    (*edges2)[i]->set_id(i);

  // 4. Split/Merge chains from touching end points
  vcl_vector<vtol_vertex_2d_sptr > vertices2;
  this->FindJunctions(edgels, // break/merge at junctions of
                      *edges2, vertices2); // distinct chains

  // 5. Copy back results into global lists
  for (unsigned int i= 0; i< edges2->size(); i++)
    edges->push_back( (*edges2)[i]);

  for (unsigned int i=0; i< vertices2.size(); i++)
    vertices->push_back( vertices2[i]);

  return true;
}


//: Return TRUE if pixel is a local maximum, and so is right on top of contour.
static bool
on_contour(const gevd_bufferxy& edgels, const int i, const int j)
{
  double pix = (1 + vnl_math::sqrt2) * floatPixel(edgels, i, j); // fuzzy threshold
  for (vxl_byte dir = 0; dir < TWOPI; dir += HALFPI) // 4-connected only
    if (floatPixel(edgels, i+DIS[dir], j+DJS[dir]) > pix)
      return false;             // should choose neighbor instead
  return true;
}


//: Delete pixel from contour, and save its location in xloc/yloc.
static void
RecordPixel(int i, int j, gevd_bufferxy& edgels,
            vcl_vector<int>& iloc, vcl_vector<int>& jloc)
{
  floatPixel(edgels, i, j) = -floatPixel(edgels, i, j); // flip sign
  iloc.push_back(i), jloc.push_back(j);
}


//:
// Find next best pixel on contour, searching for strongest response,
// and favoring 4-connected over 8-connected.
// Return 0, if no pixel is found, or direction in range [2*pi, 4*pi).
static int
NextPixel(int& i, int& j, const gevd_bufferxy& edgels)
{
  float maxpix = 0, npix;
  int maxdir = 0, dir;
  for (dir = 0; dir < TWOPI; dir += HALFPI) // 4-connected first
    if ((npix = floatPixel(edgels, i+DIS[dir], j+DJS[dir])) > maxpix) {
      maxpix = npix;
      maxdir = dir+TWOPI;
    }
  if (!maxdir) {
    for (dir = 1; dir < TWOPI; dir += HALFPI) // 8-connected next
      if ((npix = floatPixel(edgels, i+DIS[dir], j+DJS[dir])) > maxpix) {
        maxpix = npix;
        maxdir = dir+TWOPI;
      }
  }
  if (maxdir)                   // update next strongest pixel
    i += DIS[maxdir], j += DJS[maxdir];
  return maxdir;
}


//:
// Find next best pixel on contour, searching for strongest response,
// and favoring 4-connected over 8-connected.
// Return 0, if no pixel is found, or direction in range [2*pi, 4*pi).
static int
next_pixel(int& i, int& j, const vbl_array_2d<vtol_vertex_2d_sptr>& vertexMap)
{
  int maxdir = 0, dir;
  for (dir = 0; dir < TWOPI; dir += HALFPI) // 4-connected first
    if (vertexMap.get(i+DIS[dir], j+DJS[dir])) {
      maxdir = dir+TWOPI;
      break;
    }
  if (!maxdir) {
    for (dir = 1; dir < TWOPI; dir += HALFPI) // 8-connected next
      if (vertexMap.get(i+DIS[dir], j+DJS[dir])) {
        maxdir = dir+TWOPI;
        break;
      }
  }
  if (maxdir)                   // update next strongest pixel
    i += DIS[maxdir], j += DJS[maxdir];
  return maxdir;
}


//:
// Trace and collect pixels on thin contours, stronger pixels first,
// and favoring 4-connected over 8-connected. Thinning is not used,
// and so will avoid errors because of square grid tessellation.
// A chain can not cross itself. It can only touch itself or another
// chain, in which case a junction will be found later.
// The pixels of a chain include the 2 end points.
// End points and junctions are created in gevd_contour::FindJunctions.
// Return the number of chains found.  Protected.
int
gevd_contour::FindChains(gevd_bufferxy& edgels, const int njunction,
                         const int* junctionx, const int* junctiony,
                         vcl_vector<vtol_edge_2d_sptr>& edges)
{
#ifdef DEBUG
  vul_timer t;
#endif

  // 1. Save away detected junctions from extending at end points of
  // contours, without linking these contours up. This avoids random
  // order in the traversal of the contours.
  vtol_vertex_2d_sptr mark = new vtol_vertex_2d;         // dummy non zero pointer
  for (int k = 0; k < njunction; k++) {
    vertexMap->put(junctionx[k], junctiony[k], mark);
  }

  // 2. Trace elongated & thinned chains, stronger pixels first.
  // Virtual border of image should be inserted last.
  const int rmax = FRAME;
  const int xmax = edgels.GetSizeX()-rmax-1;
  const int ymax = edgels.GetSizeY()-rmax-1;
  vcl_vector<int> xloc(xmax+ymax), yloc(xmax+ymax); // work space for

  for (int j = rmax; j <= ymax; j++)
    for (int i = rmax; i <= xmax; i++)
    {
      // 2.0. Start from better pixels above noise+hysteresis
      if (floatPixel(edgels, i, j) > minStrength &&
          on_contour(edgels, i, j)) { // right on the contour
        int x = i, y = j;

        // 2.1. Prune isolated pixels
        if (!NextPixel(x, y, edgels)) {// prune isolated pixels
          floatPixel(edgels, i, j) = 0;
          continue;
        }

        // 2.2. Start collecting first 3 pixels
        xloc.clear(), yloc.clear(); // collect pixels on contour
        RecordPixel(i, j, edgels, xloc, yloc);  // first pixel
        int ii = x, jj = y;
        RecordPixel(ii, jj, edgels, xloc, yloc); // second pixel
        if (NextPixel(x, y, edgels))
          RecordPixel(x, y, edgels, xloc, yloc); // third pixel
        else {                  // reach end point
          x = i, y = j;         // revert back to start pt
          if (NextPixel(x, y, edgels)) { // reverse collection
            xloc.clear(), yloc.clear();
            RecordPixel(ii, jj, edgels, xloc, yloc); // second pixel
            RecordPixel(i, j, edgels, xloc, yloc); // first pixel
            RecordPixel(x, y, edgels, xloc, yloc); // third pixel
            ii = i, jj = j;
          } else  {             // reach other end point
            floatPixel(edgels, i, j) = 0; // prune isolated pixel-pairs
            floatPixel(edgels, ii, jj) = 0;
            continue;
          }
        }

        // 2.3. Watch out for zig-zag at 2nd pixel, from LR-TD scans
        if ((x - ii)*(ii - xloc[0]) +
            (y - jj)*(jj - yloc[0]) < 0) {
          xloc[1] = xloc[0], yloc[1] = yloc[0]; // swap first 2 points
          xloc[0] = ii, yloc[0] = jj; // to eliminate zig-zag
        }

        // 2.4. Collect both directions & extension points if 1-chain
        while (NextPixel(x, y, edgels))                // trace along first dir, 4-connected
          RecordPixel(x, y, edgels, xloc, yloc);       // and stronger first

        if (vcl_abs(xloc[0]-x) > 1 ||                  // disjoint first/last pixel
            vcl_abs(yloc[0]-y) > 1) {                  // so must be a 1-chain with end points
          if (next_pixel(x, y, *vertexMap))            // search for extra links to
            xloc.push_back(x), yloc.push_back(y);      // detected junctions
          x = xloc[0], y = yloc[0];                    // start again from first pixel

          vcl_vector<int> xloctemp( xloc.size()), yloctemp( yloc.size());
          for (unsigned int iii=0; iii< xloc.size(); iii++)
            xloctemp[iii]= xloc[xloc.size()-1-iii];
          for (unsigned int jjj=0; jjj< yloc.size(); jjj++)
            yloctemp[jjj]= yloc[yloc.size()-1-jjj];

          while (NextPixel(x, y, edgels)) // trace along other dir
            RecordPixel(x, y, edgels, xloc, yloc);
          if (next_pixel(x, y, *vertexMap)) // search for extra links to
            xloc.push_back(x), yloc.push_back(y); // detected junctions
        }
        int len = xloc.size();

        // 2.5. Check for isolated contours that are too short
        if (len < minLength) {  // zero or too few internal pixels
          for (int k = 0; k < len; k++) // zero or too few internal pixels
            floatPixel(edgels, xloc[k], yloc[k]) = 0; // prune short chains
          continue;
        }

        // Thin zig-zags on the contours? Zig-zags happen at the 2
        // end points because of extension, or inside the contour
        // because of 4/8-connected tracing through noisy chain pixels,
        // and large shifts for subpixel locations.
        // Defer the elimination of zig-zags to gevd_contour::SubPixelAccuracy()
        // or gevd_contour::EqualizeSpacing()

        // 2.6. Create network of chains, touching, possibly ending
        // at same junction, but never crossing one another
        vtol_edge_2d_sptr edge = new vtol_edge_2d();
        vdgl_edgel_chain * ec = new vdgl_edgel_chain;
        vdgl_interpolator * it = new vdgl_interpolator_linear(ec);
        vdgl_digital_curve * dc = new vdgl_digital_curve(it); // include end points

        for (int k=0; k< len; k++)
        {
          x= xloc[k];
          y= yloc[k];
          ec->add_edgel( vdgl_edgel( x, y));
          edgeMap->put(x, y, edge);
        }
        edge->set_curve(*dc);
        LookupTableInsert(edges, edge);
      }
  }

  // 3. Restore cache to original state
  for (int k = 0; k < njunction; k++)  // clear all void*/float labels
    vertexMap->put(junctionx[k], junctiony[k],NULL);
  for (int j = rmax; j <= ymax; j++)
    for (int i = rmax; i <= xmax; i++)
      if (floatPixel(edgels, i, j) < 0) // undo mark put by RecordPixel
        floatPixel(edgels, i, j) = - floatPixel(edgels, i, j);

  if (talkative)
    vcl_cout << "Find " << edges.size()
             << " chains/cycles, with pixels > " << minLength
             << " and strength > " << minStrength
#ifdef DEBUG
             << ", in " << t.real() << " msecs."
#endif
             << vcl_endl;
  return edges.size();  // number of chains found so far
}


//:
// Check that end point of a weak contour touches another stronger
// contour at an internal pixel. Localize the junction to pixel accuracy
// by searching for shortest distance from end point to chain.
// Gaussian smoothing can put local maximum change in filter response
// 1 pixel away from this junction location.
// Update junction map.
bool
DetectJunction(vtol_vertex_2d& end, int& index,
               vtol_edge_2d_sptr& weaker, vtol_edge_2d_sptr& stronger,
               const int maxSpiral,
               const gevd_bufferxy& edgels, vbl_array_2d<vtol_edge_2d_sptr>& edgeMap)
{
  // 0. Must be an end point of a dangling 1-chain
  if (end.numsup() > 1)         // avoid junction and 1-cycle
    return false;
  vcl_vector<vtol_edge_sptr> edges; end.edges(edges);
  weaker = edges[0]->cast_to_edge_2d();      // dangling edge must be a weaker contour
  vdgl_digital_curve_sptr dc = weaker->curve()->cast_to_vdgl_digital_curve();

  const int len = dc->get_interpolator()->get_edgel_chain()->size();

  // 1. Mark off pixels at end pt to find junction of a contour to itself

  const int rfuzz = vcl_min(len, 3*MINLENGTH);
  vtol_edge_2d_sptr* labels = new vtol_edge_2d_sptr[rfuzz];
  if (&end == weaker->v1()->cast_to_vertex_2d())
    for (int r = 0; r < rfuzz; r++) {
      vdgl_edgel edgel= dc->get_interpolator()->get_edgel_chain()->edgel( r);
      labels[r] = edgeMap.get( int(edgel.get_x()), int(edgel.get_y()));
      edgeMap.put(int(edgel.get_x()), int(edgel.get_y()), NULL);
    }
  else
    for (int r = 0; r < rfuzz; r++) {
      vdgl_edgel edgel= dc->get_interpolator()->get_edgel_chain()->edgel(len-1-r);
      labels[r] = edgeMap.get( int( edgel.get_x()), int( edgel.get_y()));
      edgeMap.put(int(edgel.get_x()), int(edgel.get_y()), NULL);
    }

  // 2. Find another stronger contour touched by this end point < gap.
  stronger = NULL;              // contour can join with itself
  int jx = int(end.x()), jy = int(end.y());
  for (int l = 0, n = 0; l < maxSpiral; l++) {  // increasing radius of spiral
    float maxpix = 0; int maxn = 0;     // strongest strength at this radius
    for (; n < RNS[l]; n++) {
      int x = jx+RIS[n], y = jy+RJS[n];
      if (edgeMap.get(x, y) && // find another contour or itself
          floatPixel(edgels, x, y) > maxpix) {
        maxpix = floatPixel(edgels, x, y);
        maxn = n;               // better neighbor
      }
    }
    if (maxpix) {               // location of junction on contour
      stronger = edgeMap.get(jx+RIS[maxn], jy+RJS[maxn]);
      jx += RIS[maxn], jy += RJS[maxn];
      break;
    }
  }
  // restore edgeMap around end point (undo step 1)
  if (&end == weaker->v1()->cast_to_vertex_2d())
    for (int r=0; r< rfuzz; r++) {
      vdgl_edgel edge= dc->get_interpolator()->get_edgel_chain()->edgel(r);
      edgeMap.put(int( edge.get_x()), int( edge.get_y()), labels[r]);
    }
  else
    for (int r=0; r< rfuzz; r++) {
      vdgl_edgel edgel= dc->get_interpolator()->get_edgel_chain()->edgel(len-1-r);
      edgeMap.put(int( edgel.get_x()), int( edgel.get_y()),labels[r]);
    }
  delete[] labels;

  if (!stronger)                // do not find any edge in search region
    return false;

  // 3. Find index location of junction on this contour
  index = int(INVALID);

  vdgl_digital_curve_sptr dc2 =(stronger->curve()->cast_to_vdgl_digital_curve());

  // find corresponding index on contour
  for (unsigned int n = 0; n < dc2->get_interpolator()->get_edgel_chain()->size(); ++n)
  {
    vdgl_edgel edgel= dc2->get_interpolator()->get_edgel_chain()->edgel(n);

    if ( int( edgel.get_x())== jx && int( edgel.get_y())== jy)
    {
      index = n;
      break;
    }
  }
  return true;
}


//: Confirm there is a strong jump in response near a junction.
// The location of this jump is however inaccurate, and so junctions
// can not be localized accurately along the stronger cycle.
static bool
ConfirmJunctionOnCycle(int index, float threshold,
                       vtol_edge_2d& cycle, const gevd_bufferxy& edgels)
{
  vdgl_digital_curve_sptr dc =(cycle.curve()->cast_to_vdgl_digital_curve());
  const int len = dc->get_interpolator()->get_edgel_chain()->size();
  const int wrap = 10*len;      // for positive index
  const int radius = 3;         // gap < 3, around junction pixel

  for (int n = index-radius; n <= index+radius; n++)
  {
    int nm = (n-1+wrap)%len;    // modulo operations to wrap at borders
    int np = (n+1+wrap)%len;

    vdgl_edgel edgel_m= dc->get_interpolator()->get_edgel_chain()->edgel( nm);
    vdgl_edgel edgel_p= dc->get_interpolator()->get_edgel_chain()->edgel( np);

    if (vcl_fabs(floatPixel(edgels, int( edgel_p.x()), int( edgel_p.y())) -
                 floatPixel(edgels, int( edgel_m.x()), int( edgel_m.y())))
        > threshold)
      return true;
  }
  return false;
}


//:
// Break the cycle at given index, and create new cycle from/to
// and not including index pixel. Update the chain map accordingly.
void
BreakCycle(vtol_vertex_2d& junction, const int index,
           vtol_edge_2d& stronger, vtol_edge_2d_sptr& split,
           vbl_array_2d<vtol_edge_2d_sptr>& edgeMap, vbl_array_2d<vtol_vertex_2d_sptr>& vertexMap)
{
  vdgl_digital_curve_sptr dc = (stronger.curve()->cast_to_vdgl_digital_curve());
  const int len = dc->get_interpolator()->get_edgel_chain()->size();

  // 1. Move location of junction
  int jx = int(junction.x()), jy = int(junction.y());
  vertexMap.put(jx, jy, NULL); // erase old location

  vdgl_edgel tempedgel= dc->get_interpolator()->get_edgel_chain()->edgel( index);
  jx = int( tempedgel.x()), jy = int( tempedgel.y());
  junction.set_x(jx), junction.set_y(jy); // update new location
  vertexMap.put(jx, jy, &junction);
  edgeMap.put(jx, jy, NULL);

  // 2. Create 1-cycle, including junction pixel
  split = new vtol_edge_2d();

  vdgl_edgel_chain *es= new vdgl_edgel_chain;
  vdgl_interpolator *it= new vdgl_interpolator_linear( vdgl_edgel_chain_sptr( es));
  vdgl_digital_curve *ds = new vdgl_digital_curve( vdgl_interpolator_sptr( it));

  split->set_curve(*ds->cast_to_curve());

  int i=0;
  for (int k = index; k < len; i++,k++) {
    vdgl_edgel e= dc->get_interpolator()->get_edgel_chain()->edgel( k);
    es->add_edgel( e);
    edgeMap.put(int(e.x()), int(e.y()), split);
  }
  for (int k = 0; i < len; i++,k++) {
    vdgl_edgel c= dc->get_interpolator()->get_edgel_chain()->edgel( k);
    es->add_edgel( c);
    edgeMap.put(int(c.x()), int(c.y()), split);
  }

  split->set_v1(&junction);     // link both directions v-e
  split->set_v2(&junction);
}


//: Confirm there is a strong jump in response near a junction.
// The location of this jump is however inaccurate, and so junctions
// can not be localized accurately along the stronger chain.
static bool
ConfirmJunctionOnChain(int index, float threshold,
                       vtol_edge_2d& chain, const gevd_bufferxy& edgels)
{
  vdgl_digital_curve_sptr dc = chain.curve()->cast_to_vdgl_digital_curve();
  const int len = dc->get_interpolator()->get_edgel_chain()->size()-1;

  if (len < 2*MINLENGTH-1) // will merge vertices instead of
    return false;               // breaking up chains

  const int fuzz = MINLENGTH-1; // from min length of broken chains
  const int radius = 3;         // gap < 3, around junction pixel

  for (int n = vcl_max(index-radius, fuzz); n <= vcl_min(index+radius,len-1-fuzz); n++)
  {
    vdgl_edgel cp1= dc->get_interpolator()->get_edgel_chain()->edgel(n+1);
    vdgl_edgel cm1= dc->get_interpolator()->get_edgel_chain()->edgel(n-1);

    if (vcl_fabs(floatPixel(edgels, int(cp1.x()), int(cp1.y())) -
                 floatPixel(edgels, int(cp1.x()), int(cm1.y())))
        > threshold)
    {
      return true;
    }
  }
  return false;
}


//: Break the edge at given index, and create two subchains from it.
// Update the chain map accordingly.
void
BreakChain(vtol_vertex_2d& junction, const int index,
           vtol_edge_2d& stronger,
           vtol_edge_2d_sptr& longer, vtol_edge_2d_sptr& shorter,
           vbl_array_2d<vtol_edge_2d_sptr>& edgeMap, vbl_array_2d<vtol_vertex_2d_sptr>& vertexMap)
{
  vdgl_digital_curve_sptr dc = stronger.curve()->cast_to_vdgl_digital_curve();
  vdgl_edgel_chain_sptr cxy= dc->get_interpolator()->get_edgel_chain();

  const int l0 = dc->get_interpolator()->get_edgel_chain()->size();
  const int l1 = index+1, l2 = l0-index;

  // 1. Move location of junction
  int jx = int(junction.x()), jy = int(junction.y());
  vertexMap.put(jx, jy, NULL); // erase old location
  vdgl_edgel c= dc->get_interpolator()->get_edgel_chain()->edgel( index);
  jx = int(c.x()), jy = int(c.y());
  junction.set_x(jx), junction.set_y(jy);       // update new location
  vertexMap.put(jx, jy, &junction);
  edgeMap.put( jx, jy, NULL);

  // 2. Create first subchain up to and including junction pixel.
  vtol_edge_2d_sptr edge1 = new vtol_edge_2d();    // create subchains, broken at junction.

  vdgl_edgel_chain *ec= new vdgl_edgel_chain;
  vdgl_interpolator *it= new vdgl_interpolator_linear( ec);
  vdgl_digital_curve *dc1 = new vdgl_digital_curve( it);

  edge1->set_curve(*dc1);

  vdgl_edgel_chain *cxy1= ec;

  for (int k = 0; k < l1; k++)
  {
    cxy1->add_edgel ( (*cxy)[k] );
    (*cxy1)[k] = (*cxy)[k];
    edgeMap.put( int((*cxy1)[k].x()), int((*cxy1)[k].y()),  edge1);
  }

  //vtol_vertex_2d_sptr  v1 = stronger.v1().ptr()->cast_to_vertex_2d();
  vtol_vertex_sptr  v1 = stronger.v1();

  if (v1->numsup() == 1)        // dangling chain with end pt at v1
    edgeMap.put( int((*cxy1)[0].x()), int((*cxy1)[0].y()), NULL);
  edge1->set_v1(v1.ptr());            // link both directions v-e
  edge1->set_v2(&junction);     // unlink when stronger.UnProtect()

  // 3. Create second subchain from and including junction pixel.
  vtol_edge_2d_sptr edge2 = new vtol_edge_2d();    // create second subchain

  vdgl_edgel_chain *ec2= new vdgl_edgel_chain;
  vdgl_interpolator *it2= new vdgl_interpolator_linear( ec2);
  vdgl_digital_curve *dc2= new vdgl_digital_curve( it2);

  edge2->set_curve(*dc2);

  vdgl_edgel_chain *cxy2= ec2;


  for (int k = 0; k < l2; k++)
  {
    cxy2->add_edgel( cxy->edgel( k+index));
    edgeMap.put( int((*cxy2)[k].x()), int((*cxy2)[k].y()), edge2);
  }

  //vtol_vertex_2d_sptr  v2 = stronger.v2().ptr()->cast_to_vertex_2d();
  vtol_vertex_sptr  v2 = stronger.v2();

  if (v2->numsup() == 1)        // dangling chain with end pt at v2
    edgeMap.put( int((*cxy2)[l2-1].x()), int((*cxy2)[l2-1].y()), NULL);

  edge2->set_v1(&junction);     // link both directions v-e
  edge2->set_v2(v2.ptr());            // unlink when stronger.UnProtect()

  if (l1 >= l2)                 // sort longer/shorter chains
    longer = edge1, shorter = edge2;
  else
    longer = edge2, shorter = edge1;
}


//: Break the chain at given index, and create a loop.
// Update the chain map accordingly.
void
LoopChain(vtol_vertex_2d& junction, const int index,
          vtol_edge_2d& chain,
          vtol_edge_2d_sptr& straight, vtol_edge_2d_sptr& curled,
          vbl_array_2d<vtol_edge_2d_sptr>& edgeMap, vbl_array_2d<vtol_vertex_2d_sptr>& vertexMap)
{
  vdgl_digital_curve_sptr dc = chain.curve()->cast_to_vdgl_digital_curve();

  vdgl_edgel_chain_sptr cxy= dc->get_interpolator()->get_edgel_chain();
  const int l0 = cxy->size();

  // 1. Move location of junction
  int jx = int(junction.x()), jy = int(junction.y());
  vertexMap.put(jx, jy, NULL); // erase old location
  jx = int((*cxy)[index].x()), jy = int((*cxy)[index].y());
  junction.set_x(jx), junction.set_y(jy);       // update new location
  vertexMap.put(jx, jy, &junction);
  edgeMap.put( jx, jy, NULL);

  // 1. Find straight/curled chains
  straight = new vtol_edge_2d(), curled = new vtol_edge_2d();
  const int l1 = index+1, l2 = l0-index;

  if (&junction == chain.v1()->cast_to_vertex_2d())
  { // first subchain is curled
    vdgl_edgel_chain *ec= new vdgl_edgel_chain;
    vdgl_interpolator *it= new vdgl_interpolator_linear( ec);
    vdgl_digital_curve *c= new vdgl_digital_curve( it);

    curled->set_curve(*c);

    vdgl_edgel_chain *xy= ec;
    for (int k = 0; k < l1; k++)
    {
    xy->add_edgel ( (*cxy)[k] );
    (*xy)[k] = (*cxy)[k];//, y[k] = cy[k];
      edgeMap.put( int((*xy)[k].x()), int((*xy)[k].y()), curled);
    }

    curled->set_v1(&junction);
    curled->set_v2(&junction);


    ec= new vdgl_edgel_chain;
    it= new vdgl_interpolator_linear( ec);
    c = new vdgl_digital_curve( it);    // second subchain is straight

    straight->set_curve(*c);

    xy= ec;

    for (int k = 0; k < l2; k++)
    {
      xy->add_edgel ( (*cxy)[k+index] );
      (*xy)[k] = (*cxy)[k+index];//, y[k] = dy[k];
      edgeMap.put( int((*xy)[k].x()), int((*xy)[k].y()), straight);
    }

    if (chain.v2()->numsup()==1)
    {
      edgeMap.put( int((*xy)[l2-1].x()), int((*xy)[l2-1].y()), NULL);
    }

    straight->set_v1(&junction);
    straight->set_v2(chain.v2().ptr());
  }
  else
  {                     // first subchain is straight
    vdgl_edgel_chain *ec= new vdgl_edgel_chain;
    vdgl_interpolator *it= new vdgl_interpolator_linear( ec);
    vdgl_digital_curve *c= new vdgl_digital_curve( it);

    straight->set_curve(*c);


    vdgl_edgel_chain *xy= ec;

    for (int k = 0; k < l1; k++)
    {
    xy->add_edgel ( (*cxy)[k] );
      (*xy)[k] = (*cxy)[k];//, y[k] = cy[k];
      edgeMap.put( int((*xy)[k].x()), int((*xy)[k].y()), straight);
    }

    if (chain.v1()->numsup()==1)
    {
      edgeMap.put( int((*xy)[0].x()), int((*xy)[0].y()), NULL);
    }

    straight->set_v1(chain.v1().ptr());
    straight->set_v2(&junction);


    ec= new vdgl_edgel_chain;
    it= new vdgl_interpolator_linear( ec);
    c = new vdgl_digital_curve( it);    // second subchain is curled

    curled->set_curve(*c);

    xy = ec; // ->GetX(), y = c->GetY();
    for (int k = 0; k < l2; k++)
    {
      xy->add_edgel ( (*cxy)[k+index] );
      (*xy)[k] = (*cxy)[k+index];//, y[k] = dy[k];
      edgeMap.put( int((*xy)[k].x()), int((*xy)[k].y()), curled);
    }

    curled->set_v1(&junction);
    curled->set_v2(&junction);
  }
}

//: Find number of rays connected to a vertex.
int
NumConnectedRays(vtol_vertex_2d& v)
{
  int nray = 0;
  vcl_vector<vtol_edge_sptr> segs; v.edges(segs);
  for (unsigned int i=0; i< segs.size(); i++)
  {
    if (segs[i]->v1()->cast_to_vertex_2d() == &v) nray++; // 1 for 1-chain
    if (segs[i]->v2()->cast_to_vertex_2d() == &v) nray++; // 2 for 1-cycle
  }
  return nray;
}


//: Detect touching another junction or end point, from an end point of a dangling chain.
vtol_vertex_2d_sptr
DetectTouch(const vtol_vertex_2d& end, const int maxSpiral,
            vbl_array_2d<vtol_vertex_2d_sptr>& vertexMap)
{
  const int jx = int(end.x()), jy = int(end.y());
  for (int l = 0, n = 0; l < maxSpiral; l++) {  // increasing radius of spiral
    vtol_vertex_2d_sptr  other = NULL;      // prefer junction over endpt
    int maxray = 0;            // largest number of rays
    for (; n < RNS[l]; n++) {  // 4- then 8-connected
      vtol_vertex_2d_sptr nbr = vertexMap.get(jx+RIS[n], jy+RJS[n]);
      int nray = (nbr ? NumConnectedRays(*nbr) : 0);
      if (nray > maxray) {
        maxray = nray;         // number of rays connected to it
        other = nbr;           // better neighbor
      }
    }
    if (maxray)                // find larger/other junction
      return other;
  }
  return NULL;
}


//: Find dangling edges connected to vertex
vtol_edge_2d_sptr
DanglingEdge(vtol_vertex_2d& v)
{
  vcl_vector<vtol_edge_sptr> segs; v.edges(segs);

  if (segs.size()==1)
    return segs[0]->cast_to_edge_2d();
  else
    return 0;
}


//: Merge 2 end points of a same chain.
// Update global maps.
void
MergeEndPtsOfChain(vtol_vertex_2d& endpt, vtol_vertex_2d& other, vtol_edge_2d& common,
                   vbl_array_2d<vtol_edge_2d_sptr>& edgeMap, vbl_array_2d<vtol_vertex_2d_sptr>&vertexMap)
{
  int px = int(other.x()), py = int(other.y());
  vertexMap.put(px, py, NULL); // erase old location
  edgeMap.put( px, py, &common);
  if (common.v1() == &other)  // remove links to other endpt
    common.set_v1(&endpt);
  else
    common.set_v2(&endpt);
}


//: Merge 2 touching chains into 1, deleting the 2 touching end points and their chains.
// Smooth away short kinks is delayed for later.  Update global maps.
void
MergeEndPtTouchingEndPt(vtol_vertex_2d& end1, vtol_vertex_2d