diff --git a/src/admesh/connect.cpp b/src/admesh/connect.cpp
index 5ae03597e..cbd4f1d33 100644
--- a/src/admesh/connect.cpp
+++ b/src/admesh/connect.cpp
@@ -196,13 +196,13 @@ struct HashTableEdges {
 
 	// Hash table on edges
 	std::vector<HashEdge*> 	heads;
-	HashEdge* 					tail;
-	int           					M;
+	HashEdge* 				tail;
+	int           			M;
 
 #ifndef NDEBUG
-	size_t 							malloced   	= 0;
-	size_t 							freed 	  	= 0;
-	size_t 							collisions 	= 0;
+	size_t 					malloced   	= 0;
+	size_t 					freed 	  	= 0;
+	size_t 					collisions 	= 0;
 #endif /* NDEBUG */
 
 private:
diff --git a/src/admesh/normals.cpp b/src/admesh/normals.cpp
index f20d9e68e..464f9e6d8 100644
--- a/src/admesh/normals.cpp
+++ b/src/admesh/normals.cpp
@@ -27,19 +27,16 @@
 
 #include "stl.h"
 
-static int stl_check_normal_vector(stl_file *stl, int facet_num, int normal_fix_flag);
-
 static void reverse_facet(stl_file *stl, int facet_num)
 {
-	stl->stats.facets_reversed += 1;
+	++ stl->stats.facets_reversed;
 
 	int neighbor[3] = { stl->neighbors_start[facet_num].neighbor[0], stl->neighbors_start[facet_num].neighbor[1], stl->neighbors_start[facet_num].neighbor[2] };
 	int vnot[3] = { stl->neighbors_start[facet_num].which_vertex_not[0], stl->neighbors_start[facet_num].which_vertex_not[1], stl->neighbors_start[facet_num].which_vertex_not[2] };
 
 	// reverse the facet
 	stl_vertex tmp_vertex = stl->facet_start[facet_num].vertex[0];
-	stl->facet_start[facet_num].vertex[0] =
-	stl->facet_start[facet_num].vertex[1];
+	stl->facet_start[facet_num].vertex[0] = stl->facet_start[facet_num].vertex[1];
 	stl->facet_start[facet_num].vertex[1] = tmp_vertex;
 
 	// fix the vnots of the neighboring facets
@@ -64,187 +61,164 @@ static void reverse_facet(stl_file *stl, int facet_num)
 	stl->neighbors_start[facet_num].which_vertex_not[2] = (stl->neighbors_start[facet_num].which_vertex_not[2] + 3) % 6;
 }
 
-void stl_fix_normal_directions(stl_file *stl)
+// Returns true if the normal was flipped.
+static bool check_normal_vector(stl_file *stl, int facet_num, int normal_fix_flag)
 {
-  /*  int edge_num;*/
-  /*  int vnot;*/
-  int checked = 0;
-  int facet_num;
-  /*  int next_facet;*/
-  int i;
-  int j;
-  struct stl_normal {
-    int               facet_num;
-    struct stl_normal *next;
-  };
-  struct stl_normal *head;
-  struct stl_normal *tail;
-  struct stl_normal *newn;
-  struct stl_normal *temp;
+	stl_facet *facet = &stl->facet_start[facet_num];
 
-  int reversed_count = 0;
-  int id;
-  int force_exit = 0;
+	stl_normal normal;
+	stl_calculate_normal(normal, facet);
+	stl_normalize_vector(normal);
+	stl_normal normal_dif = (normal - facet->normal).cwiseAbs();
 
-  // this may happen for malformed models, see: https://github.com/prusa3d/PrusaSlicer/issues/2209
-  if (stl->stats.number_of_facets == 0) return;
+	const float eps = 0.001f;
+	if (normal_dif(0) < eps && normal_dif(1) < eps && normal_dif(2) < eps) {
+		// Normal is within tolerance. It is not really necessary to change the values here, but just for consistency, I will.
+		facet->normal = normal;
+		return false;
+	}
 
-  /* Initialize linked list. */
-  head = new stl_normal;
-  tail = new stl_normal;
-  head->next = tail;
-  tail->next = tail;
+	stl_normal test_norm = facet->normal;
+	stl_normalize_vector(test_norm);
+	normal_dif = (normal - test_norm).cwiseAbs();
+	if (normal_dif(0) < eps && normal_dif(1) < eps && normal_dif(2) < eps) {
+		// The normal is not within tolerance, but direction is OK.
+		if (normal_fix_flag) {
+	  		facet->normal = normal;
+	  		++ stl->stats.normals_fixed;
+		}
+		return false;
+	}
 
-  /* Initialize list that keeps track of already fixed facets. */
-  std::vector<char> norm_sw(stl->stats.number_of_facets, 0);
-  /* Initialize list that keeps track of reversed facets. */
-  std::vector<int> reversed_ids(stl->stats.number_of_facets, 0);
-
-  facet_num = 0;
-  /* If normal vector is not within tolerance and backwards:
-     Arbitrarily starts at face 0.  If this one is wrong, we're screwed.  Thankfully, the chances
-     of it being wrong randomly are low if most of the triangles are right: */
-  if (stl_check_normal_vector(stl, 0, 0) == 2) {
-      reverse_facet(stl, 0);
-      reversed_ids[reversed_count++] = 0;
-  }
-
-  /* Say that we've fixed this facet: */
-  norm_sw[facet_num] = 1;
-  checked++;
-
-  for(;;) {
-    /* Add neighbors_to_list.
-       Add unconnected neighbors to the list:a  */
-    for(j = 0; j < 3; j++) {
-      /* Reverse the neighboring facets if necessary. */
-      if(stl->neighbors_start[facet_num].which_vertex_not[j] > 2) {
-        /* If the facet has a neighbor that is -1, it means that edge isn't shared by another facet */
-        if(stl->neighbors_start[facet_num].neighbor[j] != -1) {
-            if (norm_sw[stl->neighbors_start[facet_num].neighbor[j]] == 1) {
-                /* trying to modify a facet already marked as fixed, revert all changes made until now and exit (fixes: #716, #574, #413, #269, #262, #259, #230, #228, #206) */
-                for (id = reversed_count - 1; id >= 0; --id) {
-                    reverse_facet(stl, reversed_ids[id]);
-                }
-                force_exit = 1;
-                break;
-            } else {
-                reverse_facet(stl, stl->neighbors_start[facet_num].neighbor[j]);
-                reversed_ids[reversed_count++] = stl->neighbors_start[facet_num].neighbor[j];
-            }
-        }
-      }
-      /* If this edge of the facet is connected: */
-      if(stl->neighbors_start[facet_num].neighbor[j] != -1) {
-        /* If we haven't fixed this facet yet, add it to the list: */
-        if(norm_sw[stl->neighbors_start[facet_num].neighbor[j]] != 1) {
-          /* Add node to beginning of list. */
-          newn = new stl_normal;
-          newn->facet_num = stl->neighbors_start[facet_num].neighbor[j];
-          newn->next = head->next;
-          head->next = newn;
-        }
-      }
-    }
-
-    /* an error occourred, quit the for loop and exit */
-    if (force_exit) break;
-
-    /* Get next facet to fix from top of list. */
-    if(head->next != tail) {
-      facet_num = head->next->facet_num;
-      if(norm_sw[facet_num] != 1) { /* If facet is in list mutiple times */
-        norm_sw[facet_num] = 1; /* Record this one as being fixed. */
-        checked++;
-      }
-      temp = head->next;	/* Delete this facet from the list. */
-      head->next = head->next->next;
-      delete temp;
-    } else { /* if we ran out of facets to fix: */
-      /* All of the facets in this part have been fixed. */
-      stl->stats.number_of_parts += 1;
-      if(checked >= stl->stats.number_of_facets) {
-        /* All of the facets have been checked.  Bail out. */
-        break;
-      } else {
-        /* There is another part here.  Find it and continue. */
-        for(i = 0; i < stl->stats.number_of_facets; i++) {
-          if(norm_sw[i] == 0) {
-            /* This is the first facet of the next part. */
-            facet_num = i;
-            if(stl_check_normal_vector(stl, i, 0) == 2) {
-                reverse_facet(stl, i);
-                reversed_ids[reversed_count++] = i;
-            }
-
-            norm_sw[facet_num] = 1;
-            checked++;
-            break;
-          }
-        }
-      }
-    }
-  }
-  delete head;
-  delete tail;
+	test_norm *= -1.f;
+	normal_dif = (normal - test_norm).cwiseAbs();
+	if (normal_dif(0) < eps && normal_dif(1) < eps && normal_dif(2) < eps) {
+		// The normal is not within tolerance and backwards.
+		if (normal_fix_flag) {
+	  		facet->normal = normal;
+	  		++ stl->stats.normals_fixed;
+		}
+		return true;
+	}
+	if (normal_fix_flag) {
+		facet->normal = normal;
+		++ stl->stats.normals_fixed;
+	}
+	// Status is unknown.
+	return false;
 }
 
-static int stl_check_normal_vector(stl_file *stl, int facet_num, int normal_fix_flag)
+void stl_fix_normal_directions(stl_file *stl)
 {
-  /* Returns 0 if the normal is within tolerance */
-  /* Returns 1 if the normal is not within tolerance, but direction is OK */
-  /* Returns 2 if the normal is not within tolerance and backwards */
-  /* Returns 4 if the status is unknown. */
+ 	// This may happen for malformed models, see: https://github.com/prusa3d/PrusaSlicer/issues/2209
+  	if (stl->stats.number_of_facets == 0)
+  		return;
 
-  stl_facet *facet;
+	struct stl_normal {
+    	int         facet_num;
+    	stl_normal *next;
+  	};
 
-  facet = &stl->facet_start[facet_num];
+  	// Initialize linked list.
+   	stl_normal *head = new stl_normal;
+  	stl_normal *tail = new stl_normal;
+	head->next = tail;
+	tail->next = tail;
 
-  stl_normal normal;
-  stl_calculate_normal(normal, facet);
-  stl_normalize_vector(normal);
-  stl_normal normal_dif = (normal - facet->normal).cwiseAbs();
+	// Initialize list that keeps track of already fixed facets.
+	std::vector<char> norm_sw(stl->stats.number_of_facets, 0);
+	// Initialize list that keeps track of reversed facets.
+	std::vector<int>  reversed_ids(stl->stats.number_of_facets, 0);
 
-  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 = normal;
-    return 0;
-  }
+  	int facet_num = 0;
+  	int reversed_count = 0;
+  	// If normal vector is not within tolerance and backwards:
+    // Arbitrarily starts at face 0.  If this one is wrong, we're screwed. Thankfully, the chances
+    // of it being wrong randomly are low if most of the triangles are right:
+  	if (check_normal_vector(stl, 0, 0)) {
+    	reverse_facet(stl, 0);
+      	reversed_ids[reversed_count ++] = 0;
+  	}
 
-  stl_normal test_norm = facet->normal;
-  stl_normalize_vector(test_norm);
-  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 = normal;
-      stl->stats.normals_fixed += 1;
-    }
-    return 1;
-  }
+  	// Say that we've fixed this facet:
+  	norm_sw[facet_num] = 1;
+	int checked = 1;
 
-  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 = normal;
-      stl->stats.normals_fixed += 1;
-    }
-    return 2;
-  }
-  if(normal_fix_flag) {
-    facet->normal = normal;
-    stl->stats.normals_fixed += 1;
-  }
-  return 4;
+  	for (;;) {
+    	// Add neighbors_to_list. Add unconnected neighbors to the list.
+    	bool force_exit = false;
+    	for (int j = 0; j < 3; ++ j) {
+      		// Reverse the neighboring facets if necessary.
+      		if (stl->neighbors_start[facet_num].which_vertex_not[j] > 2) {
+        		// If the facet has a neighbor that is -1, it means that edge isn't shared by another facet
+        		if (stl->neighbors_start[facet_num].neighbor[j] != -1) {
+            		if (norm_sw[stl->neighbors_start[facet_num].neighbor[j]] == 1) {
+                		// trying to modify a facet already marked as fixed, revert all changes made until now and exit (fixes: #716, #574, #413, #269, #262, #259, #230, #228, #206)
+                		for (int id = reversed_count - 1; id >= 0; -- id)
+                    		reverse_facet(stl, reversed_ids[id]);
+                		force_exit = true;
+                		break;
+            		}
+            		reverse_facet(stl, stl->neighbors_start[facet_num].neighbor[j]);
+            		reversed_ids[reversed_count ++] = stl->neighbors_start[facet_num].neighbor[j];
+        		}
+      		}
+      		// If this edge of the facet is connected:
+      		if (stl->neighbors_start[facet_num].neighbor[j] != -1) {
+        		// If we haven't fixed this facet yet, add it to the list:
+        		if (norm_sw[stl->neighbors_start[facet_num].neighbor[j]] != 1) {
+	          		// Add node to beginning of list.
+	          		stl_normal *newn = new stl_normal;
+	          		newn->facet_num = stl->neighbors_start[facet_num].neighbor[j];
+	          		newn->next = head->next;
+	          		head->next = newn;
+	        	}
+	      	}
+	    }
+
+    	// an error occourred, quit the for loop and exit
+    	if (force_exit)
+    		break;
+
+    	// Get next facet to fix from top of list.
+    	if (head->next != tail) {
+      		facet_num = head->next->facet_num;
+      		if (norm_sw[facet_num] != 1) { // If facet is in list mutiple times
+        		norm_sw[facet_num] = 1; // Record this one as being fixed.
+        		++ checked;
+      		}
+      		stl_normal *temp = head->next;	// Delete this facet from the list.
+      		head->next = head->next->next;
+      		delete temp;
+    	} else { // If we ran out of facets to fix: All of the facets in this part have been fixed.
+      		++ stl->stats.number_of_parts;
+      		if (checked >= stl->stats.number_of_facets)
+        		// All of the facets have been checked.  Bail out.
+        		break;
+    		// There is another part here.  Find it and continue.
+    		for (uint32_t i = 0; i < stl->stats.number_of_facets; ++ i)
+      			if (norm_sw[i] == 0) {
+        			// This is the first facet of the next part.
+        			facet_num = i;
+        			if (check_normal_vector(stl, i, 0)) {
+            			reverse_facet(stl, i);
+            			reversed_ids[reversed_count++] = i;
+        			}
+        			norm_sw[facet_num] = 1;
+        			++ checked;
+        			break;
+      			}
+    	}
+  	}
+
+	delete head;
+	delete tail;
 }
 
 void stl_fix_normal_values(stl_file *stl)
 {
 	for (uint32_t i = 0; i < stl->stats.number_of_facets; ++ i)
-    	stl_check_normal_vector(stl, i, 1);
+    	check_normal_vector(stl, i, 1);
 }
 
 void stl_reverse_all_facets(stl_file *stl)
diff --git a/src/admesh/shared.cpp b/src/admesh/shared.cpp
index fe6d5e656..78179e863 100644
--- a/src/admesh/shared.cpp
+++ b/src/admesh/shared.cpp
@@ -33,7 +33,7 @@
 void stl_generate_shared_vertices(stl_file *stl, indexed_triangle_set &its)
 {
 	// 3 indices to vertex per face
-	its.indices.assign(stl->stats.number_of_facets, v_indices_struct());
+	its.indices.assign(stl->stats.number_of_facets, stl_triangle_vertex_indices(-1, -1, -1));
 	// Shared vertices (3D coordinates)
 	its.vertices.clear();
 	its.vertices.reserve(stl->stats.number_of_facets / 2);
@@ -46,7 +46,7 @@ void stl_generate_shared_vertices(stl_file *stl, indexed_triangle_set &its)
 
 	for (uint32_t facet_idx = 0; facet_idx < stl->stats.number_of_facets; ++ facet_idx) {
 		for (int j = 0; j < 3; ++ j) {
-			if (its.indices[facet_idx].vertex[j] != -1)
+			if (its.indices[facet_idx][j] != -1)
 				// Shared vertex was already assigned.
 				continue;
 			// Create a new shared vertex.
@@ -84,7 +84,7 @@ void stl_generate_shared_vertices(stl_file *stl, indexed_triangle_set &its)
 			    		next_edge    = pivot_vertex;
 			  		}
 				}
-				its.indices[facet_in_fan_idx].vertex[pivot_vertex] = its.vertices.size() - 1;
+				its.indices[facet_in_fan_idx][pivot_vertex] = its.vertices.size() - 1;
 				fan_traversal_facet_visited[facet_in_fan_idx] = fan_traversal_stamp;
 
 				// next_edge is an index of the starting vertex of the edge, not an index of the opposite vertex to the edge!
@@ -139,7 +139,7 @@ bool its_write_off(const indexed_triangle_set &its, const char *file)
 	for (int i = 0; i < its.vertices.size(); ++ i)
 		fprintf(fp, "\t%f %f %f\n", its.vertices[i](0), its.vertices[i](1), its.vertices[i](2));
 	for (uint32_t i = 0; i < its.indices.size(); ++ i)
-		fprintf(fp, "\t3 %d %d %d\n", its.indices[i].vertex[0], its.indices[i].vertex[1], its.indices[i].vertex[2]);
+		fprintf(fp, "\t3 %d %d %d\n", its.indices[i][0], its.indices[i][1], its.indices[i][2]);
 	fclose(fp);
 	return true;
 }
@@ -177,8 +177,8 @@ bool its_write_vrml(const indexed_triangle_set &its, const char *file)
 	fprintf(fp, "\t\t\tcoordIndex [\n");
 
 	for (size_t i = 0; i + 1 < its.indices.size(); ++ i)
-		fprintf(fp, "\t\t\t\t%d, %d, %d, -1,\n", its.indices[i].vertex[0], its.indices[i].vertex[1], its.indices[i].vertex[2]);
-	fprintf(fp, "\t\t\t\t%d, %d, %d, -1]\n", its.indices[i].vertex[0], its.indices[i].vertex[1], its.indices[i].vertex[2]);
+		fprintf(fp, "\t\t\t\t%d, %d, %d, -1,\n", its.indices[i][0], its.indices[i][1], its.indices[i][2]);
+	fprintf(fp, "\t\t\t\t%d, %d, %d, -1]\n", its.indices[i][0], its.indices[i][1], its.indices[i][2]);
 	fprintf(fp, "\t\t}\n");
 	fprintf(fp, "\t}\n");
 	fprintf(fp, "}\n");
@@ -198,7 +198,7 @@ bool its_write_obj(const indexed_triangle_set &its, const char *file)
 	for (size_t i = 0; i < its.vertices.size(); ++ i)
     	fprintf(fp, "v %f %f %f\n", its.vertices[i](0), its.vertices[i](1), its.vertices[i](2));
   	for (size_t i = 0; i < its.indices.size(); ++ i)
-    	fprintf(fp, "f %d %d %d\n", its.indices[i].vertex[0]+1, its.indices[i].vertex[1]+1, its.indices[i].vertex[2]+1);
+    	fprintf(fp, "f %d %d %d\n", its.indices[i][0]+1, its.indices[i][1]+1, its.indices[i][2]+1);
   	fclose(fp);
   	return true;
 }
@@ -220,7 +220,7 @@ bool stl_validate(const stl_file *stl, const indexed_triangle_set &its)
     // Verify validity of neighborship data.
     for (int facet_idx = 0; facet_idx < (int)stl->stats.number_of_facets; ++ facet_idx) {
         const stl_neighbors &nbr 		= stl->neighbors_start[facet_idx];
-        const int 			*vertices 	= (its.indices.empty()) ? nullptr : its.indices[facet_idx].vertex;
+        const int 			*vertices 	= (its.indices.empty()) ? nullptr : its.indices[facet_idx];
         for (int nbr_idx = 0; nbr_idx < 3; ++ nbr_idx) {
             int nbr_face = stl->neighbors_start[facet_idx].neighbor[nbr_idx];
             assert(nbr_face < (int)stl->stats.number_of_facets);
@@ -237,10 +237,10 @@ bool stl_validate(const stl_file *stl, const indexed_triangle_set &its)
 					// Has shared vertices.
 	            	if (nbr_vnot < 3) {
 	            		// Faces facet_idx and nbr_face share two vertices accross the common edge. Faces are correctly oriented.
-						assert((its.indices[nbr_face].vertex[(nbr_vnot + 1) % 3] == vertices[(nbr_idx + 1) % 3] && its.indices[nbr_face].vertex[(nbr_vnot + 2) % 3] == vertices[nbr_idx]));
+						assert((its.indices[nbr_face][(nbr_vnot + 1) % 3] == vertices[(nbr_idx + 1) % 3] && its.indices[nbr_face][(nbr_vnot + 2) % 3] == vertices[nbr_idx]));
 					} else {
 	            		// Faces facet_idx and nbr_face share two vertices accross the common edge. Faces are incorrectly oriented, one of them is flipped.
-						assert((its.indices[nbr_face].vertex[(nbr_vnot + 2) % 3] == vertices[(nbr_idx + 1) % 3] && its.indices[nbr_face].vertex[(nbr_vnot + 1) % 3] == vertices[nbr_idx]));
+						assert((its.indices[nbr_face][(nbr_vnot + 2) % 3] == vertices[(nbr_idx + 1) % 3] && its.indices[nbr_face][(nbr_vnot + 1) % 3] == vertices[nbr_idx]));
 					}
 				}
             }
diff --git a/src/admesh/stl.h b/src/admesh/stl.h
index 500d6bfdb..fce23eb3f 100644
--- a/src/admesh/stl.h
+++ b/src/admesh/stl.h
@@ -41,6 +41,7 @@
 
 typedef Eigen::Matrix<float, 3, 1, Eigen::DontAlign> stl_vertex;
 typedef Eigen::Matrix<float, 3, 1, Eigen::DontAlign> stl_normal;
+typedef Eigen::Matrix<int,   3, 1, Eigen::DontAlign> stl_triangle_vertex_indices;
 static_assert(sizeof(stl_vertex) == 12, "size of stl_vertex incorrect");
 static_assert(sizeof(stl_normal) == 12, "size of stl_normal incorrect");
 
@@ -68,12 +69,6 @@ static_assert(sizeof(stl_facet) >= SIZEOF_STL_FACET, "size of stl_facet incorrec
 
 typedef enum {binary, ascii, inmemory} stl_type;
 
-struct stl_edge {
-	stl_vertex p1;
-	stl_vertex p2;
-	int        facet_number;
-};
-
 struct stl_neighbors {
   	stl_neighbors() { reset(); }
   	void reset() {
@@ -93,12 +88,6 @@ struct stl_neighbors {
   	char  which_vertex_not[3];
 };
 
-struct v_indices_struct {
-	// -1 means no vertex index has been assigned yet
-	v_indices_struct() { vertex[0] = -1; vertex[1] = -1; vertex[2] = -1; }
-  	int   vertex[3];
-};
-
 struct stl_stats {
 	char          header[81];
 	stl_type      type;
@@ -137,8 +126,8 @@ struct stl_file {
 struct indexed_triangle_set
 {
 	void clear() { indices.clear(); vertices.clear(); }
-	std::vector<v_indices_struct> 	indices;
-	std::vector<stl_vertex>       	vertices;
+	std::vector<stl_triangle_vertex_indices> 	indices;
+	std::vector<stl_vertex>       				vertices;
 };
 
 extern bool stl_open(stl_file *stl, const char *file);
diff --git a/src/admesh/stl_io.cpp b/src/admesh/stl_io.cpp
index dc4e4a7db..7e181716c 100644
--- a/src/admesh/stl_io.cpp
+++ b/src/admesh/stl_io.cpp
@@ -22,94 +22,55 @@
 
 #include <stdlib.h>
 #include <string.h>
-#include "stl.h"
 
 #include <boost/log/trivial.hpp>
 #include <boost/nowide/cstdio.hpp>
 #include <boost/detail/endian.hpp>
 
-#if !defined(SEEK_SET)
-#define SEEK_SET 0
-#define SEEK_CUR 1
-#define SEEK_END 2
-#endif
+#include "stl.h"
 
 void stl_stats_out(stl_file *stl, FILE *file, char *input_file)
 {
-  /* this is here for Slic3r, without our config.h
-     it won't use this part of the code anyway */
+  	// This is here for Slic3r, without our config.h it won't use this part of the code anyway.
 #ifndef VERSION
 #define VERSION "unknown"
 #endif
-  fprintf(file, "\n\
-================= Results produced by ADMesh version " VERSION " ================\n");
-  fprintf(file, "\
-Input file         : %s\n", input_file);
-  if(stl->stats.type == binary) {
-    fprintf(file, "\
-File type          : Binary STL file\n");
-  } else {
-    fprintf(file, "\
-File type          : ASCII STL file\n");
-  }
-  fprintf(file, "\
-Header             : %s\n", stl->stats.header);
-  fprintf(file, "============== Size ==============\n");
-  fprintf(file, "Min X = % f, Max X = % f\n",
-          stl->stats.min(0), stl->stats.max(0));
-  fprintf(file, "Min Y = % f, Max Y = % f\n",
-          stl->stats.min(1), stl->stats.max(1));
-  fprintf(file, "Min Z = % f, Max Z = % f\n",
-          stl->stats.min(2), stl->stats.max(2));
-
-  fprintf(file, "\
-========= Facet Status ========== Original ============ Final ====\n");
-  fprintf(file, "\
-Number of facets                 : %5d               %5d\n",
-          stl->stats.original_num_facets, stl->stats.number_of_facets);
-  fprintf(file, "\
-Facets with 1 disconnected edge  : %5d               %5d\n",
-          stl->stats.facets_w_1_bad_edge, stl->stats.connected_facets_2_edge -
-          stl->stats.connected_facets_3_edge);
-  fprintf(file, "\
-Facets with 2 disconnected edges : %5d               %5d\n",
-          stl->stats.facets_w_2_bad_edge, stl->stats.connected_facets_1_edge -
-          stl->stats.connected_facets_2_edge);
-  fprintf(file, "\
-Facets with 3 disconnected edges : %5d               %5d\n",
-          stl->stats.facets_w_3_bad_edge, stl->stats.number_of_facets -
-          stl->stats.connected_facets_1_edge);
-  fprintf(file, "\
-Total disconnected facets        : %5d               %5d\n",
-          stl->stats.facets_w_1_bad_edge + stl->stats.facets_w_2_bad_edge +
-          stl->stats.facets_w_3_bad_edge, stl->stats.number_of_facets -
-          stl->stats.connected_facets_3_edge);
-
-  fprintf(file,
-          "=== Processing Statistics ===     ===== Other Statistics =====\n");
-  fprintf(file, "\
-Number of parts       : %5d        Volume   : % f\n",
-          stl->stats.number_of_parts, stl->stats.volume);
-  fprintf(file, "\
-Degenerate facets     : %5d\n", stl->stats.degenerate_facets);
-  fprintf(file, "\
-Edges fixed           : %5d\n", stl->stats.edges_fixed);
-  fprintf(file, "\
-Facets removed        : %5d\n", stl->stats.facets_removed);
-  fprintf(file, "\
-Facets added          : %5d\n", stl->stats.facets_added);
-  fprintf(file, "\
-Facets reversed       : %5d\n", stl->stats.facets_reversed);
-  fprintf(file, "\
-Backwards edges       : %5d\n", stl->stats.backwards_edges);
-  fprintf(file, "\
-Normals fixed         : %5d\n", stl->stats.normals_fixed);
+  	fprintf(file, "\n================= Results produced by ADMesh version " VERSION " ================\n");
+  	fprintf(file, "Input file         : %s\n", input_file);
+  	if (stl->stats.type == binary)
+    	fprintf(file, "File type          : Binary STL file\n");
+  	else
+    	fprintf(file, "File type          : ASCII STL file\n");
+  	fprintf(file, "Header             : %s\n", stl->stats.header);
+  	fprintf(file, "============== Size ==============\n");
+  	fprintf(file, "Min X = % f, Max X = % f\n", stl->stats.min(0), stl->stats.max(0));
+  	fprintf(file, "Min Y = % f, Max Y = % f\n", stl->stats.min(1), stl->stats.max(1));
+  	fprintf(file, "Min Z = % f, Max Z = % f\n", stl->stats.min(2), stl->stats.max(2));
+  	fprintf(file, "========= Facet Status ========== Original ============ Final ====\n");
+  	fprintf(file, "Number of facets                 : %5d               %5d\n", stl->stats.original_num_facets, stl->stats.number_of_facets);
+  	fprintf(file, "Facets with 1 disconnected edge  : %5d               %5d\n", 
+  		stl->stats.facets_w_1_bad_edge, stl->stats.connected_facets_2_edge - stl->stats.connected_facets_3_edge);
+  	fprintf(file, "Facets with 2 disconnected edges : %5d               %5d\n",
+    	stl->stats.facets_w_2_bad_edge, stl->stats.connected_facets_1_edge - stl->stats.connected_facets_2_edge);
+  	fprintf(file, "Facets with 3 disconnected edges : %5d               %5d\n",
+        stl->stats.facets_w_3_bad_edge, stl->stats.number_of_facets - stl->stats.connected_facets_1_edge);
+  	fprintf(file, "Total disconnected facets        : %5d               %5d\n",
+		stl->stats.facets_w_1_bad_edge + stl->stats.facets_w_2_bad_edge + stl->stats.facets_w_3_bad_edge, stl->stats.number_of_facets - stl->stats.connected_facets_3_edge);
+  	fprintf(file, "=== Processing Statistics ===     ===== Other Statistics =====\n");
+  	fprintf(file, "Number of parts       : %5d        Volume   : %f\n", stl->stats.number_of_parts, stl->stats.volume);
+  	fprintf(file, "Degenerate facets     : %5d\n", stl->stats.degenerate_facets);
+  	fprintf(file, "Edges fixed           : %5d\n", stl->stats.edges_fixed);
+  	fprintf(file, "Facets removed        : %5d\n", stl->stats.facets_removed);
+  	fprintf(file, "Facets added          : %5d\n", stl->stats.facets_added);
+  	fprintf(file, "Facets reversed       : %5d\n", stl->stats.facets_reversed);
+  	fprintf(file, "Backwards edges       : %5d\n", stl->stats.backwards_edges);
+  	fprintf(file, "Normals fixed         : %5d\n", stl->stats.normals_fixed);
 }
 
 bool stl_write_ascii(stl_file *stl, const char *file, const char *label)
 {
 	FILE *fp = boost::nowide::fopen(file, "w");
-  	if (fp == NULL) {
+  	if (fp == nullptr) {
 		BOOST_LOG_TRIVIAL(error) << "stl_write_ascii: Couldn't open " << file << " for writing";
     	return false;
   	}
@@ -117,19 +78,11 @@ bool stl_write_ascii(stl_file *stl, const char *file, const char *label)
 	fprintf(fp, "solid  %s\n", label);
 
 	for (uint32_t i = 0; i < stl->stats.number_of_facets; ++ i) {
-		fprintf(fp, "  facet normal % .8E % .8E % .8E\n",
-	        stl->facet_start[i].normal(0), stl->facet_start[i].normal(1),
-	        stl->facet_start[i].normal(2));
+		fprintf(fp, "  facet normal % .8E % .8E % .8E\n", 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](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](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](0), stl->facet_start[i].vertex[2](1),
-	        stl->facet_start[i].vertex[2](2));
+		fprintf(fp, "      vertex % .8E % .8E % .8E\n", 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](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](0), stl->facet_start[i].vertex[2](1), stl->facet_start[i].vertex[2](2));
 		fprintf(fp, "    endloop\n");
 		fprintf(fp, "  endfacet\n");
 	}
@@ -142,7 +95,7 @@ bool stl_write_ascii(stl_file *stl, const char *file, const char *label)
 bool stl_print_neighbors(stl_file *stl, char *file)
 {
 	FILE *fp = boost::nowide::fopen(file, "w");
-	if (fp == NULL) {
+	if (fp == nullptr) {
 		BOOST_LOG_TRIVIAL(error) << "stl_print_neighbors: Couldn't open " << file << " for writing";
     	return false;
   	}
@@ -175,7 +128,7 @@ void stl_internal_reverse_quads(char *buf, size_t cnt)
 bool stl_write_binary(stl_file *stl, const char *file, const char *label)
 {
 	FILE *fp = boost::nowide::fopen(file, "wb");
-	if (fp == NULL) {
+	if (fp == nullptr) {
 		BOOST_LOG_TRIVIAL(error) << "stl_write_binary: Couldn't open " << file << " for writing";
     	return false;
   	}
@@ -184,6 +137,9 @@ bool stl_write_binary(stl_file *stl, const char *file, const char *label)
 	for (size_t i = strlen(label); i < LABEL_SIZE; ++ i)
 		putc(0, fp);
 
+#if !defined(SEEK_SET)
+	#define SEEK_SET 0
+#endif
 	fseek(fp, LABEL_SIZE, SEEK_SET);
 #ifdef BOOST_LITTLE_ENDIAN
 	fwrite(&stl->stats.number_of_facets, 4, 1, fp);
@@ -242,7 +198,7 @@ bool stl_write_quad_object(stl_file *stl, char *file)
 	stl_vertex color;
 
 	FILE *fp = boost::nowide::fopen(file, "w");
-	if (fp == NULL) {
+	if (fp == nullptr) {
 		BOOST_LOG_TRIVIAL(error) << "stl_write_quad_object: Couldn't open " << file << " for writing";
 		return false;
 	}
@@ -255,22 +211,10 @@ bool stl_write_quad_object(stl_file *stl, char *file)
     	case 2: color = uncon_2_color; break;
     	default: color = uncon_3_color;
 	    }
-	    fprintf(fp, "%f %f %f    %1.1f %1.1f %1.1f 1\n",
-            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](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](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](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[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](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](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](0), stl->facet_start[i].vertex[2](1), stl->facet_start[i].vertex[2](2), color(0), color(1), color(2));
   }
   fclose(fp);
   return true;
@@ -279,7 +223,7 @@ bool stl_write_quad_object(stl_file *stl, char *file)
 bool stl_write_dxf(stl_file *stl, const char *file, char *label) 
 {
 	FILE *fp = boost::nowide::fopen(file, "w");
-	if (fp == NULL) {
+	if (fp == nullptr) {
 		BOOST_LOG_TRIVIAL(error) << "stl_write_quad_object: Couldn't open " << file << " for writing";
     	return false;
   	}
@@ -292,20 +236,12 @@ bool stl_write_dxf(stl_file *stl, const char *file, char *label)
 
 	fprintf(fp, "0\nSECTION\n2\nENTITIES\n");
 
-	for (uint32_t i = 0; i < stl->stats.number_of_facets; i++) {
+	for (uint32_t 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](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](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](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](0), stl->facet_start[i].vertex[2](1),
-	        stl->facet_start[i].vertex[2](2));
+		fprintf(fp, "10\n%f\n20\n%f\n30\n%f\n", 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](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](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](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/src/admesh/util.cpp b/src/admesh/util.cpp
index bb135db95..6fff8a8ed 100644
--- a/src/admesh/util.cpp
+++ b/src/admesh/util.cpp
@@ -29,16 +29,17 @@
 
 #include "stl.h"
 
-static void stl_rotate(float *x, float *y, const double c, const double s);
-static float get_area(stl_facet *facet);
-static float get_volume(stl_file *stl);
-
 void stl_verify_neighbors(stl_file *stl)
 {
 	stl->stats.backwards_edges = 0;
 
 	for (uint32_t i = 0; i < stl->stats.number_of_facets; ++ i) {
 		for (int j = 0; j < 3; ++ j) {
+			struct stl_edge {
+				stl_vertex p1;
+				stl_vertex p2;
+				int        facet_number;
+			};
 			stl_edge edge_a;
 			edge_a.p1 = stl->facet_start[i].vertex[j];
 			edge_a.p2 = stl->facet_start[i].vertex[(j + 1) % 3];
@@ -67,164 +68,140 @@ void stl_verify_neighbors(stl_file *stl)
 
 void stl_translate(stl_file *stl, float x, float y, float z)
 {
-  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_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;
 }
 
 /* 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)
 {
-  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_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;
 }
 
 void stl_scale_versor(stl_file *stl, const stl_vertex &versor)
 {
-  // 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;
+	// 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;
 }
 
 static void calculate_normals(stl_file *stl) 
 {
-  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 = normal;
-  }
+	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 = normal;
+	}
 }
 
-void
-stl_rotate_x(stl_file *stl, float angle) {
-  int i;
-  int j;
-  double radian_angle = (angle / 180.0) * M_PI;
-  double c = cos(radian_angle);
-  double s = sin(radian_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](1),
-                 &stl->facet_start[i].vertex[j](2), c, s);
-    }
-  }
-  stl_get_size(stl);
-  calculate_normals(stl);
+static void rotate_point_2d(float *x, float *y, const double c, const double s)
+{
+	double xold = *x;
+	double yold = *y;
+	*x = float(c * xold - s * yold);
+	*y = float(s * xold + c * yold);
 }
 
-void
-stl_rotate_y(stl_file *stl, float angle) {
-  int i;
-  int j;
-  double radian_angle = (angle / 180.0) * M_PI;
-  double c = cos(radian_angle);
-  double s = sin(radian_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](2),
-                 &stl->facet_start[i].vertex[j](0), c, s);
-    }
-  }
-  stl_get_size(stl);
-  calculate_normals(stl);
+void stl_rotate_x(stl_file *stl, float angle)
+{
+	double radian_angle = (angle / 180.0) * M_PI;
+	double c = cos(radian_angle);
+	double s = sin(radian_angle);
+  	for (uint32_t i = 0; i < stl->stats.number_of_facets; ++ i)
+    	for (int j = 0; j < 3; ++ j)
+      		rotate_point_2d(&stl->facet_start[i].vertex[j](1), &stl->facet_start[i].vertex[j](2), c, s);
+  	stl_get_size(stl);
+  	calculate_normals(stl);
 }
 
-void
-stl_rotate_z(stl_file *stl, float angle) {
-  int i;
-  int j;
-  double radian_angle = (angle / 180.0) * M_PI;
-  double c = cos(radian_angle);
-  double s = sin(radian_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](0),
-                 &stl->facet_start[i].vertex[j](1), c, s);
-    }
-  }
-  stl_get_size(stl);
-  calculate_normals(stl);
+void stl_rotate_y(stl_file *stl, float angle)
+{
+	double radian_angle = (angle / 180.0) * M_PI;
+	double c = cos(radian_angle);
+	double s = sin(radian_angle);
+  	for (uint32_t i = 0; i < stl->stats.number_of_facets; ++ i)
+    	for (int j = 0; j < 3; ++ j)
+			rotate_point_2d(&stl->facet_start[i].vertex[j](2), &stl->facet_start[i].vertex[j](0), c, s);
+  	stl_get_size(stl);
+  	calculate_normals(stl);
 }
 
-
-
-static void
-stl_rotate(float *x, float *y, const double c, const double s) {
-  double xold = *x;
-  double yold = *y;
-  *x = float(c * xold - s * yold);
-  *y = float(s * xold + c * yold);
+void stl_rotate_z(stl_file *stl, float angle)
+{
+	double radian_angle = (angle / 180.0) * M_PI;
+	double c = cos(radian_angle);
+	double s = sin(radian_angle);
+  	for (uint32_t i = 0; i < stl->stats.number_of_facets; ++ i)
+    	for (int j = 0; j < 3; ++ j)
+      		rotate_point_2d(&stl->facet_start[i].vertex[j](0), &stl->facet_start[i].vertex[j](1), c, s);
+  	stl_get_size(stl);
+  	calculate_normals(stl);
 }
 
 void stl_get_size(stl_file *stl)
 {
-  if (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 = stl->stats.max - stl->stats.min;
-  stl->stats.bounding_diameter = stl->stats.size.norm();
+  	if (stl->stats.number_of_facets == 0)
+  		return;
+  	stl->stats.min = stl->facet_start[0].vertex[0];
+  	stl->stats.max = stl->stats.min;
+  	for (uint32_t 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 = stl->stats.max - stl->stats.min;
+  	stl->stats.bounding_diameter = stl->stats.size.norm();
 }
 
 void stl_mirror_xy(stl_file *stl)
 {
-  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;
-    }
-  }
-  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 */
+  	for (uint32_t 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;
+	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)
 {
-  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;
-    }
-  }
-  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 */
+  	for (uint32_t 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;
+	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)
@@ -241,55 +218,55 @@ void stl_mirror_xz(stl_file *stl)
 	stl->stats.facets_reversed -= stl->stats.number_of_facets;  // for not altering stats
 }
 
+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) */
+	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)));
+	}
+
+	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.
+	//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);
+	return 0.5f * n.dot(sum);
+}
+
 static float get_volume(stl_file *stl)
 {
-  // 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;
+  	// 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)
 {
-  stl->stats.volume = get_volume(stl);
-  if(stl->stats.volume < 0.0) {
-    stl_reverse_all_facets(stl);
-    stl->stats.volume = -stl->stats.volume;
-  }
-}
-
-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) */
-  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)));
-  }
-
-  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.
-  //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);
-  return 0.5f * n.dot(sum);
+  	stl->stats.volume = get_volume(stl);
+  	if (stl->stats.volume < 0.0) {
+    	stl_reverse_all_facets(stl);
+    	stl->stats.volume = -stl->stats.volume;
+  	}
 }
 
 void stl_repair(
diff --git a/src/libslic3r/Format/3mf.cpp b/src/libslic3r/Format/3mf.cpp
index 8298fe222..654426528 100644
--- a/src/libslic3r/Format/3mf.cpp
+++ b/src/libslic3r/Format/3mf.cpp
@@ -1928,7 +1928,7 @@ namespace Slic3r {
                 stream << "     <" << TRIANGLE_TAG << " ";
                 for (int j = 0; j < 3; ++j)
                 {
-                    stream << "v" << j + 1 << "=\"" << its.indices[i].vertex[j] + volume_it->second.first_vertex_id << "\" ";
+                    stream << "v" << j + 1 << "=\"" << its.indices[i][j] + volume_it->second.first_vertex_id << "\" ";
                 }
                 stream << "/>\n";
             }
diff --git a/src/libslic3r/Format/AMF.cpp b/src/libslic3r/Format/AMF.cpp
index dcd913864..48887bc78 100644
--- a/src/libslic3r/Format/AMF.cpp
+++ b/src/libslic3r/Format/AMF.cpp
@@ -958,7 +958,7 @@ bool store_amf(const char *path, Model *model, const DynamicPrintConfig *config)
             for (size_t i = 0; i < (int)volume->mesh.its.indices.size(); ++i) {
                 stream << "        <triangle>\n";
                 for (int j = 0; j < 3; ++j)
-                stream << "          <v" << j + 1 << ">" << volume->mesh.its.indices[i].vertex[j] + vertices_offset << "</v" << j + 1 << ">\n";
+                stream << "          <v" << j + 1 << ">" << volume->mesh.its.indices[i][j] + vertices_offset << "</v" << j + 1 << ">\n";
                 stream << "        </triangle>\n";
             }
             stream << "      </volume>\n";
diff --git a/src/libslic3r/TriangleMesh.cpp b/src/libslic3r/TriangleMesh.cpp
index ae35c8a5b..5b0ecbd00 100644
--- a/src/libslic3r/TriangleMesh.cpp
+++ b/src/libslic3r/TriangleMesh.cpp
@@ -614,8 +614,8 @@ void TriangleMeshSlicer::init(const TriangleMesh *_mesh, throw_on_cancel_callbac
     for (uint32_t facet_idx = 0; facet_idx < this->mesh->stl.stats.number_of_facets; ++ facet_idx)
         for (int i = 0; i < 3; ++ i) {
             EdgeToFace &e2f = edges_map[facet_idx*3+i];
-            e2f.vertex_low  = this->mesh->its.indices[facet_idx].vertex[i];
-            e2f.vertex_high = this->mesh->its.indices[facet_idx].vertex[(i + 1) % 3];
+            e2f.vertex_low  = this->mesh->its.indices[facet_idx][i];
+            e2f.vertex_high = this->mesh->its.indices[facet_idx][(i + 1) % 3];
             e2f.face        = facet_idx;
             // 1 based indexing, to be always strictly positive.
             e2f.face_edge   = i + 1;
@@ -852,7 +852,7 @@ TriangleMeshSlicer::FacetSliceType 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)
-    const int *vertices = this->mesh->its.indices[facet_idx].vertex;
+    const stl_triangle_vertex_indices &vertices = this->mesh->its.indices[facet_idx];
     int i = (facet.vertex[1].z() == min_z) ? 1 : ((facet.vertex[2].z() == min_z) ? 2 : 0);
 
     // These are used only if the cut plane is tilted:
diff --git a/xs/xsp/TriangleMesh.xsp b/xs/xsp/TriangleMesh.xsp
index 3e90bfefd..f3153665c 100644
--- a/xs/xsp/TriangleMesh.xsp
+++ b/xs/xsp/TriangleMesh.xsp
@@ -129,9 +129,9 @@ TriangleMesh::facets()
             AV* facet = newAV();
             av_store(facets, i, newRV_noinc((SV*)facet));
             av_extend(facet, 2);
-            av_store(facet, 0, newSVnv(THIS->its.indices[i].vertex[0]));
-            av_store(facet, 1, newSVnv(THIS->its.indices[i].vertex[1]));
-            av_store(facet, 2, newSVnv(THIS->its.indices[i].vertex[2]));
+            av_store(facet, 0, newSVnv(THIS->its.indices[i][0]));
+            av_store(facet, 1, newSVnv(THIS->its.indices[i][1]));
+            av_store(facet, 2, newSVnv(THIS->its.indices[i][2]));
         }
         
         RETVAL = newRV_noinc((SV*)facets);