diff --git a/example/voronoi_basic_tutorial.cpp b/example/voronoi_basic_tutorial.cpp
index 41497be..69ebccd 100644
--- a/example/voronoi_basic_tutorial.cpp
+++ b/example/voronoi_basic_tutorial.cpp
@@ -11,9 +11,10 @@
 #include <vector>
 
 #include <boost/polygon/voronoi.hpp>
-#include <boost/polygon/voronoi_utils.hpp>
 using namespace boost::polygon;
 
+#include "voronoi_visual_utils.hpp"
+
 struct Point {
   int a;
   int b;
@@ -57,7 +58,7 @@ struct segment_traits<Segment> {
 }  // boost
 
 // Traversing Voronoi edges using edge iterator.
-int iterate_primary_edges1(const voronoi_diagram<double> &vd) {
+int iterate_primary_edges1(const voronoi_diagram<double>& vd) {
   int result = 0;
   for (voronoi_diagram<double>::const_edge_iterator it = vd.edges().begin();
        it != vd.edges().end(); ++it) {
@@ -72,8 +73,8 @@ int iterate_primary_edges2(const voronoi_diagram<double> &vd) {
   int result = 0;
   for (voronoi_diagram<double>::const_cell_iterator it = vd.cells().begin();
        it != vd.cells().end(); ++it) {
-    const voronoi_diagram<double>::cell_type &cell = *it;
-    const voronoi_diagram<double>::edge_type *edge = cell.incident_edge();
+    const voronoi_diagram<double>::cell_type& cell = *it;
+    const voronoi_diagram<double>::edge_type* edge = cell.incident_edge();
     // This is convenient way to iterate edges around Voronoi cell.
     do {
       if (edge->is_primary())
@@ -92,8 +93,8 @@ int iterate_primary_edges3(const voronoi_diagram<double> &vd) {
   int result = 0;
   for (voronoi_diagram<double>::const_vertex_iterator it = vd.vertices().begin();
        it != vd.vertices().end(); ++it) {
-    const voronoi_diagram<double>::vertex_type &vertex = *it;
-    const voronoi_diagram<double>::edge_type *edge = vertex.incident_edge();
+    const voronoi_diagram<double>::vertex_type& vertex = *it;
+    const voronoi_diagram<double>::edge_type* edge = vertex.incident_edge();
     // This is convenient way to iterate edges around Voronoi vertex.
     do {
       if (edge->is_primary())
@@ -104,35 +105,6 @@ int iterate_primary_edges3(const voronoi_diagram<double> &vd) {
   return result;
 }
 
-// Prototype of a function that renders segments.
-void draw_segment(double x1, double y1, double x2, double y2) {
-  printf("Rendering segment: ");
-  printf("%7.3f %7.3f %7.3f %7.3f\n", x1, y1, x2, y2);
-}
-
-void render_diagram(const voronoi_diagram<double> &vd,
-                    const voronoi_utils<double>::brect_type &bbox) {
-  const std::size_t VISITED = 1;
-  for (voronoi_diagram<double>::const_edge_iterator it = vd.edges().begin();
-       it != vd.edges().end(); ++it) {
-    // We use color to mark visited edges.
-    it->color(VISITED);
-    // Don't render the same edge twice.
-    if (it->twin()->color()) continue;
-    voronoi_utils<double>::point_set_type polyline;
-    if (it->is_linear())
-      voronoi_utils<double>::clip(*it, bbox, polyline);
-    else
-      // Parabolic edges are always finite.
-      voronoi_utils<double>::discretize(*it, 1E-1, polyline);
-    // Note: discretized edges may also lie outside of the bbox.
-    // So user might do additional clipping before rendering each such edge.  
-    for (std::size_t i = 1; i < polyline.size(); ++i)
-      draw_segment(polyline[i-1].x(), polyline[i-1].y(),
-                   polyline[i].x(), polyline[i].y());
-  }
-}
-
 int main() {
   // Preparing Input Geometries.
   std::vector<Point> points;
@@ -155,46 +127,21 @@ int main() {
     printf("\n");
   }
 
-  // Using color member of the Voronoi primitives..
+  // Using color member of the Voronoi primitives to store the average number of edges
+  // around each cell (including secondary edges).
   {
+    printf("Number of edges (including secondary edges) around the Voronoi cells:\n");
+    for (voronoi_diagram<double>::const_edge_iterator it = vd.edges().begin();
+         it != vd.edges().end(); ++it) {
+      std::size_t cnt = it->cell()->color();
+      it->cell()->color(cnt + 1);
+    }
     for (voronoi_diagram<double>::const_cell_iterator it = vd.cells().begin();
          it != vd.cells().end(); ++it) {
-      const voronoi_diagram<double>::cell_type &cell = *it;
-      const voronoi_diagram<double>::edge_type *edge = cell.incident_edge();
-      std::size_t count = 0;
-      do {
-        ++count;
-        edge = edge->next();
-      } while (edge != cell.incident_edge());
-      cell.color(count);
+      printf("%d ", it->color());
     }
-
-    // Count the average number of edges.
-    double total = 0;
-    for (voronoi_diagram<double>::const_cell_iterator it = vd.cells().begin();
-         it != vd.cells().end(); ++it) {
-      total += it->color();
-    }
-    total /= vd.cells().size();
-    printf("The average number of edges per Voronoi cell is equal to: %3.1f\n", total);
+    printf("\n");
     printf("\n");
   } 
-
-  // Rendering Voronoi diagram.
-  {
-    // Construct clipping bounding rectangle.
-    bounding_rectangle<double> bbox;
-    for (std::vector<Point>::iterator it = points.begin(); it != points.end(); ++it)
-      bbox.update(it->a, it->b);
-    for (std::vector<Segment>::iterator it = segments.begin(); it != segments.end(); ++it) {
-      bbox.update(it->p0.a, it->p0.b);
-      bbox.update(it->p1.a, it->p1.b);
-    }
-    // Add 10% offset to the bounding rectangle.
-    bbox = voronoi_utils<double>::scale(bbox, 1.1);
-    // Render Voronoi diagram.
-    render_diagram(vd, bbox);
-  }
-
   return 0;
 }
diff --git a/include/boost/polygon/detail/voronoi_predicates.hpp b/include/boost/polygon/detail/voronoi_predicates.hpp
index 988bd4e..c505db5 100644
--- a/include/boost/polygon/detail/voronoi_predicates.hpp
+++ b/include/boost/polygon/detail/voronoi_predicates.hpp
@@ -37,12 +37,12 @@ public:
   };
 
   template <typename Point>
-  static bool is_vertical(const Point &point1, const Point &point2) {
+  static bool is_vertical(const Point& point1, const Point& point2) {
     return point1.x() == point2.x();
   }
 
   template <typename Site>
-  static bool is_vertical(const Site &site) {
+  static bool is_vertical(const Site& site) {
     return is_vertical(site.point0(), site.point1());
   }
 
@@ -100,9 +100,9 @@ public:
     }
 
     template <typename Point>
-    static Orientation eval(const Point &point1,
-                            const Point &point2,
-                            const Point &point3) {
+    static Orientation eval(const Point& point1,
+                            const Point& point2,
+                            const Point& point3) {
       int_x2_type dx1 = static_cast<int_x2_type>(point1.x()) -
                         static_cast<int_x2_type>(point2.x());
       int_x2_type dx2 = static_cast<int_x2_type>(point2.x()) -
@@ -120,7 +120,7 @@ public:
   public:
     typedef Point point_type;
 
-    bool operator()(const point_type &lhs, const point_type &rhs) const {
+    bool operator()(const point_type& lhs, const point_type& rhs) const {
       if (lhs.x() == rhs.x())
         return lhs.y() < rhs.y();
       return lhs.x() < rhs.x();
@@ -133,7 +133,7 @@ public:
     typedef Site site_type;
     typedef Circle circle_type;
 
-    bool operator()(const site_type &lhs, const site_type &rhs) const {
+    bool operator()(const site_type& lhs, const site_type& rhs) const {
       if (lhs.x0() != rhs.x0())
         return lhs.x0() < rhs.x0();
       if (!lhs.is_segment()) {
@@ -156,7 +156,7 @@ public:
       }
     }
 
-    bool operator()(const site_type &lhs, const circle_type &rhs) const {
+    bool operator()(const site_type& lhs, const circle_type& rhs) const {
       typename ulp_cmp_type::Result xCmp =
           ulp_cmp(to_fpt(lhs.x()), to_fpt(rhs.lower_x()), ULPS);
       if (xCmp != ulp_cmp_type::EQUAL)
@@ -166,7 +166,7 @@ public:
       return yCmp == ulp_cmp_type::LESS;
     }
 
-    bool operator()(const circle_type &lhs, const site_type &rhs) const {
+    bool operator()(const circle_type& lhs, const site_type& rhs) const {
       typename ulp_cmp_type::Result xCmp =
           ulp_cmp(to_fpt(lhs.lower_x()), to_fpt(rhs.x()), ULPS);
       if (xCmp != ulp_cmp_type::EQUAL)
@@ -176,7 +176,7 @@ public:
       return yCmp == ulp_cmp_type::LESS;
     }
 
-    bool operator()(const circle_type &lhs, const circle_type &rhs) const {
+    bool operator()(const circle_type& lhs, const circle_type& rhs) const {
       typename ulp_cmp_type::Result xCmp =
           ulp_cmp(to_fpt(lhs.lower_x()), to_fpt(rhs.lower_x()), ULPSx2);
       if (xCmp != ulp_cmp_type::EQUAL)
@@ -199,9 +199,9 @@ public:
     // Returns true if a horizontal line going through a new site intersects
     // right arc at first, else returns false. If horizontal line goes
     // through intersection point of the given two arcs returns false also.
-    bool operator()(const site_type &left_site,
-                    const site_type &right_site,
-                    const site_type &new_site) const {
+    bool operator()(const site_type& left_site,
+                    const site_type& right_site,
+                    const site_type& new_site) const {
       if (!left_site.is_segment()) {
         if (!right_site.is_segment()) {
           return pp(left_site, right_site, new_site);
@@ -232,12 +232,12 @@ public:
     // Returns true if a horizontal line going through the new point site
     // intersects right arc at first, else returns false. If horizontal line
     // goes through intersection point of the given two arcs returns false.
-    bool pp(const site_type &left_site,
-            const site_type &right_site,
-            const site_type &new_site) const {
-      const point_type &left_point = left_site.point0();
-      const point_type &right_point = right_site.point0();
-      const point_type &new_point = new_site.point0();
+    bool pp(const site_type& left_site,
+            const site_type& right_site,
+            const site_type& new_site) const {
+      const point_type& left_point = left_site.point0();
+      const point_type& right_point = right_site.point0();
+      const point_type& new_point = new_site.point0();
       if (left_point.x() > right_point.x()) {
         if (new_point.y() <= left_point.y())
           return false;
@@ -257,8 +257,8 @@ public:
       return dist1 < dist2;
     }
 
-    bool ps(const site_type &left_site, const site_type &right_site,
-            const site_type &new_site, bool reverse_order) const {
+    bool ps(const site_type& left_site, const site_type& right_site,
+            const site_type& new_site, bool reverse_order) const {
       kPredicateResult fast_res = fast_ps(
         left_site, right_site, new_site, reverse_order);
       if (fast_res != UNDEFINED)
@@ -273,9 +273,9 @@ public:
       return reverse_order ^ (dist1 < dist2);
     }
 
-    bool ss(const site_type &left_site,
-            const site_type &right_site,
-            const site_type &new_site) const {
+    bool ss(const site_type& left_site,
+            const site_type& right_site,
+            const site_type& new_site) const {
       // Handle temporary segment sites.
       if (left_site.point0() == right_site.point0() &&
           left_site.point1() == right_site.point1()) {
@@ -294,7 +294,7 @@ public:
     }
 
     fpt_type find_distance_to_point_arc(
-        const site_type &site, const point_type &point) const {
+        const site_type& site, const point_type& point) const {
       fpt_type dx = to_fpt(site.x()) - to_fpt(point.x());
       fpt_type dy = to_fpt(site.y()) - to_fpt(point.y());
       // The relative error is at most 3EPS.
@@ -302,12 +302,12 @@ public:
     }
 
     fpt_type find_distance_to_segment_arc(
-        const site_type &site, const point_type &point) const {
+        const site_type& site, const point_type& point) const {
       if (is_vertical(site)) {
         return (to_fpt(site.x()) - to_fpt(point.x())) * to_fpt(0.5);
       } else {
-        const point_type &segment0 = site.point0(true);
-        const point_type &segment1 = site.point1(true);
+        const point_type& segment0 = site.point0(true);
+        const point_type& segment1 = site.point1(true);
         fpt_type a1 = to_fpt(segment1.x()) - to_fpt(segment0.x());
         fpt_type b1 = to_fpt(segment1.y()) - to_fpt(segment0.y());
         fpt_type k = get_sqrt(a1 * a1 + b1 * b1);
@@ -327,12 +327,12 @@ public:
     }
 
     kPredicateResult fast_ps(
-        const site_type &left_site, const site_type &right_site,
-        const site_type &new_site, bool reverse_order) const {
-      const point_type &site_point = left_site.point0();
-      const point_type &segment_start = right_site.point0(true);
-      const point_type &segment_end = right_site.point1(true);
-      const point_type &new_point = new_site.point0();
+        const site_type& left_site, const site_type& right_site,
+        const site_type& new_site, bool reverse_order) const {
+      const point_type& site_point = left_site.point0();
+      const point_type& segment_start = right_site.point0(true);
+      const point_type& segment_end = right_site.point1(true);
+      const point_type& new_point = new_site.point0();
 
       if (ot::eval(segment_start, segment_end, new_point) != ot::RIGHT)
         return (!right_site.is_inverse()) ? LESS : MORE;
@@ -392,11 +392,11 @@ public:
     // Comparison is only called during the new site events processing.
     // That's why one of the nodes will always lie on the sweepline and may
     // be represented as a straight horizontal line.
-    bool operator() (const node_type &node1,
-                     const node_type &node2) const {
+    bool operator() (const node_type& node1,
+                     const node_type& node2) const {
       // Get x coordinate of the rightmost site from both nodes.
-      const site_type &site1 = get_comparison_site(node1);
-      const site_type &site2 = get_comparison_site(node2);
+      const site_type& site1 = get_comparison_site(node1);
+      const site_type& site2 = get_comparison_site(node2);
 
       if (site1.x() < site2.x()) {
         // The second node contains a new site.
@@ -425,7 +425,7 @@ public:
 
   private:
     // Get the newer site.
-    const site_type &get_comparison_site(const node_type &node) const {
+    const site_type& get_comparison_site(const node_type& node) const {
       if (node.left_site().sorted_index() > node.right_site().sorted_index()) {
         return node.left_site();
       }
@@ -434,7 +434,7 @@ public:
 
     // Get comparison pair: y coordinate and direction of the newer site.
     std::pair<coordinate_type, int> get_comparison_y(
-      const node_type &node, bool is_new_node = true) const {
+      const node_type& node, bool is_new_node = true) const {
       if (node.left_site().sorted_index() ==
           node.right_site().sorted_index()) {
         return std::make_pair(node.left_site().y(), 0);
@@ -459,16 +459,16 @@ public:
     typedef typename Site::point_type point_type;
     typedef Site site_type;
 
-    bool ppp(const site_type &site1,
-             const site_type &site2,
-             const site_type &site3) const {
+    bool ppp(const site_type& site1,
+             const site_type& site2,
+             const site_type& site3) const {
       return ot::eval(site1.point0(), site2.point0(), site3.point0()) ==
              ot::RIGHT;
     }
 
-    bool pps(const site_type &site1,
-             const site_type &site2,
-             const site_type &site3,
+    bool pps(const site_type& site1,
+             const site_type& site2,
+             const site_type& site3,
              int segment_index) const {
       if (segment_index != 2) {
         typename ot::Orientation orient1 = ot::eval(site1.point0(),
@@ -492,9 +492,9 @@ public:
       return true;
     }
 
-    bool pss(const site_type &site1,
-             const site_type &site2,
-             const site_type &site3,
+    bool pss(const site_type& site1,
+             const site_type& site2,
+             const site_type& site3,
              int point_index) const {
       if (site2.point0() == site3.point0() &&
           site2.point1() == site3.point1()) {
@@ -512,9 +512,9 @@ public:
       return true;
     }
 
-    bool sss(const site_type &site1,
-             const site_type &site2,
-             const site_type &site3) const {
+    bool sss(const site_type& site1,
+             const site_type& site2,
+             const site_type& site3) const {
       if (site1.point0() == site2.point0() && site1.point1() == site2.point1())
         return false;
       if (site2.point0() == site3.point0() && site2.point1() == site3.point1())
@@ -532,10 +532,10 @@ public:
     typedef robust_sqrt_expr<big_int_type, efpt_type, to_efpt_converter>
         robust_sqrt_expr_type;
 
-    void ppp(const site_type &site1,
-             const site_type &site2,
-             const site_type &site3,
-             circle_type &circle,
+    void ppp(const site_type& site1,
+             const site_type& site2,
+             const site_type& site3,
+             circle_type& circle,
              bool recompute_c_x = true,
              bool recompute_c_y = true,
              bool recompute_lower_x = true) {
@@ -600,11 +600,11 @@ public:
     }
 
     // Recompute parameters of the circle event using high-precision library.
-    void pps(const site_type &site1,
-             const site_type &site2,
-             const site_type &site3,
+    void pps(const site_type& site1,
+             const site_type& site2,
+             const site_type& site3,
              int segment_index,
-             circle_type &c_event,
+             circle_type& c_event,
              bool recompute_c_x = true,
              bool recompute_c_y = true,
              bool recompute_lower_x = true) {
@@ -696,19 +696,19 @@ public:
     }
 
     // Recompute parameters of the circle event using high-precision library.
-    void pss(const site_type &site1,
-             const site_type &site2,
-             const site_type &site3,
+    void pss(const site_type& site1,
+             const site_type& site2,
+             const site_type& site3,
              int point_index,
-             circle_type &c_event,
+             circle_type& c_event,
              bool recompute_c_x = true,
              bool recompute_c_y = true,
              bool recompute_lower_x = true) {
       big_int_type a[2], b[2], c[2], cA[4], cB[4];
-      const point_type &segm_start1 = site2.point1(true);
-      const point_type &segm_end1 = site2.point0(true);
-      const point_type &segm_start2 = site3.point0(true);
-      const point_type &segm_end2 = site3.point1(true);
+      const point_type& segm_start1 = site2.point1(true);
+      const point_type& segm_end1 = site2.point0(true);
+      const point_type& segm_start2 = site3.point0(true);
+      const point_type& segm_end2 = site3.point1(true);
       a[0] = static_cast<int_x2_type>(segm_end1.x()) -
              static_cast<int_x2_type>(segm_start1.x());
       b[0] = static_cast<int_x2_type>(segm_end1.y()) -
@@ -829,10 +829,10 @@ public:
     }
 
     // Recompute parameters of the circle event using high-precision library.
-    void sss(const site_type &site1,
-             const site_type &site2,
-             const site_type &site3,
-             circle_type &c_event,
+    void sss(const site_type& site1,
+             const site_type& site2,
+             const site_type& site3,
+             circle_type& c_event,
              bool recompute_c_x = true,
              bool recompute_c_y = true,
              bool recompute_lower_x = true) {
@@ -992,10 +992,10 @@ public:
     typedef mp_circle_formation_functor<site_type, circle_type>
         exact_circle_formation_functor_type;
 
-    void ppp(const site_type &site1,
-             const site_type &site2,
-             const site_type &site3,
-             circle_type &c_event) {
+    void ppp(const site_type& site1,
+             const site_type& site2,
+             const site_type& site3,
+             circle_type& c_event) {
       fpt_type dif_x1 = to_fpt(site1.x()) - to_fpt(site2.x());
       fpt_type dif_x2 = to_fpt(site2.x()) - to_fpt(site3.x());
       fpt_type dif_y1 = to_fpt(site1.y()) - to_fpt(site2.y());
@@ -1040,11 +1040,11 @@ public:
       }
     }
 
-    void pps(const site_type &site1,
-             const site_type &site2,
-             const site_type &site3,
+    void pps(const site_type& site1,
+             const site_type& site2,
+             const site_type& site3,
              int segment_index,
-             circle_type &c_event) {
+             circle_type& c_event) {
       fpt_type line_a = to_fpt(site3.point1(true).y()) -
                         to_fpt(site3.point0(true).y());
       fpt_type line_b = to_fpt(site3.point0(true).x()) -
@@ -1113,15 +1113,15 @@ public:
       }
     }
 
-    void pss(const site_type &site1,
-             const site_type &site2,
-             const site_type &site3,
+    void pss(const site_type& site1,
+             const site_type& site2,
+             const site_type& site3,
              int point_index,
-             circle_type &c_event) {
-      const point_type &segm_start1 = site2.point1(true);
-      const point_type &segm_end1 = site2.point0(true);
-      const point_type &segm_start2 = site3.point0(true);
-      const point_type &segm_end2 = site3.point1(true);
+             circle_type& c_event) {
+      const point_type& segm_start1 = site2.point1(true);
+      const point_type& segm_end1 = site2.point0(true);
+      const point_type& segm_start2 = site3.point0(true);
+      const point_type& segm_end2 = site3.point1(true);
       fpt_type a1 = to_fpt(segm_end1.x()) - to_fpt(segm_start1.x());
       fpt_type b1 = to_fpt(segm_end1.y()) - to_fpt(segm_start1.y());
       fpt_type a2 = to_fpt(segm_end2.x()) - to_fpt(segm_start2.x());
@@ -1271,10 +1271,10 @@ public:
       }
     }
 
-    void sss(const site_type &site1,
-             const site_type &site2,
-             const site_type &site3,
-             circle_type &c_event) {
+    void sss(const site_type& site1,
+             const site_type& site2,
+             const site_type& site3,
+             circle_type& c_event) {
       robust_fpt_type a1(to_fpt(site1.x1(true)) - to_fpt(site1.x0(true)));
       robust_fpt_type b1(to_fpt(site1.y1(true)) - to_fpt(site1.y0(true)));
       robust_fpt_type c1(robust_cross_product(
@@ -1371,8 +1371,8 @@ public:
     // Create a circle event from the given three sites.
     // Returns true if the circle event exists, else false.
     // If exists circle event is saved into the c_event variable.
-    bool operator()(const site_type &site1, const site_type &site2,
-                    const site_type &site3, circle_type &circle) {
+    bool operator()(const site_type& site1, const site_type& site2,
+                    const site_type& site3, circle_type& circle) {
       if (!site1.is_segment()) {
         if (!site2.is_segment()) {
           if (!site3.is_segment()) {
diff --git a/include/boost/polygon/detail/voronoi_robust_fpt.hpp b/include/boost/polygon/detail/voronoi_robust_fpt.hpp
index a0123d6..478fa7b 100644
--- a/include/boost/polygon/detail/voronoi_robust_fpt.hpp
+++ b/include/boost/polygon/detail/voronoi_robust_fpt.hpp
@@ -94,7 +94,7 @@ public:
   relative_error_type re() const { return re_; }
   relative_error_type ulp() const { return re_; }
 
-  robust_fpt& operator=(const robust_fpt &that) {
+  robust_fpt& operator=(const robust_fpt& that) {
     this->fpv_ = that.fpv_;
     this->re_ = that.re_;
     return *this;
@@ -116,7 +116,7 @@ public:
     return robust_fpt(-fpv_, re_);
   }
 
-  robust_fpt& operator+=(const robust_fpt &that) {
+  robust_fpt& operator+=(const robust_fpt& that) {
     floating_point_type fpv = this->fpv_ + that.fpv_;
     if ((!is_neg(this->fpv_) && !is_neg(that.fpv_)) ||
         (!is_pos(this->fpv_) && !is_pos(that.fpv_)))
@@ -132,7 +132,7 @@ public:
     return *this;
   }
 
-  robust_fpt& operator-=(const robust_fpt &that) {
+  robust_fpt& operator-=(const robust_fpt& that) {
     floating_point_type fpv = this->fpv_ - that.fpv_;
     if ((!is_neg(this->fpv_) && !is_pos(that.fpv_)) ||
         (!is_pos(this->fpv_) && !is_neg(that.fpv_)))
@@ -148,19 +148,19 @@ public:
     return *this;
   }
 
-  robust_fpt& operator*=(const robust_fpt &that) {
+  robust_fpt& operator*=(const robust_fpt& that) {
     this->re_ += that.re_ + ROUNDING_ERROR;
     this->fpv_ *= that.fpv_;
     return *this;
   }
 
-  robust_fpt& operator/=(const robust_fpt &that) {
+  robust_fpt& operator/=(const robust_fpt& that) {
     this->re_ += that.re_ + ROUNDING_ERROR;
     this->fpv_ /= that.fpv_;
     return *this;
   }
 
-  robust_fpt operator+(const robust_fpt &that) const {
+  robust_fpt operator+(const robust_fpt& that) const {
     floating_point_type fpv = this->fpv_ + that.fpv_;
     relative_error_type re;
     if ((!is_neg(this->fpv_) && !is_neg(that.fpv_)) ||
@@ -176,7 +176,7 @@ public:
     return robust_fpt(fpv, re);
   }
 
-  robust_fpt operator-(const robust_fpt &that) const {
+  robust_fpt operator-(const robust_fpt& that) const {
     floating_point_type fpv = this->fpv_ - that.fpv_;
     relative_error_type re;
     if ((!is_neg(this->fpv_) && !is_pos(that.fpv_)) ||
@@ -192,13 +192,13 @@ public:
     return robust_fpt(fpv, re);
   }
 
-  robust_fpt operator*(const robust_fpt &that) const {
+  robust_fpt operator*(const robust_fpt& that) const {
     floating_point_type fpv = this->fpv_ * that.fpv_;
     relative_error_type re = this->re_ + that.re_ + ROUNDING_ERROR;
     return robust_fpt(fpv, re);
   }
 
-  robust_fpt operator/(const robust_fpt &that) const {
+  robust_fpt operator/(const robust_fpt& that) const {
     floating_point_type fpv = this->fpv_ / that.fpv_;
     relative_error_type re = this->re_ + that.re_ + ROUNDING_ERROR;
     return robust_fpt(fpv, re);
@@ -251,11 +251,11 @@ public:
       positive_sum_(0),
       negative_sum_(0) {}
 
-  robust_dif(const T &value) :
+  robust_dif(const T& value) :
       positive_sum_((value>0)?value:0),
       negative_sum_((value<0)?-value:0) {}
 
-  robust_dif(const T &pos, const T &neg) :
+  robust_dif(const T& pos, const T& neg) :
       positive_sum_(pos),
       negative_sum_(neg) {}
 
@@ -275,7 +275,7 @@ public:
     return robust_dif(negative_sum_, positive_sum_);
   }
 
-  robust_dif<T> &operator+=(const T &val) {
+  robust_dif<T>& operator+=(const T& val) {
     if (!is_neg(val))
       positive_sum_ += val;
     else
@@ -283,13 +283,13 @@ public:
     return *this;
   }
 
-  robust_dif<T> &operator+=(const robust_dif<T> &that) {
+  robust_dif<T>& operator+=(const robust_dif<T>& that) {
     positive_sum_ += that.positive_sum_;
     negative_sum_ += that.negative_sum_;
     return *this;
   }
 
-  robust_dif<T> &operator-=(const T &val) {
+  robust_dif<T>& operator-=(const T& val) {
     if (!is_neg(val))
       negative_sum_ += val;
     else
@@ -297,13 +297,13 @@ public:
     return *this;
   }
 
-  robust_dif<T> &operator-=(const robust_dif<T> &that) {
+  robust_dif<T>& operator-=(const robust_dif<T>& that) {
     positive_sum_ += that.negative_sum_;
     negative_sum_ += that.positive_sum_;
     return *this;
   }
 
-  robust_dif<T> &operator*=(const T &val) {
+  robust_dif<T>& operator*=(const T& val) {
     if (!is_neg(val)) {
       positive_sum_ *= val;
       negative_sum_ *= val;
@@ -315,7 +315,7 @@ public:
     return *this;
   }
 
-  robust_dif<T> &operator*=(const robust_dif<T> &that) {
+  robust_dif<T>& operator*=(const robust_dif<T>& that) {
     T positive_sum = this->positive_sum_ * that.positive_sum_ +
                      this->negative_sum_ * that.negative_sum_;
     T negative_sum = this->positive_sum_ * that.negative_sum_ +
@@ -325,7 +325,7 @@ public:
     return *this;
   }
 
-  robust_dif<T> &operator/=(const T &val) {
+  robust_dif<T>& operator/=(const T& val) {
     if (!is_neg(val)) {
       positive_sum_ /= val;
       negative_sum_ /= val;
@@ -442,7 +442,7 @@ public:
 
   // Evaluates expression (re = 4 EPS):
   // A[0] * sqrt(B[0]).
-  _fpt eval1(_int *A, _int *B) {
+  _fpt eval1(_int* A, _int* B) {
     _fpt a = convert(A[0]);
     _fpt b = convert(B[0]);
     return a * get_sqrt(b);
@@ -450,7 +450,7 @@ public:
 
   // Evaluates expression (re = 7 EPS):
   // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]).
-  _fpt eval2(_int *A, _int *B) {
+  _fpt eval2(_int* A, _int* B) {
     _fpt a = eval1(A, B);
     _fpt b = eval1(A + 1, B + 1);
     if ((!is_neg(a) && !is_neg(b)) ||
@@ -461,7 +461,7 @@ public:
 
   // Evaluates expression (re = 16 EPS):
   // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2]).
-  _fpt eval3(_int *A, _int *B) {
+  _fpt eval3(_int* A, _int* B) {
     _fpt a = eval2(A, B);
     _fpt b = eval1(A + 2, B + 2);
     if ((!is_neg(a) && !is_neg(b)) ||
@@ -478,7 +478,7 @@ public:
   // Evaluates expression (re = 25 EPS):
   // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) +
   // A[2] * sqrt(B[2]) + A[3] * sqrt(B[3]).
-  _fpt eval4(_int *A, _int *B) {
+  _fpt eval4(_int* A, _int* B) {
     _fpt a = eval2(A, B);
     _fpt b = eval2(A + 2, B + 2);
     if ((!is_neg(a) && !is_neg(b)) ||
diff --git a/include/boost/polygon/detail/voronoi_structures.hpp b/include/boost/polygon/detail/voronoi_structures.hpp
index 0fd34fb..5e8140d 100644
--- a/include/boost/polygon/detail/voronoi_structures.hpp
+++ b/include/boost/polygon/detail/voronoi_structures.hpp
@@ -30,11 +30,11 @@ public:
       x_(x),
       y_(y) {}
 
-  bool operator==(const point_2d &that) const {
+  bool operator==(const point_2d& that) const {
     return (this->x_ == that.x()) && (this->y_ == that.y());
   }
 
-  bool operator!=(const point_2d &that) const {
+  bool operator!=(const point_2d& that) const {
     return (this->x_ != that.x()) || (this->y_ != that.y());
   }
 
@@ -128,7 +128,7 @@ public:
       sorted_index_(0),
       flags_(0) {}
 
-  site_event(const point_type &point) :
+  site_event(const point_type& point) :
       point0_(point),
       point1_(point),
       sorted_index_(0),
@@ -141,18 +141,18 @@ public:
       sorted_index_(0),
       flags_(0) {}
 
-  site_event(const point_type &point1, const point_type &point2) :
+  site_event(const point_type& point1, const point_type& point2) :
       point0_(point1),
       point1_(point2),
       sorted_index_(0),
       flags_(0) {}
 
-  bool operator==(const site_event &that) const {
+  bool operator==(const site_event& that) const {
     return (this->point0_ == that.point0_) &&
            (this->point1_ == that.point1_);
   }
 
-  bool operator!=(const site_event &that) const {
+  bool operator!=(const site_event& that) const {
     return (this->point0_ != that.point0_) ||
            (this->point1_ != that.point1_);
   }
@@ -189,13 +189,13 @@ public:
     return is_inverse() ? point0_.y() : point1_.y();
   }
 
-  const point_type &point0(bool oriented = false) const {
+  const point_type& point0(bool oriented = false) const {
     if (!oriented)
       return point0_;
     return is_inverse() ? point1_ : point0_;
   }
 
-  const point_type &point1(bool oriented = false) const {
+  const point_type& point1(bool oriented = false) const {
     if (!oriented)
       return point1_;
     return is_inverse() ? point0_ : point1_;
@@ -220,7 +220,7 @@ public:
   }
 
   bool is_inverse() const {
-    return (flags_ & IS_INVERSE) ? true : false;
+    return flags_ & IS_INVERSE;
   }
 
   site_event& inverse() {
@@ -456,31 +456,31 @@ private:
 template <typename Edge, typename Circle>
 class beach_line_node_data {
 public:
-  explicit beach_line_node_data(Edge *new_edge) :
+  explicit beach_line_node_data(Edge* new_edge) :
       circle_event_(NULL),
       edge_(new_edge) {}
 
-  Circle *circle_event() const {
+  Circle* circle_event() const {
     return circle_event_;
   }
 
-  beach_line_node_data& circle_event(Circle *circle_event) {
+  beach_line_node_data& circle_event(Circle* circle_event) {
     circle_event_ = circle_event;
     return *this;
   }
 
-  Edge *edge() const {
+  Edge* edge() const {
     return edge_;
   }
 
-  beach_line_node_data& edge(Edge *new_edge) {
+  beach_line_node_data& edge(Edge* new_edge) {
     edge_ = new_edge;
     return *this;
   }
 
 private:
-  Circle *circle_event_;
-  Edge *edge_;
+  Circle* circle_event_;
+  Edge* edge_;
 };
 }  // detail
 }  // polygon
diff --git a/include/boost/polygon/voronoi.hpp b/include/boost/polygon/voronoi.hpp
index 79f4a55..9785dec 100644
--- a/include/boost/polygon/voronoi.hpp
+++ b/include/boost/polygon/voronoi.hpp
@@ -41,7 +41,7 @@ typename enable_if<
   >::type,
   void
 >::type
-insert(const Point &point, VB *vb) {
+insert(const Point& point, VB* vb) {
   vb->insert_point(x(point), y(point));
 }
 
@@ -56,7 +56,7 @@ typename enable_if<
   >::type,
   void
 >::type
-insert(PointIterator first, const PointIterator last, VB *vb) {
+insert(PointIterator first, const PointIterator last, VB* vb) {
   for (PointIterator it = first; it != last; ++it) {
     insert(*it, vb);
   }
@@ -71,7 +71,7 @@ typename enable_if<
   >::type,
   void
 >::type
-insert(const Segment &segment, VB *vb) {
+insert(const Segment& segment, VB* vb) {
   vb->insert_segment(x(low(segment)), y(low(segment)), x(high(segment)), y(high(segment)));
 }
 
@@ -86,7 +86,7 @@ typename enable_if<
   >::type,
   void
 >::type
-insert(SegmentIterator first, SegmentIterator last, VB *vb) {
+insert(SegmentIterator first, SegmentIterator last, VB* vb) {
   for (SegmentIterator it = first; it != last; ++it) {
     insert(*it, vb);
   }
@@ -103,7 +103,7 @@ typename enable_if<
   >::type,
   void
 >::type
-construct_voronoi(PointIterator first, PointIterator last, VD *vd) {
+construct_voronoi(PointIterator first, PointIterator last, VD* vd) {
   default_voronoi_builder builder;
   insert(first, last, &builder);
   builder.construct(vd);
@@ -120,7 +120,7 @@ typename enable_if<
   >::type,
   void
 >::type
-construct_voronoi(SegmentIterator first, SegmentIterator last, VD *vd) {
+construct_voronoi(SegmentIterator first, SegmentIterator last, VD* vd) {
   default_voronoi_builder builder;
   insert(first, last, &builder);
   builder.construct(vd);
@@ -147,7 +147,7 @@ typename enable_if<
   void
 >::type
 construct_voronoi(PointIterator p_first, PointIterator p_last,
-    SegmentIterator s_first, SegmentIterator s_last, VD *vd) {
+    SegmentIterator s_first, SegmentIterator s_last, VD* vd) {
   default_voronoi_builder builder;
   insert(p_first, p_last, &builder);
   insert(s_first, s_last, &builder);
diff --git a/include/boost/polygon/voronoi_builder.hpp b/include/boost/polygon/voronoi_builder.hpp
index a9050a8..6690bd1 100644
--- a/include/boost/polygon/voronoi_builder.hpp
+++ b/include/boost/polygon/voronoi_builder.hpp
@@ -94,7 +94,7 @@ public:
   template <typename OUTPUT>
   void construct(OUTPUT *output) {
     // Init structures.
-    output->builder()->reserve(site_events_.size());
+    output->_reserve(site_events_.size());
     init_sites_queue();
     init_beach_line(output);
 
@@ -107,7 +107,8 @@ public:
       } else if (site_event_iterator_ == site_events_.end()) {
         process_circle_event(output);
       } else {
-        if (event_comparison(*site_event_iterator_, circle_events_.top().first)) {
+        if (event_comparison(*site_event_iterator_,
+                             circle_events_.top().first)) {
           process_site_event(output);
         } else {
           process_circle_event(output);
@@ -121,7 +122,7 @@ public:
     beach_line_.clear();
 
     // Finish construction.
-    output->builder()->build();
+    output->_build();
   }
 
   void clear() {
@@ -172,8 +173,9 @@ private:
         site_events_.begin(), site_events_.end()), site_events_.end());
 
     // Index sites.
-    for (std::size_t cur = 0; cur < site_events_.size(); ++cur)
+    for (std::size_t cur = 0; cur < site_events_.size(); ++cur) {
       site_events_[cur].sorted_index(cur);
+    }
 
     // Init site iterator.
     site_event_iterator_ = site_events_.begin();
@@ -185,7 +187,7 @@ private:
       return;
     if (site_events_.size() == 1) {
       // Handle single site event case.
-      output->builder()->process_single_site(site_events_[0]);
+      output->_process_single_site(site_events_[0]);
       ++site_event_iterator_;
     } else {
       int skip = 0;
@@ -234,7 +236,7 @@ private:
 
       // Update the output.
       edge_type *edge =
-          output->builder()->insert_new_edge(*it_first, *it_second).first;
+          output->_insert_new_edge(*it_first, *it_second).first;
 
       // Insert a new bisector into the beach line.
       beach_line_.insert(beach_line_.end(),
@@ -395,7 +397,7 @@ private:
     const_cast<key_type &>(it_first->first).right_site(site3);
 
     // Insert the new bisector into the beach line.
-    it_first->second.edge(output->builder()->insert_new_edge(
+    it_first->second.edge(output->_insert_new_edge(
         site1, site3, circle_event, bisector1, bisector2).first);
 
     // Remove the (B, C) bisector node from the beach line.
@@ -441,7 +443,7 @@ private:
 
     // Update the output.
     std::pair<edge_type*, edge_type*> edges =
-        output->builder()->insert_new_edge(site_arc2, site_event);
+        output->_insert_new_edge(site_arc2, site_event);
     position = beach_line_.insert(position,
         typename beach_line_type::value_type(
             new_right_node, value_type(edges.second)));
diff --git a/include/boost/polygon/voronoi_diagram.hpp b/include/boost/polygon/voronoi_diagram.hpp
index eda669f..3fc4dd7 100644
--- a/include/boost/polygon/voronoi_diagram.hpp
+++ b/include/boost/polygon/voronoi_diagram.hpp
@@ -222,25 +222,25 @@ public:
   // Returns true if the edge is linear (segment, ray, line).
   // Returns false if the edge is curved (parabolic arc).
   bool is_linear() const {
-    return (color_ & BIT_IS_LINEAR) ? true : false;
+    return color_ & BIT_IS_LINEAR;
   }
 
   // Returns true if the edge is curved (parabolic arc).
   // Returns false if the edge is linear (segment, ray, line).
   bool is_curved() const {
-    return (color_ & BIT_IS_LINEAR) ? false : true;
+    return !(color_ & BIT_IS_LINEAR);
   }
 
   // Returns false if edge goes through the endpoint of the segment.
   // Returns true else.
   bool is_primary() const {
-    return (color_ & BIT_IS_PRIMARY) ? true : false;
+    return color_ & BIT_IS_PRIMARY;
   }
 
   // Returns true if edge goes through the endpoint of the segment.
   // Returns false else.
   bool is_secondary() const {
-    return (color_ & BIT_IS_PRIMARY) ? false : true;
+    return !(color_ & BIT_IS_PRIMARY);
   }
 
   color_type color() const { return color_ >> BITS_SHIFT; }
@@ -309,58 +309,12 @@ public:
   typedef typename edge_container_type::iterator edge_iterator;
   typedef typename edge_container_type::const_iterator const_edge_iterator;
 
-  // This builder class is mainly used to hide from the user methods that
-  // construct the Voronoi diagram.
-  class voronoi_diagram_builder {
-  public:
-    void vd(voronoi_diagram* vd) {
-      vd_ = vd;
-    }
-
-    bool done() {
-      return vd_ == NULL;
-    }
-
-    void reserve(int num_sites) {
-      vd_->reserve(num_sites);
-    }
-
-    template <typename SEvent>
-    void process_single_site(const SEvent& site) {
-      vd_->process_single_site(site);
-    }
-
-    template <typename SEvent>
-    std::pair<void*, void*> insert_new_edge(
-        const SEvent& site1, const SEvent& site2) {
-      return vd_->insert_new_edge(site1, site2);
-    }
-
-    template <typename SEvent, typename CEvent>
-    std::pair<void*, void*> insert_new_edge(
-        const SEvent& site1, const SEvent& site3, const CEvent& circle,
-        void* data12, void* data23) {
-      return vd_->insert_new_edge(site1, site3, circle, data12, data23);
-    }
-
-    void build() {
-      vd_->build();
-      vd_ = NULL;
-    }
-
-  private:
-    voronoi_diagram* vd_;
-  };
-
-  voronoi_diagram() {
-    builder_.vd(&(*this));
-  }
+  voronoi_diagram() {}
 
   void clear() {
     cells_.clear();
     vertices_.clear();
     edges_.clear();
-    builder_.vd(&(*this));
   }
 
   const cell_container_type& cells() const {
@@ -391,51 +345,14 @@ public:
     return vertices_.size();
   }
 
-  voronoi_diagram_builder* builder() {
-    if (builder_.done()) {
-      return NULL;
-    } else {
-      return &builder_;
-    }
-  }
-
-private:
-  typedef typename TRAITS::vertex_equality_predicate_type
-    vertex_equality_predicate_type;
-
-  friend class voronoi_diagram_builder;
-
-  void reserve(int num_sites) {
+  void _reserve(int num_sites) {
     cells_.reserve(num_sites);
     vertices_.reserve(num_sites << 1);
     edges_.reserve((num_sites << 2) + (num_sites << 1));
   }
 
-  template <typename SEvent>
-  bool is_primary_edge(const SEvent& site1, const SEvent& site2) const {
-    bool flag1 = site1.is_segment();
-    bool flag2 = site2.is_segment();
-    if (flag1 && !flag2) {
-      return (site1.point0() != site2.point0()) &&
-             (site1.point1() != site2.point0());
-    }
-    if (!flag1 && flag2) {
-      return (site2.point0() != site1.point0()) &&
-             (site2.point1() != site1.point0());
-    }
-    return true;
-  }
-
-  template <typename SEvent>
-  bool is_linear_edge(const SEvent& site1, const SEvent& site2) const {
-    if (!is_primary_edge(site1, site2)) {
-      return true;
-    }
-    return !(site1.is_segment() ^ site2.is_segment());
-  }
-
-  template <typename SEvent>
-  void process_single_site(const SEvent& site) {
+template <typename SEvent>
+  void _process_single_site(const SEvent& site) {
     cells_.push_back(cell_type(
          site.initial_index(), site.source_category(), NULL));
   }
@@ -444,7 +361,7 @@ private:
   // Takes as input left and right sites that form a new bisector.
   // Returns a pair of pointers to a new half-edges.
   template <typename SEvent>
-  std::pair<void*, void*> insert_new_edge(
+  std::pair<void*, void*> _insert_new_edge(
       const SEvent& site1, const SEvent& site2) {
     // Get sites' indexes.
     int site_index1 = site1.sorted_index();
@@ -491,7 +408,7 @@ private:
   // pointers to those half-edges. Half-edges' direction goes out of the
   // new Voronoi vertex point. Returns a pair of pointers to a new half-edges.
   template <typename SEvent, typename CEvent>
-  std::pair<void*, void*> insert_new_edge(
+  std::pair<void*, void*> _insert_new_edge(
       const SEvent& site1, const SEvent& site3, const CEvent& circle,
       void* data12, void* data23) {
     edge_type* edge12 = static_cast<edge_type*>(data12);
@@ -537,7 +454,7 @@ private:
     return std::make_pair(&new_edge1, &new_edge2);
   }
 
-  void build() {
+  void _build() {
     // Remove degenerate edges.
     edge_iterator last_edge = edges_.begin();
     for (edge_iterator it = edges_.begin(); it != edges_.end(); it += 2) {
@@ -647,6 +564,33 @@ private:
     }
   }
 
+private:
+  typedef typename TRAITS::vertex_equality_predicate_type
+    vertex_equality_predicate_type;
+
+  template <typename SEvent>
+  bool is_primary_edge(const SEvent& site1, const SEvent& site2) const {
+    bool flag1 = site1.is_segment();
+    bool flag2 = site2.is_segment();
+    if (flag1 && !flag2) {
+      return (site1.point0() != site2.point0()) &&
+             (site1.point1() != site2.point0());
+    }
+    if (!flag1 && flag2) {
+      return (site2.point0() != site1.point0()) &&
+             (site2.point1() != site1.point0());
+    }
+    return true;
+  }
+
+  template <typename SEvent>
+  bool is_linear_edge(const SEvent& site1, const SEvent& site2) const {
+    if (!is_primary_edge(site1, site2)) {
+      return true;
+    }
+    return !(site1.is_segment() ^ site2.is_segment());
+  }
+
   // Remove degenerate edge.
   void remove_edge(edge_type* edge) {
     // Update the endpoints of the incident edges to the second vertex.
@@ -676,8 +620,6 @@ private:
   cell_container_type cells_;
   vertex_container_type vertices_;
   edge_container_type edges_;
-
-  voronoi_diagram_builder builder_;
   vertex_equality_predicate_type vertex_equality_predicate_;
 
   // Disallow copy constructor and operator=