diff --git a/Triangle.NET/MeshRenderer.Core/RenderData.cs b/Triangle.NET/MeshRenderer.Core/RenderData.cs
index 7fc60fb..d2f803a 100644
--- a/Triangle.NET/MeshRenderer.Core/RenderData.cs
+++ b/Triangle.NET/MeshRenderer.Core/RenderData.cs
@@ -69,10 +69,10 @@ namespace MeshRenderer.Core
             }
 
             this.Bounds = new BoundingBox(
-                (float)data.Bounds.Xmin,
-                (float)data.Bounds.Xmax,
-                (float)data.Bounds.Ymin,
-                (float)data.Bounds.Ymax);
+                (float)data.Bounds.MinX,
+                (float)data.Bounds.MaxX,
+                (float)data.Bounds.MinY,
+                (float)data.Bounds.MaxY);
         }
 
         /// <summary>
@@ -150,10 +150,10 @@ namespace MeshRenderer.Core
             this.Triangles = triangles.ToArray();
 
             this.Bounds = new BoundingBox(
-                (float)mesh.Bounds.Xmin,
-                (float)mesh.Bounds.Xmax,
-                (float)mesh.Bounds.Ymin,
-                (float)mesh.Bounds.Ymax);
+                (float)mesh.Bounds.MinX,
+                (float)mesh.Bounds.MaxX,
+                (float)mesh.Bounds.MinY,
+                (float)mesh.Bounds.MaxY);
         }
 
         /// <summary>
diff --git a/Triangle.NET/TestApp/FormMain.cs b/Triangle.NET/TestApp/FormMain.cs
index 2681602..0bb0c18 100644
--- a/Triangle.NET/TestApp/FormMain.cs
+++ b/Triangle.NET/TestApp/FormMain.cs
@@ -8,6 +8,7 @@ using MeshExplorer.IO;
 using MeshRenderer.Core;
 using TriangleNet;
 using TriangleNet.Geometry;
+using TriangleNet.Log;
 using TriangleNet.Tools;
 
 namespace MeshExplorer
@@ -340,31 +341,38 @@ namespace MeshExplorer
 
         private bool Open(string filename)
         {
-            if (FileProcessor.ContainsMeshData(filename))
+            if (!FileProcessor.CanHandleFile(filename))
             {
-                if (DarkMessageBox.Show("Import mesh", Settings.ImportString,
-                    "Do you want to import the mesh?", MessageBoxButtons.YesNo) == DialogResult.OK)
-                {
-                    input = null;
-                    mesh = FileProcessor.Import(filename);
-
-                    if (mesh != null)
-                    {
-                        statisticView.UpdateStatistic(mesh);
-
-                        // Update settings
-                        settings.CurrentFile = Path.GetFileName(filename);
-
-                        HandleMeshImport();
-                        btnSmooth.Enabled = true; // TODO: Remove
-                    }
-                    // else Message
-
-                    return true;
-                }
+                // TODO: show message.
             }
+            else
+            {
+                if (FileProcessor.ContainsMeshData(filename))
+                {
+                    if (DarkMessageBox.Show("Import mesh", Settings.ImportString,
+                        "Do you want to import the mesh?", MessageBoxButtons.YesNo) == DialogResult.OK)
+                    {
+                        input = null;
+                        mesh = FileProcessor.Import(filename);
 
-            input = FileProcessor.Read(filename);
+                        if (mesh != null)
+                        {
+                            statisticView.UpdateStatistic(mesh);
+
+                            // Update settings
+                            settings.CurrentFile = Path.GetFileName(filename);
+
+                            HandleMeshImport();
+                            btnSmooth.Enabled = true; // TODO: Remove
+                        }
+                        // else Message
+
+                        return true;
+                    }
+                }
+
+                input = FileProcessor.Read(filename);
+            }
 
             if (input != null)
             {
@@ -737,6 +745,16 @@ namespace MeshExplorer
 
                 Behavior.Verbose = save;
 
+                if (isConsistent)
+                {
+                    SimpleLog.Instance.Info("Mesh topology appears to be consistent.");
+                }
+
+                if (isDelaunay)
+                {
+                    SimpleLog.Instance.Info("Mesh is (conforming) Delaunay.");
+                }
+
                 ShowLog();
             }
         }
diff --git a/Triangle.NET/TestApp/IO/EpsImage.cs b/Triangle.NET/TestApp/IO/EpsImage.cs
index 86450af..4218848 100644
--- a/Triangle.NET/TestApp/IO/EpsImage.cs
+++ b/Triangle.NET/TestApp/IO/EpsImage.cs
@@ -321,10 +321,10 @@ namespace MeshExplorer.IO
 
         private void UpdateMetrics(BoundingBox bounds)
         {
-            x_max = bounds.Xmax;
-            x_min = bounds.Xmin;
-            y_max = bounds.Ymax;
-            y_min = bounds.Ymin;
+            x_max = bounds.MaxX;
+            x_min = bounds.MinX;
+            y_max = bounds.MaxY;
+            y_min = bounds.MinY;
 
             // Enlarge width 5% on each side
             x_scale = x_max - x_min;
diff --git a/Triangle.NET/TestApp/IO/SvgImage.cs b/Triangle.NET/TestApp/IO/SvgImage.cs
index ec285f6..e0d4943 100644
--- a/Triangle.NET/TestApp/IO/SvgImage.cs
+++ b/Triangle.NET/TestApp/IO/SvgImage.cs
@@ -49,8 +49,8 @@ namespace MeshExplorer.IO
 
             scale = width / ((float)bounds.Width + 2 * margin);
 
-            int x_offset = -(int)((bounds.Xmin - margin) * scale);
-            int y_offset = (int)((bounds.Ymax + margin) * scale);
+            int x_offset = -(int)((bounds.MinX - margin) * scale);
+            int y_offset = (int)((bounds.MaxY + margin) * scale);
 
             int height = (int)((bounds.Height + 2 * margin) * scale);
 
diff --git a/Triangle.NET/TestApp/Topology/TopologyRenderControl.cs b/Triangle.NET/TestApp/Topology/TopologyRenderControl.cs
index 78b2bc3..fa59b5a 100644
--- a/Triangle.NET/TestApp/Topology/TopologyRenderControl.cs
+++ b/Triangle.NET/TestApp/Topology/TopologyRenderControl.cs
@@ -48,8 +48,8 @@ namespace MeshExplorer.Topology
             //zoom.ClipMargin = 10.0f;
 
             var b = mesh.Bounds;
-            zoom.Update(new BoundingBox((float)b.Xmin, (float)b.Xmax,
-                (float)b.Ymin, (float)b.Ymax));
+            zoom.Update(new BoundingBox((float)b.MinX, (float)b.MaxX,
+                (float)b.MinY, (float)b.MaxY));
 
             InitializeBuffer();
 
diff --git a/Triangle.NET/Triangle/Algorithm/Incremental.cs b/Triangle.NET/Triangle/Algorithm/Incremental.cs
index 788f560..f77f962 100644
--- a/Triangle.NET/Triangle/Algorithm/Incremental.cs
+++ b/Triangle.NET/Triangle/Algorithm/Incremental.cs
@@ -42,9 +42,9 @@ namespace TriangleNet.Algorithm
                 width = 1.0;
             }
             // Create the vertices of the bounding box.
-            mesh.infvertex1 = new Vertex(box.Xmin - 50.0 * width, box.Ymin - 40.0 * width);
-            mesh.infvertex2 = new Vertex(box.Xmax + 50.0 * width, box.Ymin - 40.0 * width);
-            mesh.infvertex3 = new Vertex(0.5 * (box.Xmin + box.Xmax), box.Ymax + 60.0 * width);
+            mesh.infvertex1 = new Vertex(box.MinX - 50.0 * width, box.MinY - 40.0 * width);
+            mesh.infvertex2 = new Vertex(box.MaxX + 50.0 * width, box.MinY - 40.0 * width);
+            mesh.infvertex3 = new Vertex(0.5 * (box.MinX + box.MaxX), box.MaxY + 60.0 * width);
 
             // Create the bounding box.
             mesh.MakeTriangle(ref inftri);
diff --git a/Triangle.NET/Triangle/Algorithm/SweepLine.cs b/Triangle.NET/Triangle/Algorithm/SweepLine.cs
index 47e58a0..49aa268 100644
--- a/Triangle.NET/Triangle/Algorithm/SweepLine.cs
+++ b/Triangle.NET/Triangle/Algorithm/SweepLine.cs
@@ -525,7 +525,7 @@ namespace TriangleNet.Algorithm
 
             // Nonexistent x value used as a flag to mark circle events in sweepline
             // Delaunay algorithm.
-            xminextreme = 10 * mesh.bounds.Xmin - 9 * mesh.bounds.Xmax;
+            xminextreme = 10 * mesh.bounds.MinX - 9 * mesh.bounds.MaxX;
 
             SweepEvent[] eventheap;
 
@@ -603,7 +603,7 @@ namespace TriangleNet.Algorithm
                 HeapDelete(eventheap, heapsize, 0);
                 heapsize--;
                 check4events = true;
-                if (nextevent.xkey < mesh.bounds.Xmin)
+                if (nextevent.xkey < mesh.bounds.MinX)
                 {
                     fliptri = nextevent.otriEvent;
                     fliptri.Oprev(ref farlefttri);
diff --git a/Triangle.NET/Triangle/Geometry/BoundingBox.cs b/Triangle.NET/Triangle/Geometry/BoundingBox.cs
index b3f30ae..946154e 100644
--- a/Triangle.NET/Triangle/Geometry/BoundingBox.cs
+++ b/Triangle.NET/Triangle/Geometry/BoundingBox.cs
@@ -19,11 +19,13 @@ namespace TriangleNet.Geometry
         /// Initializes a new instance of the <see cref="BoundingBox" /> class.
         /// </summary>
         public BoundingBox()
+            : this(double.MaxValue, double.MaxValue, -double.MaxValue, -double.MaxValue)
+        {
+        }
+
+        public BoundingBox(BoundingBox other)
+            : this(other.MinX, other.MinY, other.MaxX, other.MaxY)
         {
-            xmin = double.MaxValue;
-            ymin = double.MaxValue;
-            xmax = -double.MaxValue;
-            ymax = -double.MaxValue;
         }
 
         /// <summary>
@@ -37,39 +39,39 @@ namespace TriangleNet.Geometry
         public BoundingBox(double xmin, double ymin, double xmax, double ymax)
         {
             this.xmin = xmin;
-            this.ymin = ymin;
             this.xmax = xmax;
+            this.ymin = ymin;
             this.ymax = ymax;
         }
 
         /// <summary>
         /// Gets the minimum x value (left boundary).
         /// </summary>
-        public double Xmin
+        public double MinX
         {
             get { return xmin; }
         }
 
-        /// <summary>
-        /// Gets the minimum y value (bottom boundary).
-        /// </summary>
-        public double Ymin
-        {
-            get { return ymin; }
-        }
-
         /// <summary>
         /// Gets the maximum x value (right boundary).
         /// </summary>
-        public double Xmax
+        public double MaxX
         {
             get { return xmax; }
         }
 
+        /// <summary>
+        /// Gets the minimum y value (bottom boundary).
+        /// </summary>
+        public double MinY
+        {
+            get { return ymin; }
+        }
+
         /// <summary>
         /// Gets the maximum y value (top boundary).
         /// </summary>
-        public double Ymax
+        public double MaxY
         {
             get { return ymax; }
         }
@@ -91,11 +93,24 @@ namespace TriangleNet.Geometry
         }
 
         /// <summary>
-        /// Update bounds.
+        /// Scale bounds.
+        /// </summary>
+        /// <param name="dx">Add dx to left and right bounds.</param>
+        /// <param name="dy">Add dy to top and bottom bounds.</param>
+        public void Resize(double dx, double dy)
+        {
+            xmin -= dx;
+            xmax += dx;
+            ymin -= dy;
+            ymax += dy;
+        }
+
+        /// <summary>
+        /// Expand rectangle to include given point.
         /// </summary>
         /// <param name="x">X coordinate.</param>
         /// <param name="y">Y coordinate.</param>
-        public void Update(double x, double y)
+        public void Expand(double x, double y)
         {
             xmin = Math.Min(xmin, x);
             ymin = Math.Min(ymin, y);
@@ -104,16 +119,16 @@ namespace TriangleNet.Geometry
         }
 
         /// <summary>
-        /// Scale bounds.
+        /// Expand rectangle to include given rectangle.
         /// </summary>
-        /// <param name="dx">Add dx to left and right bounds.</param>
-        /// <param name="dy">Add dy to top and bottom bounds.</param>
-        public void Scale(double dx, double dy)
+        /// <param name="x">X coordinate.</param>
+        /// <param name="y">Y coordinate.</param>
+        public void Expand(BoundingBox other)
         {
-            xmin -= dx;
-            xmax += dx;
-            ymin -= dy;
-            ymax += dy;
+            xmin = Math.Min(xmin, other.xmin);
+            ymin = Math.Min(ymin, other.ymin);
+            xmax = Math.Max(xmax, other.xmax);
+            ymax = Math.Max(ymax, other.ymax);
         }
 
         /// <summary>
@@ -125,5 +140,17 @@ namespace TriangleNet.Geometry
         {
             return ((pt.x >= xmin) && (pt.x <= xmax) && (pt.y >= ymin) && (pt.y <= ymax));
         }
+
+        public bool Contains(BoundingBox other)
+        {
+            return (xmin <= other.MinX && other.MaxX <= xmax
+                && ymin <= other.MinY && other.MaxY <= ymax);
+        }
+
+        public bool Intersects(BoundingBox other)
+        {
+            return (other.MinX < xmax && xmin < other.MaxX
+                && other.MinY < ymax && ymin < other.MaxY);
+        }
     }
 }
diff --git a/Triangle.NET/Triangle/Geometry/InputGeometry.cs b/Triangle.NET/Triangle/Geometry/InputGeometry.cs
index 9063cc8..55b717b 100644
--- a/Triangle.NET/Triangle/Geometry/InputGeometry.cs
+++ b/Triangle.NET/Triangle/Geometry/InputGeometry.cs
@@ -139,7 +139,7 @@ namespace TriangleNet.Geometry
         public void AddPoint(double x, double y, int boundary)
         {
             points.Add(new Vertex(x, y, boundary));
-            bounds.Update(x, y);
+            bounds.Expand(x, y);
         }
 
         /// <summary>
@@ -177,7 +177,7 @@ namespace TriangleNet.Geometry
             }
 
             points.Add(new Vertex(x, y, boundary) { attributes = attribs });
-            bounds.Update(x, y);
+            bounds.Expand(x, y);
         }
 
         /// <summary>
@@ -202,7 +202,7 @@ namespace TriangleNet.Geometry
             }
 
             points.Add(v);
-            bounds.Update(v.x, v.y);
+            bounds.Expand(v.x, v.y);
         }
 
         /// <summary>
diff --git a/Triangle.NET/Triangle/MeshValidator.cs b/Triangle.NET/Triangle/MeshValidator.cs
index eda3ead..ebf8408 100644
--- a/Triangle.NET/Triangle/MeshValidator.cs
+++ b/Triangle.NET/Triangle/MeshValidator.cs
@@ -43,16 +43,21 @@ namespace TriangleNet
                     triorg = tri.Org();
                     tridest = tri.Dest();
                     if (tri.orient == 0)
-                    {   // Only test for inversion once.
+                    {
+                        // Only test for inversion once.
                         // Test if the triangle is flat or inverted.
                         triapex = tri.Apex();
                         if (Primitives.CounterClockwise(triorg, tridest, triapex) <= 0.0)
                         {
-                            logger.Warning("Triangle is flat or inverted.",
-                                "Quality.CheckMesh()");
+                            if (Behavior.Verbose)
+                            {
+                                logger.Warning("Triangle is flat or inverted.", "Quality.CheckMesh()");
+                            }
+
                             horrors++;
                         }
                     }
+
                     // Find the neighboring triangle on this edge.
                     tri.Sym(ref oppotri);
                     if (oppotri.triangle != Mesh.dummytri)
@@ -61,7 +66,7 @@ namespace TriangleNet
                         oppotri.Sym(ref oppooppotri);
                         if ((tri.triangle != oppooppotri.triangle) || (tri.orient != oppooppotri.orient))
                         {
-                            if (tri.triangle == oppooppotri.triangle)
+                            if (tri.triangle == oppooppotri.triangle && Behavior.Verbose)
                             {
                                 logger.Warning("Asymmetric triangle-triangle bond: (Right triangle, wrong orientation)",
                                     "Quality.CheckMesh()");
@@ -75,8 +80,11 @@ namespace TriangleNet
                         oppodest = oppotri.Dest();
                         if ((triorg != oppodest) || (tridest != oppoorg))
                         {
-                            logger.Warning("Mismatched edge coordinates between two triangles.",
-                                "Quality.CheckMesh()");
+                            if (Behavior.Verbose)
+                            {
+                                logger.Warning("Mismatched edge coordinates between two triangles.",
+                                    "Quality.CheckMesh()");
+                            }
 
                             horrors++;
                         }
@@ -88,18 +96,13 @@ namespace TriangleNet
             mesh.MakeVertexMap();
             foreach (var v in mesh.vertices.Values)
             {
-                if (v.tri.triangle == null)
+                if (v.tri.triangle == null && Behavior.Verbose)
                 {
                     logger.Warning("Vertex (ID " + v.id + ") not connected to mesh (duplicate input vertex?)",
                                 "Quality.CheckMesh()");
                 }
             }
 
-            if (horrors == 0) // && Behavior.Verbose
-            {
-                logger.Info("Mesh topology appears to be consistent.");
-            }
-
             // Restore the status of exact arithmetic.
             Behavior.NoExact = saveexact;
 
@@ -107,15 +110,32 @@ namespace TriangleNet
         }
 
         /// <summary>
-        /// Ensure that the mesh is (constrained) Delaunay.
+        /// Check if the mesh is (conforming) Delaunay.
         /// </summary>
         public static bool IsDelaunay(Mesh mesh)
+        {
+            return IsDelaunay(mesh, false);
+        }
+
+        /// <summary>
+        /// Check if that the mesh is (constrained) Delaunay.
+        /// </summary>
+        public static bool IsConstrainedDelaunay(Mesh mesh)
+        {
+            return IsDelaunay(mesh, true);
+        }
+
+        /// <summary>
+        /// Ensure that the mesh is (constrained) Delaunay.
+        /// </summary>
+        private static bool IsDelaunay(Mesh mesh, bool constrained)
         {
             Otri loop = default(Otri);
             Otri oppotri = default(Otri);
             Osub opposubseg = default(Osub);
-            Vertex triorg, tridest, triapex;
+            Vertex org, dest, apex;
             Vertex oppoapex;
+            Vertex inf1, inf2, inf3;
             bool shouldbedelaunay;
             int horrors;
             bool saveexact;
@@ -127,49 +147,57 @@ namespace TriangleNet
             Behavior.NoExact = false;
             horrors = 0;
 
+            inf1 = mesh.infvertex1;
+            inf2 = mesh.infvertex2;
+            inf3 = mesh.infvertex3;
+
             // Run through the list of triangles, checking each one.
             foreach (var tri in mesh.triangles.Values)
             {
                 loop.triangle = tri;
 
                 // Check all three edges of the triangle.
-                for (loop.orient = 0; loop.orient < 3;
-                     loop.orient++)
+                for (loop.orient = 0; loop.orient < 3; loop.orient++)
                 {
-                    triorg = loop.Org();
-                    tridest = loop.Dest();
-                    triapex = loop.Apex();
+                    org = loop.Org();
+                    dest = loop.Dest();
+                    apex = loop.Apex();
+
                     loop.Sym(ref oppotri);
                     oppoapex = oppotri.Apex();
+
                     // Only test that the edge is locally Delaunay if there is an
                     // adjoining triangle whose pointer is larger (to ensure that
                     // each pair isn't tested twice).
-                    shouldbedelaunay = (oppotri.triangle != Mesh.dummytri) &&
-                          !Otri.IsDead(oppotri.triangle) && loop.triangle.id < oppotri.triangle.id &&
-                          (triorg != mesh.infvertex1) && (triorg != mesh.infvertex2) &&
-                          (triorg != mesh.infvertex3) &&
-                          (tridest != mesh.infvertex1) && (tridest != mesh.infvertex2) &&
-                          (tridest != mesh.infvertex3) &&
-                          (triapex != mesh.infvertex1) && (triapex != mesh.infvertex2) &&
-                          (triapex != mesh.infvertex3) &&
-                          (oppoapex != mesh.infvertex1) && (oppoapex != mesh.infvertex2) &&
-                          (oppoapex != mesh.infvertex3);
-                    if (mesh.checksegments && shouldbedelaunay)
+                    shouldbedelaunay = (loop.triangle.id < oppotri.triangle.id) &&
+                           !Otri.IsDead(oppotri.triangle) && (oppotri.triangle != Mesh.dummytri) &&
+                          (org != inf1) && (org != inf2) && (org != inf3) &&
+                          (dest != inf1) && (dest != inf2) && (dest != inf3) &&
+                          (apex != inf1) && (apex != inf2) && (apex != inf3) &&
+                          (oppoapex != inf1) && (oppoapex != inf2) && (oppoapex != inf3);
+
+                    if (constrained && mesh.checksegments && shouldbedelaunay)
                     {
                         // If a subsegment separates the triangles, then the edge is
                         // constrained, so no local Delaunay test should be done.
                         loop.SegPivot(ref opposubseg);
+
                         if (opposubseg.seg != Mesh.dummysub)
                         {
                             shouldbedelaunay = false;
                         }
                     }
+
                     if (shouldbedelaunay)
                     {
-                        if (Primitives.NonRegular(triorg, tridest, triapex, oppoapex) > 0.0)
+                        if (Primitives.NonRegular(org, dest, apex, oppoapex) > 0.0)
                         {
-                            logger.Warning(String.Format("Non-regular pair of triangles found (IDs {0}/{1}).",
-                                loop.triangle.id, oppotri.triangle.id), "Quality.CheckDelaunay()");
+                            if (Behavior.Verbose)
+                            {
+                                logger.Warning(String.Format("Non-regular pair of triangles found (IDs {0}/{1}).",
+                                    loop.triangle.id, oppotri.triangle.id), "Quality.CheckDelaunay()");
+                            }
+
                             horrors++;
                         }
                     }
@@ -177,11 +205,6 @@ namespace TriangleNet
 
             }
 
-            if (horrors == 0) // && Behavior.Verbose
-            {
-                logger.Info("Mesh is Delaunay.");
-            }
-
             // Restore the status of exact arithmetic.
             Behavior.NoExact = saveexact;
 
diff --git a/Triangle.NET/Triangle/Tools/AdjacencyMatrix.cs b/Triangle.NET/Triangle/Tools/AdjacencyMatrix.cs
index f773033..77e23a4 100644
--- a/Triangle.NET/Triangle/Tools/AdjacencyMatrix.cs
+++ b/Triangle.NET/Triangle/Tools/AdjacencyMatrix.cs
@@ -1,5 +1,6 @@
 // -----------------------------------------------------------------------
 // <copyright file="AdjacencyMatrix.cs" company="">
+// Original Matlab code by John Burkardt, Florida State University
 // Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
 // </copyright>
 // -----------------------------------------------------------------------
diff --git a/Triangle.NET/Triangle/Tools/BoundedVoronoi.cs b/Triangle.NET/Triangle/Tools/BoundedVoronoi.cs
index 8ffd857..e851801 100644
--- a/Triangle.NET/Triangle/Tools/BoundedVoronoi.cs
+++ b/Triangle.NET/Triangle/Tools/BoundedVoronoi.cs
@@ -67,7 +67,7 @@ namespace TriangleNet.Tools
         /// <summary>
         /// Gets the list of Voronoi regions.
         /// </summary>
-        public List<VoronoiRegion> Regions
+        public ICollection<VoronoiRegion> Regions
         {
             get { return regions; }
         }
@@ -95,11 +95,11 @@ namespace TriangleNet.Tools
                 // TODO: Need a reliable way to check if a vertex is on a segment
                 if (v.type == VertexType.FreeVertex || v.Boundary == 0)
                 {
-                    ConstructBvdCell(v);
+                    ConstructCell(v);
                 }
                 else if (includeBoundary)
                 {
-                    ConstructBoundaryBvdCell(v);
+                    ConstructBoundaryCell(v);
                 }
             }
 
@@ -262,7 +262,7 @@ namespace TriangleNet.Tools
             return false;
         }
 
-        private void ConstructBvdCell(Vertex vertex)
+        private void ConstructCell(Vertex vertex)
         {
             VoronoiRegion region = new VoronoiRegion(vertex);
             regions.Add(region);
@@ -374,7 +374,7 @@ namespace TriangleNet.Tools
             region.Add(vpoints);
         }
 
-        private void ConstructBoundaryBvdCell(Vertex vertex)
+        private void ConstructBoundaryCell(Vertex vertex)
         {
             VoronoiRegion region = new VoronoiRegion(vertex);
             regions.Add(region);
diff --git a/Triangle.NET/Triangle/Tools/IVoronoi.cs b/Triangle.NET/Triangle/Tools/IVoronoi.cs
index 591efad..3b3c835 100644
--- a/Triangle.NET/Triangle/Tools/IVoronoi.cs
+++ b/Triangle.NET/Triangle/Tools/IVoronoi.cs
@@ -10,7 +10,7 @@ namespace TriangleNet.Tools
     using TriangleNet.Geometry;
 
     /// <summary>
-    /// TODO: Update summary.
+    /// Voronoi diagram interface.
     /// </summary>
     public interface IVoronoi
     {
@@ -22,6 +22,6 @@ namespace TriangleNet.Tools
         /// <summary>
         /// Gets the list of Voronoi regions.
         /// </summary>
-        List<VoronoiRegion> Regions { get; }
+        ICollection<VoronoiRegion> Regions { get; }
     }
 }
diff --git a/Triangle.NET/Triangle/Tools/QuadTree.cs b/Triangle.NET/Triangle/Tools/QuadTree.cs
index 57537b3..58dd60e 100644
--- a/Triangle.NET/Triangle/Tools/QuadTree.cs
+++ b/Triangle.NET/Triangle/Tools/QuadTree.cs
@@ -150,8 +150,8 @@ namespace TriangleNet.Tools
         {
             this.tree = tree;
 
-            this.bounds = new BoundingBox(box.Xmin, box.Ymin, box.Xmax, box.Ymax);
-            this.pivot = new Point((box.Xmin + box.Xmax) / 2, (box.Ymin + box.Ymax) / 2);
+            this.bounds = new BoundingBox(box.MinX, box.MinY, box.MaxX, box.MaxY);
+            this.pivot = new Point((box.MinX + box.MaxX) / 2, (box.MinY + box.MaxY) / 2);
 
             this.bitRegions = 0;
 
@@ -193,19 +193,19 @@ namespace TriangleNet.Tools
             BoundingBox box;
 
             // 1. region south west
-            box = new BoundingBox(bounds.Xmin, bounds.Ymin, pivot.X, pivot.Y);
+            box = new BoundingBox(bounds.MinX, bounds.MinY, pivot.X, pivot.Y);
             regions[0] = new QuadNode(box, tree);
 
             // 2. region south east
-            box = new BoundingBox(pivot.X, bounds.Ymin, bounds.Xmax, pivot.Y);
+            box = new BoundingBox(pivot.X, bounds.MinY, bounds.MaxX, pivot.Y);
             regions[1] = new QuadNode(box, tree);
 
             // 3. region north west
-            box = new BoundingBox(bounds.Xmin, pivot.Y, pivot.X, bounds.Ymax);
+            box = new BoundingBox(bounds.MinX, pivot.Y, pivot.X, bounds.MaxY);
             regions[2] = new QuadNode(box, tree);
 
             // 4. region north east
-            box = new BoundingBox(pivot.X, pivot.Y, bounds.Xmax, bounds.Ymax);
+            box = new BoundingBox(pivot.X, pivot.Y, bounds.MaxX, bounds.MaxY);
             regions[3] = new QuadNode(box, tree);
 
             Point[] triangle = new Point[3];
@@ -295,12 +295,12 @@ namespace TriangleNet.Tools
                 // we have an intersection
                 double yComponent = triangle[k].Y + t * dy;
 
-                if (yComponent < pivot.Y && yComponent >= bounds.Ymin)
+                if (yComponent < pivot.Y && yComponent >= bounds.MinY)
                 {
                     AddToRegion(index, SW);
                     AddToRegion(index, SE);
                 }
-                else if (yComponent <= bounds.Ymax)
+                else if (yComponent <= bounds.MaxY)
                 {
                     AddToRegion(index, NW);
                     AddToRegion(index, NE);
@@ -308,34 +308,34 @@ namespace TriangleNet.Tools
             }
 
             // find intersection with plane x = m_boundingBox[0].dX
-            t = (bounds.Xmin - triangle[k].X) / dx;
+            t = (bounds.MinX - triangle[k].X) / dx;
             if (t < (1 + EPS) && t > -EPS)
             {
                 // we have an intersection
                 double yComponent = triangle[k].Y + t * dy;
 
-                if (yComponent < pivot.Y && yComponent >= bounds.Ymin)
+                if (yComponent < pivot.Y && yComponent >= bounds.MinY)
                 {
                     AddToRegion(index, SW);
                 }
-                else if (yComponent <= bounds.Ymax) // TODO: check && yComponent >= pivot.Y
+                else if (yComponent <= bounds.MaxY) // TODO: check && yComponent >= pivot.Y
                 {
                     AddToRegion(index, NW);
                 }
             }
 
             // find intersection with plane x = m_boundingBox[1].dX
-            t = (bounds.Xmax - triangle[k].X) / dx;
+            t = (bounds.MaxX - triangle[k].X) / dx;
             if (t < (1 + EPS) && t > -EPS)
             {
                 // we have an intersection
                 double yComponent = triangle[k].Y + t * dy;
 
-                if (yComponent < pivot.Y && yComponent >= bounds.Ymin)
+                if (yComponent < pivot.Y && yComponent >= bounds.MinY)
                 {
                     AddToRegion(index, SE);
                 }
-                else if (yComponent <= bounds.Ymax)
+                else if (yComponent <= bounds.MaxY)
                 {
                     AddToRegion(index, NE);
                 }
@@ -353,12 +353,12 @@ namespace TriangleNet.Tools
                 // we have an intersection
                 xComponent = triangle[k].X + t * dx;
 
-                if (xComponent > pivot.X && xComponent <= bounds.Xmax)
+                if (xComponent > pivot.X && xComponent <= bounds.MaxX)
                 {
                     AddToRegion(index, SE);
                     AddToRegion(index, NE);
                 }
-                else if (xComponent >= bounds.Xmin)
+                else if (xComponent >= bounds.MinX)
                 {
                     AddToRegion(index, SW);
                     AddToRegion(index, NW);
@@ -366,34 +366,34 @@ namespace TriangleNet.Tools
             }
 
             // find intersection with plane y = m_boundingBox[0].dY
-            t = (bounds.Ymin - triangle[k].Y) / dy;
+            t = (bounds.MinY - triangle[k].Y) / dy;
             if (t < (1 + EPS) && t > -EPS)
             {
                 // we have an intersection
                 xComponent = triangle[k].X + t * dx;
 
-                if (xComponent > pivot.X && xComponent <= bounds.Xmax)
+                if (xComponent > pivot.X && xComponent <= bounds.MaxX)
                 {
                     AddToRegion(index, SE);
                 }
-                else if (xComponent >= bounds.Xmin)
+                else if (xComponent >= bounds.MinX)
                 {
                     AddToRegion(index, SW);
                 }
             }
 
             // find intersection with plane y = m_boundingBox[1].dY
-            t = (bounds.Ymax - triangle[k].Y) / dy;
+            t = (bounds.MaxY - triangle[k].Y) / dy;
             if (t < (1 + EPS) && t > -EPS)
             {
                 // we have an intersection
                 xComponent = triangle[k].X + t * dx;
 
-                if (xComponent > pivot.X && xComponent <= bounds.Xmax)
+                if (xComponent > pivot.X && xComponent <= bounds.MaxX)
                 {
                     AddToRegion(index, NE);
                 }
-                else if (xComponent >= bounds.Xmin)
+                else if (xComponent >= bounds.MinX)
                 {
                     AddToRegion(index, NW);
                 }
diff --git a/Triangle.NET/Triangle/Tools/Voronoi.cs b/Triangle.NET/Triangle/Tools/Voronoi.cs
index 6f0679f..27a8235 100644
--- a/Triangle.NET/Triangle/Tools/Voronoi.cs
+++ b/Triangle.NET/Triangle/Tools/Voronoi.cs
@@ -20,9 +20,9 @@ namespace TriangleNet.Tools
         Mesh mesh;
 
         Point[] points;
-        List<VoronoiRegion> regions;
+        Dictionary<int, VoronoiRegion> regions;
 
-        // Stores the endpoints of rays of infinite Voronoi cells
+        // Stores the endpoints of rays of unbounded Voronoi cells
         Dictionary<int, Point> rayPoints;
         int rayIndex;
 
@@ -54,9 +54,9 @@ namespace TriangleNet.Tools
         /// <summary>
         /// Gets the list of Voronoi regions.
         /// </summary>
-        public List<VoronoiRegion> Regions
+        public ICollection<VoronoiRegion> Regions
         {
-            get { return regions; }
+            get { return regions.Values; }
         }
 
         /// <summary>
@@ -77,7 +77,7 @@ namespace TriangleNet.Tools
 
             // Allocate space for voronoi diagram
             this.points = new Point[mesh.triangles.Count + mesh.hullsize];
-            this.regions = new List<VoronoiRegion>(mesh.vertices.Count);
+            this.regions = new Dictionary<int, VoronoiRegion>(mesh.vertices.Count);
 
             rayPoints = new Dictionary<int, Point>();
             rayIndex = 0;
@@ -87,12 +87,18 @@ namespace TriangleNet.Tools
             // Compute triangles circumcenters and setup bounding box
             ComputeCircumCenters();
 
+            // Add all Voronoi regions to the map.
+            foreach (var vertex in mesh.vertices.Values)
+            {
+                regions.Add(vertex.id, new VoronoiRegion(vertex));
+            }
+
             // Loop over the mesh vertices (Voronoi generators).
-            foreach (var item in mesh.vertices.Values)
+            foreach (var region in regions.Values)
             {
                 //if (item.Boundary == 0)
                 {
-                    ConstructVoronoiRegion(item);
+                    ConstructVoronoiRegion(region);
                 }
             }
         }
@@ -113,11 +119,11 @@ namespace TriangleNet.Tools
 
                 points[item.id] = pt;
 
-                bounds.Update(pt.x, pt.y);
+                bounds.Expand(pt.x, pt.y);
             }
 
             double ds = Math.Max(bounds.Width, bounds.Height);
-            bounds.Scale(ds, ds);
+            bounds.Resize(ds, ds);
         }
 
         /// <summary>
@@ -125,12 +131,11 @@ namespace TriangleNet.Tools
         /// </summary>
         /// <param name="vertex"></param>
         /// <returns>The circumcenter indices which make up the cell.</returns>
-        private void ConstructVoronoiRegion(Vertex vertex)
+        private void ConstructVoronoiRegion(VoronoiRegion region)
         {
-            VoronoiRegion region = new VoronoiRegion(vertex);
-            regions.Add(region);
+            var vertex = region.Generator as Vertex;
 
-            List<Point> vpoints = new List<Point>();
+            var vpoints = new List<Point>();
 
             Otri f = default(Otri);
             Otri f_init = default(Otri);
@@ -165,6 +170,8 @@ namespace TriangleNet.Tools
                 // Add circumcenter of current triangle
                 vpoints.Add(points[f.triangle.id]);
 
+                region.AddNeighbor(f.triangle.id, regions[f.Apex().id]);
+
                 if (f_next.Equal(f_init))
                 {
                     // Voronoi cell is complete (bounded case).
@@ -179,23 +186,21 @@ namespace TriangleNet.Tools
             // Voronoi cell is unbounded
             region.Bounded = false;
 
-            Vertex torg, tdest, tapex, intersection;
+            Vertex torg, tdest, tapex;
+            Point intersection;
             int sid, n = mesh.triangles.Count;
 
-            // Find the boundary segment id.
+            // Find the boundary segment id (we use this id to number the endpoints of infinit rays).
             f.Lprev(ref f_next);
             f_next.SegPivot(ref sub);
             sid = sub.seg.hash;
 
             // Last valid f lies at the boundary. Add the circumcenter.
             vpoints.Add(points[f.triangle.id]);
+            region.AddNeighbor(f.triangle.id, regions[f.Apex().id]);
 
             // Check if the intersection with the bounding box has already been computed.
-            if (rayPoints.ContainsKey(sid))
-            {
-                vpoints.Add(rayPoints[sid]);
-            }
-            else
+            if (!rayPoints.TryGetValue(sid, out intersection))
             {
                 torg = f.Org();
                 tapex = f.Apex();
@@ -205,13 +210,13 @@ namespace TriangleNet.Tools
                 intersection.id = n + rayIndex;
 
                 points[n + rayIndex] = intersection;
-
                 rayIndex++;
 
-                vpoints.Add(intersection);
                 rayPoints.Add(sid, intersection);
             }
 
+            vpoints.Add(intersection);
+
             // Now walk from f_init clockwise till we reach the boundary.
             vpoints.Reverse();
 
@@ -221,6 +226,7 @@ namespace TriangleNet.Tools
             while (f_prev.triangle != Mesh.dummytri)
             {
                 vpoints.Add(points[f_prev.triangle.id]);
+                region.AddNeighbor(f_prev.triangle.id, regions[f_prev.Apex().id]);
 
                 f_prev.Copy(ref f);
                 f_prev.OprevSelf();
@@ -229,12 +235,8 @@ namespace TriangleNet.Tools
             // Find the boundary segment id.
             f.SegPivot(ref sub);
             sid = sub.seg.hash;
-            
-            if (rayPoints.ContainsKey(sid))
-            {
-                vpoints.Add(rayPoints[sid]);
-            }
-            else
+
+            if (!rayPoints.TryGetValue(sid, out intersection))
             {
                 // Intersection has not been computed yet.
                 torg = f.Org();
@@ -245,20 +247,21 @@ namespace TriangleNet.Tools
                 // Set the correct id for the vertex
                 intersection.id = n + rayIndex;
 
-                points[n + rayIndex] = intersection;
-
-                rayIndex++;
-
-                vpoints.Add(intersection);
                 rayPoints.Add(sid, intersection);
+
+                points[n + rayIndex] = intersection;
+                rayIndex++;
             }
 
+            vpoints.Add(intersection);
+            region.AddNeighbor(intersection.id, regions[f.Dest().id]);
+
             // Add the new points to the region (in counter-clockwise order)
             vpoints.Reverse();
             region.Add(vpoints);
         }
 
-        private bool BoxRayIntersection(Point pt, double dx, double dy, out Vertex intersect)
+        private bool BoxRayIntersection(Point pt, double dx, double dy, out Point intersect)
         {
             double x = pt.X;
             double y = pt.Y;
@@ -266,10 +269,10 @@ namespace TriangleNet.Tools
             double t1, x1, y1, t2, x2, y2;
 
             // Bounding box
-            double minX = bounds.Xmin;
-            double maxX = bounds.Xmax;
-            double minY = bounds.Ymin;
-            double maxY = bounds.Ymax;
+            double minX = bounds.MinX;
+            double maxX = bounds.MaxX;
+            double minY = bounds.MinY;
+            double maxY = bounds.MaxY;
 
             // Check if point is inside the bounds
             if (x < minX || x > maxX || y < minY || y > maxY)
@@ -308,7 +311,7 @@ namespace TriangleNet.Tools
                 x2 = x + t2 * dx;
                 y2 = minY;
             }
-            else if (dx > 0)
+            else if (dy > 0)
             {
                 // Line going upwards: intersect with y = maxY
                 t2 = (maxY - y) / dy;
@@ -324,11 +327,11 @@ namespace TriangleNet.Tools
 
             if (t1 < t2)
             {
-                intersect = new Vertex(x1, y1, -1);
+                intersect = new Point(x1, y1);
             }
             else
             {
-                intersect = new Vertex(x2, y2, -1);
+                intersect = new Point(x2, y2);
             }
 
             return true;
diff --git a/Triangle.NET/Triangle/Tools/VoronoiRegion.cs b/Triangle.NET/Triangle/Tools/VoronoiRegion.cs
index 0d00471..433169a 100644
--- a/Triangle.NET/Triangle/Tools/VoronoiRegion.cs
+++ b/Triangle.NET/Triangle/Tools/VoronoiRegion.cs
@@ -1,6 +1,6 @@
 // -----------------------------------------------------------------------
 // <copyright file="VoronoiRegion.cs" company="">
-// TODO: Update copyright text.
+// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
 // </copyright>
 // -----------------------------------------------------------------------
 
@@ -8,10 +8,8 @@ namespace TriangleNet.Tools
 {
     using System;
     using System.Collections.Generic;
-    using System.Linq;
-    using System.Text;
-    using TriangleNet.Geometry;
     using TriangleNet.Data;
+    using TriangleNet.Geometry;
 
     /// <summary>
     /// Represents a region in the Voronoi diagram.
@@ -23,6 +21,9 @@ namespace TriangleNet.Tools
         List<Point> vertices;
         bool bounded;
 
+        // A map (vertex id) -> (neighbor across adjacent edge)
+        Dictionary<int, VoronoiRegion> neighbors;
+
         /// <summary>
         /// Gets the Voronoi region id (which is the same as the generators vertex id).
         /// </summary>
@@ -62,6 +63,8 @@ namespace TriangleNet.Tools
             this.generator = generator;
             this.vertices = new List<Point>();
             this.bounded = true;
+
+            this.neighbors = new Dictionary<int, VoronoiRegion>();
         }
 
         public void Add(Point point)
@@ -74,6 +77,32 @@ namespace TriangleNet.Tools
             this.vertices.AddRange(points);
         }
 
+        /// <summary>
+        /// Returns the neighbouring Voronoi region, that lies across the edge starting at
+        /// given vertex.
+        /// </summary>
+        /// <param name="p">Vertex defining an edge of the region.</param>
+        /// <returns>Neighbouring Voronoi region</returns>
+        /// <remarks>
+        /// The edge starting at p is well defined (vertices are ordered counterclockwise).
+        /// </remarks>
+        public VoronoiRegion GetNeighbor(Point p)
+        {
+            VoronoiRegion neighbor;
+
+            if (neighbors.TryGetValue(p.id, out neighbor))
+            {
+                return neighbor;
+            }
+
+            return null;
+        }
+
+        internal void AddNeighbor(int id, VoronoiRegion neighbor)
+        {
+            this.neighbors.Add(id, neighbor);
+        }
+
         public override string ToString()
         {
             return String.Format("R-ID {0}", id);