diff --git a/xs/src/admesh/connect.cpp b/xs/src/admesh/connect.cpp
index e9129d007..da5b66720 100644
--- a/xs/src/admesh/connect.cpp
+++ b/xs/src/admesh/connect.cpp
@@ -25,11 +25,11 @@
 #include <string.h>
 #include <math.h>
 
+#include <boost/detail/endian.hpp>
+
 #include "stl.h"
 
 
-static void stl_match_neighbors_exact(stl_file *stl,
-                                      stl_hash_edge *edge_a, stl_hash_edge *edge_b);
 static void stl_match_neighbors_nearby(stl_file *stl,
                                        stl_hash_edge *edge_a, stl_hash_edge *edge_b);
 static void stl_record_neighbors(stl_file *stl,
@@ -43,7 +43,6 @@ static int stl_load_edge_nearby(stl_file *stl, stl_hash_edge *edge,
 static void insert_hash_edge(stl_file *stl, stl_hash_edge edge,
                              void (*match_neighbors)(stl_file *stl,
                                  stl_hash_edge *edge_a, stl_hash_edge *edge_b));
-static int stl_get_hash_for_edge(int M, stl_hash_edge *edge);
 static int stl_compare_function(stl_hash_edge *edge_a, stl_hash_edge *edge_b);
 static void stl_free_edges(stl_file *stl);
 static void stl_remove_facet(stl_file *stl, int facet_number);
@@ -82,37 +81,20 @@ stl_check_facets_exact(stl_file *stl) {
 
   for(i = 0; i < stl->stats.number_of_facets; i++) {
     facet = stl->facet_start[i];
-    // Positive and negative zeros are possible in the floats, which are considered equal by the FP unit.
-    // When using a memcmp on raw floats, those numbers report to be different.
-    // Unify all +0 and -0 to +0 to make the floats equal under memcmp.
-    {
-      uint32_t *f = (uint32_t*)&facet;
-      for (int j = 0; j < 12; ++ j, ++ f) // 3x vertex + normal: 4x3 = 12 floats
-        if (*f == 0x80000000)
-            // Negative zero, switch to positive zero.
-            *f = 0;
-    }
-
-    /* If any two of the three vertices are found to be exactally the same, call them degenerate and remove the facet. */
-    if(   !memcmp(&facet.vertex[0], &facet.vertex[1],
-                  sizeof(stl_vertex))
-          || !memcmp(&facet.vertex[1], &facet.vertex[2],
-                     sizeof(stl_vertex))
-          || !memcmp(&facet.vertex[0], &facet.vertex[2],
-                     sizeof(stl_vertex))) {
+    // If any two of the three vertices are found to be exactally the same, call them degenerate and remove the facet.
+    if (facet.vertex[0] == facet.vertex[1] ||
+        facet.vertex[1] == facet.vertex[2] ||
+        facet.vertex[0] == facet.vertex[2]) {
       stl->stats.degenerate_facets += 1;
       stl_remove_facet(stl, i);
-      i--;
+      -- i;
       continue;
-
     }
     for(j = 0; j < 3; j++) {
       edge.facet_number = i;
       edge.which_edge = j;
-      stl_load_edge_exact(stl, &edge, &facet.vertex[j],
-                          &facet.vertex[(j + 1) % 3]);
-
-      insert_hash_edge(stl, edge, stl_match_neighbors_exact);
+      stl_load_edge_exact(stl, &edge, &facet.vertex[j], &facet.vertex[(j + 1) % 3]);
+      insert_hash_edge(stl, edge, stl_record_neighbors);
     }
   }
   stl_free_edges(stl);
@@ -131,28 +113,33 @@ stl_load_edge_exact(stl_file *stl, stl_hash_edge *edge,
   if (stl->error) return;
 
   {
-    float diff_x = ABS(a->x - b->x);
-    float diff_y = ABS(a->y - b->y);
-    float diff_z = ABS(a->z - b->z);
-    float max_diff = STL_MAX(diff_x, diff_y);
-    max_diff = STL_MAX(diff_z, max_diff);
-    stl->stats.shortest_edge = STL_MIN(max_diff, stl->stats.shortest_edge);
+    stl_vertex diff = (*a - *b).cwiseAbs();
+    float max_diff = std::max(diff(0), std::max(diff(1), diff(2)));
+    stl->stats.shortest_edge = std::min(max_diff, stl->stats.shortest_edge);
   }
 
   // Ensure identical vertex ordering of equal edges.
   // This method is numerically robust.
-  if ((a->x != b->x) ? 
-        (a->x < b->x) : 
-        ((a->y != b->y) ? 
-            (a->y < b->y) : 
-            (a->z < b->z))) {
-    memcpy(&edge->key[0], a, sizeof(stl_vertex));
-    memcpy(&edge->key[3], b, sizeof(stl_vertex));
+  if (stl_vertex_lower(*a, *b)) {
   } else {
-    memcpy(&edge->key[0], b, sizeof(stl_vertex));
-    memcpy(&edge->key[3], a, sizeof(stl_vertex));
+    std::swap(a, b);
     edge->which_edge += 3; /* this edge is loaded backwards */
   }
+  memcpy(&edge->key[0],                  a->data(), sizeof(stl_vertex));
+  memcpy(&edge->key[sizeof(stl_vertex)], b->data(), sizeof(stl_vertex));
+  // Switch negative zeros to positive zeros, so memcmp will consider them to be equal.
+  for (size_t i = 0; i < 6; ++ i) {
+    unsigned char *p = edge->key + i * 4;
+#ifdef BOOST_LITTLE_ENDIAN
+    if (p[0] == 0 && p[1] == 0 && p[2] == 0 && p[3] == 0x80)
+      // Negative zero, switch to positive zero.
+      p[3] = 0;
+#else /* BOOST_LITTLE_ENDIAN */
+    if (p[0] == 0x80 && p[1] == 0 && p[2] == 0 && p[3] == 0)
+      // Negative zero, switch to positive zero.
+      p[0] = 0;
+#endif /* BOOST_LITTLE_ENDIAN */
+  }
 }
 
 static void
@@ -188,21 +175,17 @@ stl_initialize_facet_check_exact(stl_file *stl) {
   }
 }
 
-static void
-insert_hash_edge(stl_file *stl, stl_hash_edge edge,
+static void insert_hash_edge(stl_file *stl, stl_hash_edge edge,
                  void (*match_neighbors)(stl_file *stl,
-                     stl_hash_edge *edge_a, stl_hash_edge *edge_b)) {
-  stl_hash_edge *link;
-  stl_hash_edge *new_edge;
-  stl_hash_edge *temp;
-  int            chain_number;
-
+                     stl_hash_edge *edge_a, stl_hash_edge *edge_b))
+{
   if (stl->error) return;
 
-  chain_number = stl_get_hash_for_edge(stl->M, &edge);
-
-  link = stl->heads[chain_number];
+  int            chain_number = edge.hash(stl->M);
+  stl_hash_edge *link = stl->heads[chain_number];
 
+  stl_hash_edge *new_edge;
+  stl_hash_edge *temp;
   if(link == stl->tail) {
     /* This list doesn't have any edges currently in it.  Add this one. */
     new_edge = (stl_hash_edge*)malloc(sizeof(stl_hash_edge));
@@ -252,30 +235,17 @@ insert_hash_edge(stl_file *stl, stl_hash_edge edge,
   }
 }
 
-
-static int
-stl_get_hash_for_edge(int M, stl_hash_edge *edge) {
-  return ((edge->key[0] / 23 + edge->key[1] / 19 + edge->key[2] / 17
-           + edge->key[3] /13  + edge->key[4] / 11 + edge->key[5] / 7 ) % M);
+// Return 1 if the edges are not matched.
+static inline int stl_compare_function(stl_hash_edge *edge_a, stl_hash_edge *edge_b)
+{
+    // Don't match edges of the same facet
+    return (edge_a->facet_number == edge_b->facet_number) || (*edge_a != *edge_b);
 }
 
-static int
-stl_compare_function(stl_hash_edge *edge_a, stl_hash_edge *edge_b) {
-  if(edge_a->facet_number == edge_b->facet_number) {
-    return 1;			/* Don't match edges of the same facet */
-  } else {
-    return memcmp(edge_a, edge_b, SIZEOF_EDGE_SORT);
-  }
-}
-
-void
-stl_check_facets_nearby(stl_file *stl, float tolerance) {
-  stl_hash_edge  edge[3];
-  stl_facet      facet;
-  int            i;
-  int            j;
-
-  if (stl->error) return;
+void stl_check_facets_nearby(stl_file *stl, float tolerance)
+{
+  if (stl->error)
+    return;
 
   if(   (stl->stats.connected_facets_1_edge == stl->stats.number_of_facets)
         && (stl->stats.connected_facets_2_edge == stl->stats.number_of_facets)
@@ -286,27 +256,19 @@ stl_check_facets_nearby(stl_file *stl, float tolerance) {
 
   stl_initialize_facet_check_nearby(stl);
 
-  for(i = 0; i < stl->stats.number_of_facets; i++) {
-    facet = stl->facet_start[i];
-    // Positive and negative zeros are possible in the floats, which are considered equal by the FP unit.
-    // When using a memcmp on raw floats, those numbers report to be different.
-    // Unify all +0 and -0 to +0 to make the floats equal under memcmp.
-    {
-      uint32_t *f = (uint32_t*)&facet;
-      for (int j = 0; j < 12; ++ j, ++ f) // 3x vertex + normal: 4x3 = 12 floats
-        if (*f == 0x80000000)
-            // Negative zero, switch to positive zero.
-            *f = 0;
-    }
-    for(j = 0; j < 3; j++) {
+  for (int i = 0; i < stl->stats.number_of_facets; ++ i) {
+    //FIXME is the copy necessary?
+    stl_facet facet = stl->facet_start[i];
+    for (int j = 0; j < 3; j++) {
       if(stl->neighbors_start[i].neighbor[j] == -1) {
-        edge[j].facet_number = i;
-        edge[j].which_edge = j;
-        if(stl_load_edge_nearby(stl, &edge[j], &facet.vertex[j],
+        stl_hash_edge edge;
+        edge.facet_number = i;
+        edge.which_edge = j;
+        if(stl_load_edge_nearby(stl, &edge, &facet.vertex[j],
                                 &facet.vertex[(j + 1) % 3],
                                 tolerance)) {
           /* only insert edges that have different keys */
-          insert_hash_edge(stl, edge[j], stl_match_neighbors_nearby);
+          insert_hash_edge(stl, edge, stl_match_neighbors_nearby);
         }
       }
     }
@@ -315,27 +277,17 @@ stl_check_facets_nearby(stl_file *stl, float tolerance) {
   stl_free_edges(stl);
 }
 
-static int
-stl_load_edge_nearby(stl_file *stl, stl_hash_edge *edge,
-                     stl_vertex *a, stl_vertex *b, float tolerance) {
+static int stl_load_edge_nearby(stl_file *stl, stl_hash_edge *edge, stl_vertex *a, stl_vertex *b, float tolerance)
+{
   // Index of a grid cell spaced by tolerance.
-  uint32_t vertex1[3] = {
-    (uint32_t)((a->x - stl->stats.min.x) / tolerance),
-    (uint32_t)((a->y - stl->stats.min.y) / tolerance),
-    (uint32_t)((a->z - stl->stats.min.z) / tolerance)
-  };
-  uint32_t vertex2[3] = {
-    (uint32_t)((b->x - stl->stats.min.x) / tolerance),
-    (uint32_t)((b->y - stl->stats.min.y) / tolerance),
-    (uint32_t)((b->z - stl->stats.min.z) / tolerance)
-  };
+  typedef Eigen::Matrix<int32_t,  3, 1, Eigen::DontAlign> Vec3i;
+  Vec3i vertex1 = (*a / tolerance).cast<int32_t>();
+  Vec3i vertex2 = (*b / tolerance).cast<int32_t>();
+  static_assert(sizeof(Vec3i) == 12, "size of Vec3i incorrect");
 
-  if(   (vertex1[0] == vertex2[0])
-        && (vertex1[1] == vertex2[1])
-        && (vertex1[2] == vertex2[2])) {
-    /* Both vertices hash to the same value */
+  if (vertex1 == vertex2)
+    // Both vertices hash to the same value
     return 0;
-  }
 
   // Ensure identical vertex ordering of edges, which vertices land into equal grid cells.
   // This method is numerically robust.
@@ -344,30 +296,27 @@ stl_load_edge_nearby(stl_file *stl, stl_hash_edge *edge,
         ((vertex1[1] != vertex2[1]) ? 
             (vertex1[1] < vertex2[1]) : 
             (vertex1[2] < vertex2[2]))) {
-    memcpy(&edge->key[0], vertex1, sizeof(stl_vertex));
-    memcpy(&edge->key[3], vertex2, sizeof(stl_vertex));
+    memcpy(&edge->key[0],                  vertex1.data(), sizeof(stl_vertex));
+    memcpy(&edge->key[sizeof(stl_vertex)], vertex2.data(), sizeof(stl_vertex));
   } else {
-    memcpy(&edge->key[0], vertex2, sizeof(stl_vertex));
-    memcpy(&edge->key[3], vertex1, sizeof(stl_vertex));
+    memcpy(&edge->key[0],                  vertex2.data(), sizeof(stl_vertex));
+    memcpy(&edge->key[sizeof(stl_vertex)], vertex1.data(), sizeof(stl_vertex));
     edge->which_edge += 3; /* this edge is loaded backwards */
   }
   return 1;
 }
 
-static void
-stl_free_edges(stl_file *stl) {
-  int i;
-  stl_hash_edge *temp;
-
-  if (stl->error) return;
+static void stl_free_edges(stl_file *stl)
+{
+  if (stl->error)
+    return;
 
   if(stl->stats.malloced != stl->stats.freed) {
-    for(i = 0; i < stl->M; i++) {
-      for(temp = stl->heads[i]; stl->heads[i] != stl->tail;
-          temp = stl->heads[i]) {
+    for (int i = 0; i < stl->M; i++) {
+      for (stl_hash_edge *temp = stl->heads[i]; stl->heads[i] != stl->tail; temp = stl->heads[i]) {
         stl->heads[i] = stl->heads[i]->next;
         free(temp);
-        stl->stats.freed++;
+        ++ stl->stats.freed;
       }
     }
   }
@@ -375,8 +324,8 @@ stl_free_edges(stl_file *stl) {
   free(stl->tail);
 }
 
-static void
-stl_initialize_facet_check_nearby(stl_file *stl) {
+static void stl_initialize_facet_check_nearby(stl_file *stl)
+{
   int i;
 
   if (stl->error) return;
@@ -467,16 +416,8 @@ stl_record_neighbors(stl_file *stl,
   }
 }
 
-static void
-stl_match_neighbors_exact(stl_file *stl,
-                          stl_hash_edge *edge_a, stl_hash_edge *edge_b) {
-  if (stl->error) return;
-  stl_record_neighbors(stl, edge_a, edge_b);
-}
-
-static void
-stl_match_neighbors_nearby(stl_file *stl,
-                           stl_hash_edge *edge_a, stl_hash_edge *edge_b) {
+static void stl_match_neighbors_nearby(stl_file *stl, stl_hash_edge *edge_a, stl_hash_edge *edge_b)
+{
   int facet1;
   int facet2;
   int vertex1;
@@ -517,9 +458,7 @@ stl_match_neighbors_nearby(stl_file *stl,
 }
 
 
-static void
-stl_change_vertices(stl_file *stl, int facet_num, int vnot,
-                    stl_vertex new_vertex) {
+static void stl_change_vertices(stl_file *stl, int facet_num, int vnot, stl_vertex new_vertex) {
   int first_facet;
   int direction;
   int next_edge;
@@ -551,30 +490,30 @@ stl_change_vertices(stl_file *stl, int facet_num, int vnot,
       }
     }
 #if 0
-    if (stl->facet_start[facet_num].vertex[pivot_vertex].x == new_vertex.x &&
-        stl->facet_start[facet_num].vertex[pivot_vertex].y == new_vertex.y &&
-        stl->facet_start[facet_num].vertex[pivot_vertex].z == new_vertex.z)
+    if (stl->facet_start[facet_num].vertex[pivot_vertex](0) == new_vertex(0) &&
+        stl->facet_start[facet_num].vertex[pivot_vertex](1) == new_vertex(1) &&
+        stl->facet_start[facet_num].vertex[pivot_vertex](2) == new_vertex(2))
       printf("Changing vertex %f,%f,%f: Same !!!\r\n", 
-        new_vertex.x, new_vertex.y, new_vertex.z);
+        new_vertex(0), new_vertex(1), new_vertex(2));
     else {
-      if (stl->facet_start[facet_num].vertex[pivot_vertex].x != new_vertex.x)
+      if (stl->facet_start[facet_num].vertex[pivot_vertex](0) != new_vertex(0))
         printf("Changing coordinate x, vertex %e (0x%08x) to %e(0x%08x)\r\n", 
-          stl->facet_start[facet_num].vertex[pivot_vertex].x,
-          *reinterpret_cast<const int*>(&stl->facet_start[facet_num].vertex[pivot_vertex].x),
-          new_vertex.x,
-          *reinterpret_cast<const int*>(&new_vertex.x));
-      if (stl->facet_start[facet_num].vertex[pivot_vertex].y != new_vertex.y)
+          stl->facet_start[facet_num].vertex[pivot_vertex](0),
+          *reinterpret_cast<const int*>(&stl->facet_start[facet_num].vertex[pivot_vertex](0)),
+          new_vertex(0),
+          *reinterpret_cast<const int*>(&new_vertex(0)));
+      if (stl->facet_start[facet_num].vertex[pivot_vertex](1) != new_vertex(1))
         printf("Changing coordinate x, vertex %e (0x%08x) to %e(0x%08x)\r\n", 
-          stl->facet_start[facet_num].vertex[pivot_vertex].y,
-          *reinterpret_cast<const int*>(&stl->facet_start[facet_num].vertex[pivot_vertex].y),
-          new_vertex.y,
-          *reinterpret_cast<const int*>(&new_vertex.y));
-      if (stl->facet_start[facet_num].vertex[pivot_vertex].z != new_vertex.z)
+          stl->facet_start[facet_num].vertex[pivot_vertex](1),
+          *reinterpret_cast<const int*>(&stl->facet_start[facet_num].vertex[pivot_vertex](1)),
+          new_vertex(1),
+          *reinterpret_cast<const int*>(&new_vertex(1)));
+      if (stl->facet_start[facet_num].vertex[pivot_vertex](2) != new_vertex(2))
         printf("Changing coordinate x, vertex %e (0x%08x) to %e(0x%08x)\r\n", 
-          stl->facet_start[facet_num].vertex[pivot_vertex].z,
-          *reinterpret_cast<const int*>(&stl->facet_start[facet_num].vertex[pivot_vertex].z),
-          new_vertex.z,
-          *reinterpret_cast<const int*>(&new_vertex.z));
+          stl->facet_start[facet_num].vertex[pivot_vertex](2),
+          *reinterpret_cast<const int*>(&stl->facet_start[facet_num].vertex[pivot_vertex](2)),
+          new_vertex(2),
+          *reinterpret_cast<const int*>(&new_vertex(2)));
     }
 #endif
     stl->facet_start[facet_num].vertex[pivot_vertex] = new_vertex;
@@ -595,7 +534,6 @@ Try using a smaller tolerance or don't do a nearby check\n");
   }
 }
 
-
 static void
 stl_which_vertices_to_change(stl_file *stl, stl_hash_edge *edge_a,
                              stl_hash_edge *edge_b, int *facet1, int *vertex1,
@@ -622,11 +560,10 @@ stl_which_vertices_to_change(stl_file *stl, stl_hash_edge *edge_a,
     v1b = (edge_b->which_edge + 1) % 3;
   }
 
-  /* Of the first pair, which vertex, if any, should be changed */
-  if(!memcmp(&stl->facet_start[edge_a->facet_number].vertex[v1a],
-             &stl->facet_start[edge_b->facet_number].vertex[v1b],
-             sizeof(stl_vertex))) {
-    /* These facets are already equal.  No need to change. */
+  // Of the first pair, which vertex, if any, should be changed
+  if(stl->facet_start[edge_a->facet_number].vertex[v1a] == 
+     stl->facet_start[edge_b->facet_number].vertex[v1b]) {
+    // These facets are already equal.  No need to change.
     *facet1 = -1;
   } else {
     if(   (stl->neighbors_start[edge_a->facet_number].neighbor[v1a] == -1)
@@ -644,10 +581,9 @@ stl_which_vertices_to_change(stl_file *stl, stl_hash_edge *edge_a,
   }
 
   /* Of the second pair, which vertex, if any, should be changed */
-  if(!memcmp(&stl->facet_start[edge_a->facet_number].vertex[v2a],
-             &stl->facet_start[edge_b->facet_number].vertex[v2b],
-             sizeof(stl_vertex))) {
-    /* These facets are already equal.  No need to change. */
+  if(stl->facet_start[edge_a->facet_number].vertex[v2a] == 
+     stl->facet_start[edge_b->facet_number].vertex[v2b]) {
+    // These facets are already equal.  No need to change.
     *facet2 = -1;
   } else {
     if(   (stl->neighbors_start[edge_a->facet_number].neighbor[v2a] == -1)
@@ -718,40 +654,35 @@ in stl_remove_facet: neighbor = %d numfacets = %d this is wrong\n",
   }
 }
 
-void
-stl_remove_unconnected_facets(stl_file *stl) {
+void stl_remove_unconnected_facets(stl_file *stl)
+{
   /* A couple of things need to be done here.  One is to remove any */
   /* completely unconnected facets (0 edges connected) since these are */
   /* useless and could be completely wrong.   The second thing that needs to */
   /* be done is to remove any degenerate facets that were created during */
   /* stl_check_facets_nearby(). */
+  if (stl->error)
+    return;
 
-  int i;
-
-  if (stl->error) return;
-
-  /* remove degenerate facets */
-  for(i = 0; i < stl->stats.number_of_facets; i++) {
-    if(   !memcmp(&stl->facet_start[i].vertex[0],
-                  &stl->facet_start[i].vertex[1], sizeof(stl_vertex))
-          || !memcmp(&stl->facet_start[i].vertex[1],
-                     &stl->facet_start[i].vertex[2], sizeof(stl_vertex))
-          || !memcmp(&stl->facet_start[i].vertex[0],
-                     &stl->facet_start[i].vertex[2], sizeof(stl_vertex))) {
+  // remove degenerate facets
+  for (int i = 0; i < stl->stats.number_of_facets; ++ i) {
+    if(stl->facet_start[i].vertex[0] == stl->facet_start[i].vertex[1] ||
+       stl->facet_start[i].vertex[0] == stl->facet_start[i].vertex[2] ||
+       stl->facet_start[i].vertex[1] == stl->facet_start[i].vertex[2]) {
       stl_remove_degenerate(stl, i);
       i--;
     }
   }
 
   if(stl->stats.connected_facets_1_edge < stl->stats.number_of_facets) {
-    /* remove completely unconnected facets */
-    for(i = 0; i < stl->stats.number_of_facets; i++) {
-      if(   (stl->neighbors_start[i].neighbor[0] == -1)
-            && (stl->neighbors_start[i].neighbor[1] == -1)
-            && (stl->neighbors_start[i].neighbor[2] == -1)) {
-        /* This facet is completely unconnected.  Remove it. */
+    // remove completely unconnected facets
+    for (int i = 0; i < stl->stats.number_of_facets; i++) {
+      if (stl->neighbors_start[i].neighbor[0] == -1 &&
+          stl->neighbors_start[i].neighbor[1] == -1 &&
+          stl->neighbors_start[i].neighbor[2] == -1) {
+        // This facet is completely unconnected.  Remove it.
         stl_remove_facet(stl, i);
-        i--;
+        -- i;
       }
     }
   }
@@ -771,30 +702,24 @@ stl_remove_degenerate(stl_file *stl, int facet) {
 
   if (stl->error) return;
 
-  if(   !memcmp(&stl->facet_start[facet].vertex[0],
-                &stl->facet_start[facet].vertex[1], sizeof(stl_vertex))
-        && !memcmp(&stl->facet_start[facet].vertex[1],
-                   &stl->facet_start[facet].vertex[2], sizeof(stl_vertex))) {
+  if (stl->facet_start[facet].vertex[0] == stl->facet_start[facet].vertex[1] &&
+      stl->facet_start[facet].vertex[1] == stl->facet_start[facet].vertex[2]) {
     /* all 3 vertices are equal.  Just remove the facet.  I don't think*/
     /* this is really possible, but just in case... */
     printf("removing a facet in stl_remove_degenerate\n");
-
     stl_remove_facet(stl, facet);
     return;
   }
 
-  if(!memcmp(&stl->facet_start[facet].vertex[0],
-             &stl->facet_start[facet].vertex[1], sizeof(stl_vertex))) {
+  if (stl->facet_start[facet].vertex[0] == stl->facet_start[facet].vertex[1]) {
     edge1 = 1;
     edge2 = 2;
     edge3 = 0;
-  } else if(!memcmp(&stl->facet_start[facet].vertex[1],
-                    &stl->facet_start[facet].vertex[2], sizeof(stl_vertex))) {
+  } else if (stl->facet_start[facet].vertex[1] == stl->facet_start[facet].vertex[2]) {
     edge1 = 0;
     edge2 = 2;
     edge3 = 1;
-  } else if(!memcmp(&stl->facet_start[facet].vertex[2],
-                    &stl->facet_start[facet].vertex[0], sizeof(stl_vertex))) {
+  } else if (stl->facet_start[facet].vertex[2] == stl->facet_start[facet].vertex[0]) {
     edge1 = 0;
     edge2 = 1;
     edge3 = 2;
@@ -883,7 +808,7 @@ stl_fill_holes(stl_file *stl) {
       stl_load_edge_exact(stl, &edge, &facet.vertex[j],
                           &facet.vertex[(j + 1) % 3]);
 
-      insert_hash_edge(stl, edge, stl_match_neighbors_exact);
+      insert_hash_edge(stl, edge, stl_record_neighbors);
     }
   }
 
@@ -939,7 +864,7 @@ stl_fill_holes(stl_file *stl) {
             stl_load_edge_exact(stl, &edge, &new_facet.vertex[k],
                                 &new_facet.vertex[(k + 1) % 3]);
 
-            insert_hash_edge(stl, edge, stl_match_neighbors_exact);
+            insert_hash_edge(stl, edge, stl_record_neighbors);
           }
           break;
         } else {
@@ -977,9 +902,7 @@ stl_add_facet(stl_file *stl, stl_facet *new_facet) {
   stl->facet_start[stl->stats.number_of_facets] = *new_facet;
 
   /* note that the normal vector is not set here, just initialized to 0 */
-  stl->facet_start[stl->stats.number_of_facets].normal.x = 0.0;
-  stl->facet_start[stl->stats.number_of_facets].normal.y = 0.0;
-  stl->facet_start[stl->stats.number_of_facets].normal.z = 0.0;
+  stl->facet_start[stl->stats.number_of_facets].normal = stl_normal::Zero();
 
   stl->neighbors_start[stl->stats.number_of_facets].neighbor[0] = -1;
   stl->neighbors_start[stl->stats.number_of_facets].neighbor[1] = -1;
diff --git a/xs/src/admesh/normals.cpp b/xs/src/admesh/normals.cpp
index b7cf9a8ab..a8faa44bd 100644
--- a/xs/src/admesh/normals.cpp
+++ b/xs/src/admesh/normals.cpp
@@ -27,12 +27,6 @@
 
 #include "stl.h"
 
-static void stl_reverse_vector(float v[]) {
-  v[0] *= -1;
-  v[1] *= -1;
-  v[2] *= -1;
-}
-
 static int stl_check_normal_vector(stl_file *stl, int facet_num, int normal_fix_flag);
 
 static void
@@ -228,102 +222,52 @@ static int stl_check_normal_vector(stl_file *stl, int facet_num, int normal_fix_
   /* Returns 2 if the normal is not within tolerance and backwards */
   /* Returns 4 if the status is unknown. */
 
-  float normal[3];
-  float test_norm[3];
   stl_facet *facet;
 
   facet = &stl->facet_start[facet_num];
 
+  stl_normal normal;
   stl_calculate_normal(normal, facet);
   stl_normalize_vector(normal);
+  stl_normal normal_dif = (normal - facet->normal).cwiseAbs();
 
-  if(   (ABS(normal[0] - facet->normal.x) < 0.001)
-        && (ABS(normal[1] - facet->normal.y) < 0.001)
-        && (ABS(normal[2] - facet->normal.z) < 0.001)) {
+  const float eps = 0.001f;
+  if (normal_dif(0) < eps && normal_dif(1) < eps && normal_dif(2) < eps) {
     /* It is not really necessary to change the values here */
     /* but just for consistency, I will. */
-    facet->normal.x = normal[0];
-    facet->normal.y = normal[1];
-    facet->normal.z = normal[2];
+    facet->normal = normal;
     return 0;
   }
 
-  test_norm[0] = facet->normal.x;
-  test_norm[1] = facet->normal.y;
-  test_norm[2] = facet->normal.z;
-
+  stl_normal test_norm = facet->normal;
   stl_normalize_vector(test_norm);
-  if(   (ABS(normal[0] - test_norm[0]) < 0.001)
-        && (ABS(normal[1] - test_norm[1]) < 0.001)
-        && (ABS(normal[2] - test_norm[2]) < 0.001)) {
+  normal_dif = (normal - test_norm).cwiseAbs();
+  if (normal_dif(0) < eps && normal_dif(1) < eps && normal_dif(2) < eps) {
     if(normal_fix_flag) {
-      facet->normal.x = normal[0];
-      facet->normal.y = normal[1];
-      facet->normal.z = normal[2];
+      facet->normal = normal;
       stl->stats.normals_fixed += 1;
     }
     return 1;
   }
 
-  stl_reverse_vector(test_norm);
-  if(   (ABS(normal[0] - test_norm[0]) < 0.001)
-        && (ABS(normal[1] - test_norm[1]) < 0.001)
-        && (ABS(normal[2] - test_norm[2]) < 0.001)) {
-    /* Facet is backwards. */
+  test_norm *= -1.f;
+  normal_dif = (normal - test_norm).cwiseAbs();
+  if (normal_dif(0) < eps && normal_dif(1) < eps && normal_dif(2) < eps) {
+    // Facet is backwards.
     if(normal_fix_flag) {
-      facet->normal.x = normal[0];
-      facet->normal.y = normal[1];
-      facet->normal.z = normal[2];
+      facet->normal = normal;
       stl->stats.normals_fixed += 1;
     }
     return 2;
   }
   if(normal_fix_flag) {
-    facet->normal.x = normal[0];
-    facet->normal.y = normal[1];
-    facet->normal.z = normal[2];
+    facet->normal = normal;
     stl->stats.normals_fixed += 1;
   }
   return 4;
 }
 
-void stl_calculate_normal(float normal[], stl_facet *facet) {
-  float v1[3] = {
-    facet->vertex[1].x - facet->vertex[0].x,
-    facet->vertex[1].y - facet->vertex[0].y,
-    facet->vertex[1].z - facet->vertex[0].z
-  };
-  float v2[3] = {
-    facet->vertex[2].x - facet->vertex[0].x,
-    facet->vertex[2].y - facet->vertex[0].y,
-    facet->vertex[2].z - facet->vertex[0].z
-  };
-  normal[0] = (float)((double)v1[1] * (double)v2[2]) - ((double)v1[2] * (double)v2[1]);
-  normal[1] = (float)((double)v1[2] * (double)v2[0]) - ((double)v1[0] * (double)v2[2]);
-  normal[2] = (float)((double)v1[0] * (double)v2[1]) - ((double)v1[1] * (double)v2[0]);
-}
-
-void stl_normalize_vector(float v[]) {
-  double length;
-  double factor;
-  float min_normal_length;
-
-  length = sqrt((double)v[0] * (double)v[0] + (double)v[1] * (double)v[1] + (double)v[2] * (double)v[2]);
-  min_normal_length = 0.000000000001;
-  if(length < min_normal_length) {
-    v[0] = 0.0;
-    v[1] = 0.0;
-    v[2] = 0.0;
-    return;
-  }
-  factor = 1.0 / length;
-  v[0] *= factor;
-  v[1] *= factor;
-  v[2] *= factor;
-}
-
-void
-stl_fix_normal_values(stl_file *stl) {
+void stl_fix_normal_values(stl_file *stl) {
   int i;
 
   if (stl->error) return;
@@ -333,20 +277,16 @@ stl_fix_normal_values(stl_file *stl) {
   }
 }
 
-void
-stl_reverse_all_facets(stl_file *stl) {
-  int i;
-  float normal[3];
+void stl_reverse_all_facets(stl_file *stl)
+{
+  if (stl->error)
+  	return;
 
-  if (stl->error) return;
-
-  for(i = 0; i < stl->stats.number_of_facets; i++) {
+  stl_normal normal;
+  for(int i = 0; i < stl->stats.number_of_facets; i++) {
     stl_reverse_facet(stl, i);
     stl_calculate_normal(normal, &stl->facet_start[i]);
     stl_normalize_vector(normal);
-    stl->facet_start[i].normal.x = normal[0];
-    stl->facet_start[i].normal.y = normal[1];
-    stl->facet_start[i].normal.z = normal[2];
+    stl->facet_start[i].normal = normal;
   }
 }
-
diff --git a/xs/src/admesh/shared.cpp b/xs/src/admesh/shared.cpp
index 8080f3574..91bb82e00 100644
--- a/xs/src/admesh/shared.cpp
+++ b/xs/src/admesh/shared.cpp
@@ -169,7 +169,7 @@ stl_write_off(stl_file *stl, char *file) {
 
   for(i = 0; i < stl->stats.shared_vertices; i++) {
     fprintf(fp, "\t%f %f %f\n",
-            stl->v_shared[i].x, stl->v_shared[i].y, stl->v_shared[i].z);
+            stl->v_shared[i](0), stl->v_shared[i](1), stl->v_shared[i](2));
   }
   for(i = 0; i < stl->stats.number_of_facets; i++) {
     fprintf(fp, "\t3 %d %d %d\n", stl->v_indices[i].vertex[0],
@@ -216,10 +216,10 @@ stl_write_vrml(stl_file *stl, char *file) {
 
   for(i = 0; i < (stl->stats.shared_vertices - 1); i++) {
     fprintf(fp, "\t\t\t\t%f %f %f,\n",
-            stl->v_shared[i].x, stl->v_shared[i].y, stl->v_shared[i].z);
+            stl->v_shared[i](0), stl->v_shared[i](1), stl->v_shared[i](2));
   }
   fprintf(fp, "\t\t\t\t%f %f %f]\n",
-          stl->v_shared[i].x, stl->v_shared[i].y, stl->v_shared[i].z);
+          stl->v_shared[i](0), stl->v_shared[i](1), stl->v_shared[i](2));
   fprintf(fp, "\t\t}\n");
   fprintf(fp, "\t\tDEF STLTriangles IndexedFaceSet {\n");
   fprintf(fp, "\t\t\tcoordIndex [\n");
@@ -254,7 +254,7 @@ void stl_write_obj (stl_file *stl, char *file) {
   }
 
   for (i = 0; i < stl->stats.shared_vertices; i++) {
-    fprintf(fp, "v %f %f %f\n", stl->v_shared[i].x, stl->v_shared[i].y, stl->v_shared[i].z);
+    fprintf(fp, "v %f %f %f\n", stl->v_shared[i](0), stl->v_shared[i](1), stl->v_shared[i](2));
   }
   for (i = 0; i < stl->stats.number_of_facets; i++) {
     fprintf(fp, "f %d %d %d\n", stl->v_indices[i].vertex[0]+1, stl->v_indices[i].vertex[1]+1, stl->v_indices[i].vertex[2]+1);
diff --git a/xs/src/admesh/stl.h b/xs/src/admesh/stl.h
index dc0f7ae19..5f7a3c5c1 100644
--- a/xs/src/admesh/stl.h
+++ b/xs/src/admesh/stl.h
@@ -27,9 +27,7 @@
 #include <stdint.h>
 #include <stddef.h>
 
-#define STL_MAX(A,B) ((A)>(B)? (A):(B))
-#define STL_MIN(A,B) ((A)<(B)? (A):(B))
-#define ABS(X)  ((X) < 0 ? -(X) : (X))
+#include <Eigen/Geometry> 
 
 // Size of the binary STL header, free form.
 #define LABEL_SIZE             80
@@ -39,31 +37,16 @@
 #define HEADER_SIZE            84
 #define STL_MIN_FILE_SIZE      284
 #define ASCII_LINES_PER_FACET  7
-// Comparing an edge by memcmp, 2x3x4 bytes = 24
-#define SIZEOF_EDGE_SORT       24
-
-typedef struct {
-  float x;
-  float y;
-  float z;
-} stl_vertex;
 
+typedef Eigen::Matrix<float, 3, 1, Eigen::DontAlign> stl_vertex;
+typedef Eigen::Matrix<float, 3, 1, Eigen::DontAlign> stl_normal;
 static_assert(sizeof(stl_vertex) == 12, "size of stl_vertex incorrect");
-
-typedef struct {
-  float x;
-  float y;
-  float z;
-} stl_normal;
-
 static_assert(sizeof(stl_normal) == 12, "size of stl_normal incorrect");
 
-typedef char stl_extra[2];
-
 typedef struct {
   stl_normal normal;
   stl_vertex vertex[3];
-  stl_extra  extra;
+  char       extra[2];
 } stl_facet;
 #define SIZEOF_STL_FACET       50
 
@@ -81,8 +64,12 @@ typedef struct {
 } stl_edge;
 
 typedef struct stl_hash_edge {
-  // Key of a hash edge: 2x binary copy of a floating point vertex.
-  uint32_t       key[6];
+  // Key of a hash edge: sorted vertices of the edge.
+  unsigned char key[2 * sizeof(stl_vertex)];
+  // Compare two keys.
+  bool operator==(const stl_hash_edge &rhs) { return memcmp(key, rhs.key, sizeof(key)) == 0; }
+  bool operator!=(const stl_hash_edge &rhs) { return ! (*this == rhs); }
+  int  hash(int M) const { return ((key[0] / 23 + key[1] / 19 + key[2] / 17 + key[3] /13  + key[4] / 11 + key[5] / 7 ) % M); }
   // Index of a facet owning this edge.
   int            facet_number;
   // Index of this edge inside the facet with an index of facet_number.
@@ -91,8 +78,6 @@ typedef struct stl_hash_edge {
   struct stl_hash_edge  *next;
 } stl_hash_edge;
 
-static_assert(offsetof(stl_hash_edge, facet_number) == SIZEOF_EDGE_SORT, "size of stl_hash_edge.key incorrect");
-
 typedef struct {
   // Index of a neighbor facet.
   int   neighbor[3];
@@ -179,8 +164,8 @@ extern void stl_fix_normal_values(stl_file *stl);
 extern void stl_reverse_all_facets(stl_file *stl);
 extern void stl_translate(stl_file *stl, float x, float y, float z);
 extern void stl_translate_relative(stl_file *stl, float x, float y, float z);
-extern void stl_scale_versor(stl_file *stl, float versor[3]);
-extern void stl_scale(stl_file *stl, float factor);
+extern void stl_scale_versor(stl_file *stl, const stl_vertex &versor);
+inline void stl_scale(stl_file *stl, float factor) { stl_scale_versor(stl, stl_vertex(factor, factor, factor)); }
 extern void stl_rotate_x(stl_file *stl, float angle);
 extern void stl_rotate_y(stl_file *stl, float angle);
 extern void stl_rotate_z(stl_file *stl, float angle);
@@ -195,8 +180,20 @@ extern void stl_write_obj(stl_file *stl, char *file);
 extern void stl_write_off(stl_file *stl, char *file);
 extern void stl_write_dxf(stl_file *stl, char *file, char *label);
 extern void stl_write_vrml(stl_file *stl, char *file);
-extern void stl_calculate_normal(float normal[], stl_facet *facet);
-extern void stl_normalize_vector(float v[]);
+inline void stl_calculate_normal(stl_normal &normal, stl_facet *facet) {
+  normal = (facet->vertex[1] - facet->vertex[0]).cross(facet->vertex[2] - facet->vertex[0]);
+}
+inline void stl_normalize_vector(stl_normal &normal) {
+  double length = normal.cast<double>().norm();
+  if (length < 0.000000000001)
+    normal = stl_normal::Zero();
+  else
+    normal *= (1.0 / length);
+}
+inline bool stl_vertex_lower(const stl_vertex &a, const stl_vertex &b) {
+  return (a(0) != b(0)) ? (a(0) < b(0)) :
+        ((a(1) != b(1)) ? (a(1) < b(1)) : (a(2) < b(2)));
+}
 extern void stl_calculate_volume(stl_file *stl);
 
 extern void stl_repair(stl_file *stl, int fixall_flag, int exact_flag, int tolerance_flag, float tolerance, int increment_flag, float increment, int nearby_flag, int iterations, int remove_unconnected_flag, int fill_holes_flag, int normal_directions_flag, int normal_values_flag, int reverse_all_flag, int verbose_flag);
@@ -204,8 +201,8 @@ extern void stl_repair(stl_file *stl, int fixall_flag, int exact_flag, int toler
 extern void stl_initialize(stl_file *stl);
 extern void stl_count_facets(stl_file *stl, const char *file);
 extern void stl_allocate(stl_file *stl);
-extern void stl_read(stl_file *stl, int first_facet, int first);
-extern void stl_facet_stats(stl_file *stl, stl_facet facet, int first);
+extern void stl_read(stl_file *stl, int first_facet, bool first);
+extern void stl_facet_stats(stl_file *stl, stl_facet facet, bool &first);
 extern void stl_reallocate(stl_file *stl);
 extern void stl_add_facet(stl_file *stl, stl_facet *new_facet);
 extern void stl_get_size(stl_file *stl);
diff --git a/xs/src/admesh/stl_io.cpp b/xs/src/admesh/stl_io.cpp
index 1603981fc..81e29d286 100644
--- a/xs/src/admesh/stl_io.cpp
+++ b/xs/src/admesh/stl_io.cpp
@@ -44,9 +44,9 @@ stl_print_edges(stl_file *stl, FILE *file) {
   for(i = 0; i < edges_allocated; i++) {
     fprintf(file, "%d, %f, %f, %f, %f, %f, %f\n",
             stl->edge_start[i].facet_number,
-            stl->edge_start[i].p1.x, stl->edge_start[i].p1.y,
-            stl->edge_start[i].p1.z, stl->edge_start[i].p2.x,
-            stl->edge_start[i].p2.y, stl->edge_start[i].p2.z);
+            stl->edge_start[i].p1(0), stl->edge_start[i].p1(1),
+            stl->edge_start[i].p1(2), stl->edge_start[i].p2(0),
+            stl->edge_start[i].p2(1), stl->edge_start[i].p2(2));
   }
 }
 
@@ -75,11 +75,11 @@ File type          : ASCII STL file\n");
 Header             : %s\n", stl->stats.header);
   fprintf(file, "============== Size ==============\n");
   fprintf(file, "Min X = % f, Max X = % f\n",
-          stl->stats.min.x, stl->stats.max.x);
+          stl->stats.min(0), stl->stats.max(0));
   fprintf(file, "Min Y = % f, Max Y = % f\n",
-          stl->stats.min.y, stl->stats.max.y);
+          stl->stats.min(1), stl->stats.max(1));
   fprintf(file, "Min Z = % f, Max Z = % f\n",
-          stl->stats.min.z, stl->stats.max.z);
+          stl->stats.min(2), stl->stats.max(2));
 
   fprintf(file, "\
 ========= Facet Status ========== Original ============ Final ====\n");
@@ -149,18 +149,18 @@ stl_write_ascii(stl_file *stl, const char *file, const char *label) {
 
   for(i = 0; i < stl->stats.number_of_facets; i++) {
     fprintf(fp, "  facet normal % .8E % .8E % .8E\n",
-            stl->facet_start[i].normal.x, stl->facet_start[i].normal.y,
-            stl->facet_start[i].normal.z);
+            stl->facet_start[i].normal(0), stl->facet_start[i].normal(1),
+            stl->facet_start[i].normal(2));
     fprintf(fp, "    outer loop\n");
     fprintf(fp, "      vertex % .8E % .8E % .8E\n",
-            stl->facet_start[i].vertex[0].x, stl->facet_start[i].vertex[0].y,
-            stl->facet_start[i].vertex[0].z);
+            stl->facet_start[i].vertex[0](0), stl->facet_start[i].vertex[0](1),
+            stl->facet_start[i].vertex[0](2));
     fprintf(fp, "      vertex % .8E % .8E % .8E\n",
-            stl->facet_start[i].vertex[1].x, stl->facet_start[i].vertex[1].y,
-            stl->facet_start[i].vertex[1].z);
+            stl->facet_start[i].vertex[1](0), stl->facet_start[i].vertex[1](1),
+            stl->facet_start[i].vertex[1](2));
     fprintf(fp, "      vertex % .8E % .8E % .8E\n",
-            stl->facet_start[i].vertex[2].x, stl->facet_start[i].vertex[2].y,
-            stl->facet_start[i].vertex[2].z);
+            stl->facet_start[i].vertex[2](0), stl->facet_start[i].vertex[2](1),
+            stl->facet_start[i].vertex[2](2));
     fprintf(fp, "    endloop\n");
     fprintf(fp, "  endfacet\n");
   }
@@ -264,9 +264,9 @@ void
 stl_write_vertex(stl_file *stl, int facet, int vertex) {
   if (stl->error) return;
   printf("  vertex %d/%d % .8E % .8E % .8E\n", vertex, facet,
-         stl->facet_start[facet].vertex[vertex].x,
-         stl->facet_start[facet].vertex[vertex].y,
-         stl->facet_start[facet].vertex[vertex].z);
+         stl->facet_start[facet].vertex[vertex](0),
+         stl->facet_start[facet].vertex[vertex](1),
+         stl->facet_start[facet].vertex[vertex](2));
 }
 
 void
@@ -309,10 +309,10 @@ stl_write_quad_object(stl_file *stl, char *file) {
   int       i;
   int       j;
   char      *error_msg;
-  stl_vertex connect_color;
-  stl_vertex uncon_1_color;
-  stl_vertex uncon_2_color;
-  stl_vertex uncon_3_color;
+  stl_vertex connect_color = stl_vertex::Zero();
+  stl_vertex uncon_1_color = stl_vertex::Zero();
+  stl_vertex uncon_2_color = stl_vertex::Zero();
+  stl_vertex uncon_3_color = stl_vertex::Zero();
   stl_vertex color;
 
   if (stl->error) return;
@@ -330,19 +330,6 @@ stl_write_quad_object(stl_file *stl, char *file) {
     return;
   }
 
-  connect_color.x = 0.0;
-  connect_color.y = 0.0;
-  connect_color.z = 1.0;
-  uncon_1_color.x = 0.0;
-  uncon_1_color.y = 1.0;
-  uncon_1_color.z = 0.0;
-  uncon_2_color.x = 1.0;
-  uncon_2_color.y = 1.0;
-  uncon_2_color.z = 1.0;
-  uncon_3_color.x = 1.0;
-  uncon_3_color.y = 0.0;
-  uncon_3_color.z = 0.0;
-
   fprintf(fp, "CQUAD\n");
   for(i = 0; i < stl->stats.number_of_facets; i++) {
     j = ((stl->neighbors_start[i].neighbor[0] == -1) +
@@ -358,21 +345,21 @@ stl_write_quad_object(stl_file *stl, char *file) {
       color = uncon_3_color;
     }
     fprintf(fp, "%f %f %f    %1.1f %1.1f %1.1f 1\n",
-            stl->facet_start[i].vertex[0].x,
-            stl->facet_start[i].vertex[0].y,
-            stl->facet_start[i].vertex[0].z, color.x, color.y, color.z);
+            stl->facet_start[i].vertex[0](0),
+            stl->facet_start[i].vertex[0](1),
+            stl->facet_start[i].vertex[0](2), color(0), color(1), color(2));
     fprintf(fp, "%f %f %f    %1.1f %1.1f %1.1f 1\n",
-            stl->facet_start[i].vertex[1].x,
-            stl->facet_start[i].vertex[1].y,
-            stl->facet_start[i].vertex[1].z, color.x, color.y, color.z);
+            stl->facet_start[i].vertex[1](0),
+            stl->facet_start[i].vertex[1](1),
+            stl->facet_start[i].vertex[1](2), color(0), color(1), color(2));
     fprintf(fp, "%f %f %f    %1.1f %1.1f %1.1f 1\n",
-            stl->facet_start[i].vertex[2].x,
-            stl->facet_start[i].vertex[2].y,
-            stl->facet_start[i].vertex[2].z, color.x, color.y, color.z);
+            stl->facet_start[i].vertex[2](0),
+            stl->facet_start[i].vertex[2](1),
+            stl->facet_start[i].vertex[2](2), color(0), color(1), color(2));
     fprintf(fp, "%f %f %f    %1.1f %1.1f %1.1f 1\n",
-            stl->facet_start[i].vertex[2].x,
-            stl->facet_start[i].vertex[2].y,
-            stl->facet_start[i].vertex[2].z, color.x, color.y, color.z);
+            stl->facet_start[i].vertex[2](0),
+            stl->facet_start[i].vertex[2](1),
+            stl->facet_start[i].vertex[2](2), color(0), color(1), color(2));
   }
   fclose(fp);
 }
@@ -409,17 +396,17 @@ stl_write_dxf(stl_file *stl, char *file, char *label) {
   for(i = 0; i < stl->stats.number_of_facets; i++) {
     fprintf(fp, "0\n3DFACE\n8\n0\n");
     fprintf(fp, "10\n%f\n20\n%f\n30\n%f\n",
-            stl->facet_start[i].vertex[0].x, stl->facet_start[i].vertex[0].y,
-            stl->facet_start[i].vertex[0].z);
+            stl->facet_start[i].vertex[0](0), stl->facet_start[i].vertex[0](1),
+            stl->facet_start[i].vertex[0](2));
     fprintf(fp, "11\n%f\n21\n%f\n31\n%f\n",
-            stl->facet_start[i].vertex[1].x, stl->facet_start[i].vertex[1].y,
-            stl->facet_start[i].vertex[1].z);
+            stl->facet_start[i].vertex[1](0), stl->facet_start[i].vertex[1](1),
+            stl->facet_start[i].vertex[1](2));
     fprintf(fp, "12\n%f\n22\n%f\n32\n%f\n",
-            stl->facet_start[i].vertex[2].x, stl->facet_start[i].vertex[2].y,
-            stl->facet_start[i].vertex[2].z);
+            stl->facet_start[i].vertex[2](0), stl->facet_start[i].vertex[2](1),
+            stl->facet_start[i].vertex[2](2));
     fprintf(fp, "13\n%f\n23\n%f\n33\n%f\n",
-            stl->facet_start[i].vertex[2].x, stl->facet_start[i].vertex[2].y,
-            stl->facet_start[i].vertex[2].z);
+            stl->facet_start[i].vertex[2](0), stl->facet_start[i].vertex[2](1),
+            stl->facet_start[i].vertex[2](2));
   }
 
   fprintf(fp, "0\nENDSEC\n0\nEOF\n");
diff --git a/xs/src/admesh/stlinit.cpp b/xs/src/admesh/stlinit.cpp
index e572ce930..6bb4bf633 100644
--- a/xs/src/admesh/stlinit.cpp
+++ b/xs/src/admesh/stlinit.cpp
@@ -40,7 +40,7 @@ stl_open(stl_file *stl, const char *file) {
   stl_initialize(stl);
   stl_count_facets(stl, file);
   stl_allocate(stl);
-  stl_read(stl, 0, 1);
+  stl_read(stl, 0, true);
   if (!stl->error) fclose(stl->fp);
 }
 
@@ -227,7 +227,7 @@ stl_open_merge(stl_file *stl, char *file_to_merge) {
      Start at num_facets_so_far, the index to the first unused facet.  Also say
      that this isn't our first time so we should augment stats like min and max
      instead of erasing them. */
-  stl_read(stl, num_facets_so_far, 0);
+  stl_read(stl, num_facets_so_far, false);
 
   /* Restore the stl information we overwrote (for stl_read) so that it still accurately
      reflects the subject part: */
@@ -255,8 +255,7 @@ stl_reallocate(stl_file *stl) {
 /* Reads the contents of the file pointed to by stl->fp into the stl structure,
    starting at facet first_facet.  The second argument says if it's our first
    time running this for the stl and therefore we should reset our max and min stats. */
-void
-stl_read(stl_file *stl, int first_facet, int first) {
+void stl_read(stl_file *stl, int first_facet, bool first) {
   stl_facet facet;
   int   i;
 
@@ -294,11 +293,11 @@ stl_read(stl_file *stl, int first_facet, int first) {
       assert(res_normal == 3);
       int res_outer_loop = fscanf(stl->fp, " outer loop");
       assert(res_outer_loop == 0);
-      int res_vertex1    = fscanf(stl->fp, " vertex %f %f %f", &facet.vertex[0].x, &facet.vertex[0].y, &facet.vertex[0].z);
+      int res_vertex1    = fscanf(stl->fp, " vertex %f %f %f", &facet.vertex[0](0), &facet.vertex[0](1), &facet.vertex[0](2));
       assert(res_vertex1 == 3);
-      int res_vertex2    = fscanf(stl->fp, " vertex %f %f %f", &facet.vertex[1].x, &facet.vertex[1].y, &facet.vertex[1].z);
+      int res_vertex2    = fscanf(stl->fp, " vertex %f %f %f", &facet.vertex[1](0), &facet.vertex[1](1), &facet.vertex[1](2));
       assert(res_vertex2 == 3);
-      int res_vertex3    = fscanf(stl->fp, " vertex %f %f %f", &facet.vertex[2].x, &facet.vertex[2].y, &facet.vertex[2].z);
+      int res_vertex3    = fscanf(stl->fp, " vertex %f %f %f", &facet.vertex[2](0), &facet.vertex[2](1), &facet.vertex[2](2));
       assert(res_vertex3 == 3);
       int res_endloop    = fscanf(stl->fp, " endloop");
       assert(res_endloop == 0);
@@ -311,9 +310,9 @@ stl_read(stl_file *stl, int first_facet, int first) {
       }
 
       // The facet normal has been parsed as a single string as to workaround for not a numbers in the normal definition.
-	  if (sscanf(normal_buf[0], "%f", &facet.normal.x) != 1 ||
-		  sscanf(normal_buf[1], "%f", &facet.normal.y) != 1 ||
-		  sscanf(normal_buf[2], "%f", &facet.normal.z) != 1) {
+	  if (sscanf(normal_buf[0], "%f", &facet.normal(0)) != 1 ||
+		  sscanf(normal_buf[1], "%f", &facet.normal(1)) != 1 ||
+		  sscanf(normal_buf[2], "%f", &facet.normal(2)) != 1) {
 		  // Normal was mangled. Maybe denormals or "not a number" were stored?
 		  // Just reset the normal and silently ignore it.
 		  memset(&facet.normal, 0, sizeof(facet.normal));
@@ -326,104 +325,45 @@ stl_read(stl_file *stl, int first_facet, int first) {
       // It may be worth to round these numbers to zero during loading to reduce the number of errors reported
       // during the STL import.
       for (size_t j = 0; j < 3; ++ j) {
-        if (facet.vertex[j].x > -1e-12f && facet.vertex[j].x < 1e-12f)
-            printf("stl_read: facet %d.x = %e\r\n", j, facet.vertex[j].x);
-        if (facet.vertex[j].y > -1e-12f && facet.vertex[j].y < 1e-12f)
-            printf("stl_read: facet %d.y = %e\r\n", j, facet.vertex[j].y);
-        if (facet.vertex[j].z > -1e-12f && facet.vertex[j].z < 1e-12f)
-            printf("stl_read: facet %d.z = %e\r\n", j, facet.vertex[j].z);
+        if (facet.vertex[j](0) > -1e-12f && facet.vertex[j](0) < 1e-12f)
+            printf("stl_read: facet %d(0) = %e\r\n", j, facet.vertex[j](0));
+        if (facet.vertex[j](1) > -1e-12f && facet.vertex[j](1) < 1e-12f)
+            printf("stl_read: facet %d(1) = %e\r\n", j, facet.vertex[j](1));
+        if (facet.vertex[j](2) > -1e-12f && facet.vertex[j](2) < 1e-12f)
+            printf("stl_read: facet %d(2) = %e\r\n", j, facet.vertex[j](2));
       }
 #endif
 
-#if 1
-    {
-      // Positive and negative zeros are possible in the floats, which are considered equal by the FP unit.
-      // When using a memcmp on raw floats, those numbers report to be different.
-      // Unify all +0 and -0 to +0 to make the floats equal under memcmp.
-      uint32_t *f = (uint32_t*)&facet;
-      for (int j = 0; j < 12; ++ j, ++ f) // 3x vertex + normal: 4x3 = 12 floats
-        if (*f == 0x80000000)
-          // Negative zero, switch to positive zero.
-          *f = 0;
-    }
-#else
-    {
-      // Due to the nature of the floating point numbers, close to zero values may be represented with singificantly higher precision 
-      // than the rest of the vertices. Round them to zero.
-      float *f = (float*)&facet;
-      for (int j = 0; j < 12; ++ j, ++ f) // 3x vertex + normal: 4x3 = 12 floats
-        if (*f > -1e-12f && *f < 1e-12f)
-          // Negative zero, switch to positive zero.
-          *f = 0;
-    }
-#endif
     /* Write the facet into memory. */
-    memcpy(stl->facet_start+i, &facet, SIZEOF_STL_FACET);
+    stl->facet_start[i] = facet;
     stl_facet_stats(stl, facet, first);
-    first = 0;
   }
-  stl->stats.size.x = stl->stats.max.x - stl->stats.min.x;
-  stl->stats.size.y = stl->stats.max.y - stl->stats.min.y;
-  stl->stats.size.z = stl->stats.max.z - stl->stats.min.z;
-  stl->stats.bounding_diameter = sqrt(
-                                   stl->stats.size.x * stl->stats.size.x +
-                                   stl->stats.size.y * stl->stats.size.y +
-                                   stl->stats.size.z * stl->stats.size.z
-                                 );
+  stl->stats.size = stl->stats.max - stl->stats.min;
+  stl->stats.bounding_diameter = stl->stats.size.norm();
 }
 
-void
-stl_facet_stats(stl_file *stl, stl_facet facet, int first) {
-  float diff_x;
-  float diff_y;
-  float diff_z;
-  float max_diff;
+void stl_facet_stats(stl_file *stl, stl_facet facet, bool &first)
+{
+  if (stl->error)
+  	return;
 
-  if (stl->error) return;
+  // While we are going through all of the facets, let's find the
+  // maximum and minimum values for x, y, and z
 
-  /* while we are going through all of the facets, let's find the  */
-  /* maximum and minimum values for x, y, and z  */
-
-  /* Initialize the max and min values the first time through*/
   if (first) {
-    stl->stats.max.x = facet.vertex[0].x;
-    stl->stats.min.x = facet.vertex[0].x;
-    stl->stats.max.y = facet.vertex[0].y;
-    stl->stats.min.y = facet.vertex[0].y;
-    stl->stats.max.z = facet.vertex[0].z;
-    stl->stats.min.z = facet.vertex[0].z;
-
-    diff_x = ABS(facet.vertex[0].x - facet.vertex[1].x);
-    diff_y = ABS(facet.vertex[0].y - facet.vertex[1].y);
-    diff_z = ABS(facet.vertex[0].z - facet.vertex[1].z);
-    max_diff = STL_MAX(diff_x, diff_y);
-    max_diff = STL_MAX(diff_z, max_diff);
-    stl->stats.shortest_edge = max_diff;
-
-    first = 0;
+	// Initialize the max and min values the first time through
+    stl->stats.min = facet.vertex[0];
+    stl->stats.max = facet.vertex[0];
+    stl_vertex diff = (facet.vertex[1] - facet.vertex[0]).cwiseAbs();
+    stl->stats.shortest_edge = std::max(diff(0), std::max(diff(1), diff(2)));
+    first = false;
   }
 
-  /* now find the max and min values */
-  stl->stats.max.x = STL_MAX(stl->stats.max.x, facet.vertex[0].x);
-  stl->stats.min.x = STL_MIN(stl->stats.min.x, facet.vertex[0].x);
-  stl->stats.max.y = STL_MAX(stl->stats.max.y, facet.vertex[0].y);
-  stl->stats.min.y = STL_MIN(stl->stats.min.y, facet.vertex[0].y);
-  stl->stats.max.z = STL_MAX(stl->stats.max.z, facet.vertex[0].z);
-  stl->stats.min.z = STL_MIN(stl->stats.min.z, facet.vertex[0].z);
-
-  stl->stats.max.x = STL_MAX(stl->stats.max.x, facet.vertex[1].x);
-  stl->stats.min.x = STL_MIN(stl->stats.min.x, facet.vertex[1].x);
-  stl->stats.max.y = STL_MAX(stl->stats.max.y, facet.vertex[1].y);
-  stl->stats.min.y = STL_MIN(stl->stats.min.y, facet.vertex[1].y);
-  stl->stats.max.z = STL_MAX(stl->stats.max.z, facet.vertex[1].z);
-  stl->stats.min.z = STL_MIN(stl->stats.min.z, facet.vertex[1].z);
-
-  stl->stats.max.x = STL_MAX(stl->stats.max.x, facet.vertex[2].x);
-  stl->stats.min.x = STL_MIN(stl->stats.min.x, facet.vertex[2].x);
-  stl->stats.max.y = STL_MAX(stl->stats.max.y, facet.vertex[2].y);
-  stl->stats.min.y = STL_MIN(stl->stats.min.y, facet.vertex[2].y);
-  stl->stats.max.z = STL_MAX(stl->stats.max.z, facet.vertex[2].z);
-  stl->stats.min.z = STL_MIN(stl->stats.min.z, facet.vertex[2].z);
+  // Now find the max and min values.
+  for (size_t i = 0; i < 3; ++ i) {
+  	stl->stats.min = stl->stats.min.cwiseMin(facet.vertex[i]);
+  	stl->stats.max = stl->stats.max.cwiseMin(facet.vertex[i]);
+  }
 }
 
 void
diff --git a/xs/src/admesh/util.cpp b/xs/src/admesh/util.cpp
index f3bf59b56..fba1ee76b 100644
--- a/xs/src/admesh/util.cpp
+++ b/xs/src/admesh/util.cpp
@@ -62,7 +62,7 @@ stl_verify_neighbors(stl_file *stl) {
         edge_b.p1 = stl->facet_start[neighbor].vertex[(vnot + 1) % 3];
         edge_b.p2 = stl->facet_start[neighbor].vertex[(vnot + 2) % 3];
       }
-      if(memcmp(&edge_a, &edge_b, SIZEOF_EDGE_SORT) != 0) {
+      if (edge_a.p1 != edge_b.p1 || edge_a.p2 != edge_b.p2) {
         /* These edges should match but they don't.  Print results. */
         printf("edge %d of facet %d doesn't match edge %d of facet %d\n",
                j, i, vnot + 1, neighbor);
@@ -73,114 +73,67 @@ stl_verify_neighbors(stl_file *stl) {
   }
 }
 
-void
-stl_translate(stl_file *stl, float x, float y, float z) {
-  int i;
-  int j;
-
-  if (stl->error) return;
-
-  for(i = 0; i < stl->stats.number_of_facets; i++) {
-    for(j = 0; j < 3; j++) {
-      stl->facet_start[i].vertex[j].x -= (stl->stats.min.x - x);
-      stl->facet_start[i].vertex[j].y -= (stl->stats.min.y - y);
-      stl->facet_start[i].vertex[j].z -= (stl->stats.min.z - z);
-    }
-  }
-  stl->stats.max.x -= (stl->stats.min.x - x);
-  stl->stats.max.y -= (stl->stats.min.y - y);
-  stl->stats.max.z -= (stl->stats.min.z - z);
-  stl->stats.min.x = x;
-  stl->stats.min.y = y;
-  stl->stats.min.z = z;
+void stl_translate(stl_file *stl, float x, float y, float z)
+{
+  if (stl->error)
+  	return;
 
+  stl_vertex new_min(x, y, z);
+  stl_vertex shift = new_min - stl->stats.min;
+  for (int i = 0; i < stl->stats.number_of_facets; ++ i)
+    for (int j = 0; j < 3; ++ j)
+      stl->facet_start[i].vertex[j] += shift;
+  stl->stats.min = new_min;
+  stl->stats.max += shift;
   stl_invalidate_shared_vertices(stl);
 }
 
 /* Translates the stl by x,y,z, relatively from wherever it is currently */
-void
-stl_translate_relative(stl_file *stl, float x, float y, float z) {
-  int i;
-  int j;
-
-  if (stl->error) return;
-
-  for(i = 0; i < stl->stats.number_of_facets; i++) {
-    for(j = 0; j < 3; j++) {
-      stl->facet_start[i].vertex[j].x += x;
-      stl->facet_start[i].vertex[j].y += y;
-      stl->facet_start[i].vertex[j].z += z;
-    }
-  }
-  stl->stats.min.x += x;
-  stl->stats.min.y += y;
-  stl->stats.min.z += z;
-  stl->stats.max.x += x;
-  stl->stats.max.y += y;
-  stl->stats.max.z += z;
+void stl_translate_relative(stl_file *stl, float x, float y, float z)
+{
+  if (stl->error)
+  	return;
 
+  stl_vertex shift(x, y, z);
+  for (int i = 0; i < stl->stats.number_of_facets; ++ i)
+    for (int j = 0; j < 3; ++ j)
+      stl->facet_start[i].vertex[j] += shift;
+  stl->stats.min += shift;
+  stl->stats.max += shift;
   stl_invalidate_shared_vertices(stl);
 }
 
-void
-stl_scale_versor(stl_file *stl, float versor[3]) {
-  int i;
-  int j;
-
-  if (stl->error) return;
-
-  /* scale extents */
-  stl->stats.min.x *= versor[0];
-  stl->stats.min.y *= versor[1];
-  stl->stats.min.z *= versor[2];
-  stl->stats.max.x *= versor[0];
-  stl->stats.max.y *= versor[1];
-  stl->stats.max.z *= versor[2];
-
-  /* scale size */
-  stl->stats.size.x *= versor[0];
-  stl->stats.size.y *= versor[1];
-  stl->stats.size.z *= versor[2];
-
-  /* scale volume */
-  if (stl->stats.volume > 0.0) {
-    stl->stats.volume *= (versor[0] * versor[1] * versor[2]);
-  }
-
-  for(i = 0; i < stl->stats.number_of_facets; i++) {
-    for(j = 0; j < 3; j++) {
-      stl->facet_start[i].vertex[j].x *= versor[0];
-      stl->facet_start[i].vertex[j].y *= versor[1];
-      stl->facet_start[i].vertex[j].z *= versor[2];
-    }
-  }
+void stl_scale_versor(stl_file *stl, const stl_vertex &versor)
+{
+  if (stl->error)
+  	return;
 
+  // Scale extents.
+  auto s = versor.array();
+  stl->stats.min.array() *= s;
+  stl->stats.max.array() *= s;
+  // Scale size.
+  stl->stats.size.array() *= s;
+  // Scale volume.
+  if (stl->stats.volume > 0.0)
+    stl->stats.volume *= versor(0) * versor(1) * versor(2);
+  // Scale the mesh.
+  for (int i = 0; i < stl->stats.number_of_facets; ++ i)
+    for (int j = 0; j < 3; ++ j)
+      stl->facet_start[i].vertex[j].array() *= s;
   stl_invalidate_shared_vertices(stl);
 }
 
-void
-stl_scale(stl_file *stl, float factor) {
-  float versor[3];
-
-  if (stl->error) return;
-
-  versor[0] = factor;
-  versor[1] = factor;
-  versor[2] = factor;
-  stl_scale_versor(stl, versor);
-}
-
-static void calculate_normals(stl_file *stl) {
-  float normal[3];
-
-  if (stl->error) return;
+static void calculate_normals(stl_file *stl) 
+{
+  if (stl->error)
+  	return;
 
+  stl_normal normal;
   for(uint32_t i = 0; i < stl->stats.number_of_facets; i++) {
     stl_calculate_normal(normal, &stl->facet_start[i]);
     stl_normalize_vector(normal);
-    stl->facet_start[i].normal.x = normal[0];
-    stl->facet_start[i].normal.y = normal[1];
-    stl->facet_start[i].normal.z = normal[2];
+    stl->facet_start[i].normal = normal;
   }
 }
 
@@ -193,9 +146,9 @@ void stl_transform(stl_file *stl, float *trafo3x4) {
     for (i_vertex = 0; i_vertex < 3; ++ i_vertex) {
       stl_vertex &v_dst = vertices[i_vertex];
       stl_vertex  v_src = v_dst;
-      v_dst.x = trafo3x4[0] * v_src.x + trafo3x4[1] * v_src.y + trafo3x4[2]  * v_src.z + trafo3x4[3];
-      v_dst.y = trafo3x4[4] * v_src.x + trafo3x4[5] * v_src.y + trafo3x4[6]  * v_src.z + trafo3x4[7];
-      v_dst.z = trafo3x4[8] * v_src.x + trafo3x4[9] * v_src.y + trafo3x4[10] * v_src.z + trafo3x4[11];
+      v_dst(0) = trafo3x4[0] * v_src(0) + trafo3x4[1] * v_src(1) + trafo3x4[2]  * v_src(2) + trafo3x4[3];
+      v_dst(1) = trafo3x4[4] * v_src(0) + trafo3x4[5] * v_src(1) + trafo3x4[6]  * v_src(2) + trafo3x4[7];
+      v_dst(2) = trafo3x4[8] * v_src(0) + trafo3x4[9] * v_src(1) + trafo3x4[10] * v_src(2) + trafo3x4[11];
     }
   }
   stl_get_size(stl);
@@ -214,8 +167,8 @@ stl_rotate_x(stl_file *stl, float angle) {
 
   for(i = 0; i < stl->stats.number_of_facets; i++) {
     for(j = 0; j < 3; j++) {
-      stl_rotate(&stl->facet_start[i].vertex[j].y,
-                 &stl->facet_start[i].vertex[j].z, c, s);
+      stl_rotate(&stl->facet_start[i].vertex[j](1),
+                 &stl->facet_start[i].vertex[j](2), c, s);
     }
   }
   stl_get_size(stl);
@@ -234,8 +187,8 @@ stl_rotate_y(stl_file *stl, float angle) {
 
   for(i = 0; i < stl->stats.number_of_facets; i++) {
     for(j = 0; j < 3; j++) {
-      stl_rotate(&stl->facet_start[i].vertex[j].z,
-                 &stl->facet_start[i].vertex[j].x, c, s);
+      stl_rotate(&stl->facet_start[i].vertex[j](2),
+                 &stl->facet_start[i].vertex[j](0), c, s);
     }
   }
   stl_get_size(stl);
@@ -254,8 +207,8 @@ stl_rotate_z(stl_file *stl, float angle) {
 
   for(i = 0; i < stl->stats.number_of_facets; i++) {
     for(j = 0; j < 3; j++) {
-      stl_rotate(&stl->facet_start[i].vertex[j].x,
-                 &stl->facet_start[i].vertex[j].y, c, s);
+      stl_rotate(&stl->facet_start[i].vertex[j](0),
+                 &stl->facet_start[i].vertex[j](1), c, s);
     }
   }
   stl_get_size(stl);
@@ -272,142 +225,98 @@ stl_rotate(float *x, float *y, const double c, const double s) {
   *y = float(s * xold + c * yold);
 }
 
-extern void
-stl_get_size(stl_file *stl) {
-  int i;
-  int j;
-
-  if (stl->error) return;
-  if (stl->stats.number_of_facets == 0) return;
-
-  stl->stats.min.x = stl->facet_start[0].vertex[0].x;
-  stl->stats.min.y = stl->facet_start[0].vertex[0].y;
-  stl->stats.min.z = stl->facet_start[0].vertex[0].z;
-  stl->stats.max.x = stl->facet_start[0].vertex[0].x;
-  stl->stats.max.y = stl->facet_start[0].vertex[0].y;
-  stl->stats.max.z = stl->facet_start[0].vertex[0].z;
-
-  for(i = 0; i < stl->stats.number_of_facets; i++) {
-    for(j = 0; j < 3; j++) {
-      stl->stats.min.x = STL_MIN(stl->stats.min.x,
-                                 stl->facet_start[i].vertex[j].x);
-      stl->stats.min.y = STL_MIN(stl->stats.min.y,
-                                 stl->facet_start[i].vertex[j].y);
-      stl->stats.min.z = STL_MIN(stl->stats.min.z,
-                                 stl->facet_start[i].vertex[j].z);
-      stl->stats.max.x = STL_MAX(stl->stats.max.x,
-                                 stl->facet_start[i].vertex[j].x);
-      stl->stats.max.y = STL_MAX(stl->stats.max.y,
-                                 stl->facet_start[i].vertex[j].y);
-      stl->stats.max.z = STL_MAX(stl->stats.max.z,
-                                 stl->facet_start[i].vertex[j].z);
+void stl_get_size(stl_file *stl)
+{
+  if (stl->error || stl->stats.number_of_facets == 0)
+  	return;
+  stl->stats.min = stl->facet_start[0].vertex[0];
+  stl->stats.max = stl->stats.min;
+  for (int i = 0; i < stl->stats.number_of_facets; ++ i) {
+  	const stl_facet &face = stl->facet_start[i];
+    for (int j = 0; j < 3; ++ j) {
+      stl->stats.min = stl->stats.min.cwiseMin(face.vertex[j]);
+      stl->stats.max = stl->stats.max.cwiseMax(face.vertex[j]);
     }
   }
-  stl->stats.size.x = stl->stats.max.x - stl->stats.min.x;
-  stl->stats.size.y = stl->stats.max.y - stl->stats.min.y;
-  stl->stats.size.z = stl->stats.max.z - stl->stats.min.z;
-  stl->stats.bounding_diameter = sqrt(
-                                   stl->stats.size.x * stl->stats.size.x +
-                                   stl->stats.size.y * stl->stats.size.y +
-                                   stl->stats.size.z * stl->stats.size.z
-                                 );
+  stl->stats.size = stl->stats.max - stl->stats.min;
+  stl->stats.bounding_diameter = stl->stats.size.norm();
 }
 
-void
-stl_mirror_xy(stl_file *stl) {
-  int i;
-  int j;
-  float temp_size;
+void stl_mirror_xy(stl_file *stl)
+{
+  if (stl->error) 
+  	return;
 
-  if (stl->error) return;
-
-  for(i = 0; i < stl->stats.number_of_facets; i++) {
-    for(j = 0; j < 3; j++) {
-      stl->facet_start[i].vertex[j].z *= -1.0;
+  for(int i = 0; i < stl->stats.number_of_facets; i++) {
+    for(int j = 0; j < 3; j++) {
+      stl->facet_start[i].vertex[j](2) *= -1.0;
     }
   }
-  temp_size = stl->stats.min.z;
-  stl->stats.min.z = stl->stats.max.z;
-  stl->stats.max.z = temp_size;
-  stl->stats.min.z *= -1.0;
-  stl->stats.max.z *= -1.0;
+  float temp_size = stl->stats.min(2);
+  stl->stats.min(2) = stl->stats.max(2);
+  stl->stats.max(2) = temp_size;
+  stl->stats.min(2) *= -1.0;
+  stl->stats.max(2) *= -1.0;
   stl_reverse_all_facets(stl);
   stl->stats.facets_reversed -= stl->stats.number_of_facets;  /* for not altering stats */
 }
 
-void
-stl_mirror_yz(stl_file *stl) {
-  int i;
-  int j;
-  float temp_size;
-
+void stl_mirror_yz(stl_file *stl)
+{
   if (stl->error) return;
 
-  for(i = 0; i < stl->stats.number_of_facets; i++) {
-    for(j = 0; j < 3; j++) {
-      stl->facet_start[i].vertex[j].x *= -1.0;
+  for (int i = 0; i < stl->stats.number_of_facets; i++) {
+    for (int j = 0; j < 3; j++) {
+      stl->facet_start[i].vertex[j](0) *= -1.0;
     }
   }
-  temp_size = stl->stats.min.x;
-  stl->stats.min.x = stl->stats.max.x;
-  stl->stats.max.x = temp_size;
-  stl->stats.min.x *= -1.0;
-  stl->stats.max.x *= -1.0;
+  float temp_size = stl->stats.min(0);
+  stl->stats.min(0) = stl->stats.max(0);
+  stl->stats.max(0) = temp_size;
+  stl->stats.min(0) *= -1.0;
+  stl->stats.max(0) *= -1.0;
   stl_reverse_all_facets(stl);
   stl->stats.facets_reversed -= stl->stats.number_of_facets;  /* for not altering stats */
 }
 
-void
-stl_mirror_xz(stl_file *stl) {
-  int i;
-  int j;
-  float temp_size;
+void stl_mirror_xz(stl_file *stl)
+{
+  if (stl->error)
+  	return;
 
-  if (stl->error) return;
-
-  for(i = 0; i < stl->stats.number_of_facets; i++) {
-    for(j = 0; j < 3; j++) {
-      stl->facet_start[i].vertex[j].y *= -1.0;
+  for (int i = 0; i < stl->stats.number_of_facets; i++) {
+    for (int j = 0; j < 3; j++) {
+      stl->facet_start[i].vertex[j](1) *= -1.0;
     }
   }
-  temp_size = stl->stats.min.y;
-  stl->stats.min.y = stl->stats.max.y;
-  stl->stats.max.y = temp_size;
-  stl->stats.min.y *= -1.0;
-  stl->stats.max.y *= -1.0;
+  float temp_size = stl->stats.min(1);
+  stl->stats.min(1) = stl->stats.max(1);
+  stl->stats.max(1) = temp_size;
+  stl->stats.min(1) *= -1.0;
+  stl->stats.max(1) *= -1.0;
   stl_reverse_all_facets(stl);
   stl->stats.facets_reversed -= stl->stats.number_of_facets;  /* for not altering stats */
 }
 
-static float get_volume(stl_file *stl) {
-  stl_vertex p0;
-  stl_vertex p;
-  stl_normal n;
-  float height;
-  float area;
-  float volume = 0.0;
+static float get_volume(stl_file *stl)
+{
+  if (stl->error)
+  	return 0;
 
-  if (stl->error) return 0;
-
-  /* Choose a point, any point as the reference */
-  p0.x = stl->facet_start[0].vertex[0].x;
-  p0.y = stl->facet_start[0].vertex[0].y;
-  p0.z = stl->facet_start[0].vertex[0].z;
-
-  for(uint32_t i = 0; i < stl->stats.number_of_facets; i++) {
-    p.x = stl->facet_start[i].vertex[0].x - p0.x;
-    p.y = stl->facet_start[i].vertex[0].y - p0.y;
-    p.z = stl->facet_start[i].vertex[0].z - p0.z;
-    /* Do dot product to get distance from point to plane */
-    n = stl->facet_start[i].normal;
-    height = (n.x * p.x) + (n.y * p.y) + (n.z * p.z);
-    area = get_area(&stl->facet_start[i]);
+  // Choose a point, any point as the reference.
+  stl_vertex p0 = stl->facet_start[0].vertex[0];
+  float volume = 0.f;
+  for(uint32_t i = 0; i < stl->stats.number_of_facets; ++ i) {
+    // Do dot product to get distance from point to plane.
+    float height = stl->facet_start[i].normal.dot(stl->facet_start[i].vertex[0] - p0);
+    float area   = get_area(&stl->facet_start[i]);
     volume += (area * height) / 3.0f;
   }
   return volume;
 }
 
-void stl_calculate_volume(stl_file *stl) {
+void stl_calculate_volume(stl_file *stl)
+{
   if (stl->error) return;
   stl->stats.volume = get_volume(stl);
   if(stl->stats.volume < 0.0) {
@@ -416,35 +325,32 @@ void stl_calculate_volume(stl_file *stl) {
   }
 }
 
-static float get_area(stl_facet *facet) {
-  double cross[3][3];
-  float sum[3];
-  float n[3];
-  float area;
-  int i;
-
+static float get_area(stl_facet *facet)
+{
   /* cast to double before calculating cross product because large coordinates
      can result in overflowing product
     (bad area is responsible for bad volume and bad facets reversal) */
-  for(i = 0; i < 3; i++) {
-    cross[i][0]=(((double)facet->vertex[i].y * (double)facet->vertex[(i + 1) % 3].z) -
-                 ((double)facet->vertex[i].z * (double)facet->vertex[(i + 1) % 3].y));
-    cross[i][1]=(((double)facet->vertex[i].z * (double)facet->vertex[(i + 1) % 3].x) -
-                 ((double)facet->vertex[i].x * (double)facet->vertex[(i + 1) % 3].z));
-    cross[i][2]=(((double)facet->vertex[i].x * (double)facet->vertex[(i + 1) % 3].y) -
-                 ((double)facet->vertex[i].y * (double)facet->vertex[(i + 1) % 3].x));
+  double cross[3][3];
+  for (int i = 0; i < 3; i++) {
+    cross[i][0]=(((double)facet->vertex[i](1) * (double)facet->vertex[(i + 1) % 3](2)) -
+                 ((double)facet->vertex[i](2) * (double)facet->vertex[(i + 1) % 3](1)));
+    cross[i][1]=(((double)facet->vertex[i](2) * (double)facet->vertex[(i + 1) % 3](0)) -
+                 ((double)facet->vertex[i](0) * (double)facet->vertex[(i + 1) % 3](2)));
+    cross[i][2]=(((double)facet->vertex[i](0) * (double)facet->vertex[(i + 1) % 3](1)) -
+                 ((double)facet->vertex[i](1) * (double)facet->vertex[(i + 1) % 3](0)));
   }
 
-  sum[0] = cross[0][0] + cross[1][0] + cross[2][0];
-  sum[1] = cross[0][1] + cross[1][1] + cross[2][1];
-  sum[2] = cross[0][2] + cross[1][2] + cross[2][2];
+  stl_normal sum;
+  sum(0) = cross[0][0] + cross[1][0] + cross[2][0];
+  sum(1) = cross[0][1] + cross[1][1] + cross[2][1];
+  sum(2) = cross[0][2] + cross[1][2] + cross[2][2];
 
-  /* This should already be done.  But just in case, let's do it again */
+  // This should already be done.  But just in case, let's do it again.
+  //FIXME this is questionable. the "sum" normal should be accurate, while the normal "n" may be calculated with a low accuracy.
+  stl_normal n;
   stl_calculate_normal(n, facet);
   stl_normalize_vector(n);
-
-  area = 0.5 * (n[0] * sum[0] + n[1] * sum[1] + n[2] * sum[2]);
-  return area;
+  return 0.5f * n.dot(sum);
 }
 
 void stl_repair(stl_file *stl,
diff --git a/xs/src/libslic3r/Format/3mf.cpp b/xs/src/libslic3r/Format/3mf.cpp
index cd6f45c70..6337601cd 100644
--- a/xs/src/libslic3r/Format/3mf.cpp
+++ b/xs/src/libslic3r/Format/3mf.cpp
@@ -1485,7 +1485,7 @@ namespace Slic3r {
                 stl_facet& facet = stl.facet_start[i];
                 for (unsigned int v = 0; v < 3; ++v)
                 {
-                    ::memcpy((void*)&facet.vertex[v].x, (const void*)&geometry.vertices[geometry.triangles[src_start_id + ii + v] * 3], 3 * sizeof(float));
+                    ::memcpy(facet.vertex[v].data(), (const void*)&geometry.vertices[geometry.triangles[src_start_id + ii + v] * 3], 3 * sizeof(float));
                 }
             }
 
@@ -1844,9 +1844,9 @@ namespace Slic3r {
             for (int i = 0; i < stl.stats.shared_vertices; ++i)
             {
                 stream << "     <" << VERTEX_TAG << " ";
-                stream << "x=\"" << stl.v_shared[i].x << "\" ";
-                stream << "y=\"" << stl.v_shared[i].y << "\" ";
-                stream << "z=\"" << stl.v_shared[i].z << "\" />\n";
+                stream << "x=\"" << stl.v_shared[i](0) << "\" ";
+                stream << "y=\"" << stl.v_shared[i](1) << "\" ";
+                stream << "z=\"" << stl.v_shared[i](2) << "\" />\n";
             }
         }
 
diff --git a/xs/src/libslic3r/Format/AMF.cpp b/xs/src/libslic3r/Format/AMF.cpp
index 1ad0c082a..136b23060 100644
--- a/xs/src/libslic3r/Format/AMF.cpp
+++ b/xs/src/libslic3r/Format/AMF.cpp
@@ -402,7 +402,7 @@ void AMFParserContext::endElement(const char * /* name */)
         for (size_t i = 0; i < m_volume_facets.size();) {
             stl_facet &facet = stl.facet_start[i/3];
             for (unsigned int v = 0; v < 3; ++ v)
-                memcpy(&facet.vertex[v].x, &m_object_vertices[m_volume_facets[i ++] * 3], 3 * sizeof(float));
+                memcpy(facet.vertex[v].data(), &m_object_vertices[m_volume_facets[i ++] * 3], 3 * sizeof(float));
         }
         stl_get_size(&stl);
         m_volume->mesh.repair();
@@ -760,9 +760,9 @@ bool store_amf(const char *path, Model *model, Print* print, bool export_print_c
             for (size_t i = 0; i < stl.stats.shared_vertices; ++ i) {
                 stream << "         <vertex>\n";
                 stream << "           <coordinates>\n";
-                stream << "             <x>" << stl.v_shared[i].x << "</x>\n";
-                stream << "             <y>" << stl.v_shared[i].y << "</y>\n";
-                stream << "             <z>" << stl.v_shared[i].z << "</z>\n";
+                stream << "             <x>" << stl.v_shared[i](0) << "</x>\n";
+                stream << "             <y>" << stl.v_shared[i](1) << "</y>\n";
+                stream << "             <z>" << stl.v_shared[i](2) << "</z>\n";
                 stream << "           </coordinates>\n";
                 stream << "         </vertex>\n";
             }
diff --git a/xs/src/libslic3r/Format/OBJ.cpp b/xs/src/libslic3r/Format/OBJ.cpp
index ea6b5604c..ee5756083 100644
--- a/xs/src/libslic3r/Format/OBJ.cpp
+++ b/xs/src/libslic3r/Format/OBJ.cpp
@@ -57,14 +57,14 @@ bool load_obj(const char *path, Model *model, const char *object_name_in)
             continue;
         stl_facet &facet = stl.facet_start[i_face ++];
         size_t     num_normals = 0;
-        stl_normal normal = { 0.f };
+        stl_normal normal(stl_normal::Zero());
         for (unsigned int v = 0; v < 3; ++ v) {
             const ObjParser::ObjVertex &vertex = data.vertices[i++];
-            memcpy(&facet.vertex[v].x, &data.coordinates[vertex.coordIdx*4], 3 * sizeof(float));
+            memcpy(facet.vertex[v].data(), &data.coordinates[vertex.coordIdx*4], 3 * sizeof(float));
             if (vertex.normalIdx != -1) {
-                normal.x += data.normals[vertex.normalIdx*3];
-                normal.y += data.normals[vertex.normalIdx*3+1];
-                normal.z += data.normals[vertex.normalIdx*3+2];
+                normal(0) += data.normals[vertex.normalIdx*3];
+                normal(1) += data.normals[vertex.normalIdx*3+1];
+                normal(2) += data.normals[vertex.normalIdx*3+2];
                 ++ num_normals;
             }
         }
@@ -74,33 +74,27 @@ bool load_obj(const char *path, Model *model, const char *object_name_in)
             facet2.vertex[0] = facet.vertex[0];
             facet2.vertex[1] = facet.vertex[2];
 			const ObjParser::ObjVertex &vertex = data.vertices[i++];
-			memcpy(&facet2.vertex[2].x, &data.coordinates[vertex.coordIdx * 4], 3 * sizeof(float));
+			memcpy(facet2.vertex[2].data(), &data.coordinates[vertex.coordIdx * 4], 3 * sizeof(float));
 			if (vertex.normalIdx != -1) {
-                normal.x += data.normals[vertex.normalIdx*3];
-                normal.y += data.normals[vertex.normalIdx*3+1];
-                normal.z += data.normals[vertex.normalIdx*3+2];
+                normal(0) += data.normals[vertex.normalIdx*3];
+                normal(1) += data.normals[vertex.normalIdx*3+1];
+                normal(2) += data.normals[vertex.normalIdx*3+2];
                 ++ num_normals;
             }
             if (num_normals == 4) {
                 // Normalize an average normal of a quad.
-                float len = sqrt(facet.normal.x*facet.normal.x + facet.normal.y*facet.normal.y + facet.normal.z*facet.normal.z);
+                float len = facet.normal.norm();
                 if (len > EPSILON) {
-                    normal.x /= len;
-                    normal.y /= len;
-                    normal.z /= len;
+                    normal /= len;
                     facet.normal = normal;
                     facet2.normal = normal;
                 }
             }
         } else if (num_normals == 3) {
             // Normalize an average normal of a triangle.
-            float len = sqrt(facet.normal.x*facet.normal.x + facet.normal.y*facet.normal.y + facet.normal.z*facet.normal.z);
-            if (len > EPSILON) {
-                normal.x /= len;
-                normal.y /= len;
-                normal.z /= len;
-                facet.normal = normal;
-            }
+            float len = facet.normal.norm();
+            if (len > EPSILON)
+                facet.normal = normal / len;
         }
 	}
     stl_get_size(&stl);
diff --git a/xs/src/libslic3r/Format/PRUS.cpp b/xs/src/libslic3r/Format/PRUS.cpp
index 9f42fc151..3cf3fc075 100644
--- a/xs/src/libslic3r/Format/PRUS.cpp
+++ b/xs/src/libslic3r/Format/PRUS.cpp
@@ -260,8 +260,8 @@ bool load_prus(const char *path, Model *model)
 							mesh.repair();
 							// Transform the model.
 							stl_transform(&stl, &trafo[0][0]);
-							if (std::abs(stl.stats.min.z) < EPSILON)
-								stl.stats.min.z = 0.;
+							if (std::abs(stl.stats.min(2)) < EPSILON)
+								stl.stats.min(2) = 0.;
 							// Add a mesh to a model.
 							if (mesh.facets_count() > 0)
                                 mesh_valid = true;
@@ -309,11 +309,11 @@ bool load_prus(const char *path, Model *model)
 						assert(res_normal == 3);
                         int res_outer_loop	= line_reader.next_line_scanf(" outer loop");
 						assert(res_outer_loop == 0);
-						int res_vertex1 = line_reader.next_line_scanf(" vertex %f %f %f", &facet.vertex[0].x, &facet.vertex[0].y, &facet.vertex[0].z);
+						int res_vertex1 = line_reader.next_line_scanf(" vertex %f %f %f", &facet.vertex[0](0), &facet.vertex[0](1), &facet.vertex[0](2));
 						assert(res_vertex1 == 3);
-						int res_vertex2 = line_reader.next_line_scanf(" vertex %f %f %f", &facet.vertex[1].x, &facet.vertex[1].y, &facet.vertex[1].z);
+						int res_vertex2 = line_reader.next_line_scanf(" vertex %f %f %f", &facet.vertex[1](0), &facet.vertex[1](1), &facet.vertex[1](2));
 						assert(res_vertex2 == 3);
-						int res_vertex3 = line_reader.next_line_scanf(" vertex %f %f %f", &facet.vertex[2].x, &facet.vertex[2].y, &facet.vertex[2].z);
+						int res_vertex3 = line_reader.next_line_scanf(" vertex %f %f %f", &facet.vertex[2](0), &facet.vertex[2](1), &facet.vertex[2](2));
 						assert(res_vertex3 == 3);
 						int res_endloop = line_reader.next_line_scanf(" endloop");
 						assert(res_endloop == 0);
@@ -324,9 +324,9 @@ bool load_prus(const char *path, Model *model)
                             break;
                         }
                         // The facet normal has been parsed as a single string as to workaround for not a numbers in the normal definition.
-                        if (sscanf(normal_buf[0], "%f", &facet.normal.x) != 1 ||
-                            sscanf(normal_buf[1], "%f", &facet.normal.y) != 1 ||
-                            sscanf(normal_buf[2], "%f", &facet.normal.z) != 1) {
+                        if (sscanf(normal_buf[0], "%f", &facet.normal(0)) != 1 ||
+                            sscanf(normal_buf[1], "%f", &facet.normal(1)) != 1 ||
+                            sscanf(normal_buf[2], "%f", &facet.normal(2)) != 1) {
                             // Normal was mangled. Maybe denormals or "not a number" were stored?
                             // Just reset the normal and silently ignore it.
                             memset(&facet.normal, 0, sizeof(facet.normal));
diff --git a/xs/src/libslic3r/Model.cpp b/xs/src/libslic3r/Model.cpp
index e0f226476..60f949837 100644
--- a/xs/src/libslic3r/Model.cpp
+++ b/xs/src/libslic3r/Model.cpp
@@ -642,24 +642,13 @@ BoundingBoxf3 ModelObject::tight_bounding_box(bool include_modifiers) const
                     {
                         // original point
                         const stl_vertex& v = facet.vertex[i];
-                        Vec3d p((double)v.x, (double)v.y, (double)v.z);
-
                         // scale
-                        p(0) *= inst->scaling_factor;
-                        p(1) *= inst->scaling_factor;
-                        p(2) *= inst->scaling_factor;
-
-                        // rotate Z
-                        double x = p(0);
-                        double y = p(1);
-                        p(0) = c * x - s * y;
-                        p(1) = s * x + c * y;
-
-                        // translate
-                        p(0) += inst->offset(0);
-                        p(1) += inst->offset(1);
-
-                        bb.merge(p);
+                        Vec3d p = v.cast<double>() * inst->scaling_factor;
+                        // rotate Z, translate
+                        Vec3d q(c * p(0) - s * p(1) + inst->offset(0),
+                                s * p(0) + c * p(1) + inst->offset(1), 
+                                p(2));
+                        bb.merge(q);
                     }
                 }
             }
@@ -770,7 +759,7 @@ void ModelObject::rotate(float angle, const Axis &axis)
     for (ModelVolume *v : this->volumes)
     {
         v->mesh.rotate(angle, axis);
-        min_z = std::min(min_z, v->mesh.stl.stats.min.z);
+        min_z = std::min(min_z, v->mesh.stl.stats.min(2));
     }
 
     if (min_z != 0.0f)
@@ -927,24 +916,13 @@ void ModelObject::check_instances_print_volume_state(const BoundingBoxf3& print_
                     {
                         // original point
                         const stl_vertex& v = facet.vertex[i];
-                        Vec3d p((double)v.x, (double)v.y, (double)v.z);
-
                         // scale
-                        p(0) *= inst->scaling_factor;
-                        p(1) *= inst->scaling_factor;
-                        p(2) *= inst->scaling_factor;
-
+                        Vec3d p = v.cast<double>() * inst->scaling_factor;
                         // rotate Z
-                        double x = p(0);
-                        double y = p(1);
-                        p(0) = c * x - s * y;
-                        p(1) = s * x + c * y;
-
-                        // translate
-                        p(0) += inst->offset(0);
-                        p(1) += inst->offset(1);
-
-                        bb.merge(p);
+                        bb.merge(Vec3d(
+                            c * p(0) - s * p(1) + inst->offset(0),
+                            s * p(0) + c * p(1) + inst->offset(1), 
+                            p(2)));
                     }
                 }
 
@@ -1081,7 +1059,7 @@ BoundingBoxf3 ModelInstance::transform_mesh_bounding_box(const TriangleMesh* mes
         const stl_facet &facet = mesh->stl.facet_start[i];
         for (int j = 0; j < 3; ++ j) {
             const stl_vertex &v = facet.vertex[j];
-			bbox.merge(Vec3d(c * v.x - s * v.y, s * v.x + c * v.y, v.z));
+			bbox.merge(Vec3d(c * v(0) - s * v(1), s * v(0) + c * v(1), v(2)));
         }
     }
     if (! empty(bbox)) {
diff --git a/xs/src/libslic3r/SlicingAdaptive.cpp b/xs/src/libslic3r/SlicingAdaptive.cpp
index ff0da7636..2ef4aec8c 100644
--- a/xs/src/libslic3r/SlicingAdaptive.cpp
+++ b/xs/src/libslic3r/SlicingAdaptive.cpp
@@ -15,8 +15,8 @@ void SlicingAdaptive::clear()
 std::pair<float, float> face_z_span(const stl_facet *f)
 {
 	return std::pair<float, float>(
-		std::min(std::min(f->vertex[0].z, f->vertex[1].z), f->vertex[2].z),
-		std::max(std::max(f->vertex[0].z, f->vertex[1].z), f->vertex[2].z));
+		std::min(std::min(f->vertex[0](2), f->vertex[1](2)), f->vertex[2](2)),
+		std::max(std::max(f->vertex[0](2), f->vertex[1](2)), f->vertex[2](2)));
 }
 
 void SlicingAdaptive::prepare()
@@ -40,7 +40,7 @@ void SlicingAdaptive::prepare()
 	// 3) Generate Z components of the facet normals.
 	m_face_normal_z.assign(m_faces.size(), 0.f);
     for (size_t iface = 0; iface < m_faces.size(); ++ iface)
-    	m_face_normal_z[iface] = m_faces[iface]->normal.z;
+    	m_face_normal_z[iface] = m_faces[iface]->normal(2);
 }
 
 float SlicingAdaptive::cusp_height(float z, float cusp_value, int &current_facet)
diff --git a/xs/src/libslic3r/TriangleMesh.cpp b/xs/src/libslic3r/TriangleMesh.cpp
index 782c50cf5..3d090bb7b 100644
--- a/xs/src/libslic3r/TriangleMesh.cpp
+++ b/xs/src/libslic3r/TriangleMesh.cpp
@@ -51,31 +51,16 @@ TriangleMesh::TriangleMesh(const Pointf3s &points, const std::vector<Vec3crd>& f
 
     for (int i = 0; i < stl.stats.number_of_facets; i++) {
         stl_facet facet;
-
-        const Vec3d& ref_f1 = points[facets[i](0)];
-        facet.vertex[0].x = ref_f1(0);
-        facet.vertex[0].y = ref_f1(1);
-        facet.vertex[0].z = ref_f1(2);
-
-        const Vec3d& ref_f2 = points[facets[i](1)];
-        facet.vertex[1].x = ref_f2(0);
-        facet.vertex[1].y = ref_f2(1);
-        facet.vertex[1].z = ref_f2(2);
-
-        const Vec3d& ref_f3 = points[facets[i](2)];
-        facet.vertex[2].x = ref_f3(0);
-        facet.vertex[2].y = ref_f3(1);
-        facet.vertex[2].z = ref_f3(2);
-        
+        facet.vertex[0] = points[facets[i](0)].cast<float>();
+        facet.vertex[1] = points[facets[i](1)].cast<float>();
+        facet.vertex[2] = points[facets[i](2)].cast<float>();
         facet.extra[0] = 0;
         facet.extra[1] = 0;
 
-        float normal[3];
+        stl_normal normal;
         stl_calculate_normal(normal, &facet);
         stl_normalize_vector(normal);
-        facet.normal.x = normal[0];
-        facet.normal.y = normal[1];
-        facet.normal.z = normal[2];
+        facet.normal = normal;
 
         stl.facet_start[i] = facet;
     }
@@ -302,11 +287,7 @@ void TriangleMesh::scale(float factor)
 
 void TriangleMesh::scale(const Vec3d &versor)
 {
-    float fversor[3];
-    fversor[0] = versor(0);
-    fversor[1] = versor(1);
-    fversor[2] = versor(2);
-    stl_scale_versor(&this->stl, fversor);
+    stl_scale_versor(&this->stl, versor.cast<float>());
     stl_invalidate_shared_vertices(&this->stl);
 }
 
@@ -390,10 +371,9 @@ void TriangleMesh::transform(const float* matrix3x4)
 void TriangleMesh::align_to_origin()
 {
     this->translate(
-        -(this->stl.stats.min.x),
-        -(this->stl.stats.min.y),
-        -(this->stl.stats.min.z)
-    );
+        - this->stl.stats.min(0),
+        - this->stl.stats.min(1),
+        - this->stl.stats.min(2));
 }
 
 void TriangleMesh::rotate(double angle, Point* center)
@@ -518,7 +498,7 @@ TriangleMesh::split() const
         stl_clear_error(&mesh->stl);
         stl_allocate(&mesh->stl);
         
-        int first = 1;
+        bool first = true;
         for (std::deque<int>::const_iterator facet = facets.begin(); facet != facets.end(); ++facet) {
             mesh->stl.facet_start[facet - facets.begin()] = this->stl.facet_start[*facet];
             stl_facet_stats(&mesh->stl, this->stl.facet_start[*facet], first);
@@ -561,9 +541,9 @@ ExPolygons TriangleMesh::horizontal_projection() const
         stl_facet* facet = &this->stl.facet_start[i];
         Polygon p;
         p.points.resize(3);
-        p.points[0] = Point::new_scale(facet->vertex[0].x, facet->vertex[0].y);
-        p.points[1] = Point::new_scale(facet->vertex[1].x, facet->vertex[1].y);
-        p.points[2] = Point::new_scale(facet->vertex[2].x, facet->vertex[2].y);
+        p.points[0] = Point::new_scale(facet->vertex[0](0), facet->vertex[0](1));
+        p.points[1] = Point::new_scale(facet->vertex[1](0), facet->vertex[1](1));
+        p.points[2] = Point::new_scale(facet->vertex[2](0), facet->vertex[2](1));
         p.make_counter_clockwise();  // do this after scaling, as winding order might change while doing that
         pp.emplace_back(p);
     }
@@ -578,8 +558,8 @@ Polygon TriangleMesh::convex_hull()
     Points pp;
     pp.reserve(this->stl.stats.shared_vertices);
     for (int i = 0; i < this->stl.stats.shared_vertices; ++ i) {
-        stl_vertex* v = &this->stl.v_shared[i];
-        pp.emplace_back(Point::new_scale(v->x, v->y));
+        const stl_vertex &v = this->stl.v_shared[i];
+        pp.emplace_back(Point::new_scale(v(0), v(1)));
     }
     return Slic3r::Geometry::convex_hull(pp);
 }
@@ -589,12 +569,8 @@ TriangleMesh::bounding_box() const
 {
     BoundingBoxf3 bb;
     bb.defined = true;
-    bb.min(0) = this->stl.stats.min.x;
-    bb.min(1) = this->stl.stats.min.y;
-    bb.min(2) = this->stl.stats.min.z;
-    bb.max(0) = this->stl.stats.max.x;
-    bb.max(1) = this->stl.stats.max.y;
-    bb.max(2) = this->stl.stats.max.z;
+    bb.min = this->stl.stats.min.cast<double>();
+    bb.max = this->stl.stats.max.cast<double>();
     return bb;
 }
 
@@ -619,11 +595,8 @@ TriangleMeshSlicer::TriangleMeshSlicer(TriangleMesh* _mesh) :
     facets_edges.assign(_mesh->stl.stats.number_of_facets * 3, -1);
     v_scaled_shared.assign(_mesh->stl.v_shared, _mesh->stl.v_shared + _mesh->stl.stats.shared_vertices);
     // Scale the copied vertices.
-    for (int i = 0; i < this->mesh->stl.stats.shared_vertices; ++ i) {
-        this->v_scaled_shared[i].x /= float(SCALING_FACTOR);
-        this->v_scaled_shared[i].y /= float(SCALING_FACTOR);
-        this->v_scaled_shared[i].z /= float(SCALING_FACTOR);
-    }
+    for (int i = 0; i < this->mesh->stl.stats.shared_vertices; ++ i)
+        this->v_scaled_shared[i] *= float(1. / SCALING_FACTOR);
 
     // Create a mapping from triangle edge into face.
     struct EdgeToFace {
@@ -779,14 +752,14 @@ void TriangleMeshSlicer::_slice_do(size_t facet_idx, std::vector<IntersectionLin
     const stl_facet &facet = this->mesh->stl.facet_start[facet_idx];
     
     // find facet extents
-    const float min_z = fminf(facet.vertex[0].z, fminf(facet.vertex[1].z, facet.vertex[2].z));
-    const float max_z = fmaxf(facet.vertex[0].z, fmaxf(facet.vertex[1].z, facet.vertex[2].z));
+    const float min_z = fminf(facet.vertex[0](2), fminf(facet.vertex[1](2), facet.vertex[2](2)));
+    const float max_z = fmaxf(facet.vertex[0](2), fmaxf(facet.vertex[1](2), facet.vertex[2](2)));
     
     #ifdef SLIC3R_DEBUG
     printf("\n==> FACET %d (%f,%f,%f - %f,%f,%f - %f,%f,%f):\n", facet_idx,
-        facet.vertex[0].x, facet.vertex[0].y, facet.vertex[0].z,
-        facet.vertex[1].x, facet.vertex[1].y, facet.vertex[1].z,
-        facet.vertex[2].x, facet.vertex[2].y, facet.vertex[2].z);
+        facet.vertex[0].x, facet.vertex[0].y, facet.vertex[0](2),
+        facet.vertex[1].x, facet.vertex[1].y, facet.vertex[1](2),
+        facet.vertex[2].x, facet.vertex[2].y, facet.vertex[2](2));
     printf("z: min = %.2f, max = %.2f\n", min_z, max_z);
     #endif
     
@@ -806,18 +779,18 @@ void TriangleMeshSlicer::_slice_do(size_t facet_idx, std::vector<IntersectionLin
             if (il.edge_type == feHorizontal) {
                 // Insert all three edges of the face.
                 const int *vertices = this->mesh->stl.v_indices[facet_idx].vertex;
-                const bool reverse  = this->mesh->stl.facet_start[facet_idx].normal.z < 0;
+                const bool reverse  = this->mesh->stl.facet_start[facet_idx].normal(2) < 0;
                 for (int j = 0; j < 3; ++ j) {
                     int               a_id     = vertices[j % 3];
                     int               b_id     = vertices[(j+1) % 3];
                     if (reverse)
                         std::swap(a_id, b_id);
-                    const stl_vertex *a = &this->v_scaled_shared[a_id];
-                    const stl_vertex *b = &this->v_scaled_shared[b_id];
-                    il.a(0)    = a->x;
-                    il.a(1)    = a->y;
-                    il.b(0)    = b->x;
-                    il.b(1)    = b->y;
+                    const stl_vertex &a = this->v_scaled_shared[a_id];
+                    const stl_vertex &b = this->v_scaled_shared[b_id];
+                    il.a(0)    = a(0);
+                    il.a(1)    = a(1);
+                    il.b(0)    = b(0);
+                    il.b(1)    = b(1);
                     il.a_id   = a_id;
                     il.b_id   = b_id;
                     (*lines)[layer_idx].emplace_back(il);
@@ -863,66 +836,63 @@ bool TriangleMeshSlicer::slice_facet(
     // Reorder vertices so that the first one is the one with lowest Z.
     // This is needed to get all intersection lines in a consistent order
     // (external on the right of the line)
-    int i = (facet.vertex[1].z == min_z) ? 1 : ((facet.vertex[2].z == min_z) ? 2 : 0);
+    int i = (facet.vertex[1](2) == min_z) ? 1 : ((facet.vertex[2](2) == min_z) ? 2 : 0);
     for (int j = i; j - i < 3; ++ j) {  // loop through facet edges
         int               edge_id  = this->facets_edges[facet_idx * 3 + (j % 3)];
         const int        *vertices = this->mesh->stl.v_indices[facet_idx].vertex;
         int               a_id     = vertices[j % 3];
         int               b_id     = vertices[(j+1) % 3];
-        const stl_vertex *a = &this->v_scaled_shared[a_id];
-        const stl_vertex *b = &this->v_scaled_shared[b_id];
+        const stl_vertex &a = this->v_scaled_shared[a_id];
+        const stl_vertex &b = this->v_scaled_shared[b_id];
         
         // Is edge or face aligned with the cutting plane?
-        if (a->z == slice_z && b->z == slice_z) {
+        if (a(2) == slice_z && b(2) == slice_z) {
             // Edge is horizontal and belongs to the current layer.
             const stl_vertex &v0 = this->v_scaled_shared[vertices[0]];
             const stl_vertex &v1 = this->v_scaled_shared[vertices[1]];
             const stl_vertex &v2 = this->v_scaled_shared[vertices[2]];
+            bool              swap = false;
             if (min_z == max_z) {
                 // All three vertices are aligned with slice_z.
                 line_out->edge_type = feHorizontal;
-                if (this->mesh->stl.facet_start[facet_idx].normal.z < 0) {
+                if (this->mesh->stl.facet_start[facet_idx].normal(2) < 0) {
                     // If normal points downwards this is a bottom horizontal facet so we reverse its point order.
-                    std::swap(a, b);
-                    std::swap(a_id, b_id);
+                    swap = true;
                 }
-            } else if (v0.z < slice_z || v1.z < slice_z || v2.z < slice_z) {
+            } else if (v0(2) < slice_z || v1(2) < slice_z || v2(2) < slice_z) {
                 // Two vertices are aligned with the cutting plane, the third vertex is below the cutting plane.
                 line_out->edge_type = feTop;
-                std::swap(a, b);
-                std::swap(a_id, b_id);
+                swap = true;
             } else {
                 // Two vertices are aligned with the cutting plane, the third vertex is above the cutting plane.
                 line_out->edge_type = feBottom;
             }
-            line_out->a(0)  = a->x;
-            line_out->a(1)  = a->y;
-            line_out->b(0)  = b->x;
-            line_out->b(1)  = b->y;
-            line_out->a_id   = a_id;
-            line_out->b_id   = b_id;
+            line_out->a = to_2d(swap ? b : a).cast<coord_t>();
+            line_out->b = to_2d(swap ? a : b).cast<coord_t>();
+            line_out->a_id = swap ? b_id : a_id;
+            line_out->b_id = swap ? a_id : b_id;
             return true;
         }
 
-        if (a->z == slice_z) {
+        if (a(2) == slice_z) {
             // Only point a alings with the cutting plane.
             points_on_layer[num_points_on_layer ++] = num_points;
             IntersectionPoint &point = points[num_points ++];
-            point(0)       = a->x;
-            point(1)       = a->y;
+            point(0)       = a(0);
+            point(1)       = a(1);
             point.point_id  = a_id;
-        } else if (b->z == slice_z) {
+        } else if (b(2) == slice_z) {
             // Only point b alings with the cutting plane.
             points_on_layer[num_points_on_layer ++] = num_points;
             IntersectionPoint &point = points[num_points ++];
-            point(0)       = b->x;
-            point(1)       = b->y;
+            point(0)       = b(0);
+            point(1)       = b(1);
             point.point_id  = b_id;
-        } else if ((a->z < slice_z && b->z > slice_z) || (b->z < slice_z && a->z > slice_z)) {
+        } else if ((a(2) < slice_z && b(2) > slice_z) || (b(2) < slice_z && a(2) > slice_z)) {
             // A general case. The face edge intersects the cutting plane. Calculate the intersection point.
             IntersectionPoint &point = points[num_points ++];
-            point(0)       = b->x + (a->x - b->x) * (slice_z - b->z) / (a->z - b->z);
-            point(1)       = b->y + (a->y - b->y) * (slice_z - b->z) / (a->z - b->z);
+            point(0)       = b(0) + (a(0) - b(0)) * (slice_z - b(2)) / (a(2) - b(2));
+            point(1)       = b(1) + (a(1) - b(1)) * (slice_z - b(2)) / (a(2) - b(2));
             point.edge_id   = edge_id;
         }
     }
@@ -1389,8 +1359,8 @@ void TriangleMeshSlicer::cut(float z, TriangleMesh* upper, TriangleMesh* lower)
         stl_facet* facet = &this->mesh->stl.facet_start[facet_idx];
         
         // find facet extents
-        float min_z = std::min(facet->vertex[0].z, std::min(facet->vertex[1].z, facet->vertex[2].z));
-        float max_z = std::max(facet->vertex[0].z, std::max(facet->vertex[1].z, facet->vertex[2].z));
+        float min_z = std::min(facet->vertex[0](2), std::min(facet->vertex[1](2), facet->vertex[2](2)));
+        float max_z = std::max(facet->vertex[0](2), std::max(facet->vertex[1](2), facet->vertex[2](2)));
         
         // intersect facet with cutting plane
         IntersectionLine line;
@@ -1417,47 +1387,47 @@ void TriangleMeshSlicer::cut(float z, TriangleMesh* upper, TriangleMesh* lower)
 
             // look for the vertex on whose side of the slicing plane there are no other vertices
             int isolated_vertex;
-            if ( (facet->vertex[0].z > z) == (facet->vertex[1].z > z) ) {
+            if ( (facet->vertex[0](2) > z) == (facet->vertex[1](2) > z) ) {
                 isolated_vertex = 2;
-            } else if ( (facet->vertex[1].z > z) == (facet->vertex[2].z > z) ) {
+            } else if ( (facet->vertex[1](2) > z) == (facet->vertex[2](2) > z) ) {
                 isolated_vertex = 0;
             } else {
                 isolated_vertex = 1;
             }
             
             // get vertices starting from the isolated one
-            stl_vertex* v0 = &facet->vertex[isolated_vertex];
-            stl_vertex* v1 = &facet->vertex[(isolated_vertex+1) % 3];
-            stl_vertex* v2 = &facet->vertex[(isolated_vertex+2) % 3];
+            const stl_vertex &v0 = facet->vertex[isolated_vertex];
+            const stl_vertex &v1 = facet->vertex[(isolated_vertex+1) % 3];
+            const stl_vertex &v2 = facet->vertex[(isolated_vertex+2) % 3];
             
             // intersect v0-v1 and v2-v0 with cutting plane and make new vertices
             stl_vertex v0v1, v2v0;
-            v0v1.x = v1->x + (v0->x - v1->x) * (z - v1->z) / (v0->z - v1->z);
-            v0v1.y = v1->y + (v0->y - v1->y) * (z - v1->z) / (v0->z - v1->z);
-            v0v1.z = z;
-            v2v0.x = v2->x + (v0->x - v2->x) * (z - v2->z) / (v0->z - v2->z);
-            v2v0.y = v2->y + (v0->y - v2->y) * (z - v2->z) / (v0->z - v2->z);
-            v2v0.z = z;
+            v0v1(0) = v1(0) + (v0(0) - v1(0)) * (z - v1(2)) / (v0(2) - v1(2));
+            v0v1(1) = v1(1) + (v0(1) - v1(1)) * (z - v1(2)) / (v0(2) - v1(2));
+            v0v1(2) = z;
+            v2v0(0) = v2(0) + (v0(0) - v2(0)) * (z - v2(2)) / (v0(2) - v2(2));
+            v2v0(1) = v2(1) + (v0(1) - v2(1)) * (z - v2(2)) / (v0(2) - v2(2));
+            v2v0(2) = z;
             
             // build the triangular facet
             stl_facet triangle;
             triangle.normal = facet->normal;
-            triangle.vertex[0] = *v0;
+            triangle.vertex[0] = v0;
             triangle.vertex[1] = v0v1;
             triangle.vertex[2] = v2v0;
             
             // build the facets forming a quadrilateral on the other side
             stl_facet quadrilateral[2];
             quadrilateral[0].normal = facet->normal;
-            quadrilateral[0].vertex[0] = *v1;
-            quadrilateral[0].vertex[1] = *v2;
+            quadrilateral[0].vertex[0] = v1;
+            quadrilateral[0].vertex[1] = v2;
             quadrilateral[0].vertex[2] = v0v1;
             quadrilateral[1].normal = facet->normal;
-            quadrilateral[1].vertex[0] = *v2;
+            quadrilateral[1].vertex[0] = v2;
             quadrilateral[1].vertex[1] = v2v0;
             quadrilateral[1].vertex[2] = v0v1;
             
-            if (v0->z > z) {
+            if (v0(2) > z) {
                 if (upper != NULL) stl_add_facet(&upper->stl, &triangle);
                 if (lower != NULL) {
                     stl_add_facet(&lower->stl, &quadrilateral[0]);
@@ -1489,13 +1459,11 @@ void TriangleMeshSlicer::cut(float z, TriangleMesh* upper, TriangleMesh* lower)
             Polygon p = *polygon;
             p.reverse();
             stl_facet facet;
-            facet.normal.x = 0;
-            facet.normal.y = 0;
-            facet.normal.z = -1;
+            facet.normal = stl_normal(0, 0, -1.f);
             for (size_t i = 0; i <= 2; ++i) {
-                facet.vertex[i].x = unscale<float>(p.points[i](0));
-                facet.vertex[i].y = unscale<float>(p.points[i](1));
-                facet.vertex[i].z = z;
+                facet.vertex[i](0) = unscale<float>(p.points[i](0));
+                facet.vertex[i](1) = unscale<float>(p.points[i](1));
+                facet.vertex[i](2) = z;
             }
             stl_add_facet(&upper->stl, &facet);
         }
@@ -1515,13 +1483,11 @@ void TriangleMeshSlicer::cut(float z, TriangleMesh* upper, TriangleMesh* lower)
         // convert triangles to facets and append them to mesh
         for (Polygons::const_iterator polygon = triangles.begin(); polygon != triangles.end(); ++polygon) {
             stl_facet facet;
-            facet.normal.x = 0;
-            facet.normal.y = 0;
-            facet.normal.z = 1;
+            facet.normal = stl_normal(0, 0, 1.f);
             for (size_t i = 0; i <= 2; ++i) {
-                facet.vertex[i].x = unscale<float>(polygon->points[i](0));
-                facet.vertex[i].y = unscale<float>(polygon->points[i](1));
-                facet.vertex[i].z = z;
+                facet.vertex[i](0) = unscale<float>(polygon->points[i](0));
+                facet.vertex[i](1) = unscale<float>(polygon->points[i](1));
+                facet.vertex[i](2) = z;
             }
             stl_add_facet(&lower->stl, &facet);
         }
diff --git a/xs/src/slic3r/GUI/3DScene.cpp b/xs/src/slic3r/GUI/3DScene.cpp
index 93b9a27ab..d02a9c082 100644
--- a/xs/src/slic3r/GUI/3DScene.cpp
+++ b/xs/src/slic3r/GUI/3DScene.cpp
@@ -39,7 +39,7 @@ void GLIndexedVertexArray::load_mesh_flat_shading(const TriangleMesh &mesh)
     for (int i = 0; i < mesh.stl.stats.number_of_facets; ++ i) {
         const stl_facet &facet = mesh.stl.facet_start[i];
         for (int j = 0; j < 3; ++ j)
-            this->push_geometry(facet.vertex[j].x, facet.vertex[j].y, facet.vertex[j].z, facet.normal.x, facet.normal.y, facet.normal.z);
+            this->push_geometry(facet.vertex[j](0), facet.vertex[j](1), facet.vertex[j](2), facet.normal(0), facet.normal(1), facet.normal(2));
     }
 }
 
@@ -55,7 +55,7 @@ void GLIndexedVertexArray::load_mesh_full_shading(const TriangleMesh &mesh)
     for (int i = 0; i < mesh.stl.stats.number_of_facets; ++i) {
         const stl_facet &facet = mesh.stl.facet_start[i];
         for (int j = 0; j < 3; ++j)
-            this->push_geometry(facet.vertex[j].x, facet.vertex[j].y, facet.vertex[j].z, facet.normal.x, facet.normal.y, facet.normal.z);
+            this->push_geometry(facet.vertex[j](0), facet.vertex[j](1), facet.vertex[j](2), facet.normal(0), facet.normal(1), facet.normal(2));
 
         this->push_triangle(vertices_count, vertices_count + 1, vertices_count + 2);
         vertices_count += 3;
diff --git a/xs/xsp/TriangleMesh.xsp b/xs/xsp/TriangleMesh.xsp
index 4a269bd84..2bc20da3a 100644
--- a/xs/xsp/TriangleMesh.xsp
+++ b/xs/xsp/TriangleMesh.xsp
@@ -60,14 +60,14 @@ TriangleMesh::ReadFromPerl(vertices, facets)
         for (int i = 0; i < stl.stats.number_of_facets; i++) {
             AV* facet_av = (AV*)SvRV(*av_fetch(facets_av, i, 0));
             stl_facet facet;
-            facet.normal.x = 0;
-            facet.normal.y = 0;
-            facet.normal.z = 0;
+            facet.normal(0) = 0;
+            facet.normal(1) = 0;
+            facet.normal(2) = 0;
             for (unsigned int v = 0; v <= 2; v++) {
                 AV* vertex_av = (AV*)SvRV(*av_fetch(vertices_av, SvIV(*av_fetch(facet_av, v, 0)), 0));
-                facet.vertex[v].x = SvNV(*av_fetch(vertex_av, 0, 0));
-                facet.vertex[v].y = SvNV(*av_fetch(vertex_av, 1, 0));
-                facet.vertex[v].z = SvNV(*av_fetch(vertex_av, 2, 0));
+                facet.vertex[v](0) = SvNV(*av_fetch(vertex_av, 0, 0));
+                facet.vertex[v](1) = SvNV(*av_fetch(vertex_av, 1, 0));
+                facet.vertex[v](2) = SvNV(*av_fetch(vertex_av, 2, 0));
             }
             facet.extra[0] = 0;
             facet.extra[1] = 0;
@@ -110,9 +110,9 @@ TriangleMesh::vertices()
             AV* vertex = newAV();
             av_store(vertices, i, newRV_noinc((SV*)vertex));
             av_extend(vertex, 2);
-            av_store(vertex, 0, newSVnv(THIS->stl.v_shared[i].x));
-            av_store(vertex, 1, newSVnv(THIS->stl.v_shared[i].y));
-            av_store(vertex, 2, newSVnv(THIS->stl.v_shared[i].z));
+            av_store(vertex, 0, newSVnv(THIS->stl.v_shared[i](0)));
+            av_store(vertex, 1, newSVnv(THIS->stl.v_shared[i](1)));
+            av_store(vertex, 2, newSVnv(THIS->stl.v_shared[i](2)));
         }
         
         RETVAL = newRV_noinc((SV*)vertices);
@@ -155,9 +155,9 @@ TriangleMesh::normals()
             AV* facet = newAV();
             av_store(normals, i, newRV_noinc((SV*)facet));
             av_extend(facet, 2);
-            av_store(facet, 0, newSVnv(THIS->stl.facet_start[i].normal.x));
-            av_store(facet, 1, newSVnv(THIS->stl.facet_start[i].normal.y));
-            av_store(facet, 2, newSVnv(THIS->stl.facet_start[i].normal.z));
+            av_store(facet, 0, newSVnv(THIS->stl.facet_start[i].normal(0)));
+            av_store(facet, 1, newSVnv(THIS->stl.facet_start[i].normal(1)));
+            av_store(facet, 2, newSVnv(THIS->stl.facet_start[i].normal(2)));
         }
         
         RETVAL = newRV_noinc((SV*)normals);
@@ -169,9 +169,9 @@ TriangleMesh::size()
     CODE:
         AV* size = newAV();
         av_extend(size, 2);
-        av_store(size, 0, newSVnv(THIS->stl.stats.size.x));
-        av_store(size, 1, newSVnv(THIS->stl.stats.size.y));
-        av_store(size, 2, newSVnv(THIS->stl.stats.size.z));
+        av_store(size, 0, newSVnv(THIS->stl.stats.size(0)));
+        av_store(size, 1, newSVnv(THIS->stl.stats.size(1)));
+        av_store(size, 2, newSVnv(THIS->stl.stats.size(2)));
         RETVAL = newRV_noinc((SV*)size);
     OUTPUT:
         RETVAL
@@ -216,12 +216,12 @@ TriangleMesh::cut(z, upper, lower)
 std::vector<double>
 TriangleMesh::bb3()
     CODE:
-        RETVAL.push_back(THIS->stl.stats.min.x);
-        RETVAL.push_back(THIS->stl.stats.min.y);
-        RETVAL.push_back(THIS->stl.stats.max.x);
-        RETVAL.push_back(THIS->stl.stats.max.y);
-        RETVAL.push_back(THIS->stl.stats.min.z);
-        RETVAL.push_back(THIS->stl.stats.max.z);
+        RETVAL.push_back(THIS->stl.stats.min(0));
+        RETVAL.push_back(THIS->stl.stats.min(1));
+        RETVAL.push_back(THIS->stl.stats.max(0));
+        RETVAL.push_back(THIS->stl.stats.max(1));
+        RETVAL.push_back(THIS->stl.stats.min(2));
+        RETVAL.push_back(THIS->stl.stats.max(2));
     OUTPUT:
         RETVAL