diff --git a/Triangle.NET/TestApp/FormMain.cs b/Triangle.NET/TestApp/FormMain.cs
index e2f5bb3..d466434 100644
--- a/Triangle.NET/TestApp/FormMain.cs
+++ b/Triangle.NET/TestApp/FormMain.cs
@@ -366,7 +366,7 @@ namespace MeshExplorer
             {
                 if (FileProcessor.ContainsMeshData(filename))
                 {
-                    if (DarkMessageBox.Show("Import mesh", Settings.ImportString,
+                    if (filename.EndsWith(".ele") || DarkMessageBox.Show("Import mesh", Settings.ImportString,
                         "Do you want to import the mesh?", MessageBoxButtons.YesNo) == DialogResult.OK)
                     {
                         input = null;
@@ -619,16 +619,35 @@ namespace MeshExplorer
 
             if (mesh.IsPolygon)
             {
-                this.voronoi = new BoundedVoronoi(mesh);
+                try
+                {
+                    this.voronoi = new BoundedVoronoi(mesh);
+                }
+                catch (Exception ex)
+                {
+                    if (!meshControlView.ParamConformDelChecked)
+                    {
+                        DarkMessageBox.Show("Exception - Bounded Voronoi", Settings.VoronoiString, MessageBoxButtons.OK);
+                    }
+                    else
+                    {
+                        DarkMessageBox.Show("Exception - Bounded Voronoi", ex.Message, MessageBoxButtons.OK);
+                    }
+
+                    this.voronoi = null;
+                }
             }
             else
             {
                 this.voronoi = new StandardVoronoi(mesh);
             }
 
-            // HACK: List<Vertex> -> ICollection<Point> ? Nope, no way.
-            //           Vertex[] -> ICollection<Point> ? Well, ok.
-            renderManager.Set(voronoi.Vertices.ToArray(), voronoi.Edges, false);
+            if (this.voronoi != null)
+            {
+                // HACK: List<Vertex> -> ICollection<Point> ? Nope, no way.
+                //           Vertex[] -> ICollection<Point> ? Well, ok.
+                renderManager.Set(voronoi.Vertices.ToArray(), voronoi.Edges, false);
+            }
         }
 
         private void ShowLog()
diff --git a/Triangle.NET/TestApp/Settings.cs b/Triangle.NET/TestApp/Settings.cs
index df8d601..e21d9b0 100644
--- a/Triangle.NET/TestApp/Settings.cs
+++ b/Triangle.NET/TestApp/Settings.cs
@@ -22,6 +22,9 @@ namespace MeshExplorer
         public static string ImportString = "The selected file has associated mesh information. " +
             "You can choose to import the mesh or just read the geometry.";
 
+        public static string VoronoiString = "Make sure you use the \"Confoming Delaunay\" option " +
+            "when building the Voronoi diagram from a constrained mesh.";
+
         // Open file dialog
         public string OfdDirectory { get; set; }
         public string OfdFilter { get; set; }
diff --git a/Triangle.NET/Triangle/Geometry/Edge.cs b/Triangle.NET/Triangle/Geometry/Edge.cs
index 25ae78c..f8daeb6 100644
--- a/Triangle.NET/Triangle/Geometry/Edge.cs
+++ b/Triangle.NET/Triangle/Geometry/Edge.cs
@@ -6,12 +6,6 @@
 
 namespace TriangleNet.Geometry
 {
-    using System;
-    using System.Collections.Generic;
-    using System.Linq;
-    using System.Text;
-    using TriangleNet.Topology;
-
     /// <summary>
     /// Represents a straight line segment in 2D space.
     /// </summary>
@@ -54,11 +48,11 @@ namespace TriangleNet.Geometry
         /// <summary>
         /// Initializes a new instance of the <see cref="Edge" /> class.
         /// </summary>
-        public Edge(int p0, int p1, int boundary)
+        public Edge(int p0, int p1, int label)
         {
             this.P0 = p0;
             this.P1 = p1;
-            this.Label = boundary;
+            this.Label = label;
         }
     }
 }
diff --git a/Triangle.NET/Triangle/Geometry/IEdge.cs b/Triangle.NET/Triangle/Geometry/IEdge.cs
index 66bca45..a7edf5f 100644
--- a/Triangle.NET/Triangle/Geometry/IEdge.cs
+++ b/Triangle.NET/Triangle/Geometry/IEdge.cs
@@ -21,6 +21,9 @@ namespace TriangleNet.Geometry
         /// <summary>
         /// Gets or sets a general-purpose label.
         /// </summary>
+        /// <remarks>
+        /// This is used for the segments boundary mark.
+        /// </remarks>
         int Label { get; }
     }
 }
diff --git a/Triangle.NET/Triangle/Geometry/ISegment.cs b/Triangle.NET/Triangle/Geometry/ISegment.cs
index 86a7d48..5a99983 100644
--- a/Triangle.NET/Triangle/Geometry/ISegment.cs
+++ b/Triangle.NET/Triangle/Geometry/ISegment.cs
@@ -12,9 +12,9 @@ namespace TriangleNet.Geometry
     public interface ISegment : IEdge
     {
         /// <summary>
-        /// Gets the segments endpoint.
+        /// Gets the vertex at given index.
         /// </summary>
-        /// <param name="index">The vertex index (0 or 1).</param>
+        /// <param name="index">The local index (0 or 1).</param>
         Vertex GetVertex(int index);
 
         /// <summary>
diff --git a/Triangle.NET/Triangle/Geometry/Point.cs b/Triangle.NET/Triangle/Geometry/Point.cs
index 0ad4eb1..f71061e 100644
--- a/Triangle.NET/Triangle/Geometry/Point.cs
+++ b/Triangle.NET/Triangle/Geometry/Point.cs
@@ -14,9 +14,13 @@ namespace TriangleNet.Geometry
     public class Point : IComparable<Point>, IEquatable<Point>
     {
         internal int id;
+        internal int label;
+
         internal double x;
         internal double y;
-        internal int label;
+#if USE_Z
+        internal double z;
+#endif
 
         public Point()
             : this(0, 0, 0)
@@ -64,6 +68,17 @@ namespace TriangleNet.Geometry
             set { this.y = value; }
         }
 
+#if USE_Z
+        /// <summary>
+        /// Gets or sets the vertex z coordinate.
+        /// </summary>
+        public double Z
+        {
+            get { return this.z; }
+            set { this.z = value; }
+        }
+#endif
+
         /// <summary>
         /// Gets or sets a general-purpose label.
         /// </summary>
diff --git a/Triangle.NET/Triangle/Geometry/Vertex.cs b/Triangle.NET/Triangle/Geometry/Vertex.cs
index 43325f0..49f1b58 100644
--- a/Triangle.NET/Triangle/Geometry/Vertex.cs
+++ b/Triangle.NET/Triangle/Geometry/Vertex.cs
@@ -18,8 +18,9 @@ namespace TriangleNet.Geometry
         // Hash for dictionary. Will be set by mesh instance.
         internal int hash;
 
+#if USE_ATTRIBS
         internal double[] attributes;
-
+#endif
         internal VertexType type;
         internal Otri tri;
 
@@ -27,7 +28,7 @@ namespace TriangleNet.Geometry
         /// Initializes a new instance of the <see cref="Vertex" /> class.
         /// </summary>
         public Vertex()
-            : this(0, 0, 0, 0)
+            : this(0, 0, 0)
         {
         }
 
@@ -37,7 +38,7 @@ namespace TriangleNet.Geometry
         /// <param name="x">The x coordinate of the vertex.</param>
         /// <param name="y">The y coordinate of the vertex.</param>
         public Vertex(double x, double y)
-            : this(x, y, 0, 0)
+            : this(x, y, 0)
         {
         }
 
@@ -48,10 +49,12 @@ namespace TriangleNet.Geometry
         /// <param name="y">The y coordinate of the vertex.</param>
         /// <param name="mark">The boundary mark.</param>
         public Vertex(double x, double y, int mark)
-            : this(x, y, mark, 0)
+            : base(x, y, mark)
         {
+            this.type = VertexType.InputVertex;
         }
 
+#if USE_ATTRIBS
         /// <summary>
         /// Initializes a new instance of the <see cref="Vertex" /> class.
         /// </summary>
@@ -60,18 +63,18 @@ namespace TriangleNet.Geometry
         /// <param name="mark">The boundary mark.</param>
         /// <param name="attribs">The number of point attributes.</param>
         public Vertex(double x, double y, int mark, int attribs)
-            : base(x, y, mark)
+            : this(x, y, mark)
         {
-            this.type = VertexType.InputVertex;
-
             if (attribs > 0)
             {
                 this.attributes = new double[attribs];
             }
         }
+#endif
 
         #region Public properties
 
+#if USE_ATTRIBS
         /// <summary>
         /// Gets the vertex attributes (may be null).
         /// </summary>
@@ -79,6 +82,7 @@ namespace TriangleNet.Geometry
         {
             get { return this.attributes; }
         }
+#endif
 
         /// <summary>
         /// Gets the vertex type.
diff --git a/Triangle.NET/Triangle/IO/TriangleReader.cs b/Triangle.NET/Triangle/IO/TriangleReader.cs
index e3d6587..43a1dec 100644
--- a/Triangle.NET/Triangle/IO/TriangleReader.cs
+++ b/Triangle.NET/Triangle/IO/TriangleReader.cs
@@ -8,12 +8,10 @@
 namespace TriangleNet.IO
 {
     using System;
-    using System.IO;
-    using System.Globalization;
-    using TriangleNet.Topology;
-    using TriangleNet.Logging;
-    using TriangleNet.Geometry;
     using System.Collections.Generic;
+    using System.Globalization;
+    using System.IO;
+    using TriangleNet.Geometry;
 
     /// <summary>
     /// Helper methods for reading Triangle file formats.
@@ -74,6 +72,7 @@ namespace TriangleNet.IO
 
             if (attributes > 0)
             {
+#if USE_ATTRIBS
                 var attribs = new double[attributes];
 
                 // Read the vertex attributes.
@@ -86,6 +85,7 @@ namespace TriangleNet.IO
                 }
 
                 v.attributes = attribs;
+#endif
             }
 
             data.Add(v);
diff --git a/Triangle.NET/Triangle/IO/TriangleWriter.cs b/Triangle.NET/Triangle/IO/TriangleWriter.cs
index 4c07bae..976ab9c 100644
--- a/Triangle.NET/Triangle/IO/TriangleWriter.cs
+++ b/Triangle.NET/Triangle/IO/TriangleWriter.cs
@@ -7,12 +7,11 @@
 
 namespace TriangleNet.IO
 {
-    using System;
-    using System.IO;
-    using System.Globalization;
-    using TriangleNet.Topology;
-    using TriangleNet.Geometry;
     using System.Collections.Generic;
+    using System.Globalization;
+    using System.IO;
+    using TriangleNet.Geometry;
+    using TriangleNet.Topology;
 
     /// <summary>
     /// Helper methods for writing Triangle file formats.
@@ -51,6 +50,7 @@ namespace TriangleNet.IO
         private static void WriteNodes(StreamWriter writer, Mesh mesh)
         {
             int outvertices = mesh.vertices.Count;
+            int nextras = mesh.nextras;
 
             Behavior behavior = mesh.behavior;
 
@@ -63,7 +63,7 @@ namespace TriangleNet.IO
             {
                 // Number of vertices, number of dimensions, number of vertex attributes,
                 // and number of boundary markers (zero or one).
-                writer.WriteLine("{0} {1} {2} {3}", outvertices, mesh.mesh_dim, mesh.nextras,
+                writer.WriteLine("{0} {1} {2} {3}", outvertices, mesh.mesh_dim, nextras,
                     behavior.UseBoundaryMarkers ? "1" : "0");
 
                 if (mesh.numbering == NodeNumbering.None)
@@ -76,7 +76,7 @@ namespace TriangleNet.IO
                 {
                     // If numbering is linear, just use the dictionary values.
                     WriteNodes(writer, mesh.vertices.Values, behavior.UseBoundaryMarkers,
-                        mesh.nextras, behavior.Jettison);
+                        nextras, behavior.Jettison);
                 }
                 else
                 {
@@ -92,7 +92,7 @@ namespace TriangleNet.IO
                     }
 
                     WriteNodes(writer, nodes, behavior.UseBoundaryMarkers,
-                        mesh.nextras, behavior.Jettison);
+                        nextras, behavior.Jettison);
                 }
             }
         }
@@ -114,11 +114,13 @@ namespace TriangleNet.IO
                     // Vertex number, x and y coordinates.
                     writer.Write("{0} {1} {2}", index, vertex.x.ToString(nfi), vertex.y.ToString(nfi));
 
+#if USE_ATTRIBS
                     // Write attributes.
                     for (int j = 0; j < attribs; j++)
                     {
                         writer.Write(" {0}", vertex.attributes[j].ToString(nfi));
                     }
+#endif
 
                     if (markers)
                     {
diff --git a/Triangle.NET/Triangle/Mesh.cs b/Triangle.NET/Triangle/Mesh.cs
index 800f7a5..f4b0808 100644
--- a/Triangle.NET/Triangle/Mesh.cs
+++ b/Triangle.NET/Triangle/Mesh.cs
@@ -48,10 +48,11 @@ namespace TriangleNet
         internal List<Point> holes;
         internal List<RegionPointer> regions;
 
+        // TODO: remove mesh_dim, invertices and insegments
+
         // Other variables.
         internal Rectangle bounds; // x and y bounds.
         internal int invertices;     // Number of input vertices.
-        internal int inelements;     // Number of input triangles.
         internal int insegments;     // Number of input segments.
         internal int undeads;        // Number of input vertices that don't appear in the mesh.
         internal int mesh_dim;       // Dimension (ought to be 2).
@@ -257,7 +258,6 @@ namespace TriangleNet
 
         public void Refine(QualityOptions quality)
         {
-            inelements = triangles.Count;
             invertices = vertices.Count;
 
             if (behavior.Poly)
@@ -416,8 +416,10 @@ namespace TriangleNet
 
             var v = points[0];
 
+#if USE_ATTRIBS
             // Check attributes.
             this.nextras = v.attributes == null ? 0 : v.attributes.Length;
+#endif
 
             // Simple heuristic to check if ids are already set.  We assume that if the
             // first two vertex ids are distinct, then all input vertices have pairwise
diff --git a/Triangle.NET/Triangle/Meshing/Algorithm/SweepLine.cs b/Triangle.NET/Triangle/Meshing/Algorithm/SweepLine.cs
index e979467..31c1ba3 100644
--- a/Triangle.NET/Triangle/Meshing/Algorithm/SweepLine.cs
+++ b/Triangle.NET/Triangle/Meshing/Algorithm/SweepLine.cs
@@ -76,8 +76,8 @@ namespace TriangleNet.Meshing.Algorithm
             splaynodes = new List<SplayNode>();
             splayroot = null;
 
-            CreateHeap(out eventheap);//, out events, out freeevents);
-            heapsize = mesh.invertices;
+            heapsize = points.Count;
+            CreateHeap(out eventheap, heapsize);//, out events, out freeevents);
 
             mesh.MakeTriangle(ref lefttri);
             mesh.MakeTriangle(ref righttri);
@@ -387,14 +387,14 @@ namespace TriangleNet.Meshing.Algorithm
             Heapify(heap, heapsize - 1, eventnum);
         }
 
-        void CreateHeap(out SweepEvent[] eventheap)
+        void CreateHeap(out SweepEvent[] eventheap, int size)
         {
             Vertex thisvertex;
             int maxevents;
             int i;
             SweepEvent evt;
 
-            maxevents = (3 * mesh.invertices) / 2;
+            maxevents = (3 * size) / 2;
             eventheap = new SweepEvent[maxevents];
 
             i = 0;
diff --git a/Triangle.NET/Triangle/Meshing/ConstraintMesher.cs b/Triangle.NET/Triangle/Meshing/ConstraintMesher.cs
index 8aa944b..c13d98d 100644
--- a/Triangle.NET/Triangle/Meshing/ConstraintMesher.cs
+++ b/Triangle.NET/Triangle/Meshing/ConstraintMesher.cs
@@ -710,17 +710,22 @@ namespace TriangleNet.Meshing
             newvertex = new Vertex(
                 torg.x + split * (tdest.x - torg.x),
                 torg.y + split * (tdest.y - torg.y),
-                splitsubseg.seg.boundary,
-                mesh.nextras);
+                splitsubseg.seg.boundary
+#if USE_ATTRIBS
+                , mesh.nextras
+#endif
+                );
 
             newvertex.hash = mesh.hash_vtx++;
             newvertex.id = newvertex.hash;
 
+#if USE_ATTRIBS
             // Interpolate its attributes.
             for (int i = 0; i < mesh.nextras; i++)
             {
                 newvertex.attributes[i] = torg.attributes[i] + split * (tdest.attributes[i] - torg.attributes[i]);
             }
+#endif
 
             mesh.vertices.Add(newvertex.hash, newvertex);
 
diff --git a/Triangle.NET/Triangle/Meshing/Converter.cs b/Triangle.NET/Triangle/Meshing/Converter.cs
index ce5f5c7..af7b024 100644
--- a/Triangle.NET/Triangle/Meshing/Converter.cs
+++ b/Triangle.NET/Triangle/Meshing/Converter.cs
@@ -48,7 +48,6 @@ namespace TriangleNet.Meshing
 
             mesh.TransferNodes(polygon.Points);
 
-            mesh.inelements = elements;
             mesh.regions.AddRange(polygon.Regions);
             mesh.behavior.useRegions = polygon.Regions.Count > 0;
 
diff --git a/Triangle.NET/Triangle/Meshing/QualityMesher.cs b/Triangle.NET/Triangle/Meshing/QualityMesher.cs
index 4cb967b..3d39216 100644
--- a/Triangle.NET/Triangle/Meshing/QualityMesher.cs
+++ b/Triangle.NET/Triangle/Meshing/QualityMesher.cs
@@ -12,6 +12,7 @@ namespace TriangleNet.Meshing
     using TriangleNet.Geometry;
     using TriangleNet.Logging;
     using TriangleNet.Meshing.Data;
+    using Tools;
     using TriangleNet.Topology;
 
     /// <summary>
@@ -600,8 +601,11 @@ namespace TriangleNet.Meshing
                     newvertex = new Vertex(
                         eorg.x + split * (edest.x - eorg.x),
                         eorg.y + split * (edest.y - eorg.y),
-                        currentenc.seg.boundary,
-                        mesh.nextras);
+                        currentenc.seg.boundary
+#if USE_ATTRIBS
+                        , mesh.nextras
+#endif
+                    );
 
                     newvertex.type = VertexType.SegmentVertex;
 
@@ -609,14 +613,14 @@ namespace TriangleNet.Meshing
                     newvertex.id = newvertex.hash;
 
                     mesh.vertices.Add(newvertex.hash, newvertex);
-
+#if USE_ATTRIBS
                     // Interpolate attributes.
                     for (int i = 0; i < mesh.nextras; i++)
                     {
                         newvertex.attributes[i] = eorg.attributes[i]
                             + split * (edest.attributes[i] - eorg.attributes[i]);
                     }
-
+#endif
                     if (!Behavior.NoExact)
                     {
                         // Roundoff in the above calculation may yield a 'newvertex'
@@ -689,54 +693,6 @@ namespace TriangleNet.Meshing
             }
         }
 
-        /// <summary>
-        /// Linear interpolation of vertex attributes.
-        /// </summary>
-        /// <param name="newvertex"></param>
-        private void InterpolateAttribs(Vertex newvertex)
-        {
-            Vertex org = newvertex_tri.vertices[0];
-            Vertex dest = newvertex_tri.vertices[1];
-            Vertex apex = newvertex_tri.vertices[2];
-
-            double xdo, ydo, xao, yao;
-            double denominator;
-            double dx, dy;
-            double xi, eta;
-
-            int nextras = mesh.nextras;
-
-            // Compute the circumcenter of the triangle.
-            xdo = dest.x - org.x;
-            ydo = dest.y - org.y;
-            xao = apex.x - org.x;
-            yao = apex.y - org.y;
-
-            denominator = 0.5 / (xdo * yao - xao * ydo);
-
-            //dx = (yao * dodist - ydo * aodist) * denominator;
-            //dy = (xdo * aodist - xao * dodist) * denominator;
-
-            dx = newvertex.x - org.x;
-            dy = newvertex.y - org.y;
-
-            // To interpolate vertex attributes for the new vertex inserted at
-            // the circumcenter, define a coordinate system with a xi-axis,
-            // directed from the triangle's origin to its destination, and
-            // an eta-axis, directed from its origin to its apex.
-            // Calculate the xi and eta coordinates of the circumcenter.
-            xi = (yao * dx - xao * dy) * (2.0 * denominator);
-            eta = (xdo * dy - ydo * dx) * (2.0 * denominator);
-
-            for (int i = 0; i < nextras; i++)
-            {
-                // Interpolate the vertex attributes.
-                newvertex.attributes[i] = org.attributes[i]
-                    + xi * (dest.attributes[i] - org.attributes[i])
-                    + eta * (apex.attributes[i] - org.attributes[i]);
-            }
-        }
-
         /// <summary>
         /// Inserts a vertex at the circumcenter of a triangle. Deletes 
         /// the newly inserted vertex if it encroaches upon a segment.
@@ -793,7 +749,12 @@ namespace TriangleNet.Meshing
                 {
                     // The new vertex must be in the interior, and therefore is a
                     // free vertex with a marker of zero.
-                    Vertex newvertex = new Vertex(newloc.x, newloc.y, 0, mesh.nextras);
+                    Vertex newvertex = new Vertex(newloc.x, newloc.y, 0
+#if USE_ATTRIBS
+                        , mesh.nextras
+#endif
+                        );
+
                     newvertex.type = VertexType.FreeVertex;
 
                     // Ensure that the handle 'badotri' does not represent the longest
@@ -820,12 +781,12 @@ namespace TriangleNet.Meshing
                     {
                         newvertex.hash = mesh.hash_vtx++;
                         newvertex.id = newvertex.hash;
-
+#if USE_ATTRIBS
                         if (mesh.nextras > 0)
                         {
-                            InterpolateAttribs(newvertex);
+                            Interpolation.InterpolateAttributes(newvertex, newvertex.tri.tri, mesh.nextras);
                         }
-
+#endif
                         mesh.vertices.Add(newvertex.hash, newvertex);
 
                         if (mesh.steinerleft > 0)
diff --git a/Triangle.NET/Triangle/Tools/Interpolation.cs b/Triangle.NET/Triangle/Tools/Interpolation.cs
new file mode 100644
index 0000000..6fa6410
--- /dev/null
+++ b/Triangle.NET/Triangle/Tools/Interpolation.cs
@@ -0,0 +1,106 @@
+
+namespace TriangleNet.Tools
+{
+    using TriangleNet.Geometry;
+
+    public static class Interpolation
+    {
+#if USE_ATTRIBS
+        /// <summary>
+        /// Linear interpolation of vertex attributes.
+        /// </summary>
+        /// <param name="vertex">The interpolation vertex.</param>
+        /// <param name="triangle">The triangle containing the vertex.</param>
+        /// <param name="n">The number of vertex attributes.</param>
+        /// <remarks>
+        /// The vertex is expected to lie inside the triangle.
+        /// </remarks>
+        public static void InterpolateAttributes(Vertex vertex, ITriangle triangle, int n)
+        {
+            Vertex org = triangle.GetVertex(0);
+            Vertex dest = triangle.GetVertex(1);
+            Vertex apex = triangle.GetVertex(2);
+
+            double xdo, ydo, xao, yao;
+            double denominator;
+            double dx, dy;
+            double xi, eta;
+
+            // Compute the circumcenter of the triangle.
+            xdo = dest.x - org.x;
+            ydo = dest.y - org.y;
+            xao = apex.x - org.x;
+            yao = apex.y - org.y;
+
+            denominator = 0.5 / (xdo * yao - xao * ydo);
+
+            //dx = (yao * dodist - ydo * aodist) * denominator;
+            //dy = (xdo * aodist - xao * dodist) * denominator;
+
+            dx = vertex.x - org.x;
+            dy = vertex.y - org.y;
+
+            // To interpolate vertex attributes for the new vertex inserted at
+            // the circumcenter, define a coordinate system with a xi-axis,
+            // directed from the triangle's origin to its destination, and
+            // an eta-axis, directed from its origin to its apex.
+            // Calculate the xi and eta coordinates of the circumcenter.
+            xi = (yao * dx - xao * dy) * (2.0 * denominator);
+            eta = (xdo * dy - ydo * dx) * (2.0 * denominator);
+        
+            for (int i = 0; i < n; i++)
+            {
+                // Interpolate the vertex attributes.
+                vertex.attributes[i] = org.attributes[i]
+                    + xi * (dest.attributes[i] - org.attributes[i])
+                    + eta * (apex.attributes[i] - org.attributes[i]);
+            }
+        }
+#endif
+
+#if USE_Z
+        /// <summary>
+        /// Linear interpolation of a scalar value.
+        /// </summary>
+        /// <param name="p">The interpolation point.</param>
+        /// <param name="triangle">The triangle containing the point.</param>
+        /// <remarks>
+        /// The point is expected to lie inside the triangle.
+        /// </remarks>
+        public static void InterpolateZ(Point p, ITriangle triangle)
+        {
+            Vertex org = triangle.GetVertex(0);
+            Vertex dest = triangle.GetVertex(1);
+            Vertex apex = triangle.GetVertex(2);
+
+            double xdo, ydo, xao, yao;
+            double denominator;
+            double dx, dy;
+            double xi, eta;
+
+            // Compute the circumcenter of the triangle.
+            xdo = dest.x - org.x;
+            ydo = dest.y - org.y;
+            xao = apex.x - org.x;
+            yao = apex.y - org.y;
+
+            denominator = 0.5 / (xdo * yao - xao * ydo);
+
+            //dx = (yao * dodist - ydo * aodist) * denominator;
+            //dy = (xdo * aodist - xao * dodist) * denominator;
+
+            dx = p.x - org.x;
+            dy = p.y - org.y;
+
+            // To interpolate z value for the given point inserted, define a
+            // coordinate system with a xi-axis, directed from the triangle's
+            // origin to its destination, and an eta-axis, directed from its
+            // origin to its apex.
+            xi = (yao * dx - xao * dy) * (2.0 * denominator);
+            eta = (xdo * dy - ydo * dx) * (2.0 * denominator);
+
+            p.z = org.z + xi * (dest.z - org.z) + eta * (apex.z - org.z);
+        }
+#endif
+    }
+}
diff --git a/Triangle.NET/Triangle/Triangle.csproj b/Triangle.NET/Triangle/Triangle.csproj
index 60aa0cf..d451d4b 100644
--- a/Triangle.NET/Triangle/Triangle.csproj
+++ b/Triangle.NET/Triangle/Triangle.csproj
@@ -23,7 +23,7 @@
     <DebugType>full</DebugType>
     <Optimize>false</Optimize>
     <OutputPath>bin\Debug\</OutputPath>
-    <DefineConstants>DEBUG;TRACE</DefineConstants>
+    <DefineConstants>TRACE;DEBUG;USE_ATTRIBS</DefineConstants>
     <ErrorReport>prompt</ErrorReport>
     <WarningLevel>4</WarningLevel>
   </PropertyGroup>
@@ -31,7 +31,7 @@
     <DebugType>pdbonly</DebugType>
     <Optimize>true</Optimize>
     <OutputPath>bin\Release\</OutputPath>
-    <DefineConstants>TRACE</DefineConstants>
+    <DefineConstants>TRACE;USE_ATTRIBS</DefineConstants>
     <ErrorReport>prompt</ErrorReport>
     <WarningLevel>4</WarningLevel>
   </PropertyGroup>
@@ -84,6 +84,7 @@
     <Compile Include="Smoothing\VoronoiFactory.cs" />
     <Compile Include="Tools\AdjacencyMatrix.cs" />
     <Compile Include="Tools\CuthillMcKee.cs" />
+    <Compile Include="Tools\Interpolation.cs" />
     <Compile Include="Tools\IntersectionHelper.cs" />
     <Compile Include="Tools\PointSorter.cs" />
     <Compile Include="Tools\PolygonValidator.cs" />