Add forgotten libavoid files

(bzr r8856)

10 files changed, 4714 insertions, 0 deletions

diff --git a/src/libavoid/.dirstamp b/src/libavoid/.dirstamp
new file mode 100644
index 000000000..e69de29bb
--- /dev/null
+++ b/src/libavoid/.dirstamp
diff --git a/src/libavoid/assertions.h b/src/libavoid/assertions.h
new file mode 100644
index 000000000..0725c4482
--- /dev/null
+++ b/src/libavoid/assertions.h
@@ -0,0 +1,49 @@
+/*
+ * vim: ts=4 sw=4 et tw=0 wm=0
+ *
+ * libavoid - Fast, Incremental, Object-avoiding Line Router
+ *
+ * Copyright (C) 2009  Monash University
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * See the file LICENSE.LGPL distributed with the library.
+ *
+ * Licensees holding a valid commercial license may use this file in
+ * accordance with the commercial license agreement provided with the 
+ * library.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ *
+ * Author(s):   Michael Wybrow <mjwybrow@users.sourceforge.net>
+*/
+
+#ifndef AVOID_ASSERTIONS_H
+#define AVOID_ASSERTIONS_H
+
+#ifdef NDEBUG 
+
+  #define COLA_ASSERT(expr)  static_cast<void>(0)
+
+#else // Not NDEBUG
+
+  #if defined(USE_ASSERT_EXCEPTIONS)
+
+    #include "libvpsc/assertions.h"
+
+  #else
+
+    #include <cassert>
+    #define COLA_ASSERT(expr)  assert(expr)
+
+  #endif
+
+#endif
+
+
+#endif // AVOID_ASSERTIONS_H
+
diff --git a/src/libavoid/geomtypes.cpp b/src/libavoid/geomtypes.cpp
new file mode 100644
index 000000000..10bb95a7a
--- /dev/null
+++ b/src/libavoid/geomtypes.cpp
@@ -0,0 +1,548 @@
+/*
+ * vim: ts=4 sw=4 et tw=0 wm=0
+ *
+ * libavoid - Fast, Incremental, Object-avoiding Line Router
+ *
+ * Copyright (C) 2004-2009  Monash University
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * See the file LICENSE.LGPL distributed with the library.
+ *
+ * Licensees holding a valid commercial license may use this file in
+ * accordance with the commercial license agreement provided with the 
+ * library.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ *
+ * Author(s):   Michael Wybrow <mjwybrow@users.sourceforge.net>
+*/
+
+
+#include <cmath>
+#include <cfloat>
+#include <cstdlib>
+
+#include "libavoid/geomtypes.h"
+#include "libavoid/shape.h"
+#include "libavoid/router.h"
+#include "libavoid/assertions.h"
+
+
+namespace Avoid
+{
+
+    
+Point::Point() :
+    id(0),
+    vn(kUnassignedVertexNumber)
+{
+}
+
+
+Point::Point(const double xv, const double yv) :
+    x(xv),
+    y(yv),
+    id(0),
+    vn(kUnassignedVertexNumber)
+{
+}
+
+
+bool Point::operator==(const Point& rhs) const
+{
+    if ((x == rhs.x) && (y == rhs.y))
+    {
+        return true;
+    }
+    return false;
+}
+
+
+bool Point::operator!=(const Point& rhs) const
+{
+    if ((x != rhs.x) || (y != rhs.y))
+    {
+        return true;
+    }
+    return false;
+}
+
+
+// Just defined to allow std::set<Point>.  Not particularly meaningful!
+bool Point::operator<(const Point& rhs) const
+{
+    if (x == rhs.x)
+    {
+        return (y < rhs.y);
+    }
+    return (x < rhs.x);
+}
+
+
+double& Point::operator[](const unsigned int dimension)
+{
+    COLA_ASSERT((dimension == 0) || (dimension == 1));
+    return ((dimension == 0) ? x : y);
+}
+
+
+const double& Point::operator[](const unsigned int dimension) const
+{
+    COLA_ASSERT((dimension == 0) || (dimension == 1));
+    return ((dimension == 0) ? x : y);
+}
+
+
+ReferencingPolygon::ReferencingPolygon(const Polygon& poly, const Router *router)
+    : PolygonInterface(),
+      _id(poly._id),
+      ps(poly.size())
+{
+    COLA_ASSERT(router != NULL);
+    for (size_t i = 0; i < poly.size(); ++i)
+    {
+        const Polygon *polyPtr = NULL;
+        for (ShapeRefList::const_iterator sh = router->shapeRefs.begin();
+                sh != router->shapeRefs.end(); ++sh) 
+        {
+            if ((*sh)->id() == poly.ps[i].id)
+            {
+                const Polygon& poly = (*sh)->polygon();
+                polyPtr = &poly;
+                break;
+            }
+        }
+        COLA_ASSERT(polyPtr != NULL);
+        ps[i] = std::make_pair(polyPtr, poly.ps[i].vn);
+    }
+}
+
+
+ReferencingPolygon::ReferencingPolygon()
+    : PolygonInterface()
+{
+    clear();
+}
+
+
+void ReferencingPolygon::clear(void)
+{
+    ps.clear();
+}
+
+
+bool ReferencingPolygon::empty(void) const
+{
+    return ps.empty();
+}
+
+
+size_t ReferencingPolygon::size(void) const
+{
+    return ps.size();
+}
+
+
+int ReferencingPolygon::id(void) const
+{
+    return _id;
+}
+
+
+const Point& ReferencingPolygon::at(size_t index) const 
+{
+    COLA_ASSERT(index < size());
+    const Polygon& poly = *(ps[index].first);
+    unsigned short poly_index = ps[index].second;
+    COLA_ASSERT(poly_index < poly.size());
+
+    return poly.ps[poly_index];
+}
+
+
+void PolygonInterface::getBoundingRect(double *minX, double *minY,
+        double *maxX, double *maxY) const
+{
+    double progressiveMinX = DBL_MAX;
+    double progressiveMinY = DBL_MAX;
+    double progressiveMaxX = -DBL_MAX;
+    double progressiveMaxY = -DBL_MAX;
+
+    for (size_t i = 0; i < size(); ++i)
+    {
+        progressiveMinX = std::min(progressiveMinX, at(i).x);
+        progressiveMinY = std::min(progressiveMinY, at(i).y);
+        progressiveMaxX = std::max(progressiveMaxX, at(i).x);
+        progressiveMaxY = std::max(progressiveMaxY, at(i).y);
+    }
+
+    if (minX)
+    {
+        *minX = progressiveMinX;
+    }
+    if (maxX)
+    {
+        *maxX = progressiveMaxX;
+    }
+    if (minY)
+    {
+        *minY = progressiveMinY;
+    }
+    if (maxY)
+    {
+        *maxY = progressiveMaxY;
+    }
+}
+
+
+Polygon::Polygon()
+    : PolygonInterface()
+{
+    clear();
+}
+
+
+Polygon::Polygon(const int pn)
+    : PolygonInterface(),
+      ps(pn)
+{
+}
+
+
+Polygon::Polygon(const PolygonInterface& poly)
+    : PolygonInterface(),
+      _id(poly.id()),
+      ps(poly.size())
+{
+    for (size_t i = 0; i < poly.size(); ++i)
+    {
+        ps[i] = poly.at(i);
+    }
+}
+
+
+void Polygon::clear(void)
+{
+    ps.clear();
+    ts.clear();
+}
+
+
+bool Polygon::empty(void) const
+{
+    return ps.empty();
+}
+
+
+size_t Polygon::size(void) const
+{
+    return ps.size();
+}
+
+
+int Polygon::id(void) const
+{
+    return _id;
+}
+
+
+const Point& Polygon::at(size_t index) const
+{
+    COLA_ASSERT(index < size());
+
+    return ps[index];
+}
+
+
+static const unsigned int SHORTEN_NONE  = 0;
+static const unsigned int SHORTEN_START = 1;
+static const unsigned int SHORTEN_END   = 2;
+static const unsigned int SHORTEN_BOTH  = SHORTEN_START | SHORTEN_END;
+
+// shorten_line():
+//     Given the two endpoints of a line segment, this function adjusts the
+//     endpoints of the line to shorten the line by shorten_length at either
+//     or both ends.
+//
+static void shorten_line(double& x1, double& y1, double& x2, double& y2, 
+        const unsigned int mode, const double shorten_length)
+{
+    if (mode == SHORTEN_NONE)
+    {
+        return;
+    }
+    
+    double rise = y1 - y2;
+    double run = x1 - x2;
+    double disty = fabs(rise);
+    double distx = fabs(run);
+
+    // Handle case where shorten length is greater than the length of the
+    // line segment.
+    if ((mode == SHORTEN_BOTH) &&
+            (((distx > disty) && ((shorten_length * 2) > distx)) ||
+             ((disty >= distx) && ((shorten_length * 2) > disty))))
+    {
+        x1 = x2 = x1 - (run / 2); 
+        y1 = y2 = y1 - (rise / 2); 
+        return;
+    }
+    else if ((mode == SHORTEN_START) && 
+            (((distx > disty) && (shorten_length > distx)) ||
+             ((disty >= distx) && (shorten_length > disty))))
+    {
+        x1 = x2;
+        y1 = y2;
+        return;
+    }
+    else if ((mode == SHORTEN_END) && 
+            (((distx > disty) && (shorten_length > distx)) ||
+             ((disty >= distx) && (shorten_length > disty))))
+    {
+        x2 = x1;
+        y2 = y1;
+        return;
+    }
+
+    // Handle orthogonal line segments.
+    if (x1 == x2)
+    {
+        // Vertical
+        int sign = (y1 < y2) ? 1: -1;
+        
+        if (mode & SHORTEN_START)
+        {
+            y1 += (sign * shorten_length);
+        }
+        if (mode & SHORTEN_END)
+        {
+            y2 -= (sign * shorten_length);
+        }
+        return;
+    }
+    else if (y1 == y2)
+    {
+        // Horizontal
+        int sign = (x1 < x2) ? 1: -1;
+        
+        if (mode & SHORTEN_START)
+        {
+            x1 += (sign * shorten_length);
+        }
+        if (mode & SHORTEN_END)
+        {
+            x2 -= (sign * shorten_length);
+        }
+        return;
+    }
+    
+    int xpos = (x1 < x2) ? -1 : 1;
+    int ypos = (y1 < y2) ? -1 : 1;
+    
+    double tangent = rise / run;
+   
+    if (mode & SHORTEN_END)
+    {
+        if (disty > distx)
+        {
+            y2 += shorten_length * ypos;
+            x2 += shorten_length * ypos * (1 / tangent);
+        }
+        else if (disty < distx)
+        {
+            y2 += shorten_length * xpos * tangent;
+            x2 += shorten_length * xpos;
+        }
+    }
+
+    if (mode & SHORTEN_START)
+    {
+        if (disty > distx)
+        {
+            y1 -= shorten_length * ypos;
+            x1 -= shorten_length * ypos * (1 / tangent);
+        }
+        else if (disty < distx)
+        {
+            y1 -= shorten_length * xpos * tangent;
+            x1 -= shorten_length * xpos;
+        }
+    }
+}
+
+
+void Polygon::translate(const double xDist, const double yDist)
+{
+    for (size_t i = 0; i < size(); ++i)
+    {
+        ps[i].x += xDist;
+        ps[i].y += yDist;
+    }
+}
+
+
+Polygon Polygon::simplify(void) const
+{
+    Polygon simplified = *this;
+    std::vector<Point>::iterator it = simplified.ps.begin();
+    if (it != simplified.ps.end()) ++it;
+
+    // Combine collinear line segments into single segments:
+    for (size_t j = 2; j < simplified.size(); )
+    {
+        if (vecDir(simplified.ps[j - 2], simplified.ps[j - 1], 
+                simplified.ps[j]) == 0)
+        {
+            // These three points make up two collinear segments, so just
+            // compine them into a single segment.
+            it = simplified.ps.erase(it);
+        }
+        else
+        {
+            ++j;
+            ++it;
+        }
+    }
+
+    return simplified;
+}
+
+
+#define mid(a, b) ((a < b) ? a + ((b - a) / 2) : b + ((a - b) / 2))
+
+
+// curvedPolyline():
+//     Returns a curved approximation of this multi-segment PolyLine, with 
+//     the corners replaced by smooth Bezier curves.  curve_amount describes
+//     how large to make the curves.
+//     The ts value for each point in the returned Polygon describes the 
+//     drawing operation: 'M' (move) marks the first point, a line segment 
+//     is marked with an 'L' and three 'C's (along with the previous point) 
+//     describe the control points of a Bezier curve.
+//
+Polygon Polygon::curvedPolyline(const double curve_amount,
+        const bool closed) const
+{
+    Polygon simplified = this->simplify();
+
+    Polygon curved;
+    size_t num_of_points = size();
+    if (num_of_points <= 2)
+    {
+        // There is only a single segment, do nothing.
+        curved = *this;
+        curved.ts.push_back('M');
+        curved.ts.push_back('L');
+        return curved;
+    }
+
+    // Build the curved polyline:
+    curved._id = _id;
+    double last_x = 0;
+    double last_y = 0;
+    if (closed)
+    {
+        double x1 = simplified.ps[0].x;
+        double y1 = simplified.ps[0].y;
+        double x2 = simplified.ps[1].x;
+        double y2 = simplified.ps[1].y;
+        shorten_line(x1, y1, x2, y2, SHORTEN_START, curve_amount);
+        curved.ps.push_back(Point(x1, y1));
+        curved.ts.push_back('M');
+    }
+    else
+    {
+        curved.ps.push_back(ps[0]);
+        curved.ts.push_back('M');
+    }
+   
+    size_t simpSize = simplified.size();
+    size_t finish = (closed) ? simpSize + 2 : simpSize;
+    for (size_t j = 1; j < finish; ++j)
+    {
+        double x1 = simplified.ps[(simpSize + j - 1) % simpSize].x;
+        double y1 = simplified.ps[(simpSize + j - 1) % simpSize].y;
+        double x2 = simplified.ps[j % simpSize].x;
+        double y2 = simplified.ps[j % simpSize].y;
+
+        double old_x = x1;
+        double old_y = y1;
+        
+        unsigned int mode = SHORTEN_BOTH;
+        if (!closed)
+        {
+            if (j == 1)
+            {
+                mode = SHORTEN_END;
+            }
+            else if (j == (size() - 1))
+            {
+                mode = SHORTEN_START;
+            }
+        }
+        shorten_line(x1, y1, x2, y2, mode, curve_amount);
+
+        if (j > 1)
+        {
+            curved.ts.insert(curved.ts.end(), 3, 'C');
+            curved.ps.push_back(Point(mid(last_x, old_x), mid(last_y, old_y)));
+            curved.ps.push_back(Point(mid(x1, old_x), mid(y1, old_y)));
+            curved.ps.push_back(Point(x1, y1));
+        }
+        if (closed && (j == (finish - 1)))
+        {
+            // Close the path.
+            curved.ts.push_back('Z');
+            curved.ps.push_back(Point(x1, y1));
+            break;
+        }
+        curved.ts.push_back('L');
+        curved.ps.push_back(Point(x2, y2));
+            
+        last_x = x2;
+        last_y = y2;
+    }
+    
+    return curved;
+}
+
+
+Rectangle::Rectangle(const Point& topLeft, const Point& bottomRight)
+    : Polygon(4)
+{
+    double xMin = std::min(topLeft.x, bottomRight.x);
+    double xMax = std::max(topLeft.x, bottomRight.x);
+    double yMin = std::min(topLeft.y, bottomRight.y);
+    double yMax = std::max(topLeft.y, bottomRight.y);
+
+    ps[0] = Point(xMax, yMin);
+    ps[1] = Point(xMax, yMax);
+    ps[2] = Point(xMin, yMax);
+    ps[3] = Point(xMin, yMin);
+}
+
+
+Rectangle::Rectangle(const Point& centre, const double width, 
+        const double height)
+{
+    double halfWidth  = width / 2.0;
+    double halfHeight = height / 2.0;
+    double xMin = centre.x - halfWidth;
+    double xMax = centre.x + halfWidth;
+    double yMin = centre.y - halfHeight;
+    double yMax = centre.y + halfHeight;
+
+    ps[0] = Point(xMax, yMin);
+    ps[1] = Point(xMax, yMax);
+    ps[2] = Point(xMin, yMax);
+    ps[3] = Point(xMin, yMin);
+}
+
+
+}
+
diff --git a/src/libavoid/makefile b/src/libavoid/makefile
new file mode 100644
index 000000000..e4f83a52d
--- /dev/null
+++ b/src/libavoid/makefile
@@ -0,0 +1,17 @@
+# Convenience stub makefile to call the real Makefile.
+
+
+
+OBJEXT = o
+
+# Explicit so that it's the default rule.
+all:
+	cd .. && $(MAKE) libavoid/all
+
+clean %.a %.$(OBJEXT):
+	cd .. && $(MAKE) libavoid/$@
+
+.PHONY: all clean
+
+.SUFFIXES:
+.SUFFIXES: .a .$(OBJEXT)
diff --git a/src/libavoid/orthogonal.cpp b/src/libavoid/orthogonal.cpp
new file mode 100644
index 000000000..747fd1f86
--- /dev/null
+++ b/src/libavoid/orthogonal.cpp
@@ -0,0 +1,2343 @@
+/*
+ * vim: ts=4 sw=4 et tw=0 wm=0
+ *
+ * libavoid - Fast, Incremental, Object-avoiding Line Router
+ *
+ * Copyright (C) 2009  Monash University
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * See the file LICENSE.LGPL distributed with the library.
+ *
+ * Licensees holding a valid commercial license may use this file in
+ * accordance with the commercial license agreement provided with the 
+ * library.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ *
+ * Author(s):   Michael Wybrow <mjwybrow@users.sourceforge.net>
+*/
+
+
+#include <cstdlib>
+#include <cfloat>
+#include <cmath>
+#include <set>
+#include <list>
+#include <algorithm>
+
+#include "libavoid/router.h"
+#include "libavoid/geomtypes.h"
+#include "libavoid/shape.h"
+#include "libavoid/orthogonal.h"
+#include "libavoid/connector.h"
+#include "libavoid/vpsc.h"
+#include "libavoid/assertions.h"
+
+#ifdef LIBAVOID_SDL
+  #include <SDL_gfxPrimitives.h>
+#endif
+
+
+namespace Avoid {
+
+
+static const double CHANNEL_MAX = 100000000;
+
+static const size_t XDIM = 0;
+static const size_t YDIM = 1;
+
+
+class ShiftSegment 
+{
+    public:
+        // For shiftable segments.
+        ShiftSegment(ConnRef *conn, const size_t low, const size_t high, 
+                bool isSBend, const size_t dim, double minLim, double maxLim)
+            : connRef(conn),
+              indexLow(low),
+              indexHigh(high),
+              sBend(isSBend),
+              fixed(false),
+              dimension(dim),
+              variable(NULL),
+              minSpaceLimit(minLim),
+              maxSpaceLimit(maxLim)
+        {
+        }
+        // For fixed segments.
+        ShiftSegment(ConnRef *conn, const size_t low, const size_t high, 
+                const size_t dim)
+            : connRef(conn),
+              indexLow(low),
+              indexHigh(high),
+              sBend(false),
+              fixed(true),
+              dimension(dim),
+              variable(NULL)
+        {
+            // This has no space to shift.
+            minSpaceLimit = lowPoint()[dim];
+            maxSpaceLimit = lowPoint()[dim];
+        }
+        Point& lowPoint(void)
+        {
+            return connRef->displayRoute().ps[indexLow];
+        }
+        Point& highPoint(void)
+        {
+            return connRef->displayRoute().ps[indexHigh];
+        }
+        const Point& lowPoint(void) const
+        {
+            return connRef->displayRoute().ps[indexLow];
+        }
+        const Point& highPoint(void) const 
+        {
+            return connRef->displayRoute().ps[indexHigh];
+        }
+        const int fixedOrder(bool& isFixed) const
+        {
+            if (fixed)
+            {
+                isFixed = true;
+                return 0;
+            }
+            if (lowC())
+            {
+                return 1;
+            }
+            else if (highC())
+            {
+                return -1;
+            }
+            return 0;
+        }
+        const int order(void) const
+        {
+            if (lowC())
+            {
+                return -1;
+            }
+            else if (highC())
+            {
+                return 1;
+            }
+            return 0;
+        }
+        bool operator<(const ShiftSegment& rhs) const
+        {
+            const Point& lowPt = lowPoint();
+            const Point& rhsLowPt = rhs.lowPoint();
+            
+            if (lowPt[dimension] != rhsLowPt[dimension])
+            {
+                return lowPt[dimension] < rhsLowPt[dimension];
+            }
+            return this < &rhs;
+        }
+        // This counts segments that are colliear and share an endpoint as
+        // overlapping.  This allows them to be nudged apart where possible.
+        bool overlapsWith(const ShiftSegment& rhs, const size_t dim) const
+        {
+            size_t altDim = (dim + 1) % 2;
+            const Point& lowPt = lowPoint();
+            const Point& highPt = highPoint();
+            const Point& rhsLowPt = rhs.lowPoint();
+            const Point& rhsHighPt = rhs.highPoint();
+            if ( (lowPt[altDim] <= rhsHighPt[altDim]) &&
+                    (rhsLowPt[altDim] <= highPt[altDim]))
+            {
+                if ( (minSpaceLimit <= rhs.maxSpaceLimit) &&
+                        (rhs.minSpaceLimit <= maxSpaceLimit))
+                {
+                    return true;
+                }
+            }
+            return false;
+        }
+
+        ConnRef *connRef;
+        size_t indexLow;
+        size_t indexHigh;
+        bool sBend;
+        bool fixed;
+        size_t dimension;
+        Variable *variable;
+        double minSpaceLimit;
+        double maxSpaceLimit;
+    private:
+        const bool lowC(void) const
+        {
+            // This is true if this is a cBend and its adjoining points
+            // are at lower positions.
+            if (!sBend && !fixed && (minSpaceLimit == lowPoint()[dimension]))
+            {
+                return true;
+            }
+            return false;
+        }
+        const bool highC(void) const
+        {
+            // This is true if this is a cBend and its adjoining points
+            // are at higher positions.
+            if (!sBend && !fixed && (maxSpaceLimit == lowPoint()[dimension]))
+            {
+                return true;
+            }
+            return false;
+        }
+};
+typedef std::list<ShiftSegment> ShiftSegmentList;
+
+bool cmpShiftSegment(const ShiftSegment& u, const ShiftSegment& v)
+{
+    return u < v;
+}
+
+
+struct Node;
+struct CmpNodePos { bool operator()(const Node* u, const Node* v) const; };
+
+typedef std::set<Node*,CmpNodePos> NodeSet;
+struct Node 
+{
+    ShapeRef *v;
+    VertInf *c;
+    ShiftSegment *ss;
+    double pos;
+    double min[2], max[2];
+    Node *firstAbove, *firstBelow;
+    NodeSet::iterator iter;
+
+    Node(ShapeRef *v, const double p)
+        : v(v),
+          c(NULL),
+          ss(NULL),
+          pos(p),
+          firstAbove(NULL),
+          firstBelow(NULL)
+    {   
+        //COLA_ASSERT(r->width()<1e40);
+        v->polygon().getBoundingRect(&min[0], &min[1], &max[0], &max[1]);
+    }   
+    Node(VertInf *c, const double p)
+        : v(NULL),
+          c(c),
+          ss(NULL),
+          pos(p),
+          firstAbove(NULL),
+          firstBelow(NULL)
+    {
+        min[0] = max[0] = c->point.x;
+        min[1] = max[1] = c->point.y;
+    }   
+    Node(ShiftSegment *ss, const double p)
+        : v(NULL),
+          c(NULL),
+          ss(ss),
+          pos(p),
+          firstAbove(NULL),
+          firstBelow(NULL)
+    {
+        // These values shouldn't ever be used, so they don't matter.
+        min[0] = max[0] = min[1] = max[1] = 0;
+    }   
+    ~Node() 
+    {
+    }
+    // Find the first Node above in the scanline that is a shape edge,
+    // and does not have an open or close event at this position (meaning
+    // it is just about to be removed).
+    double firstObstacleAbove(size_t dim)
+    {
+        Node *curr = firstAbove;
+        while (curr && (curr->ss || (curr->max[dim] > pos)))
+        {
+            curr = curr->firstAbove;
+        }
+       
+        if (curr)
+        {
+            return curr->max[dim];
+        }
+        return -DBL_MAX;
+    }
+    // Find the first Node below in the scanline that is a shape edge,
+    // and does not have an open or close event at this position (meaning
+    // it is just about to be removed).
+    double firstObstacleBelow(size_t dim)
+    {
+        Node *curr = firstBelow;
+        while (curr && (curr->ss || (curr->min[dim] < pos)))
+        {
+            curr = curr->firstBelow;
+        }
+        
+        if (curr)
+        {
+            return curr->min[dim];
+        }
+        return DBL_MAX;
+    }
+    // Mark all connector segments above in the scanline as being able 
+    // to see to this shape edge.
+    void markShiftSegmentsAbove(size_t dim)
+    {
+        Node *curr = firstAbove;
+        while (curr && (curr->ss || (curr->pos > min[dim])))
+        {
+            if (curr->ss && (curr->pos <= min[dim]))
+            {
+                curr->ss->maxSpaceLimit = 
+                        std::min(min[dim], curr->ss->maxSpaceLimit);
+            }
+            curr = curr->firstAbove;
+        }
+    }
+    // Mark all connector segments below in the scanline as being able 
+    // to see to this shape edge.
+    void markShiftSegmentsBelow(size_t dim)
+    {
+        Node *curr = firstBelow;
+        while (curr && (curr->ss || (curr->pos < max[dim])))
+        {
+            if (curr->ss && (curr->pos >= max[dim]))
+            {
+                curr->ss->minSpaceLimit = 
+                        std::max(max[dim], curr->ss->minSpaceLimit);
+            }
+            curr = curr->firstBelow;
+        }
+    }
+    bool findFirstPointAboveAndBelow(const size_t dim, double& firstAbovePos,
+            double& firstBelowPos, double& lastAbovePos, double& lastBelowPos)
+    {
+        bool clearVisibility = true;
+        firstAbovePos = -DBL_MAX;
+        firstBelowPos = DBL_MAX;
+        // We start looking left from the right side of the shape, 
+        // and vice versa. 
+        lastAbovePos = max[dim];
+        lastBelowPos = min[dim];
+
+        // Find the first blocking edge above this point.  Don't count the
+        // edges as we are travelling out of shapes we are inside, but then
+        // mark clearVisibility as false.
+        Node *curr = firstAbove;
+        while (curr && (curr->max[dim] > min[dim]))
+        {
+            lastAbovePos = std::min(curr->min[dim], lastAbovePos);
+            if ((curr->max[dim] >= min[dim]) && (curr->max[dim] <= max[dim]))
+            {
+                lastAbovePos = std::min(curr->max[dim], lastAbovePos);
+            }
+            lastBelowPos = std::max(curr->max[dim], lastBelowPos);
+            clearVisibility = false;
+            curr = curr->firstAbove;
+        }
+        if (curr)
+        {
+            firstAbovePos = curr->max[dim];
+        }
+        while (curr)
+        {
+            // There might be a larger shape after this one in the ordering.
+            if (curr->max[dim] < min[dim])
+            {
+                firstAbovePos = std::max(curr->max[dim], firstAbovePos);
+            }
+            curr = curr->firstAbove;
+        }
+        
+        // Find the first blocking edge below this point.  Don't count the
+        // edges as we are travelling out of shapes we are inside, but then
+        // mark clearVisibility as false.
+        curr = firstBelow;
+        while (curr && (curr->min[dim] < max[dim]))
+        {
+            lastBelowPos = std::max(curr->max[dim], lastBelowPos);
+            if ((curr->min[dim] >= min[dim]) && (curr->min[dim] <= max[dim]))
+            {
+                lastBelowPos = std::max(curr->min[dim], lastBelowPos);
+            }
+            lastAbovePos = std::min(curr->min[dim], lastAbovePos);
+            clearVisibility = false;
+            curr = curr->firstBelow;
+        }
+        if (curr)
+        {
+            firstBelowPos = curr->min[dim];
+        }
+        while (curr)
+        {
+            // There might be a larger shape after this one in the ordering.
+            if (curr->min[dim] > max[dim])
+            {
+                firstBelowPos = std::min(curr->min[dim], firstBelowPos);
+            }
+            curr = curr->firstBelow;
+        }
+
+        return clearVisibility;
+    }
+    double firstPointAbove(size_t dim) 
+    { 
+        Node *curr = firstAbove; 
+        while (curr && (curr->max[dim] >= pos)) 
+        { 
+            curr = curr->firstAbove; 
+        } 
+        
+        if (curr) 
+        { 
+            return curr->max[dim]; 
+        } 
+        return -DBL_MAX; 
+    } 
+    double firstPointBelow(size_t dim) 
+    { 
+        Node *curr = firstBelow; 
+        while (curr && (curr->min[dim] <= pos)) 
+        { 
+            curr = curr->firstBelow; 
+        } 
+         
+        if (curr) 
+        { 
+            return curr->min[dim]; 
+        } 
+        return DBL_MAX; 
+    } 
+    // This is a bit inefficient, but we won't need to do it once we have 
+    // connection points.
+    bool isInsideShape(size_t dimension)
+    {
+        for (Node *curr = firstBelow; curr; curr = curr->firstBelow)
+        {
+            if ((curr->min[dimension] < pos) && (pos < curr->max[dimension]))
+            {
+                return true;
+            }
+        }
+        for (Node *curr = firstAbove; curr; curr = curr->firstAbove)
+        {
+            if ((curr->min[dimension] < pos) && (pos < curr->max[dimension]))
+            {
+                return true;
+            }
+        }
+        return false;
+    }
+};
+
+
+bool CmpNodePos::operator() (const Node* u, const Node* v) const 
+{
+    if (u->pos != v->pos) 
+    {
+        return u->pos < v->pos;
+    }
+    
+    // Use the pointers to the base objects to differentiate them.
+    void *up = (u->v) ? (void *) u->v : 
+            ((u->c) ? (void *) u->c : (void *) u->ss);
+    void *vp = (v->v) ? (void *) v->v : 
+            ((v->c) ? (void *) v->c : (void *) v->ss);
+    return up < vp;
+}
+
+
+// Note: Open must come first.
+typedef enum {
+    Open = 1,
+    SegOpen = 2,
+    ConnPoint = 3, 
+    SegClose = 4,
+    Close = 5
+} EventType;
+
+
+struct Event
+{
+    Event(EventType t, Node *v, double p) 
+        : type(t),
+          v(v),
+          pos(p)
+    {};
+    EventType type;
+    Node *v;
+    double pos;
+};
+
+Event **events;
+
+
+// Used for quicksort.  Must return <0, 0, or >0.
+int compare_events(const void *a, const void *b)
+{
+	Event *ea = *(Event**) a;
+	Event *eb = *(Event**) b;
+    if (ea->pos != eb->pos)
+    {
+        return (ea->pos < eb->pos) ? -1 : 1;
+    }
+    if (ea->type != eb->type)
+    {
+        return ea->type - eb->type;
+    }
+    COLA_ASSERT(ea->v != eb->v);
+    return ea->v - eb->v;
+}
+
+
+// Returns a bitfield of the direction of visibility (in this dimension)
+// made up of ConnDirDown (for visibility towards lower position values) 
+// and ConnDirUp (for visibility towards higher position values).
+//
+static ConnDirFlags getPosVertInfDirection(VertInf *v, size_t dim)
+{
+    if (dim == XDIM) // X-dimension
+    {
+        unsigned int dirs = v->visDirections & (ConnDirLeft | ConnDirRight);
+        if (dirs == (ConnDirLeft | ConnDirRight))
+        {
+            return (ConnDirDown | ConnDirUp);
+        }
+        else if (dirs == ConnDirLeft)
+        {
+            return ConnDirDown;
+        }
+        else if (dirs == ConnDirRight)
+        {
+            return ConnDirUp;
+        }
+    }
+    else if (dim == YDIM) // Y-dimension
+    {
+        unsigned int dirs = v->visDirections & (ConnDirDown | ConnDirUp);
+        if (dirs == (ConnDirDown | ConnDirUp))
+        {
+            return (ConnDirDown | ConnDirUp);
+        }
+        else if (dirs == ConnDirDown)
+        {
+            // For libavoid the Y-axis points downwards, so in terms of 
+            // smaller or larger position values, Down is Up and vice versa.
+            return ConnDirUp;
+        }
+        else if (dirs == ConnDirUp)
+        {
+            // For libavoid the Y-axis points downwards, so in terms of 
+            // smaller or larger position values, Down is Up and vice versa.
+            return ConnDirDown;
+        }
+    }
+
+    // Can occur for ConnDirNone visibility.
+    return ConnDirNone;
+}
+
+
+struct PosVertInf
+{
+    PosVertInf(double p, VertInf *vI, ConnDirFlags d = ConnDirNone)
+        : pos(p),
+          vert(vI),
+          dir(d)
+    {
+    }
+    
+    bool operator<(const PosVertInf& rhs) const 
+    {
+        if (pos != rhs.pos)
+        {
+            return pos < rhs.pos;
+        }
+        return vert < rhs.vert;
+    }
+
+    double pos;
+    VertInf *vert;
+
+    // A bitfield marking the direction of visibility (in this dimension)
+    // made up of ConnDirDown (for visibility towards lower position values) 
+    // and ConnDirUp (for visibility towards higher position values).
+    //
+    ConnDirFlags dir;
+};
+
+
+struct CmpVertInf
+{ 
+    bool operator()(const VertInf* u, const VertInf* v) const
+    {
+        // Comparator for VertSet, an ordered set of VertInf pointers.
+        // It is assumed vertical sets of points will all have the same
+        // x position and horizontal sets all share a y position, so this
+        // method can be used to sort both these sets.
+        COLA_ASSERT((u->point.x == v->point.x) || (u->point.y == v->point.y));
+        if (u->point.x != v->point.x)
+        {
+            return u->point.x < v->point.x;
+        }
+        else if (u->point.y != v->point.y)
+        {
+            return u->point.y < v->point.y;
+        }
+        return u < v;
+    }
+};
+
+
+typedef std::set<VertInf *, CmpVertInf> VertSet;
+
+// A set of points to break the line segment, 
+// along with vertices for these points.
+typedef std::set<PosVertInf> BreakpointSet;
+
+// Temporary structure used to store the possible horizontal visibility 
+// lines arising from the vertical sweep.
+class LineSegment 
+{
+public:
+    LineSegment(const double& b, const double& f, const double& p, 
+            bool ss = false, VertInf *bvi = NULL, VertInf *fvi = NULL)
+        : begin(b),
+          finish(f),
+          pos(p),
+          shapeSide(false)
+    {
+        COLA_ASSERT(begin < finish);
+
+        if (bvi)
+        {
+            vertInfs.insert(bvi);
+        }
+        if (fvi)
+        {
+            vertInfs.insert(fvi);
+        }
+    }
+    LineSegment(const double& bf, const double& p, VertInf *bfvi = NULL)
+        : begin(bf),
+          finish(bf),
+          pos(p),
+          shapeSide(false)
+    {
+        if (bfvi)
+        {
+            vertInfs.insert(bfvi);
+        }
+    }
+ 
+    // Order by begin, pos, finish.
+    bool operator<(const LineSegment& rhs) const 
+    {
+        if (begin != rhs.begin)
+        {
+            return begin < rhs.begin;
+        }
+        if (pos != rhs.pos)
+        {
+            return pos < rhs.pos;
+        }
+        if (finish != rhs.finish)
+        {
+            return finish < rhs.finish;
+        }
+        COLA_ASSERT(shapeSide == rhs.shapeSide);
+        return false;
+    }
+
+    bool overlaps(const LineSegment& rhs) const
+    {
+        if ((begin == rhs.begin) && (pos == rhs.pos) &&
+                (finish == rhs.finish))
+        {
+            // Lines are exactly equal.
+            return true;
+        }
+        
+        if (pos == rhs.pos)
+        {
+            if (((begin >= rhs.begin) && (begin <= rhs.finish)) ||
+                ((rhs.begin >= begin) && (rhs.begin <= finish)) )
+            {
+                // They are colinear and overlap by some amount.
+                return true;
+            }
+        }
+        return false;
+    }
+
+    void mergeVertInfs(const LineSegment& segment)
+    {
+        begin = std::min(begin, segment.begin);
+        finish = std::max(finish, segment.finish);
+        vertInfs.insert(segment.vertInfs.begin(), segment.vertInfs.end());
+    }
+    
+    VertInf *beginVertInf(void) const
+    {
+        if (vertInfs.empty())
+        {
+            return NULL;
+        }
+        return *vertInfs.begin();
+    }
+    VertInf *finishVertInf(void) const
+    {
+        if (vertInfs.empty())
+        {
+            return NULL;
+        }
+        return *vertInfs.rbegin();
+    }
+
+    VertInf *commitPositionX(Router *router, double posX)
+    {
+        VertInf *found = NULL;
+        for (VertSet::iterator v = vertInfs.begin();
+                v != vertInfs.end(); ++v)
+        {
+            if ((*v)->point.x == posX)
+            {
+                found = *v;
+                break;
+            }
+        }
+        if (!found)
+        {
+            found = new VertInf(router, dummyOrthogID, Point(posX, pos));
+            vertInfs.insert(found);
+        }
+        return found;
+    }
+    // Set begin endpoint vertex if none has been assigned.
+    void commitBegin(Router *router, VertInf *vert = NULL)
+    {
+        if (vert)
+        {
+            vertInfs.insert(vert);
+        }
+
+        if (vertInfs.empty() ||
+                ((*vertInfs.begin())->point.x != begin))
+        {
+            vertInfs.insert(new
+                    VertInf(router, dummyOrthogID, Point(begin, pos)));
+        }
+    }
+
+    // Set begin endpoint vertex if none has been assigned.
+    void commitFinish(Router *router, VertInf *vert = NULL)
+    {
+        if (vert)
+        {
+            vertInfs.insert(vert);
+        }
+
+        if (vertInfs.empty() ||
+                ((*vertInfs.rbegin())->point.x != finish))
+        {
+            vertInfs.insert(new
+                    VertInf(router, dummyOrthogID, Point(finish, pos)));
+        }
+    }
+
+    // Converts a section of the points list to a set of breakPoints.  
+    // Returns the first of the intersection points occuring at finishPos.
+    VertSet::iterator addSegmentsUpTo(Router *router, double finishPos)
+    {
+        VertSet::iterator firstIntersectionPt = vertInfs.end();
+        for (VertSet::iterator vert = vertInfs.begin(); 
+                vert != vertInfs.end(); ++vert)
+        {
+            if ((*vert)->point.x > finishPos)
+            {
+                // We're done.
+                break;
+            }
+            
+            breakPoints.insert(PosVertInf((*vert)->point.x, (*vert),
+                        getPosVertInfDirection(*vert, XDIM)));
+
+            if ((firstIntersectionPt == vertInfs.end()) && 
+                    ((*vert)->point.x == finishPos))
+            {
+                firstIntersectionPt = vert;
+            }
+        }
+        // Returns the first of the intersection points at finishPos.
+        return firstIntersectionPt;
+    }
+
+    // Add visibility edge(s) for this segment.  There may be multiple if 
+    // one of the endpoints is shared by multiple connector endpoints.
+    void addEdgeHorizontal(Router *router)
+    {
+        commitBegin(router);
+        commitFinish(router);
+        
+        addSegmentsUpTo(router, finish);
+    }
+
+    // Add visibility edge(s) for this segment up until an intersection.
+    // Then, move the segment beginning to the intersection point, so we
+    // later only consider the remainder of the segment.
+    // There may be multiple segments added to the graph if the beginning 
+    // endpoint of the segment is shared by multiple connector endpoints.
+    VertSet addEdgeHorizontalTillIntersection(Router *router, 
+            LineSegment& vertLine)
+    {
+        VertSet intersectionSet;
+
+        commitBegin(router);
+
+        // Does a vertex already exist for this point.
+        commitPositionX(router, vertLine.pos);
+   
+        // Generate segments and set end iterator to the first point 
+        // at the intersection position.
+        VertSet::iterator restBegin = addSegmentsUpTo(router, vertLine.pos);
+
+        // Add the intersections points to intersectionSet.
+        VertSet::iterator restEnd = restBegin;
+        while ((restEnd != vertInfs.end()) && 
+                (*restEnd)->point.x == vertLine.pos)
+        {
+            ++restEnd;
+        }
+        intersectionSet.insert(restBegin, restEnd);
+
+        // Adjust segment to remove processed portion.
+        begin = vertLine.pos;
+        vertInfs.erase(vertInfs.begin(), restBegin);
+
+        return intersectionSet;
+    }
+                
+    // Insert vertical breakpoints.
+    void insertBreakpointsBegin(Router *router, LineSegment& vertLine)
+    {
+        VertInf *vert = NULL;
+        if (pos == vertLine.begin && vertLine.beginVertInf())
+        {
+            vert = vertLine.beginVertInf();
+        }
+        else if (pos == vertLine.finish && vertLine.finishVertInf())
+        {
+            vert = vertLine.finishVertInf();
+        }
+        commitBegin(router, vert);
+
+        for (VertSet::iterator v = vertInfs.begin();
+                v != vertInfs.end(); ++v)
+        {
+            if ((*v)->point.x == begin)
+            {
+                vertLine.breakPoints.insert(PosVertInf(pos, *v, 
+                        getPosVertInfDirection(*v, YDIM)));
+            }
+        }
+    }
+
+    // Insert vertical breakpoints.
+    void insertBreakpointsFinish(Router *router, LineSegment& vertLine)
+    {
+        VertInf *vert = NULL;
+        if (pos == vertLine.begin && vertLine.beginVertInf())
+        {
+            vert = vertLine.beginVertInf();
+        }
+        else if (pos == vertLine.finish && vertLine.finishVertInf())
+        {
+            vert = vertLine.finishVertInf();
+        }
+        commitFinish(router, vert);
+
+        for (VertSet::iterator v = vertInfs.begin();
+                v != vertInfs.end(); ++v)
+        {
+            if ((*v)->point.x == finish)
+            {
+                vertLine.breakPoints.insert(PosVertInf(pos, *v,
+                        getPosVertInfDirection(*v, YDIM)));
+            }
+        }
+    }
+    void generateVisibilityEdgesFromBreakpointSet(Router *router, size_t dim)
+    {
+        if ((breakPoints.begin())->pos != begin)
+        {
+            if (!beginVertInf())
+            {
+                Point point(pos, pos);
+                point[dim] = begin;
+                // Add begin point if it didn't intersect another line.
+                VertInf *vert = new VertInf(router, dummyOrthogID, point);
+                breakPoints.insert(PosVertInf(begin, vert));
+            }
+        }
+        if ((breakPoints.rbegin())->pos != finish)
+        {
+            if (!finishVertInf())
+            {
+                Point point(pos, pos);
+                point[dim] = finish;
+                // Add finish point if it didn't intersect another line.
+                VertInf *vert = new VertInf(router, dummyOrthogID, point);
+                breakPoints.insert(PosVertInf(finish, vert));
+            }
+        }
+
+        const bool orthogonal = true;
+        BreakpointSet::iterator vert, last;
+        for (vert = last = breakPoints.begin(); vert != breakPoints.end();)
+        {
+            BreakpointSet::iterator firstPrev = last;
+            while (last->vert->point[dim] != vert->vert->point[dim])
+            {
+                COLA_ASSERT(vert != last);
+                // Assert points are not at the same position.
+                COLA_ASSERT(vert->vert->point != last->vert->point);
+
+                if ( !(vert->vert->id.isShape || last->vert->id.isShape))
+                {
+                    // Here we have a pair of two endpoints that are both
+                    // connector endpoints and both are inside a shape.
+                    
+                    // Give vert visibility back to the first non-connector
+                    // endpoint vertex (i.e., the side of the shape).
+                    BreakpointSet::iterator side = last;
+                    while (!side->vert->id.isShape)
+                    {
+                        if (side == breakPoints.begin())
+                        {
+                            break;
+                        }
+                        --side;
+                    }
+                    bool canSeeDown = (vert->dir & ConnDirDown);
+                    if (canSeeDown && side->vert->id.isShape)
+                    {
+                        EdgeInf *edge = new 
+                                EdgeInf(side->vert, vert->vert, orthogonal);
+                        edge->setDist(vert->vert->point[dim] - 
+                                side->vert->point[dim]);
+                    }
+
+                    // Give last visibility back to the first non-connector
+                    // endpoint vertex (i.e., the side of the shape).
+                    side = vert;
+                    while ((side != breakPoints.end()) && 
+                            !side->vert->id.isShape)
+                    {
+                        ++side;
+                    }
+                    bool canSeeUp = (last->dir & ConnDirUp);
+                    if (canSeeUp && (side != breakPoints.end()))
+                    {
+                        EdgeInf *edge = new 
+                                EdgeInf(last->vert, side->vert, orthogonal);
+                        edge->setDist(side->vert->point[dim] - 
+                                last->vert->point[dim]);
+                    }
+                }
+                
+                // The normal case.
+                //
+                // Note: It's okay to give two connector endpoints visbility 
+                // here since we only consider the partner endpoint as a 
+                // candidate while searching if it is the other endpoint of
+                // the connector in question.
+                //
+                bool generateEdge = true;
+                if (!last->vert->id.isShape && !(last->dir & ConnDirUp))
+                {
+                    generateEdge = false;
+                }
+                else if (!vert->vert->id.isShape && !(vert->dir & ConnDirDown))
+                {
+                    generateEdge = false;
+                }
+                if (generateEdge)
+                {
+                    EdgeInf *edge = 
+                            new EdgeInf(last->vert, vert->vert, orthogonal);
+                    edge->setDist(vert->vert->point[dim] - 
+                            last->vert->point[dim]);
+                }
+
+                ++last;
+            }
+
+            ++vert;
+
+            if ((vert != breakPoints.end()) &&
+                    (last->vert->point[dim] == vert->vert->point[dim]))
+            {
+                // Still looking at same pair, just reset prev number pointer.
+                last = firstPrev;
+            }
+            else
+            {
+                // vert has moved to the beginning of a number number group.
+                // Last is now in the right place, so do nothing.
+            }
+        }
+    }
+
+    double begin;
+    double finish;
+    double pos;
+    bool shapeSide;
+    
+    VertSet vertInfs;
+    BreakpointSet breakPoints;
+private:
+	// MSVC wants to generate the assignment operator and the default 
+	// constructor, but fails.  Therefore we declare them private and 
+	// don't implement them.
+    LineSegment & operator=(LineSegment const &);
+    LineSegment();
+};
+
+typedef std::list<LineSegment> SegmentList;
+
+class SegmentListWrapper
+{
+    public:
+        LineSegment *insert(LineSegment segment)
+        {
+            SegmentList::iterator found = _list.end();
+            for (SegmentList::iterator curr = _list.begin();
+                    curr != _list.end(); ++curr)
+            {
+                if (curr->overlaps(segment))
+                {
+                    if (found != _list.end())
+                    {
+                        // This is not the first segment that overlaps,
+                        // so we need to merge and then delete an existing
+                        // segment.
+                        curr->mergeVertInfs(*found);
+                        _list.erase(found);
+                        found = curr;
+                    }
+                    else
+                    {
+                        // This is the first overlapping segment, so just 
+                        // merge the new segment with this one.
+                        curr->mergeVertInfs(segment);
+                        found = curr;
+                    }
+                }
+            }
+
+            if (found == _list.end())
+            {
+                // Add this line.
+                _list.push_back(segment);
+                return &(_list.back());
+            }
+
+            return &(*found);
+        }
+        SegmentList& list(void)
+        {
+            return _list;
+        }
+    private:
+        SegmentList _list;
+};
+
+
+// Given a router instance and a set of possible horizontal segments, and a
+// possible vertical visibility segment, compute and add edges to the
+// orthogonal visibility graph for all the visibility edges.
+static void intersectSegments(Router *router, SegmentList& segments, 
+        LineSegment& vertLine)
+{
+    COLA_ASSERT(vertLine.beginVertInf() == NULL);
+    COLA_ASSERT(vertLine.finishVertInf() == NULL);
+    for (SegmentList::iterator it = segments.begin(); it != segments.end(); )
+    {
+        LineSegment& horiLine = *it;
+
+        bool inVertSegRegion = ((vertLine.begin <= horiLine.pos) &&
+                                (vertLine.finish >= horiLine.pos));
+
+        if (horiLine.finish < vertLine.pos)
+        {
+            // Add horizontal visibility segment.
+            horiLine.addEdgeHorizontal(router);
+
+            size_t dim = XDIM; // x-dimension
+            horiLine.generateVisibilityEdgesFromBreakpointSet(router, dim);
+
+            // We've now swept past this horizontal segment, so delete.
+            it = segments.erase(it);
+            continue;
+        }
+        else if (horiLine.begin > vertLine.pos)
+        {
+            // We've yet to reach this segment in the sweep, so ignore.
+            ++it;
+            continue;
+        }
+        else if (horiLine.begin == vertLine.pos)
+        {
+            if (inVertSegRegion)
+            {
+                horiLine.insertBreakpointsBegin(router, vertLine);
+            }
+        }
+        else if (horiLine.finish == vertLine.pos)
+        {
+            if (inVertSegRegion)
+            {
+                // Add horizontal visibility segment.
+                horiLine.addEdgeHorizontal(router);
+            
+                horiLine.insertBreakpointsFinish(router, vertLine);
+                
+                size_t dim = XDIM; // x-dimension
+                horiLine.generateVisibilityEdgesFromBreakpointSet(router, dim);
+
+                // And we've now finished with the segment, so delete.
+                it = segments.erase(it);
+                continue;
+            }
+        }
+        else
+        {
+            COLA_ASSERT(horiLine.begin < vertLine.pos);
+            COLA_ASSERT(horiLine.finish > vertLine.pos);
+
+            if (inVertSegRegion)
+            {
+                // Add horizontal visibility segment.
+                VertSet intersectionVerts = 
+                        horiLine.addEdgeHorizontalTillIntersection(
+                            router, vertLine);
+
+                for (VertSet::iterator v = intersectionVerts.begin();
+                        v != intersectionVerts.end(); ++v)
+                {
+                    vertLine.breakPoints.insert(PosVertInf(horiLine.pos, *v,
+                            getPosVertInfDirection(*v, YDIM)));
+                }
+            }
+        }
+        ++it;
+    }
+
+    // Split breakPoints set into visibility segments.
+    size_t dimension = YDIM; // y-dimension
+    vertLine.generateVisibilityEdgesFromBreakpointSet(router, dimension);
+}
+
+
+// Processes an event for the vertical sweep used for computing the static 
+// orthogonal visibility graph.  This adds possible visibility sgments to 
+// the segments list.
+// The first pass is adding the event to the scanline, the second is for
+// processing the event and the third for removing it from the scanline.
+static void processEventVert(Router *router, NodeSet& scanline, 
+        SegmentListWrapper& segments, Event *e, unsigned int pass)
+{
+    Node *v = e->v;
+    
+    if ( ((pass == 1) && (e->type == Open)) ||
+         ((pass == 2) && (e->type == ConnPoint)) )
+    {
+        std::pair<NodeSet::iterator, bool> result = scanline.insert(v);
+        v->iter = result.first;
+        COLA_ASSERT(result.second);
+
+        NodeSet::iterator it = v->iter;
+        // Work out neighbours
+        if (it != scanline.begin()) 
+        {
+            Node *u = *(--it);
+            v->firstAbove = u;
+            u->firstBelow = v;
+        }
+        it = v->iter;
+        if (++it != scanline.end()) 
+        {
+            Node *u = *it;
+            v->firstBelow = u;
+            u->firstAbove = v;
+        }
+    }
+    
+    if (pass == 2)
+    {
+        if ((e->type == Open) || (e->type == Close))
+        {
+            // Shape edge positions.
+            double minShape = v->min[0];
+            double maxShape = v->max[0];
+            // As far as we can see.
+            double minLimit, maxLimit;
+            double minLimitMax, maxLimitMin;
+            v->findFirstPointAboveAndBelow(0, minLimit, maxLimit,
+                    minLimitMax, maxLimitMin);
+
+            // Only difference between Open and Close is whether the line
+            // segments are at the top or bottom of the shape.  Decide here.
+            double lineY = (e->type == Open) ? v->min[1] : v->max[1];
+
+            if (minLimitMax >= maxLimitMin)
+            {
+                // Insert possible visibility segments.
+                VertInf *vI1 = new VertInf(router, dummyOrthogID, 
+                            Point(minShape, lineY));
+                VertInf *vI2 = new VertInf(router, dummyOrthogID, 
+                            Point(maxShape, lineY));
+                
+                // There are no overlapping shapes, so give full visibility.
+                if (minLimit < minShape)
+                {
+                    segments.insert(LineSegment(minLimit, minShape, lineY,
+                                true, NULL, vI1));
+                }
+                segments.insert(LineSegment(minShape, maxShape, lineY, 
+                            true, vI1, vI2));
+                if (maxShape < maxLimit)
+                {
+                    segments.insert(LineSegment(maxShape, maxLimit, lineY,
+                                true, vI2, NULL));
+                }
+            }
+            else
+            {
+                if ((minLimitMax > minLimit) && (minLimitMax >= minShape))
+                {
+                    segments.insert(LineSegment(minLimit, minLimitMax, lineY,
+                                true, NULL, NULL));
+                }
+                if ((maxLimitMin < maxLimit) && (maxLimitMin <= maxShape))
+                {
+                    segments.insert(LineSegment(maxLimitMin, maxLimit, lineY,
+                                true, NULL, NULL));
+                }
+            }
+        }
+        else if (e->type == ConnPoint)
+        {
+            // Connection point.
+            VertInf *centreVert = e->v->c;
+            Point& cp = centreVert->point;
+
+            // As far as we can see.
+            double minLimit = v->firstPointAbove(0);
+            double maxLimit = v->firstPointBelow(0);
+            bool inShape = v->isInsideShape(0);
+
+            LineSegment *line1 = NULL, *line2 = NULL;
+            if (!inShape || (centreVert->visDirections & ConnDirLeft))
+            {
+                line1 = segments.insert(LineSegment(minLimit, cp.x, e->pos, 
+                        true, NULL, centreVert));
+            }
+            if (!inShape || (centreVert->visDirections & ConnDirRight))
+            {
+                line2 = segments.insert(LineSegment(cp.x, maxLimit, e->pos, 
+                        true, centreVert, NULL));
+            }
+            if (!line1 && !line2)
+            {
+                // Add a point segment for the centre point.
+                segments.insert(LineSegment(cp.x, e->pos, centreVert));
+            }
+            
+            if (!inShape)
+            {
+                // This is not contained within a shape so add a normal
+                // visibility graph point here too (since paths won't route
+                // *through* connector endpoint vertices).
+                if (line1 || line2)
+                {
+                    VertInf *cent = new VertInf(router, dummyOrthogID, cp);
+                    if (line1)
+                    {
+                        line1->vertInfs.insert(cent);
+                    }
+                    if (line2)
+                    {
+                        line2->vertInfs.insert(cent);
+                    }
+                }
+            }
+        }
+    }
+    
+    if ( ((pass == 3) && (e->type == Close)) ||
+         ((pass == 2) && (e->type == ConnPoint)) )
+    {
+        // Clean up neighbour pointers.
+        Node *l = v->firstAbove, *r = v->firstBelow;
+        if (l != NULL) 
+        {
+            l->firstBelow = v->firstBelow;
+        }
+        if (r != NULL)
+        {
+            r->firstAbove = v->firstAbove;
+        }
+
+        if (e->type == ConnPoint)
+        {
+            scanline.erase(v->iter);
+            delete v;
+        }
+        else  // if (e->type == Close)
+        {
+            size_t result;
+            result = scanline.erase(v);
+            COLA_ASSERT(result == 1);
+            delete v;
+        }
+    }
+}
+
+
+// Processes an event for the vertical sweep used for computing the static 
+// orthogonal visibility graph.  This adds possible visibility sgments to 
+// the segments list.
+// The first pass is adding the event to the scanline, the second is for
+// processing the event and the third for removing it from the scanline.
+static void processEventHori(Router *router, NodeSet& scanline, 
+        SegmentListWrapper& segments, Event *e, unsigned int pass)
+{
+    Node *v = e->v;
+    
+    if ( ((pass == 1) && (e->type == Open)) ||
+         ((pass == 2) && (e->type == ConnPoint)) )
+    {
+        std::pair<NodeSet::iterator, bool> result = scanline.insert(v);
+        v->iter = result.first;
+        COLA_ASSERT(result.second);
+
+        NodeSet::iterator it = v->iter;
+        // Work out neighbours
+        if (it != scanline.begin()) 
+        {
+            Node *u = *(--it);
+            v->firstAbove = u;
+            u->firstBelow = v;
+        }
+        it = v->iter;
+        if (++it != scanline.end()) 
+        {
+            Node *u = *it;
+            v->firstBelow = u;
+            u->firstAbove = v;
+        }
+    }
+    
+    if (pass == 2)
+    {
+        if ((e->type == Open) || (e->type == Close))
+        {
+            // Shape edge positions.
+            double minShape = v->min[1];
+            double maxShape = v->max[1];
+            // As far as we can see.
+            double minLimit, maxLimit;
+            double minLimitMax, maxLimitMin;
+            v->findFirstPointAboveAndBelow(1, minLimit, maxLimit,
+                    minLimitMax, maxLimitMin);
+
+            // Only difference between Open and Close is whether the line
+            // segments are at the left or right of the shape.  Decide here.
+            double lineX = (e->type == Open) ? v->min[0] : v->max[0];
+
+            if (minLimitMax >= maxLimitMin)
+            {
+                LineSegment vertSeg = LineSegment(minLimit, maxLimit, lineX);
+                segments.insert(vertSeg);
+            }
+            else
+            {
+                if ((minLimitMax > minLimit) && (minLimitMax >= minShape))
+                {
+                    LineSegment vertSeg = 
+                            LineSegment(minLimit, minLimitMax, lineX);
+                    segments.insert(vertSeg);
+                }
+                if ((maxLimitMin < maxLimit) && (maxLimitMin <= maxShape))
+                {
+                    LineSegment vertSeg = 
+                            LineSegment(maxLimitMin, maxLimit, lineX);
+                    segments.insert(vertSeg);
+                }
+            }
+        }
+        else if (e->type == ConnPoint)
+        {
+            // Connection point.
+            VertInf *centreVert = e->v->c;
+            Point& cp = centreVert->point;
+
+            // As far as we can see.
+            double minLimit = v->firstPointAbove(1);
+            double maxLimit = v->firstPointBelow(1);
+            bool inShape = v->isInsideShape(1);
+            
+            if (!inShape || (centreVert->visDirections & ConnDirUp))
+            {
+                segments.insert(LineSegment(minLimit, cp.y, e->pos));
+            }
+            if (!inShape || (centreVert->visDirections & ConnDirDown))
+            {
+                segments.insert(LineSegment(cp.y, maxLimit, e->pos));
+            }
+        }
+    }
+    
+    if ( ((pass == 3) && (e->type == Close)) ||
+         ((pass == 2) && (e->type == ConnPoint)) )
+    {
+        // Clean up neighbour pointers.
+        Node *l = v->firstAbove, *r = v->firstBelow;
+        if (l != NULL) 
+        {
+            l->firstBelow = v->firstBelow;
+        }
+        if (r != NULL)
+        {
+            r->firstAbove = v->firstAbove;
+        }
+
+        if (e->type == ConnPoint)
+        {
+            scanline.erase(v->iter);
+            delete v;
+        }
+        else  // if (e->type == Close)
+        {
+            size_t result;
+            result = scanline.erase(v);
+            COLA_ASSERT(result == 1);
+            delete v;
+        }
+    }
+}
+
+
+extern void generateStaticOrthogonalVisGraph(Router *router)
+{
+    const size_t n = router->shapeRefs.size();
+    const unsigned cpn = router->vertices.connsSize();
+    // Set up the events for the vertical sweep.
+    size_t totalEvents = (2 * n) + cpn;
+    events = new Event*[totalEvents];
+    unsigned ctr = 0;
+    ShapeRefList::iterator shRefIt = router->shapeRefs.begin();
+    for (unsigned i = 0; i < n; i++)
+    {
+        ShapeRef *shRef = *shRefIt;
+        double minX, minY, maxX, maxY;
+        shRef->polygon().getBoundingRect(&minX, &minY, &maxX, &maxY);
+        double midX = minX + ((maxX - minX) / 2);
+        Node *v = new Node(shRef, midX);
+        events[ctr++] = new Event(Open, v, minY);
+        events[ctr++] = new Event(Close, v, maxY);
+
+        ++shRefIt;
+    }
+    for (VertInf *curr = router->vertices.connsBegin(); 
+            curr && (curr != router->vertices.shapesBegin()); 
+            curr = curr->lstNext)
+    {
+        Point& point = curr->point;
+
+        Node *v = new Node(curr, point.x);
+        events[ctr++] = new Event(ConnPoint, v, point.y);
+    }
+    qsort((Event*)events, (size_t) totalEvents, sizeof(Event*), compare_events);
+
+    // Process the vertical sweep.
+    // We do multiple passes over sections of the list so we can add relevant
+    // entries to the scanline that might follow, before process them.
+    SegmentListWrapper segments;
+    NodeSet scanline;
+    double thisPos = (totalEvents > 0) ? events[0]->pos : 0;
+    unsigned int posStartIndex = 0;
+    unsigned int posFinishIndex = 0;
+    for (unsigned i = 0; i <= totalEvents; ++i)
+    {
+        // If we have finished the current scanline or all events, then we
+        // process the events on the current scanline in a couple of passes.
+        if ((i == totalEvents) || (events[i]->pos != thisPos))
+        {
+            posFinishIndex = i;
+            for (int pass = 2; pass <= 3; ++pass)
+            {
+                for (unsigned j = posStartIndex; j < posFinishIndex; ++j)
+                {
+                    processEventVert(router, scanline, segments, 
+                            events[j], pass);
+                }
+            }
+
+            if (i == totalEvents)
+            {
+                // We have cleaned up, so we can now break out of loop.
+                break;
+            }
+
+            thisPos = events[i]->pos;
+            posStartIndex = i;
+        }
+
+        // Do the first sweep event handling -- building the correct 
+        // structure of the scanline.
+        const int pass = 1;
+        processEventVert(router, scanline, segments, events[i], pass);
+    }
+    COLA_ASSERT(scanline.size() == 0);
+    for (unsigned i = 0; i < totalEvents; ++i)
+    {
+        delete events[i];
+    }
+
+    segments.list().sort();
+
+    // Set up the events for the horizontal sweep.
+    SegmentListWrapper vertSegments;
+    ctr = 0;
+    shRefIt = router->shapeRefs.begin();
+    for (unsigned i = 0; i < n; i++)
+    {
+        ShapeRef *shRef = *shRefIt;
+        double minX, minY, maxX, maxY;
+        shRef->polygon().getBoundingRect(&minX, &minY, &maxX, &maxY);
+        double midY = minY + ((maxY - minY) / 2);
+        Node *v = new Node(shRef, midY);
+        events[ctr++] = new Event(Open, v, minX);
+        events[ctr++] = new Event(Close, v, maxX);
+
+        ++shRefIt;
+    }
+    for (VertInf *curr = router->vertices.connsBegin(); 
+            curr && (curr != router->vertices.shapesBegin()); 
+            curr = curr->lstNext)
+    {
+        Point& point = curr->point;
+
+        Node *v = new Node(curr, point.y);
+        events[ctr++] = new Event(ConnPoint, v, point.x);
+    }
+    qsort((Event*)events, (size_t) totalEvents, sizeof(Event*), compare_events);
+
+    // Process the horizontal sweep
+    thisPos = (totalEvents > 0) ? events[0]->pos : 0;
+    posStartIndex = 0;
+    posFinishIndex = 0;
+    for (unsigned i = 0; i <= totalEvents; ++i)
+    {
+        // If we have finished the current scanline or all events, then we
+        // process the events on the current scanline in a couple of passes.
+        if ((i == totalEvents) || (events[i]->pos != thisPos))
+        {
+            posFinishIndex = i;
+            for (int pass = 2; pass <= 3; ++pass)
+            {
+                for (unsigned j = posStartIndex; j < posFinishIndex; ++j)
+                {
+                    processEventHori(router, scanline, vertSegments, 
+                            events[j], pass);
+                }
+            }
+            
+            // Process the merged line segments.
+            vertSegments.list().sort();
+            for (SegmentList::iterator curr = vertSegments.list().begin();
+                    curr != vertSegments.list().end(); ++curr)
+            {
+                intersectSegments(router, segments.list(), *curr);
+            }
+            vertSegments.list().clear();
+
+            if (i == totalEvents)
+            {
+                // We have cleaned up, so we can now break out of loop.
+                break;
+            }
+
+            thisPos = events[i]->pos;
+            posStartIndex = i;
+        }
+
+        // Do the first sweep event handling -- building the correct 
+        // structure of the scanline.
+        const int pass = 1;
+        processEventHori(router, scanline, vertSegments, events[i], pass);
+    }
+    COLA_ASSERT(scanline.size() == 0);
+    for (unsigned i = 0; i < totalEvents; ++i)
+    {
+        delete events[i];
+    }
+    delete [] events;
+
+    // Add portions of the horizontal line that are after the final vertical
+    // position we considered.
+    for (SegmentList::iterator it = segments.list().begin(); 
+            it != segments.list().end(); )
+    {
+        LineSegment& horiLine = *it;
+
+        horiLine.addEdgeHorizontal(router);
+        
+        size_t dim = XDIM; // x-dimension
+        horiLine.generateVisibilityEdgesFromBreakpointSet(router, dim);
+
+        it = segments.list().erase(it);
+    }
+}
+
+
+//============================================================================
+//                           Path Adjustment code
+//============================================================================
+
+
+
+
+// Processes sweep events used to determine each horizontal and vertical 
+// line segment in a connector's channel of visibility.  
+// Four calls to this function are made at each position by the scanline:
+//   1) Handle all Close event processing.
+//   2) Remove Close event objects from the scanline.
+//   3) Add Open event objects to the scanline.
+//   4) Handle all Open event processing.
+//
+static void processShiftEvent(Router *router, NodeSet& scanline, 
+        ShiftSegmentList& segments, Event *e, size_t dim,
+        unsigned int pass)
+{
+    Node *v = e->v;
+    
+    if ( ((pass == 3) && (e->type == Open)) ||
+         ((pass == 3) && (e->type == SegOpen)) )
+    {
+        std::pair<NodeSet::iterator, bool> result = scanline.insert(v);
+        v->iter = result.first;
+        COLA_ASSERT(result.second);
+
+        NodeSet::iterator it = v->iter;
+        // Work out neighbours
+        if (it != scanline.begin()) 
+        {
+            Node *u = *(--it);
+            v->firstAbove = u;
+            u->firstBelow = v;
+        }
+        it = v->iter;
+        if (++it != scanline.end()) 
+        {
+            Node *u = *it;
+            v->firstBelow = u;
+            u->firstAbove = v;
+        }
+    }
+    
+    if ( ((pass == 4) && (e->type == Open)) ||
+         ((pass == 4) && (e->type == SegOpen)) ||
+         ((pass == 1) && (e->type == SegClose)) ||
+         ((pass == 1) && (e->type == Close)) )
+    {
+        if (v->ss)
+        {
+            // As far as we can see.
+            double minLimit = v->firstObstacleAbove(dim);
+            double maxLimit = v->firstObstacleBelow(dim);
+
+            v->ss->minSpaceLimit = 
+                    std::max(minLimit, v->ss->minSpaceLimit);
+            v->ss->maxSpaceLimit = 
+                    std::min(maxLimit, v->ss->maxSpaceLimit);
+        }
+        else
+        {
+            v->markShiftSegmentsAbove(dim);
+            v->markShiftSegmentsBelow(dim);
+        }
+    }
+    
+    if ( ((pass == 2) && (e->type == SegClose)) ||
+         ((pass == 2) && (e->type == Close)) )
+    {
+        // Clean up neighbour pointers.
+        Node *l = v->firstAbove, *r = v->firstBelow;
+        if (l != NULL) 
+        {
+            l->firstBelow = v->firstBelow;
+        }
+        if (r != NULL)
+        {
+            r->firstAbove = v->firstAbove;
+        }
+
+        size_t result;
+        result = scanline.erase(v);
+        COLA_ASSERT(result == 1);
+        delete v;
+    }
+}
+
+
+static void buildOrthogonalChannelInfo(Router *router, 
+        const size_t dim, ShiftSegmentList& segmentList)
+{
+    if (router->routingPenalty(segmentPenalty) == 0)
+    {
+        // This code assumes the routes are pretty optimal, so we don't
+        // do this adjustment if the routes have no segment penalty.
+        return;
+    }
+
+    size_t altDim = (dim + 1) % 2;
+    // For each connector.
+    for (ConnRefList::const_iterator curr = router->connRefs.begin(); 
+            curr != router->connRefs.end(); ++curr) 
+    {
+        if ((*curr)->routingType() != ConnType_Orthogonal)
+        {
+            continue;
+        }
+        Polygon& displayRoute = (*curr)->displayRoute();
+        // Determine all line segments that we are interested in shifting. 
+        // We don't consider the first or last segment of a path.
+        for (size_t i = 1; i < displayRoute.size(); ++i)
+        {
+            if (displayRoute.ps[i - 1][dim] == displayRoute.ps[i][dim])
+            {
+                // It's a segment in the dimension we are processing,
+                size_t indexLow = i - 1;
+                size_t indexHigh = i;
+                if (displayRoute.ps[i - 1][altDim] > displayRoute.ps[i][altDim])
+                {
+                    indexLow = i;
+                    indexHigh = i - 1;
+                }
+                COLA_ASSERT(displayRoute.at(indexLow)[altDim] < 
+                        displayRoute.at(indexHigh)[altDim]);
+
+                if ((i == 1) || ((i + 1) == displayRoute.size()))
+                {
+                    // The first and last segment of a connector can't be 
+                    // shifted.  We call them fixed segments.  Note: this
+                    // will change if we later allow connection channels.
+                    segmentList.push_back(
+                            ShiftSegment(*curr, indexLow, indexHigh, dim));
+                    continue;
+                }
+
+                // The segment probably has space to be shifted.
+                double minLim = -CHANNEL_MAX;
+                double maxLim = CHANNEL_MAX;
+                bool isSBend = false;
+
+                double prevPos = displayRoute.ps[i - 2][dim];
+                double nextPos = displayRoute.ps[i + 1][dim];
+                if (((prevPos < displayRoute.ps[i][dim]) &&
+                     (nextPos > displayRoute.ps[i][dim])) 
+                     ||
+                    ((prevPos > displayRoute.ps[i][dim]) &&
+                     (nextPos < displayRoute.ps[i][dim])) )
+                {
+                    isSBend = true;
+
+                    // Determine limits if the s-bend is not due to an 
+                    // obstacle.  In this case we need to limit the channel 
+                    // to the span of the adjoining segments to this one.
+                    if ((prevPos < displayRoute.ps[i][dim]) &&
+                        (nextPos > displayRoute.ps[i][dim]))
+                    {
+                        minLim = std::max(minLim, prevPos);
+                        maxLim = std::min(maxLim, nextPos);
+                    }
+                    else
+                    {
+                        minLim = std::max(minLim, nextPos);
+                        maxLim = std::min(maxLim, prevPos);
+                    }
+                }
+                else
+                {
+                    // isCBend: Both adjoining segments are in the same
+                    // direction.  We indicate this for later by setting 
+                    // the maxLim or minLim to the segment position.
+                    if (prevPos < displayRoute.ps[i][dim])
+                    {
+                        minLim = displayRoute.ps[i][dim];
+                    }
+                    else
+                    {
+                        maxLim = displayRoute.ps[i][dim];
+                    }
+                }
+
+                segmentList.push_back(ShiftSegment(*curr, indexLow, 
+                            indexHigh, isSBend, dim, minLim, maxLim));
+            }
+        }
+    }
+    if (segmentList.empty())
+    {
+        // There are no segments, so we can just return now.
+        return;
+    }
+    
+    // Do a sweep and shift these segments.
+    const size_t n = router->shapeRefs.size();
+    const size_t cpn = segmentList.size();
+    // Set up the events for the sweep.
+    size_t totalEvents = 2 * (n + cpn);
+    events = new Event*[totalEvents];
+    unsigned ctr = 0;
+    ShapeRefList::iterator shRefIt = router->shapeRefs.begin();
+    for (unsigned i = 0; i < n; i++)
+    {
+        ShapeRef *shRef = *shRefIt;
+        Point min, max;
+        shRef->polygon().getBoundingRect(&min.x, &min.y, &max.x, &max.y);
+        double mid = min[dim] + ((max[dim] - min[dim]) / 2);
+        Node *v = new Node(shRef, mid);
+        events[ctr++] = new Event(Open, v, min[altDim]);
+        events[ctr++] = new Event(Close, v, max[altDim]);
+
+        ++shRefIt;
+    }
+    for (ShiftSegmentList::iterator curr = segmentList.begin(); 
+            curr != segmentList.end(); ++curr)
+    {
+        const Point& lowPt = curr->lowPoint();
+        const Point& highPt = curr->highPoint();
+
+        COLA_ASSERT(lowPt[dim] == highPt[dim]);
+        COLA_ASSERT(lowPt[altDim] < highPt[altDim]);
+        Node *v = new Node(&(*curr), lowPt[dim]);
+        events[ctr++] = new Event(SegOpen, v, lowPt[altDim]);
+        events[ctr++] = new Event(SegClose, v, highPt[altDim]);
+    }
+    qsort((Event*)events, (size_t) totalEvents, sizeof(Event*), compare_events);
+
+    // Process the sweep.
+    // We do multiple passes over sections of the list so we can add relevant
+    // entries to the scanline that might follow, before process them.
+    NodeSet scanline;
+    double thisPos = (totalEvents > 0) ? events[0]->pos : 0;
+    unsigned int posStartIndex = 0;
+    unsigned int posFinishIndex = 0;
+    for (unsigned i = 0; i <= totalEvents; ++i)
+    {
+        // If we have finished the current scanline or all events, then we
+        // process the events on the current scanline in a couple of passes.
+        if ((i == totalEvents) || (events[i]->pos != thisPos))
+        {
+            posFinishIndex = i;
+            for (int pass = 2; pass <= 4; ++pass)
+            {
+                for (unsigned j = posStartIndex; j < posFinishIndex; ++j)
+                {
+                    processShiftEvent(router, scanline, segmentList, events[j], 
+                            dim, pass);
+                }
+            }
+
+            if (i == totalEvents)
+            {
+                // We have cleaned up, so we can now break out of loop.
+                break;
+            }
+
+            thisPos = events[i]->pos;
+            posStartIndex = i;
+        }
+
+        // Do the first sweep event handling -- building the correct 
+        // structure of the scanline.
+        const int pass = 1;
+        processShiftEvent(router, scanline, segmentList, events[i],
+                dim, pass);
+    }
+    COLA_ASSERT(scanline.size() == 0);
+    for (unsigned i = 0; i < totalEvents; ++i)
+    {
+        delete events[i];
+    }
+    delete [] events;
+}
+
+
+static void simplifyOrthogonalRoutes(Router *router)
+{
+    // Simplify routes.
+    for (ConnRefList::const_iterator curr = router->connRefs.begin(); 
+            curr != router->connRefs.end(); ++curr) 
+    {
+        if ((*curr)->routingType() != ConnType_Orthogonal)
+        {
+            continue;
+        }
+        (*curr)->set_route((*curr)->displayRoute().simplify());
+    }
+}
+
+
+static void buildOrthogonalNudgingOrderInfo(Router *router, 
+        PtOrderMap& pointOrders)
+{
+    // Simplify routes.
+    simplifyOrthogonalRoutes(router);
+
+    int crossingsN = 0;
+
+    // Do segment splitting.
+    for (ConnRefList::const_iterator curr = router->connRefs.begin(); 
+            curr != router->connRefs.end(); ++curr) 
+    {
+        if ((*curr)->routingType() != ConnType_Orthogonal)
+        {
+            continue;
+        }
+        ConnRef *conn = *curr;
+        
+        for (ConnRefList::const_iterator curr2 = router->connRefs.begin(); 
+                curr2 != router->connRefs.end(); ++curr2) 
+        {
+            if ((*curr2)->routingType() != ConnType_Orthogonal)
+            {
+                continue;
+            }
+            ConnRef *conn2 = *curr2;
+            
+            if (conn == conn2)
+            {
+                continue;
+            }
+            
+            Avoid::Polygon& route = conn->displayRoute();
+            Avoid::Polygon& route2 = conn2->displayRoute();
+            splitBranchingSegments(route2, true, route);
+        }
+    }
+
+    for (ConnRefList::const_iterator curr = router->connRefs.begin(); 
+            curr != router->connRefs.end(); ++curr) 
+    {
+        if ((*curr)->routingType() != ConnType_Orthogonal)
+        {
+            continue;
+        }
+        ConnRef *conn = *curr;
+        
+        for (ConnRefList::const_iterator curr2 = curr; 
+                curr2 != router->connRefs.end(); ++curr2) 
+        {
+            if ((*curr2)->routingType() != ConnType_Orthogonal)
+            {
+                continue;
+            }
+            ConnRef *conn2 = *curr2;
+
+            if (conn == conn2)
+            {
+                continue;
+            }
+            
+            Avoid::Polygon& route = conn->displayRoute();
+            Avoid::Polygon& route2 = conn2->displayRoute();
+            bool checkForBranchingSegments = false;
+            int crossings = 0;
+            for (size_t i = 1; i < route.size(); ++i)
+            {
+                const bool finalSegment = ((i + 1) == route.size());
+                crossings += countRealCrossings(route2, true, route, i, 
+                        checkForBranchingSegments, finalSegment, NULL, 
+                        &pointOrders, conn2, conn).first;
+            }
+            if (crossings > 0)
+            {
+                crossingsN += crossings;
+            }
+        }
+    }
+    
+    // Sort the point orders.
+    PtOrderMap::iterator finish = pointOrders.end();
+    for (PtOrderMap::iterator it = pointOrders.begin(); it != finish; ++it)
+    {
+        //const VertID& ptID = it->first;
+        PtOrder& order = it->second;
+
+        for (size_t dim = XDIM; dim <= YDIM; ++dim)
+        {
+            order.sort(dim);
+        }
+    }
+}
+
+
+class CmpLineOrder 
+{
+    public:
+        CmpLineOrder(PtOrderMap& ord, const size_t dim)
+            : orders(ord),
+              dimension(dim)
+        {
+        }
+        bool operator()(const ShiftSegment& lhs, const ShiftSegment& rhs,
+                bool *comparable = NULL) const
+        {
+            if (comparable)
+            {
+                *comparable = true;
+            }
+            Point lhsLow  = lhs.lowPoint(); 
+            Point rhsLow  = rhs.lowPoint(); 
+#ifndef NDEBUG
+            const Point& lhsHigh = lhs.highPoint(); 
+            const Point& rhsHigh = rhs.highPoint(); 
+#endif
+            size_t altDim = (dimension + 1) % 2;
+
+            COLA_ASSERT(lhsLow[dimension] == lhsHigh[dimension]);
+            COLA_ASSERT(rhsLow[dimension] == rhsHigh[dimension]);
+
+            if (lhsLow[dimension] != rhsLow[dimension])
+            {
+                return lhsLow[dimension] < rhsLow[dimension];
+            }
+            
+            // If one of these is fixed, then determine order based on 
+            // fixed segment, that is, order so the fixed segment doesn't 
+            // block movement.
+            bool oneIsFixed = false;
+            const int lhsFixedOrder = lhs.fixedOrder(oneIsFixed);
+            const int rhsFixedOrder = rhs.fixedOrder(oneIsFixed);
+            if (oneIsFixed && (lhsFixedOrder != rhsFixedOrder))
+            {
+                return lhsFixedOrder < rhsFixedOrder;
+            }
+
+            // C-bends that did not have a clear order with s-bends might 
+            // not have a good ordering here, so compare their order in 
+            // terms of C-bend direction and S-bends and use that if it
+            // differs for the two segments.
+            const int lhsOrder = lhs.order();
+            const int rhsOrder = rhs.order();
+            if (lhsOrder != rhsOrder)
+            {
+                return lhsOrder < rhsOrder;
+            }
+
+            // Need to index using the original point into the map, so find it.
+            Point& unchanged = (lhsLow[altDim] > rhsLow[altDim]) ?
+                    lhsLow : rhsLow;
+
+            PtOrder& lowOrder = orders[unchanged];
+            int lhsPos = lowOrder.positionFor(lhs.connRef, dimension);
+            int rhsPos = lowOrder.positionFor(rhs.connRef, dimension);
+            if ((lhsPos == -1) || (rhsPos == -1))
+            {
+                // A value for rhsPos or lhsPos mean the points are not directly
+                // comparable, meaning they are at the same position but cannot
+                // overlap (they are just collinear.  The relative order for 
+                // these segments is not important since we do not constrain
+                // them against each other.
+                //COLA_ASSERT(lhs.overlapsWith(rhs, dimension) == false);
+                // We do need to be consistent though.
+                if (comparable)
+                {
+                    *comparable = false;
+                }
+                return lhsLow[altDim] < rhsLow[altDim];
+            }
+
+            return lhsPos < rhsPos;
+        }
+
+        PtOrderMap& orders;
+        const size_t dimension;
+};
+
+
+// We can use the normaal sort algorithm for lists since it is not possible 
+// to comapre all elements, but there will be an ordering defined between 
+// most of the elements.  Hence we order these, using insertion sort, and 
+// the case of them not being able to be compared is handled by not setting 
+// up any constraints between such segments when doing the nudging.
+//
+static ShiftSegmentList linesort(ShiftSegmentList origList, 
+        CmpLineOrder& comparison)
+{
+    ShiftSegmentList resultList;
+
+    while (!origList.empty())
+    {
+        // Get and remove the first element from the origList.
+        ShiftSegment segment = origList.front();
+        origList.pop_front();
+
+        // Find the insertion point in the resultList.
+        ShiftSegmentList::iterator curr;
+        for (curr = resultList.begin(); curr != resultList.end(); ++curr)
+        {
+            bool comparable = false;
+            bool lessThan = comparison(segment, *curr, &comparable);
+
+            if (comparable && lessThan)
+            {
+                // If it is comparable and lessThan, then we have found the
+                // insertion point.
+                break;
+            }
+        }
+
+        // Insert the element into the reultList at the required point.
+        resultList.insert(curr, segment);
+    }
+
+    return resultList;
+}
+
+
+typedef std::list<ShiftSegment *> ShiftSegmentPtrList;
+
+
+static void nudgeOrthogonalRoutes(Router *router, size_t dimension, 
+        PtOrderMap& pointOrders, ShiftSegmentList& segmentList)
+{
+    // Do the actual nudging.
+    ShiftSegmentList currentRegion;
+    while (!segmentList.empty())
+    {
+        // Take a reference segment
+        ShiftSegment& currentSegment = segmentList.front();
+        // Then, find the segments that overlap this one.
+        currentRegion.clear();
+        currentRegion.push_back(currentSegment);
+        segmentList.erase(segmentList.begin());
+        for (ShiftSegmentList::iterator curr = segmentList.begin();
+                curr != segmentList.end(); )
+        {
+            bool overlaps = false;
+            for (ShiftSegmentList::iterator curr2 = currentRegion.begin();
+                    curr2 != currentRegion.end(); ++curr2)
+            {
+                if (curr->overlapsWith(*curr2, dimension))
+                {
+                    overlaps = true;
+                    break;
+                }
+            }
+            if (overlaps)
+            {
+                currentRegion.push_back(*curr);
+                segmentList.erase(curr);
+                // Consider segments from the beginning, since we mave have
+                // since passed segments that overlap with the new set.
+                curr = segmentList.begin();
+            }
+            else
+            {
+                ++curr;
+            }
+        }
+        CmpLineOrder lineSortComp(pointOrders, dimension);
+        currentRegion = linesort(currentRegion, lineSortComp);
+        
+        if (currentRegion.size() == 1)
+        {
+            // Save creating the solver instance if there is just one
+            // immovable segment.
+            if (!currentRegion.front().sBend)
+            {
+                continue;
+            }
+        }
+
+        // Process these segments.
+        Variables vs;
+        Constraints cs;
+        ShiftSegmentPtrList prevVars;
+        // IDs:
+        const int freeID    = 0;
+        const int fixedID   = 1;
+        // Weights:
+        double freeWeight   = 0.00001;
+        double strongWeight = 0.001;
+        double fixedWeight  = 100000;
+        //printf("-------------------------------------------------------\n");
+        //printf("Nudge -- size: %d\n", (int) currentRegion.size());
+        for (ShiftSegmentList::iterator currSegment = currentRegion.begin();
+                currSegment != currentRegion.end(); ++currSegment)
+        {
+            Point& lowPt = currSegment->lowPoint();
+            
+            // Create a solver variable for the position of this segment.
+            int varID = freeID;
+            double idealPos = lowPt[dimension];
+            double weight = freeWeight;
+            if (currSegment->sBend)
+            {
+                COLA_ASSERT(currSegment->minSpaceLimit > -CHANNEL_MAX);
+                COLA_ASSERT(currSegment->maxSpaceLimit < CHANNEL_MAX);
+                
+                // For s-bends, take the middle as ideal.
+                idealPos = currSegment->minSpaceLimit +
+                        ((currSegment->maxSpaceLimit -
+                          currSegment->minSpaceLimit) / 2);
+            }
+            else if (currSegment->fixed)
+            {
+                // Fixed segments shouldn't get moved.
+                weight = fixedWeight;
+                varID = fixedID;
+            }
+            else
+            {
+                // Set a higher weight for c-bends to stop them sometimes 
+                // getting pushed out into channels by more-free connectors
+                // to the "inner" side of them.
+                weight = strongWeight;
+            }
+            currSegment->variable = new Variable(varID, idealPos, weight);
+            vs.push_back(currSegment->variable);
+            size_t index = vs.size() - 1;
+            //printf("line  %.15f  pos: %g   min: %g  max: %g\n",
+            //        lowPt[dimension], idealPos, currSegment->minSpaceLimit,
+            //        currSegment->maxSpaceLimit);
+
+            // Constrain position in relation to previously seen segments,
+            // if necessary (i.e. when they could overlap).
+            for (ShiftSegmentPtrList::iterator prevVarIt = prevVars.begin();
+                    prevVarIt != prevVars.end(); )
+            {
+                ShiftSegment *prevSeg = *prevVarIt;
+                Variable *prevVar = prevSeg->variable;
+
+                if (currSegment->overlapsWith(*prevSeg, dimension) &&
+                        (!(currSegment->fixed) || !(prevSeg->fixed)))
+                {
+                    // If there is a previous segment to the left that 
+                    // could overlap this in the shift direction, then 
+                    // constrain the two segments to be separated.
+                    // Though don't add the constraint if both the 
+                    // segments are fixed in place.
+                    cs.push_back(new Constraint(prevVar, vs[index],
+                            router->orthogonalNudgeDistance()));
+                    prevVarIt = prevVars.erase(prevVarIt);
+                }
+                else
+                {
+                    ++prevVarIt;
+                }
+            }
+
+            if (!currSegment->fixed)
+            {
+                // If this segment sees a channel boundary to its left, 
+                // then constrain its placement as such.
+                if (currSegment->minSpaceLimit > -CHANNEL_MAX)
+                {
+                    vs.push_back(new Variable(fixedID, 
+                                currSegment->minSpaceLimit, fixedWeight));
+                    cs.push_back(new Constraint(vs[vs.size() - 1], vs[index], 
+                                0.0));
+                }
+                
+                // If this segment sees a channel boundary to its right, 
+                // then constrain its placement as such.
+                if (currSegment->maxSpaceLimit < CHANNEL_MAX)
+                {
+                    vs.push_back(new Variable(fixedID, 
+                                currSegment->maxSpaceLimit, fixedWeight));
+                    cs.push_back(new Constraint(vs[index], vs[vs.size() - 1],
+                                0.0));
+                }
+            }
+            prevVars.push_back(&(*currSegment));
+        }
+#if 0
+        for(unsigned i=0;i<vs.size();i++) {
+            printf("-vs[%d]=%f\n",i,vs[i]->desiredPosition);
+        }
+#endif
+        IncSolver f(vs,cs);
+        f.solve();
+        bool satisfied = true;
+        for (size_t i = 0; i < vs.size(); ++i) 
+        {
+            if (vs[i]->id == fixedID)
+            {
+                if (fabs(vs[i]->finalPosition - vs[i]->desiredPosition) > 0.01)
+                {
+                    satisfied = false;
+                    break;
+                }
+            }
+        }
+        if (satisfied)
+        {
+            for (ShiftSegmentList::iterator currSegment = currentRegion.begin();
+                    currSegment != currentRegion.end(); ++currSegment)
+            {
+                Point& lowPt = currSegment->lowPoint();
+                Point& highPt = currSegment->highPoint();
+                double newPos = currSegment->variable->finalPosition;
+                //printf("Pos: %X, %g\n", (int) currSegment->connRef, newPos);
+                lowPt[dimension] = newPos;
+                highPt[dimension] = newPos;
+            }
+        }
+#if 0
+        for(unsigned i=0;i<vs.size();i++) {
+            printf("+vs[%d]=%f\n",i,vs[i]->finalPosition);
+        }
+#endif
+        for_each(vs.begin(),vs.end(),delete_object());
+        for_each(cs.begin(),cs.end(),delete_object());
+    }
+}
+
+
+extern void improveOrthogonalRoutes(Router *router)
+{
+    router->timers.Register(tmOrthogNudge, timerStart);
+    for (size_t dimension = 0; dimension < 2; ++dimension)
+    {
+        // Build nudging info.
+        // XXX: We need to build the point orders separately in each
+        //      dimension since things move.  There is probably a more 
+        //      efficient way to do this.
+        PtOrderMap pointOrders;
+        buildOrthogonalNudgingOrderInfo(router, pointOrders);
+
+        // Simplify routes.
+        simplifyOrthogonalRoutes(router);
+
+        // Do the centring and nudging.
+        ShiftSegmentList segLists;
+        buildOrthogonalChannelInfo(router, dimension, segLists);
+        nudgeOrthogonalRoutes(router, dimension, pointOrders, segLists);
+    }
+    router->timers.Stop();
+}
+
+
+}
diff --git a/src/libavoid/orthogonal.h b/src/libavoid/orthogonal.h
new file mode 100644
index 000000000..4fb974099
--- /dev/null
+++ b/src/libavoid/orthogonal.h
@@ -0,0 +1,39 @@
+/*
+ * vim: ts=4 sw=4 et tw=0 wm=0
+ *
+ * libavoid - Fast, Incremental, Object-avoiding Line Router
+ *
+ * Copyright (C) 2009  Monash University
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * See the file LICENSE.LGPL distributed with the library.
+ *
+ * Licensees holding a valid commercial license may use this file in
+ * accordance with the commercial license agreement provided with the 
+ * library.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ *
+ * Author(s):   Michael Wybrow <mjwybrow@users.sourceforge.net>
+*/
+
+
+#ifndef AVOID_ORTHOGONAL_H
+#define AVOID_ORTHOGONAL_H
+
+namespace Avoid {
+
+
+extern void generateStaticOrthogonalVisGraph(Router *router);
+
+extern void improveOrthogonalRoutes(Router *router);
+
+
+}
+
+#endif
diff --git a/src/libavoid/viscluster.cpp b/src/libavoid/viscluster.cpp
new file mode 100644
index 000000000..a127c5a17
--- /dev/null
+++ b/src/libavoid/viscluster.cpp
@@ -0,0 +1,96 @@
+/*
+ * vim: ts=4 sw=4 et tw=0 wm=0
+ *
+ * libavoid - Fast, Incremental, Object-avoiding Line Router
+ *
+ * Copyright (C) 2004-2008  Monash University
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * See the file LICENSE.LGPL distributed with the library.
+ *
+ * Licensees holding a valid commercial license may use this file in
+ * accordance with the commercial license agreement provided with the 
+ * library.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ *
+ * Author(s):   Michael Wybrow <mjwybrow@users.sourceforge.net>
+*/
+
+
+#include "libavoid/viscluster.h"
+#include "libavoid/router.h"
+#include "libavoid/assertions.h"
+
+
+namespace Avoid {
+
+
+ClusterRef::ClusterRef(Router *router, unsigned int id, Polygon& ply)
+    : _router(router)
+    , _poly(ply, router)
+    , _active(false)
+{
+    _id = router->assignId(id);
+}
+
+
+ClusterRef::~ClusterRef()
+{
+}
+
+
+void ClusterRef::makeActive(void)
+{
+    COLA_ASSERT(!_active);
+    
+    // Add to connRefs list.
+    _pos = _router->clusterRefs.insert(_router->clusterRefs.begin(), this);
+
+    _active = true;
+}
+
+
+void ClusterRef::makeInactive(void)
+{
+    COLA_ASSERT(_active);
+    
+    // Remove from connRefs list.
+    _router->clusterRefs.erase(_pos);
+
+    _active = false;
+}
+    
+
+void ClusterRef::setNewPoly(Polygon& poly)
+{
+    _poly = ReferencingPolygon(poly, _router);
+}
+
+
+unsigned int ClusterRef::id(void)
+{
+    return _id;
+}
+
+
+ReferencingPolygon& ClusterRef::polygon(void)
+{
+    return _poly;
+}
+
+
+Router *ClusterRef::router(void)
+{
+    return _router;
+}
+
+
+}
+
+
diff --git a/src/libavoid/viscluster.h b/src/libavoid/viscluster.h
new file mode 100644
index 000000000..5827e5070
--- /dev/null
+++ b/src/libavoid/viscluster.h
@@ -0,0 +1,67 @@
+/*
+ * vim: ts=4 sw=4 et tw=0 wm=0
+ *
+ * libavoid - Fast, Incremental, Object-avoiding Line Router
+ *
+ * Copyright (C) 2004-2008  Monash University
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * See the file LICENSE.LGPL distributed with the library.
+ *
+ * Licensees holding a valid commercial license may use this file in
+ * accordance with the commercial license agreement provided with the 
+ * library.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ *
+ * Author(s):   Michael Wybrow <mjwybrow@users.sourceforge.net>
+*/
+
+
+#ifndef AVOID_CLUSTER_H
+#define AVOID_CLUSTER_H
+
+#include <list>
+
+#include "libavoid/geometry.h"
+
+
+namespace Avoid {
+
+class Router;
+class ClusterRef;
+typedef std::list<ClusterRef *> ClusterRefList;
+
+
+class ClusterRef
+{
+    public:
+        ClusterRef(Router *router, unsigned int id, Polygon& poly);
+        ~ClusterRef();
+        void setNewPoly(Polygon& poly);
+        unsigned int id(void);
+        ReferencingPolygon& polygon(void);
+        Router *router(void);
+        void makeActive(void);
+        void makeInactive(void);
+
+    private:
+        Router *_router;
+        unsigned int _id;
+        ReferencingPolygon _poly;
+        bool _active;
+        ClusterRefList::iterator _pos;
+};
+
+
+}
+
+
+#endif
+
+
diff --git a/src/libavoid/vpsc.cpp b/src/libavoid/vpsc.cpp
new file mode 100644
index 000000000..c9a072375
--- /dev/null
+++ b/src/libavoid/vpsc.cpp
@@ -0,0 +1,1300 @@
+/*
+ * vim: ts=4 sw=4 et tw=0 wm=0
+ *
+ * libavoid - Fast, Incremental, Object-avoiding Line Router
+ *
+ * Copyright (C) 2005-2009  Monash University
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * See the file LICENSE.LGPL distributed with the library.
+ *
+ * Licensees holding a valid commercial license may use this file in
+ * accordance with the commercial license agreement provided with the 
+ * library.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ *
+ * Author(s):   Tim Dwyer  <Tim.Dwyer@csse.monash.edu.au>
+ *
+ * --------------
+ *
+ * This file contains a slightly modified version of Solver() from libvpsc:
+ * A solver for the problem of Variable Placement with Separation Constraints.
+ * It has the following changes from the Adaptagrams VPSC version:
+ *  -  The required VPSC code has been consolidated into a single file.
+ *  -  Unnecessary code (like Solver) has been removed.
+ *  -  The PairingHeap code has been replaced by a STL priority_queue.
+ *
+ * Modifications:  Michael Wybrow  <mjwybrow@users.sourceforge.net>
+ *
+*/
+
+#include <iostream>
+#include <math.h>
+#include <sstream>
+#include <map>
+#include <cfloat>
+#include <cstdio>
+
+#include "libavoid/vpsc.h"
+#include "libavoid/assertions.h"
+
+
+using namespace std;
+
+namespace Avoid {
+
+static const double ZERO_UPPERBOUND=-1e-10;
+static const double LAGRANGIAN_TOLERANCE=-1e-4;
+
+IncSolver::IncSolver(vector<Variable*> const &vs, vector<Constraint *> const &cs) 
+    : m(cs.size()), 
+      cs(cs), 
+      n(vs.size()), 
+      vs(vs) 
+{
+    for(unsigned i=0;i<n;++i) {
+        vs[i]->in.clear();
+        vs[i]->out.clear();
+    }
+    for(unsigned i=0;i<m;++i) {
+        Constraint *c=cs[i];
+        c->left->out.push_back(c);
+        c->right->in.push_back(c);
+    }
+    bs=new Blocks(vs);
+#ifdef LIBVPSC_LOGGING
+    printBlocks();
+    //COLA_ASSERT(!constraintGraphIsCyclic(n,vs));
+#endif
+
+    inactive=cs;
+    for(Constraints::iterator i=inactive.begin();i!=inactive.end();++i) {
+        (*i)->active=false;
+    }
+}
+IncSolver::~IncSolver() {
+    delete bs;
+}
+
+// useful in debugging
+void IncSolver::printBlocks() {
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    for(set<Block*>::iterator i=bs->begin();i!=bs->end();++i) {
+        Block *b=*i;
+        f<<"  "<<*b<<endl;
+    }
+    for(unsigned i=0;i<m;i++) {
+        f<<"  "<<*cs[i]<<endl;
+    }
+#endif
+}
+
+/*
+ * Stores the relative positions of the variables in their finalPosition
+ * field.
+ */
+void IncSolver::copyResult() {
+    for(Variables::const_iterator i=vs.begin();i!=vs.end();++i) {
+        Variable* v=*i;
+        v->finalPosition=v->position();
+        COLA_ASSERT(v->finalPosition==v->finalPosition);
+    }
+}
+
+struct node {
+    set<node*> in;
+    set<node*> out;
+};
+// useful in debugging - cycles would be BAD
+bool IncSolver::constraintGraphIsCyclic(const unsigned n, Variable* const vs[]) {
+    map<Variable*, node*> varmap;
+    vector<node*> graph;
+    for(unsigned i=0;i<n;i++) {
+        node *u=new node;
+        graph.push_back(u);
+        varmap[vs[i]]=u;
+    }
+    for(unsigned i=0;i<n;i++) {
+        for(vector<Constraint*>::iterator c=vs[i]->in.begin();c!=vs[i]->in.end();++c) {
+            Variable *l=(*c)->left;
+            varmap[vs[i]]->in.insert(varmap[l]);
+        }
+
+        for(vector<Constraint*>::iterator c=vs[i]->out.begin();c!=vs[i]->out.end();++c) {
+            Variable *r=(*c)->right;
+            varmap[vs[i]]->out.insert(varmap[r]);
+        }
+    }
+    while(graph.size()>0) {
+        node *u=NULL;
+        vector<node*>::iterator i=graph.begin();
+        for(;i!=graph.end();++i) {
+            u=*i;
+            if(u->in.size()==0) {
+                break;
+            }
+        }
+        if(i==graph.end() && graph.size()>0) {
+            //cycle found!
+            return true;
+        } else {
+            graph.erase(i);
+            for(set<node*>::iterator j=u->out.begin();j!=u->out.end();++j) {
+                node *v=*j;
+                v->in.erase(u);
+            }
+            delete u;
+        }
+    }
+    for(unsigned i=0; i<graph.size(); ++i) {
+        delete graph[i];
+    }
+    return false;
+}
+
+// useful in debugging - cycles would be BAD
+bool IncSolver::blockGraphIsCyclic() {
+    map<Block*, node*> bmap;
+    vector<node*> graph;
+    for(set<Block*>::const_iterator i=bs->begin();i!=bs->end();++i) {
+        Block *b=*i;
+        node *u=new node;
+        graph.push_back(u);
+        bmap[b]=u;
+    }
+    for(set<Block*>::const_iterator i=bs->begin();i!=bs->end();++i) {
+        Block *b=*i;
+        b->setUpInConstraints();
+        Constraint *c=b->findMinInConstraint();
+        while(c!=NULL) {
+            Block *l=c->left->block;
+            bmap[b]->in.insert(bmap[l]);
+            b->deleteMinInConstraint();
+            c=b->findMinInConstraint();
+        }
+
+        b->setUpOutConstraints();
+        c=b->findMinOutConstraint();
+        while(c!=NULL) {
+            Block *r=c->right->block;
+            bmap[b]->out.insert(bmap[r]);
+            b->deleteMinOutConstraint();
+            c=b->findMinOutConstraint();
+        }
+    }
+    while(graph.size()>0) {
+        node *u=NULL;
+        vector<node*>::iterator i=graph.begin();
+        for(;i!=graph.end();++i) {
+            u=*i;
+            if(u->in.size()==0) {
+                break;
+            }
+        }
+        if(i==graph.end() && graph.size()>0) {
+            //cycle found!
+            return true;
+        } else {
+            graph.erase(i);
+            for(set<node*>::iterator j=u->out.begin();j!=u->out.end();++j) {
+                node *v=*j;
+                v->in.erase(u);
+            }
+            delete u;
+        }
+    }
+    for(unsigned i=0; i<graph.size(); i++) {
+        delete graph[i];
+    }
+    return false;
+}
+
+bool IncSolver::solve() {
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"solve_inc()..."<<endl;
+#endif
+    satisfy();
+    double lastcost = DBL_MAX, cost = bs->cost();
+    while(fabs(lastcost-cost)>0.0001) {
+        satisfy();
+        lastcost=cost;
+        cost = bs->cost();
+#ifdef LIBVPSC_LOGGING
+        f<<"  bs->size="<<bs->size()<<", cost="<<cost<<endl;
+#endif
+    }
+    copyResult();
+    return bs->size()!=n; 
+}
+/*
+ * incremental version of satisfy that allows refinement after blocks are
+ * moved.
+ *
+ *  - move blocks to new positions
+ *  - repeatedly merge across most violated constraint until no more
+ *    violated constraints exist
+ *
+ * Note: there is a special case to handle when the most violated constraint
+ * is between two variables in the same block.  Then, we must split the block
+ * over an active constraint between the two variables.  We choose the 
+ * constraint with the most negative lagrangian multiplier. 
+ */
+bool IncSolver::satisfy() {
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"satisfy_inc()..."<<endl;
+#endif
+    splitBlocks();
+    //long splitCtr = 0;
+    Constraint* v = NULL;
+    //CBuffer buffer(inactive);
+    while((v=mostViolated(inactive))
+            &&(v->equality || v->slack() < ZERO_UPPERBOUND && !v->active)) 
+    {
+        COLA_ASSERT(!v->active);
+        Block *lb = v->left->block, *rb = v->right->block;
+        if(lb != rb) {
+            lb->merge(rb,v);
+        } else {
+            if(lb->isActiveDirectedPathBetween(v->right,v->left)) {
+                // cycle found, relax the violated, cyclic constraint
+                v->unsatisfiable=true;
+                continue;
+                //UnsatisfiableException e;
+                //lb->getActiveDirectedPathBetween(e.path,v->right,v->left);
+                //e.path.push_back(v);
+                //throw e;
+            }
+            //if(splitCtr++>10000) {
+                //throw "Cycle Error!";
+            //}
+            // constraint is within block, need to split first
+            try {
+                Constraint* splitConstraint
+                    =lb->splitBetween(v->left,v->right,lb,rb);
+                if(splitConstraint!=NULL) {
+                    COLA_ASSERT(!splitConstraint->active);
+                    inactive.push_back(splitConstraint);
+                } else {
+                    v->unsatisfiable=true;
+                    continue;
+                }
+            } catch(UnsatisfiableException e) {
+                e.path.push_back(v);
+                std::cerr << "Unsatisfiable:" << std::endl;
+                for(std::vector<Constraint*>::iterator r=e.path.begin();
+                        r!=e.path.end();++r)
+                {
+                    std::cerr << **r <<std::endl;
+                }
+                v->unsatisfiable=true;
+                continue;
+            }
+            if(v->slack()>=0) {
+                COLA_ASSERT(!v->active);
+                // v was satisfied by the above split!
+                inactive.push_back(v);
+                bs->insert(lb);
+                bs->insert(rb);
+            } else {
+                bs->insert(lb->merge(rb,v));
+            }
+        }
+        bs->cleanup();
+#ifdef LIBVPSC_LOGGING
+        f<<"...remaining blocks="<<bs->size()<<", cost="<<bs->cost()<<endl;
+#endif
+    }
+#ifdef LIBVPSC_LOGGING
+    f<<"  finished merges."<<endl;
+#endif
+    bs->cleanup();
+    bool activeConstraints=false;
+    for(unsigned i=0;i<m;i++) {
+        v=cs[i];
+        if(v->active) activeConstraints=true;
+        if(v->slack() < ZERO_UPPERBOUND) {
+            ostringstream s;
+            s<<"Unsatisfied constraint: "<<*v;
+#ifdef LIBVPSC_LOGGING
+            ofstream f(LOGFILE,ios::app);
+            f<<s.str()<<endl;
+#endif
+            throw s.str().c_str();
+        }
+    }
+#ifdef LIBVPSC_LOGGING
+    f<<"  finished cleanup."<<endl;
+    printBlocks();
+#endif
+    copyResult();
+    return activeConstraints;
+}
+void IncSolver::moveBlocks() {
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"moveBlocks()..."<<endl;
+#endif
+    for(set<Block*>::const_iterator i(bs->begin());i!=bs->end();++i) {
+        Block *b = *i;
+        b->updateWeightedPosition();
+        //b->posn = b->wposn / b->weight;
+    }
+#ifdef LIBVPSC_LOGGING
+    f<<"  moved blocks."<<endl;
+#endif
+}
+void IncSolver::splitBlocks() {
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+#endif
+    moveBlocks();
+    splitCnt=0;
+    // Split each block if necessary on min LM
+    for(set<Block*>::const_iterator i(bs->begin());i!=bs->end();++i) {
+        Block* b = *i;
+        Constraint* v=b->findMinLM();
+        if(v!=NULL && v->lm < LAGRANGIAN_TOLERANCE) {
+            COLA_ASSERT(!v->equality);
+#ifdef LIBVPSC_LOGGING
+            f<<"    found split point: "<<*v<<" lm="<<v->lm<<endl;
+#endif
+            splitCnt++;
+            Block *b = v->left->block, *l=NULL, *r=NULL;
+            COLA_ASSERT(v->left->block == v->right->block);
+            //double pos = b->posn;
+            b->split(l,r,v);
+            //l->posn=r->posn=pos;
+            //l->wposn = l->posn * l->weight;
+            //r->wposn = r->posn * r->weight;
+            l->updateWeightedPosition();
+            r->updateWeightedPosition();
+            bs->insert(l);
+            bs->insert(r);
+            b->deleted=true;
+            COLA_ASSERT(!v->active);
+            inactive.push_back(v);
+#ifdef LIBVPSC_LOGGING
+            f<<"  new blocks: "<<*l<<" and "<<*r<<endl;
+#endif
+        }
+    }
+    //if(splitCnt>0) { std::cout<<"  splits: "<<splitCnt<<endl; }
+#ifdef LIBVPSC_LOGGING
+    f<<"  finished splits."<<endl;
+#endif
+    bs->cleanup();
+}
+
+/*
+ * Scan constraint list for the most violated constraint, or the first equality
+ * constraint
+ */
+Constraint* IncSolver::mostViolated(Constraints &l) {
+    double minSlack = DBL_MAX;
+    Constraint* v=NULL;
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"Looking for most violated..."<<endl;
+#endif
+    Constraints::iterator end = l.end();
+    Constraints::iterator deletePoint = end;
+    for(Constraints::iterator i=l.begin();i!=end;++i) {
+        Constraint *c=*i;
+        double slack = c->slack();
+        if(c->equality || slack < minSlack) {
+            minSlack=slack;    
+            v=c;
+            deletePoint=i;
+            if(c->equality) break;
+        }
+    }
+    // Because the constraint list is not order dependent we just
+    // move the last element over the deletePoint and resize
+    // downwards.  There is always at least 1 element in the
+    // vector because of search.
+    if(deletePoint != end && (minSlack < ZERO_UPPERBOUND && !v->active || v->equality)) {
+        *deletePoint = l[l.size()-1];
+        l.resize(l.size()-1);
+    }
+#ifdef LIBVPSC_LOGGING
+    f<<"  most violated is: "<<*v<<endl;
+#endif
+    return v;
+}
+
+
+using std::set;
+using std::vector;
+using std::iterator;
+using std::list;
+using std::copy;
+#define __NOTNAN(p) (p)==(p)
+
+long blockTimeCtr;
+
+Blocks::Blocks(vector<Variable*> const &vs) : vs(vs),nvs(vs.size()) {
+    blockTimeCtr=0;
+    for(int i=0;i<nvs;i++) {
+        insert(new Block(vs[i]));
+    }
+}
+Blocks::~Blocks(void)
+{
+    blockTimeCtr=0;
+    for(set<Block*>::iterator i=begin();i!=end();++i) {
+        delete *i;
+    }
+    clear();
+}
+
+/*
+ * returns a list of variables with total ordering determined by the constraint 
+ * DAG
+ */
+list<Variable*> *Blocks::totalOrder() {
+    list<Variable*> *order = new list<Variable*>;
+    for(int i=0;i<nvs;i++) {
+        vs[i]->visited=false;
+    }
+    for(int i=0;i<nvs;i++) {
+        if(vs[i]->in.size()==0) {
+            dfsVisit(vs[i],order);
+        }
+    }
+    return order;
+}
+// Recursive depth first search giving total order by pushing nodes in the DAG
+// onto the front of the list when we finish searching them
+void Blocks::dfsVisit(Variable *v, list<Variable*> *order) {
+    v->visited=true;
+    vector<Constraint*>::iterator it=v->out.begin();
+    for(;it!=v->out.end();++it) {
+        Constraint *c=*it;
+        if(!c->right->visited) {
+            dfsVisit(c->right, order);
+        }
+    }    
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"  order="<<*v<<endl;
+#endif
+    order->push_front(v);
+}
+/*
+ * Processes incoming constraints, most violated to least, merging with the
+ * neighbouring (left) block until no more violated constraints are found
+ */
+void Blocks::mergeLeft(Block *r) {    
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"mergeLeft called on "<<*r<<endl;
+#endif
+    r->timeStamp=++blockTimeCtr;
+    r->setUpInConstraints();
+    Constraint *c=r->findMinInConstraint();
+    while (c != NULL && c->slack()<0) {
+#ifdef LIBVPSC_LOGGING
+        f<<"mergeLeft on constraint: "<<*c<<endl;
+#endif
+        r->deleteMinInConstraint();
+        Block *l = c->left->block;        
+        if (l->in==NULL) l->setUpInConstraints();
+        double dist = c->right->offset - c->left->offset - c->gap;
+        if (r->vars->size() < l->vars->size()) {
+            dist=-dist;
+            std::swap(l, r);
+        }
+        blockTimeCtr++;
+        r->merge(l, c, dist);
+        r->mergeIn(l);
+        r->timeStamp=blockTimeCtr;
+        removeBlock(l);
+        c=r->findMinInConstraint();
+    }        
+#ifdef LIBVPSC_LOGGING
+    f<<"merged "<<*r<<endl;
+#endif
+}    
+/*
+ * Symmetrical to mergeLeft
+ */
+void Blocks::mergeRight(Block *l) {    
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"mergeRight called on "<<*l<<endl;
+#endif    
+    l->setUpOutConstraints();
+    Constraint *c = l->findMinOutConstraint();
+    while (c != NULL && c->slack()<0) {        
+#ifdef LIBVPSC_LOGGING
+        f<<"mergeRight on constraint: "<<*c<<endl;
+#endif
+        l->deleteMinOutConstraint();
+        Block *r = c->right->block;
+        r->setUpOutConstraints();
+        double dist = c->left->offset + c->gap - c->right->offset;
+        if (l->vars->size() > r->vars->size()) {
+            dist=-dist;
+            std::swap(l, r);
+        }
+        l->merge(r, c, dist);
+        l->mergeOut(r);
+        removeBlock(r);
+        c=l->findMinOutConstraint();
+    }    
+#ifdef LIBVPSC_LOGGING
+    f<<"merged "<<*l<<endl;
+#endif
+}
+void Blocks::removeBlock(Block *doomed) {
+    doomed->deleted=true;
+    //erase(doomed);
+}
+void Blocks::cleanup() {
+    vector<Block*> bcopy(begin(),end());
+    for(vector<Block*>::iterator i=bcopy.begin();i!=bcopy.end();++i) {
+        Block *b=*i;
+        if(b->deleted) {
+            erase(b);
+            delete b;
+        }
+    }
+}
+/*
+ * Splits block b across constraint c into two new blocks, l and r (c's left
+ * and right sides respectively)
+ */
+void Blocks::split(Block *b, Block *&l, Block *&r, Constraint *c) {
+    b->split(l,r,c);
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"Split left: "<<*l<<endl;
+    f<<"Split right: "<<*r<<endl;
+#endif
+    r->posn = b->posn;
+    //COLA_ASSERT(r->weight!=0);
+    //r->wposn = r->posn * r->weight;
+    mergeLeft(l);
+    // r may have been merged!
+    r = c->right->block;
+    r->updateWeightedPosition();
+    //r->posn = r->wposn / r->weight;
+    mergeRight(r);
+    removeBlock(b);
+
+    insert(l);
+    insert(r);
+    COLA_ASSERT(__NOTNAN(l->posn));
+    COLA_ASSERT(__NOTNAN(r->posn));
+}
+/*
+ * returns the cost total squared distance of variables from their desired
+ * positions
+ */
+double Blocks::cost() {
+    double c = 0;
+    for(set<Block*>::iterator i=begin();i!=end();++i) {
+        c += (*i)->cost();
+    }
+    return c;
+}
+
+void PositionStats::addVariable(Variable* v) {
+    double ai=scale/v->scale;
+    double bi=v->offset/v->scale;
+    double wi=v->weight;
+    AB+=wi*ai*bi;
+    AD+=wi*ai*v->desiredPosition;
+    A2+=wi*ai*ai;
+    /*
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f << "adding v[" << v->id << "], blockscale=" << scale << ", despos=" 
+      << v->desiredPosition << ", ai=" << ai << ", bi=" << bi
+      << ", AB=" << AB << ", AD=" << AD << ", A2=" << A2;
+#endif
+*/
+}
+void Block::addVariable(Variable* v) {
+    v->block=this;
+    vars->push_back(v);
+    if(ps.A2==0) ps.scale=v->scale;
+    //weight+= v->weight;
+    //wposn += v->weight * (v->desiredPosition - v->offset);
+    //posn=wposn/weight;
+    ps.addVariable(v);
+    posn=(ps.AD - ps.AB) / ps.A2;
+    COLA_ASSERT(__NOTNAN(posn));
+    /*
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f << ", posn=" << posn << endl;
+#endif
+*/
+}
+Block::Block(Variable* const v)
+    : vars(new vector<Variable*>)
+    , posn(0)
+    //, weight(0)
+    //, wposn(0)
+    , deleted(false)
+    , timeStamp(0)
+    , in(NULL)
+    , out(NULL)
+{
+    if(v!=NULL) {
+        v->offset=0;
+        addVariable(v);
+    }
+}
+
+void Block::updateWeightedPosition() {
+    //wposn=0;
+    ps.AB=ps.AD=ps.A2=0;
+    for (Vit v=vars->begin();v!=vars->end();++v) {
+        //wposn += ((*v)->desiredPosition - (*v)->offset) * (*v)->weight;
+        ps.addVariable(*v);
+    }
+    posn=(ps.AD - ps.AB) / ps.A2;
+    COLA_ASSERT(__NOTNAN(posn));
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f << ", posn=" << posn << endl;
+#endif
+}
+Block::~Block(void)
+{
+    delete vars;
+    delete in;
+    delete out;
+}
+void Block::setUpInConstraints() {
+    setUpConstraintHeap(in,true);
+}
+void Block::setUpOutConstraints() {
+    setUpConstraintHeap(out,false);
+}
+void Block::setUpConstraintHeap(Heap* &h,bool in) {
+    delete h;
+    h = new Heap();
+    for (Vit i=vars->begin();i!=vars->end();++i) {
+        Variable *v=*i;
+        vector<Constraint*> *cs=in?&(v->in):&(v->out);
+        for (Cit j=cs->begin();j!=cs->end();++j) {
+            Constraint *c=*j;
+            c->timeStamp=blockTimeCtr;
+            if (c->left->block != this && in || c->right->block != this && !in) {
+                h->push(c);
+            }
+        }
+    }
+}    
+Block* Block::merge(Block* b, Constraint* c) {
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"  merging on: "<<*c<<",c->left->offset="<<c->left->offset<<",c->right->offset="<<c->right->offset<<endl;
+#endif
+    double dist = c->right->offset - c->left->offset - c->gap;
+    Block *l=c->left->block;
+    Block *r=c->right->block;
+    if (l->vars->size() < r->vars->size()) {
+        r->merge(l,c,dist);
+    } else {
+               l->merge(r,c,-dist);
+    }
+    Block* mergeBlock=b->deleted?this:b;
+#ifdef LIBVPSC_LOGGING
+    f<<"  merged block="<<*mergeBlock<<endl;
+#endif
+    return mergeBlock;
+}
+/*
+ * Merges b into this block across c.  Can be either a
+ * right merge or a left merge
+ * @param b block to merge into this
+ * @param c constraint being merged
+ * @param distance separation required to satisfy c
+ */
+void Block::merge(Block *b, Constraint *c, double dist) {
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"    merging: "<<*b<<"dist="<<dist<<endl;
+#endif
+    c->active=true;
+    //wposn+=b->wposn-dist*b->weight;
+    //weight+=b->weight;
+    for(Vit i=b->vars->begin();i!=b->vars->end();++i) {
+        Variable *v=*i;
+        //v->block=this;
+        //vars->push_back(v);
+        v->offset+=dist;
+        addVariable(v);
+    }
+#ifdef LIBVPSC_LOGGING
+    for(Vit i=vars->begin();i!=vars->end();++i) {
+        Variable *v=*i;
+        f<<"    v["<<v->id<<"]: d="<<v->desiredPosition
+            <<" a="<<v->scale<<" o="<<v->offset
+            <<endl;
+    }
+    f<<"  AD="<<ps.AD<<" AB="<<ps.AB<<" A2="<<ps.A2<<endl;
+#endif
+    //posn=wposn/weight;
+    //COLA_ASSERT(wposn==ps.AD - ps.AB);
+    posn=(ps.AD - ps.AB) / ps.A2;
+    COLA_ASSERT(__NOTNAN(posn));
+    b->deleted=true;
+}
+
+void Block::mergeIn(Block *b) {
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"  merging constraint heaps... "<<endl;
+#endif
+    // We check the top of the heaps to remove possible internal constraints
+    findMinInConstraint();
+    b->findMinInConstraint();
+    while (!b->in->empty())
+    {
+        in->push(b->in->top());
+        b->in->pop();
+    }
+#ifdef LIBVPSC_LOGGING
+    f<<"  merged heap: "<<*in<<endl;
+#endif
+}
+void Block::mergeOut(Block *b) {    
+    findMinOutConstraint();
+    b->findMinOutConstraint();
+    while (!b->out->empty())
+    {
+        out->push(b->out->top());
+        b->out->pop();
+    }
+}
+Constraint *Block::findMinInConstraint() {
+    Constraint *v = NULL;
+    vector<Constraint*> outOfDate;
+    while (!in->empty()) {
+        v = in->top();
+        Block *lb=v->left->block;
+        Block *rb=v->right->block;
+        // rb may not be this if called between merge and mergeIn
+#ifdef LIBVPSC_LOGGING
+        ofstream f(LOGFILE,ios::app);
+        f<<"  checking constraint ... "<<*v;
+        f<<"    timestamps: left="<<lb->timeStamp<<" right="<<rb->timeStamp<<" constraint="<<v->timeStamp<<endl;
+#endif
+        if(lb == rb) {
+            // constraint has been merged into the same block
+#ifdef LIBVPSC_LOGGING
+            if(v->slack()<0) {
+                f<<"  violated internal constraint found! "<<*v<<endl;
+                f<<"     lb="<<*lb<<endl;
+                f<<"     rb="<<*rb<<endl;
+            }
+#endif
+            in->pop();
+#ifdef LIBVPSC_LOGGING
+            f<<" ... skipping internal constraint"<<endl;
+#endif
+        } else if(v->timeStamp < lb->timeStamp) {
+            // block at other end of constraint has been moved since this
+            in->pop();
+            outOfDate.push_back(v);
+#ifdef LIBVPSC_LOGGING
+            f<<"    reinserting out of date (reinsert later)"<<endl;
+#endif
+        } else {
+            break;
+        }
+    }
+    for(Cit i=outOfDate.begin();i!=outOfDate.end();++i) {
+        v=*i;
+        v->timeStamp=blockTimeCtr;
+        in->push(v);
+    }
+    if(in->empty()) {
+        v=NULL;
+    } else {
+        v=in->top();
+    }
+    return v;
+}
+Constraint *Block::findMinOutConstraint() {
+    if(out->empty()) return NULL;
+    Constraint *v = out->top();
+    while (v->left->block == v->right->block) {
+        out->pop();
+        if(out->empty()) return NULL;
+        v = out->top();
+    }
+    return v;
+}
+void Block::deleteMinInConstraint() {
+    in->pop();
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"deleteMinInConstraint... "<<endl;
+    f<<"  result: "<<*in<<endl;
+#endif
+}
+void Block::deleteMinOutConstraint() {
+    out->pop();
+}
+inline bool Block::canFollowLeft(Constraint const* c, Variable const* last) const {
+    return c->left->block==this && c->active && last!=c->left;
+}
+inline bool Block::canFollowRight(Constraint const* c, Variable const* last) const {
+    return c->right->block==this && c->active && last!=c->right;
+}
+
+// computes the derivative of v and the lagrange multipliers
+// of v's out constraints (as the recursive sum of those below.
+// Does not backtrack over u.
+// also records the constraint with minimum lagrange multiplier
+// in min_lm
+double Block::compute_dfdv(Variable* const v, Variable* const u,
+               Constraint *&min_lm) {
+    double dfdv=v->dfdv();
+    for(Cit it=v->out.begin();it!=v->out.end();++it) {
+        Constraint *c=*it;
+        if(canFollowRight(c,u)) {
+            c->lm=compute_dfdv(c->right,v,min_lm);
+            dfdv+=c->lm*c->left->scale;
+            if(!c->equality&&(min_lm==NULL||c->lm<min_lm->lm)) min_lm=c;
+        }
+    }
+    for(Cit it=v->in.begin();it!=v->in.end();++it) {
+        Constraint *c=*it;
+        if(canFollowLeft(c,u)) {
+            c->lm=-compute_dfdv(c->left,v,min_lm);
+            dfdv-=c->lm*c->right->scale;
+            if(!c->equality&&(min_lm==NULL||c->lm<min_lm->lm)) min_lm=c;
+        }
+    }
+    return dfdv/v->scale;
+}
+double Block::compute_dfdv(Variable* const v, Variable* const u) {
+    double dfdv = v->dfdv();
+    for(Cit it = v->out.begin(); it != v->out.end(); ++it) {
+        Constraint *c = *it;
+        if(canFollowRight(c,u)) {
+            c->lm =   compute_dfdv(c->right,v);
+            dfdv += c->lm * c->left->scale;
+        }
+    }
+    for(Cit it=v->in.begin();it!=v->in.end();++it) {
+        Constraint *c = *it;
+        if(canFollowLeft(c,u)) {
+            c->lm = - compute_dfdv(c->left,v);
+            dfdv -= c->lm * c->right->scale;
+        }
+    }
+    return dfdv/v->scale;
+}
+
+// The top level v and r are variables between which we want to find the
+// constraint with the smallest lm.  
+// Similarly, m is initially NULL and is only assigned a value if the next
+// variable to be visited is r or if a possible min constraint is returned from
+// a nested call (rather than NULL).
+// Then, the search for the m with minimum lm occurs as we return from
+// the recursion (checking only constraints traversed left-to-right 
+// in order to avoid creating any new violations).
+// We also do not consider equality constraints as potential split points
+bool Block::split_path(
+    Variable* r, 
+    Variable* const v, 
+    Variable* const u, 
+    Constraint* &m,
+    bool desperation=false
+    ) 
+{
+    for(Cit it(v->in.begin());it!=v->in.end();++it) {
+        Constraint *c=*it;
+        if(canFollowLeft(c,u)) {
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"  left split path: "<<*c<<endl;
+#endif
+            if(c->left==r) {
+                if(desperation&&!c->equality) m=c;
+                return true;
+            } else {
+                if(split_path(r,c->left,v,m)) {
+                    if(desperation && !c->equality && (!m||c->lm<m->lm)) {
+                               m=c;
+                    }
+                    return true;
+                }
+            }
+        }
+    }
+    for(Cit it(v->out.begin());it!=v->out.end();++it) {
+        Constraint *c=*it;
+        if(canFollowRight(c,u)) {
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"  right split path: "<<*c<<endl;
+#endif
+            if(c->right==r) {
+                if(!c->equality) m=c;
+                return true;
+            } else {
+                if(split_path(r,c->right,v,m)) {
+                    if(!c->equality && (!m||c->lm<m->lm))
+                               m=c;
+                    return true;
+                }
+            }
+        }
+    }
+    return false;
+}
+/*
+Block::Pair Block::compute_dfdv_between(
+        Variable* r, Variable* const v, Variable* const u, 
+        const Direction dir = NONE, bool changedDirection = false) {
+    double dfdv=v->weight*(v->position() - v->desiredPosition);
+    Constraint *m=NULL;
+    for(Cit it(v->in.begin());it!=v->in.end();++it) {
+        Constraint *c=*it;
+        if(canFollowLeft(c,u)) {
+            if(dir==RIGHT) { 
+                changedDirection = true; 
+            }
+            if(c->left==r) {
+                       r=NULL;
+                    if(!c->equality) m=c; 
+            }
+            Pair p=compute_dfdv_between(r,c->left,v,
+                    LEFT,changedDirection);
+            dfdv -= c->lm = -p.first;
+            if(r && p.second) 
+                m = p.second;
+        }
+    }
+    for(Cit it(v->out.begin());it!=v->out.end();++it) {
+        Constraint *c=*it;
+        if(canFollowRight(c,u)) {
+            if(dir==LEFT) { 
+                changedDirection = true; 
+            }
+            if(c->right==r) {
+                       r=NULL; 
+                    if(!c->equality) m=c; 
+            }
+            Pair p=compute_dfdv_between(r,c->right,v,
+                    RIGHT,changedDirection);
+            dfdv += c->lm = p.first;
+            if(r && p.second) 
+                m = changedDirection && !c->equality && c->lm < p.second->lm 
+                    ? c 
+                    : p.second;
+        }
+    }
+    return Pair(dfdv,m);
+}
+*/
+
+// resets LMs for all active constraints to 0 by
+// traversing active constraint tree starting from v,
+// not back tracking over u
+void Block::reset_active_lm(Variable* const v, Variable* const u) {
+    for(Cit it=v->out.begin();it!=v->out.end();++it) {
+        Constraint *c=*it;
+        if(canFollowRight(c,u)) {
+            c->lm=0;
+            reset_active_lm(c->right,v);
+        }
+    }
+    for(Cit it=v->in.begin();it!=v->in.end();++it) {
+        Constraint *c=*it;
+        if(canFollowLeft(c,u)) {
+            c->lm=0;
+            reset_active_lm(c->left,v);
+        }
+    }
+}
+void Block::list_active(Variable* const v, Variable* const u) {
+    for(Cit it=v->out.begin();it!=v->out.end();++it) {
+        Constraint *c=*it;
+        if(canFollowRight(c,u)) {
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"  "<<*c<<endl;
+#endif
+            list_active(c->right,v);
+        }
+    }
+    for(Cit it=v->in.begin();it!=v->in.end();++it) {
+        Constraint *c=*it;
+        if(canFollowLeft(c,u)) {
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"  "<<*c<<endl;
+#endif
+            list_active(c->left,v);
+        }
+    }
+}
+/*
+ * finds the constraint with the minimum lagrange multiplier, that is, the constraint
+ * that most wants to split
+ */
+Constraint *Block::findMinLM() {
+    Constraint *min_lm=NULL;
+    reset_active_lm(vars->front(),NULL);
+    compute_dfdv(vars->front(),NULL,min_lm);
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"  langrangians: "<<endl;
+    list_active(vars->front(),NULL);
+#endif
+    return min_lm;
+}
+Constraint *Block::findMinLMBetween(Variable* const lv, Variable* const rv) {
+    reset_active_lm(vars->front(),NULL);
+    compute_dfdv(vars->front(),NULL);
+    Constraint *min_lm=NULL;
+    split_path(rv,lv,NULL,min_lm);
+#if 0
+    if(min_lm==NULL) {
+        split_path(rv,lv,NULL,min_lm,true);
+    }
+#else
+    if(min_lm==NULL) {
+        fprintf(stderr,"Couldn't find split point!\n");
+        UnsatisfiableException e;
+        getActivePathBetween(e.path,lv,rv,NULL);
+        throw e;
+    }
+    COLA_ASSERT(min_lm!=NULL);
+#endif
+    return min_lm;
+}
+
+// populates block b by traversing the active constraint tree adding variables as they're 
+// visited.  Starts from variable v and does not backtrack over variable u.
+void Block::populateSplitBlock(Block *b, Variable* v, Variable const* u) {
+    b->addVariable(v);
+    for (Cit c=v->in.begin();c!=v->in.end();++c) {
+        if (canFollowLeft(*c,u))
+            populateSplitBlock(b, (*c)->left, v);
+    }
+    for (Cit c=v->out.begin();c!=v->out.end();++c) {
+        if (canFollowRight(*c,u)) 
+            populateSplitBlock(b, (*c)->right, v);
+    }
+}
+/*
+ * Returns the active path between variables u and v... not back tracking over w
+ */
+bool Block::getActivePathBetween(Constraints& path, Variable const* u,
+               Variable const* v, Variable const *w) const {
+    if(u==v) return true;
+    for (Cit_const c=u->in.begin();c!=u->in.end();++c) {
+        if (canFollowLeft(*c,w)) {
+            if(getActivePathBetween(path, (*c)->left, v, u)) {
+                path.push_back(*c);
+                return true;
+            }
+        }
+    }
+    for (Cit_const c=u->out.begin();c!=u->out.end();++c) {
+        if (canFollowRight(*c,w)) {
+            if(getActivePathBetween(path, (*c)->right, v, u)) {
+                path.push_back(*c);
+                return true;
+            }
+        }
+    }
+    return false;
+}
+// Search active constraint tree from u to see if there is a directed path to v.
+// Returns true if path is found with all constraints in path having their visited flag
+// set true.
+bool Block::isActiveDirectedPathBetween(Variable const* u, Variable const* v) const {
+    if(u==v) return true;
+    for (Cit_const c=u->out.begin();c!=u->out.end();++c) {
+        if(canFollowRight(*c,NULL)) {
+            if(isActiveDirectedPathBetween((*c)->right,v)) {
+                return true;
+            }
+        }
+    }
+    return false;
+}
+bool Block::getActiveDirectedPathBetween(
+        Constraints& path, Variable const* u, Variable const* v) const {
+    if(u==v) return true;
+    for (Cit_const c=u->out.begin();c!=u->out.end();++c) {
+        if(canFollowRight(*c,NULL)) {
+            if(getActiveDirectedPathBetween(path,(*c)->right,v)) {
+                path.push_back(*c);
+                return true;
+            }
+        }
+    }
+    return false;
+}
+/*
+ * Block needs to be split because of a violated constraint between vl and vr.
+ * We need to search the active constraint tree between l and r and find the constraint
+ * with min lagrangrian multiplier and split at that point.
+ * Returns the split constraint
+ */
+Constraint* Block::splitBetween(Variable* const vl, Variable* const vr,
+               Block* &lb, Block* &rb) {
+#ifdef LIBVPSC_LOGGING
+    ofstream f(LOGFILE,ios::app);
+    f<<"  need to split between: "<<*vl<<" and "<<*vr<<endl;
+#endif
+    Constraint *c=findMinLMBetween(vl, vr);
+#ifdef LIBVPSC_LOGGING
+    f<<"  going to split on: "<<*c<<endl;
+#endif
+    if(c!=NULL) {
+        split(lb,rb,c);
+        deleted = true;
+    }
+    return c;
+}
+
+/*
+ * Creates two new blocks, l and r, and splits this block across constraint c,
+ * placing the left subtree of constraints (and associated variables) into l
+ * and the right into r.
+ */
+void Block::split(Block* &l, Block* &r, Constraint* c) {
+    c->active=false;
+    l=new Block();
+    populateSplitBlock(l,c->left,c->right);
+    //COLA_ASSERT(l->weight>0);
+    r=new Block();
+    populateSplitBlock(r,c->right,c->left);
+    //COLA_ASSERT(r->weight>0);
+}
+
+/*
+ * Computes the cost (squared euclidean distance from desired positions) of the
+ * current positions for variables in this block
+ */
+double Block::cost() {
+    double c = 0;
+    for (Vit v=vars->begin();v!=vars->end();++v) {
+        double diff = (*v)->position() - (*v)->desiredPosition;
+        c += (*v)->weight * diff * diff;
+    }
+    return c;
+}
+ostream& operator <<(ostream &os, const Block& b)
+{
+    os<<"Block(posn="<<b.posn<<"):";
+    for(Block::Vit v=b.vars->begin();v!=b.vars->end();++v) {
+        os<<" "<<**v;
+    }
+    if(b.deleted) {
+        os<<" Deleted!";
+    }
+    return os;
+}
+
+Constraint::Constraint(Variable *left, Variable *right, double gap, bool equality)
+: left(left),
+  right(right),
+  gap(gap),
+  timeStamp(0),
+  active(false),
+  equality(equality),
+  unsatisfiable(false)
+{
+    // In hindsight I think it's probably better to build the constraint DAG
+    // (by creating variable in/out lists) when needed, rather than in advance
+    //left->out.push_back(this);
+    //right->in.push_back(this);
+}
+Constraint::~Constraint() {
+    // see constructor: the following is just way too slow.  
+    // Better to create a
+    // new DAG on demand than maintain the lists dynamically.
+    //Constraints::iterator i;
+    //for(i=left->out.begin(); i!=left->out.end(); i++) {
+        //if(*i==this) break;
+    //}
+    //left->out.erase(i);
+    //for(i=right->in.begin(); i!=right->in.end(); i++) {
+        //if(*i==this) break;
+    //}
+    //right->in.erase(i);
+}
+double Constraint::slack() const { 
+    return unsatisfiable ? DBL_MAX
+           : right->scale * right->position() 
+         - gap - left->scale * left->position(); 
+}
+std::ostream& operator <<(std::ostream &os, const Constraint &c)
+{
+    if(&c==NULL) {
+        os<<"NULL";
+    } else {
+        const char *type=c.equality?"=":"<=";
+        std::ostringstream lscale, rscale;
+        if(c.left->scale!=1) {
+            lscale << c.left->scale << "*";
+        }
+        if(c.right->scale!=1) {
+            rscale << c.right->scale << "*";
+        }
+        os<<lscale.str()<<*c.left<<"+"<<c.gap<<type<<rscale.str()<<*c.right;
+        if(c.left->block&&c.right->block)
+            os<<"("<<c.slack()<<")"<<(c.active?"-active":"")
+                <<"(lm="<<c.lm<<")";
+        else
+            os<<"(vars have no position)";
+    }
+    return os;
+}
+
+bool CompareConstraints::operator() (
+    Constraint *const &l, Constraint *const &r
+) const {
+    double const sl = 
+        l->left->block->timeStamp > l->timeStamp
+        ||l->left->block==l->right->block
+        ?-DBL_MAX:l->slack();
+    double const sr = 
+        r->left->block->timeStamp > r->timeStamp
+        ||r->left->block==r->right->block
+        ?-DBL_MAX:r->slack();
+    if(sl==sr) {
+        // arbitrary choice based on id
+        if(l->left->id==r->left->id) {
+            if(l->right->id<r->right->id) return true;
+            return false;
+        }
+        if(l->left->id<r->left->id) return true;
+        return false;
+    }
+    return sl > sr;
+}
+
+std::ostream& operator <<(std::ostream &os, const Variable &v) {
+    if(v.block)
+        os << "(" << v.id << "=" << v.position() << ")";
+    else
+        os << "(" << v.id << "=" << v.desiredPosition << ")";
+    return os;
+}
+
+}
diff --git a/src/libavoid/vpsc.h b/src/libavoid/vpsc.h
new file mode 100644
index 000000000..4d6d8ce61
--- /dev/null
+++ b/src/libavoid/vpsc.h
@@ -0,0 +1,255 @@
+/*
+ * vim: ts=4 sw=4 et tw=0 wm=0
+ *
+ * libavoid - Fast, Incremental, Object-avoiding Line Router
+ *
+ * Copyright (C) 2005-2009  Monash University
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * See the file LICENSE.LGPL distributed with the library.
+ *
+ * Licensees holding a valid commercial license may use this file in
+ * accordance with the commercial license agreement provided with the 
+ * library.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ *
+ * Author(s):   Tim Dwyer  <Tim.Dwyer@csse.monash.edu.au>
+ *
+ * --------------
+ *
+ * This file contains a slightly modified version of Solver() from libvpsc:
+ * A solver for the problem of Variable Placement with Separation Constraints.
+ * It has the following changes from the Adaptagrams VPSC version:
+ *  -  The required VPSC code has been consolidated into a single file.
+ *  -  Unnecessary code (like Solver) has been removed.
+ *  -  The PairingHeap code has been replaced by a STL priority_queue.
+ *
+ * Modifications:  Michael Wybrow  <mjwybrow@users.sourceforge.net>
+ *
+*/
+
+#ifndef LIBAVOID_VPSC_H
+#define LIBAVOID_VPSC_H
+
+#include <vector>
+#include <list>
+#include <set>
+#include <queue>
+
+namespace Avoid {
+
+class Variable;
+class Constraint;
+typedef std::vector<Variable*> Variables;
+typedef std::vector<Constraint*> Constraints;
+class CompareConstraints {
+public:
+    bool operator() (Constraint *const &l, Constraint *const &r) const;
+};
+struct PositionStats {
+    PositionStats() : scale(0), AB(0), AD(0), A2(0) {}
+    void addVariable(Variable* const v);
+    double scale;
+    double AB;
+    double AD;
+    double A2;
+};
+
+typedef std::priority_queue<Constraint*,std::vector<Constraint*>,
+        CompareConstraints> Heap;
+
+class Block
+{
+    typedef Variables::iterator Vit;
+    typedef Constraints::iterator Cit;
+    typedef Constraints::const_iterator Cit_const;
+
+    friend std::ostream& operator <<(std::ostream &os,const Block &b);
+public:
+    Variables *vars;
+    double posn;
+    //double weight;
+    //double wposn;
+    PositionStats ps;
+    Block(Variable* const v=NULL);
+    ~Block(void);
+    Constraint* findMinLM();
+    Constraint* findMinLMBetween(Variable* const lv, Variable* const rv);
+    Constraint* findMinInConstraint();
+    Constraint* findMinOutConstraint();
+    void deleteMinInConstraint();
+    void deleteMinOutConstraint();
+    void updateWeightedPosition();
+    void merge(Block *b, Constraint *c, double dist);
+    Block* merge(Block *b, Constraint *c);
+    void mergeIn(Block *b);
+    void mergeOut(Block *b);
+    void split(Block *&l, Block *&r, Constraint *c);
+    Constraint* splitBetween(Variable* vl, Variable* vr, Block* &lb, Block* &rb);
+    void setUpInConstraints();
+    void setUpOutConstraints();
+    double cost();
+    bool deleted;
+    long timeStamp;
+    Heap *in;
+    Heap *out;
+    bool getActivePathBetween(Constraints& path, Variable const* u,
+               Variable const* v, Variable const *w) const;
+    bool isActiveDirectedPathBetween(
+            Variable const* u, Variable const* v) const;
+    bool getActiveDirectedPathBetween(Constraints& path, Variable const * u, Variable const * v) const;
+private:
+    typedef enum {NONE, LEFT, RIGHT} Direction;
+    typedef std::pair<double, Constraint*> Pair;
+    void reset_active_lm(Variable* const v, Variable* const u);
+    void list_active(Variable* const v, Variable* const u);
+    double compute_dfdv(Variable* const v, Variable* const u);
+    double compute_dfdv(Variable* const v, Variable* const u, Constraint *&min_lm);
+    bool split_path(Variable*, Variable* const, Variable* const, 
+            Constraint* &min_lm, bool desperation);
+    bool canFollowLeft(Constraint const* c, Variable const* last) const;
+    bool canFollowRight(Constraint const* c, Variable const* last) const;
+    void populateSplitBlock(Block *b, Variable* v, Variable const* u);
+    void addVariable(Variable* v);
+    void setUpConstraintHeap(Heap* &h,bool in);
+};
+
+
+class Constraint;
+typedef std::vector<Constraint*> Constraints;
+class Variable
+{
+    friend std::ostream& operator <<(std::ostream &os, const Variable &v);
+    friend class Block;
+    friend class Constraint;
+    friend class IncSolver;
+public:
+    int id; // useful in log files
+    double desiredPosition;
+    double finalPosition;
+    double weight; // how much the variable wants to 
+                   // be at it's desired position
+    double scale; // translates variable to another space
+    double offset;
+    Block *block;
+    bool visited;
+    bool fixedDesiredPosition;
+    Constraints in;
+    Constraints out;
+    char *toString();
+    inline Variable(const int id, const double desiredPos=-1.0, 
+            const double weight=1.0, const double scale=1.0)
+        : id(id)
+        , desiredPosition(desiredPos)
+        , weight(weight)
+        , scale(scale)
+        , offset(0)
+        , block(NULL)
+        , visited(false)
+        , fixedDesiredPosition(false)
+    {
+    }
+    double dfdv() const {
+        return 2. * weight * ( position() - desiredPosition );
+    }
+private:
+    double position() const {
+        return (block->ps.scale*block->posn+offset)/scale;
+    }
+};
+
+
+class Constraint
+{
+    friend std::ostream& operator <<(std::ostream &os,const Constraint &c);
+public:
+    Variable *left;
+    Variable *right;
+    double gap;
+    double lm;
+    Constraint(Variable *left, Variable *right, double gap, bool equality=false);
+    ~Constraint();
+    double slack() const;
+    long timeStamp;
+    bool active;
+    const bool equality;
+    bool unsatisfiable;
+};
+/*
+ * A block structure defined over the variables such that each block contains
+ * 1 or more variables, with the invariant that all constraints inside a block
+ * are satisfied by keeping the variables fixed relative to one another
+ */
+class Blocks : public std::set<Block*>
+{
+public:
+    Blocks(Variables const &vs);
+    ~Blocks(void);
+    void mergeLeft(Block *r);
+    void mergeRight(Block *l);
+    void split(Block *b, Block *&l, Block *&r, Constraint *c);
+    std::list<Variable*> *totalOrder();
+    void cleanup();
+    double cost();
+private:
+    void dfsVisit(Variable *v, std::list<Variable*> *order);
+    void removeBlock(Block *doomed);
+    Variables const &vs;
+    int nvs;
+};
+
+extern long blockTimeCtr;
+
+struct UnsatisfiableException {
+    Constraints path;
+};
+struct UnsatisfiedConstraint {
+    UnsatisfiedConstraint(Constraint& c):c(c) {}
+    Constraint& c;
+};
+/*
+ * Variable Placement with Separation Constraints problem instance
+ */
+class IncSolver {
+public:
+    unsigned splitCnt;
+    bool satisfy();
+    bool solve();
+    void moveBlocks();
+    void splitBlocks();
+    IncSolver(Variables const &vs, Constraints const &cs);
+
+    ~IncSolver();
+    Variables const & getVariables() { return vs; }
+protected:
+    Blocks *bs;
+    unsigned m;
+    Constraints const &cs;
+    unsigned n;
+    Variables const &vs;
+    void printBlocks();
+    void copyResult();
+private:
+    bool constraintGraphIsCyclic(const unsigned n, Variable* const vs[]);
+    bool blockGraphIsCyclic();
+    Constraints inactive;
+    Constraints violated;
+    Constraint* mostViolated(Constraints &l);
+};
+
+struct delete_object
+{
+    template <typename T>
+    void operator()(T *ptr){ delete ptr;}
+};
+
+
+}
+
+#endif // AVOID_VPSC_H