private void voxelizeTrianglesDirect(double voxelUnit) {
		// Triangle normals
		Vector3d[] triangleNormals = new Vector3d[triangles.size()];
		// Edge normals
		Vector3d[] edgeNormals = new Vector3d[edges.size()];
		for (int i = 0; i < edgeNormals.length; i++) {
			edgeNormals[i] = new Vector3d();
		}
		// Vertex normals
		Vector3d[] vertexNormals = new Vector3d[vertices.size()];
		for (int i = 0; i < vertexNormals.length; i++) {
			vertexNormals[i] = new Vector3d();
		}

		// Calculate triangle normals
		for (int i = 0; i < triangles.size(); i++) {
			Triangle t = triangles.get(i);
			Vector3d vertex0 = vertices.get(t.v0);
			Vector3d vertex1 = vertices.get(t.v1);
			Vector3d vertex2 = vertices.get(t.v2);
			triangleNormals[i] = new Vector3d();
			triangleNormals[i].cross(new Vector3d(vertex1.x - vertex0.x,
					vertex1.y - vertex0.y, vertex1.z - vertex0.z),
					new Vector3d(vertex2.x - vertex0.x, vertex2.y - vertex0.y,
							vertex2.z - vertex0.z));
			triangleNormals[i].normalize();
		}

		// Calculate edge normals
		for (int i = 0; i < triangles.size(); i++) {
			Triangle t = triangles.get(i);
			int i01 = edges.indexOf(new Edge(t.v0, t.v1));
			int i12 = edges.indexOf(new Edge(t.v1, t.v2));
			int i20 = edges.indexOf(new Edge(t.v2, t.v0));
			edgeNormals[i01].add(triangleNormals[i]);
			edgeNormals[i12].add(triangleNormals[i]);
			edgeNormals[i20].add(triangleNormals[i]);
		}

		// Calculate angle weighted vertex normals
		for (int i = 0; i < triangles.size(); i++) {
			Triangle t = triangles.get(i);
			Vector3d vertex0 = vertices.get(t.v0);
			Vector3d vertex1 = vertices.get(t.v1);
			Vector3d vertex2 = vertices.get(t.v2);
			vertexNormals[t.v0].scaleAdd(new Vector3d(vertex1.x - vertex0.x,
					vertex1.y - vertex0.y, vertex1.z - vertex0.z)
					.angle(new Vector3d(vertex2.x - vertex0.x, vertex2.y
							- vertex0.y, vertex2.z - vertex0.z)),
					triangleNormals[i], vertexNormals[t.v0]);
			vertexNormals[t.v1].scaleAdd(new Vector3d(vertex0.x - vertex1.x,
					vertex0.y - vertex1.y, vertex0.z - vertex1.z)
					.angle(new Vector3d(vertex2.x - vertex1.x, vertex2.y
							- vertex1.y, vertex2.z - vertex1.z)),
					triangleNormals[i], vertexNormals[t.v1]);
			vertexNormals[t.v2].scaleAdd(new Vector3d(vertex1.x - vertex2.x,
					vertex1.y - vertex2.y, vertex1.z - vertex2.z)
					.angle(new Vector3d(vertex0.x - vertex2.x, vertex0.y
							- vertex2.y, vertex0.z - vertex2.z)),
					triangleNormals[i], vertexNormals[t.v2]);
		}

		for (int it = 0; it < triangles.size(); it++) {
			Triangle t = triangles.get(it);
			Vector3d vertex0 = vertices.get(t.v0);
			Vector3d vertex1 = vertices.get(t.v1);
			Vector3d vertex2 = vertices.get(t.v2);

			Vector4d trianglePlane = new Vector4d(triangleNormals[it]);
			trianglePlane.w = -triangleNormals[it].dot(vertex0);

			// PlaneBounds (inward normal)
			Vector3d n01 = new Vector3d();
			n01.cross(triangleNormals[it], new Vector3d(vertex0.x - vertex1.x,
					vertex0.y - vertex1.y, vertex0.z - vertex1.z));
			Vector4d planeBound01 = new Vector4d(n01);
			planeBound01.w = -n01.dot(vertex0);

			Vector3d n12 = new Vector3d();
			n12.cross(triangleNormals[it], new Vector3d(vertex1.x - vertex2.x,
					vertex1.y - vertex2.y, vertex1.z - vertex2.z));
			Vector4d planeBound12 = new Vector4d(n12);
			planeBound12.w = -n12.dot(vertex1);

			Vector3d n20 = new Vector3d();
			n20.cross(triangleNormals[it], new Vector3d(vertex2.x - vertex0.x,
					vertex2.y - vertex0.y, vertex2.z - vertex0.z));
			Vector4d planeBound20 = new Vector4d(n20);
			planeBound20.w = -n20.dot(vertex2);

			// Edge Bounds
			Vector3d v01 = new Vector3d();
			v01.sub(vertex1, vertex0);
			Vector4d cylinderBound01 = new Vector4d(v01);
			cylinderBound01.w = -v01.dot(vertex0);

			Vector3d v10 = new Vector3d();
			v10.sub(vertex0, vertex1);
			Vector4d cylinderBound10 = new Vector4d(v10);
			cylinderBound10.w = -v10.dot(vertex1);

			Vector3d v12 = new Vector3d();
			v12.sub(vertex1, vertex2);
			Vector4d cylinderBound12 = new Vector4d(v12);
			cylinderBound12.w = -v12.dot(vertex1);

			Vector3d v21 = new Vector3d();
			v21.sub(vertex2, vertex1);
			Vector4d cylinderBound21 = new Vector4d(v21);
			cylinderBound21.w = -v21.dot(vertex2);

			Vector3d v20 = new Vector3d();
			v20.sub(vertex2, vertex0);
			Vector4d cylinderBound20 = new Vector4d(v20);
			cylinderBound20.w = -v20.dot(vertex2);

			Vector3d v02 = new Vector3d();
			v02.sub(vertex0, vertex2);
			Vector4d cylinderBound02 = new Vector4d(v02);
			cylinderBound02.w = -v02.dot(vertex0);

			for (int i = 0; i < xResolution; i++) {
				for (int j = 0; j < yResolution; j++) {
					for (int k = 0; k < zResolution; k++) {
						Vector3d point = new Vector3d(signedDistanceFieldMin.x
								+ i * voxelUnit, signedDistanceFieldMin.y + j
								* voxelUnit, signedDistanceFieldMin.z + k
								* voxelUnit);
						double td = distance(point, trianglePlane);
						double pd01 = distance(point, planeBound01);
						double pd12 = distance(point, planeBound12);
						double pd20 = distance(point, planeBound20);
						double cd01 = distance(point, cylinderBound01);
						double cd10 = distance(point, cylinderBound10);
						double cd12 = distance(point, cylinderBound12);
						double cd21 = distance(point, cylinderBound21);
						double cd20 = distance(point, cylinderBound20);
						double cd02 = distance(point, cylinderBound02);

						if ((pd01 <= 0.0 && pd12 <= 0.0 && pd20 <= 0.0)) {
							// inside
							if (Math.abs(td) < Math
									.abs(signedDistanceField[i][j][k])) {
								signedDistanceField[i][j][k] = td;
							}
						} else if (cd01 >= 0.0 && cd10 >= 0.0 && pd01 >= 0.0) {
							// semi outside 01
							Vector3d p = new Vector3d();
							p.sub(point, vertex0);
							p.scale(p.dot(v01) / v01.dot(v01), v01);
							p.add(vertex0);
							Vector3d separation = new Vector3d();
							separation.sub(point, p);
							double distance = separation.length();
							if (distance < Math
									.abs(signedDistanceField[i][j][k])) {
								signedDistanceField[i][j][k] = distance;
								if (Math.signum(separation.dot(edgeNormals[edges
										.indexOf(new Edge(t.v0, t.v1))])) < 0.0) {
									signedDistanceField[i][j][k] = -distance;
								}
							}
						} else if (cd12 >= 0.0 && cd21 >= 0.0 && pd12 >= 0.0) {
							// semi outside 12
							Vector3d p = new Vector3d();
							p.sub(point, vertex1);
							p.scale(p.dot(v12) / v12.dot(v12), v12);
							p.add(vertex1);
							Vector3d separation = new Vector3d();
							separation.sub(point, p);
							double distance = separation.length();
							if (distance < Math
									.abs(signedDistanceField[i][j][k])) {
								signedDistanceField[i][j][k] = distance;
								if (Math.signum(separation.dot(edgeNormals[edges
										.indexOf(new Edge(t.v1, t.v2))])) < 0.0) {
									signedDistanceField[i][j][k] = --distance;
								}
							}
						} else if (cd20 >= 0.0 && cd02 >= 0.0 && pd20 >= 0.0) {
							// semi outside 20
							Vector3d p = new Vector3d();
							p.sub(point, vertex2);
							p.scale(p.dot(v20) / v20.dot(v20), v20);
							p.add(vertex2);
							Vector3d separation = new Vector3d();
							separation.sub(point, p);
							double distance = separation.length();
							if (distance < Math
									.abs(signedDistanceField[i][j][k])) {
								signedDistanceField[i][j][k] = distance;
								if (Math.signum(separation.dot(edgeNormals[edges
										.indexOf(new Edge(t.v2, t.v0))])) < 0.0) {
									signedDistanceField[i][j][k] = -distance;
								}
							}
						} else {
							// fully outside
							Vector3d sep0 = new Vector3d();
							sep0.sub(point, vertex0);
							double d0 = sep0.length();
							Vector3d sep1 = new Vector3d();
							sep1.sub(point, vertex1);
							double d1 = sep1.length();
							Vector3d sep2 = new Vector3d();
							sep2.sub(point, vertex2);
							double d2 = sep2.length();
							if (d0 < d1 && d0 < d2) {
								if (d0 < Math.abs(signedDistanceField[i][j][k])) {
									signedDistanceField[i][j][k] = d0;
									if (Math.signum(sep0
											.dot(vertexNormals[t.v0])) < 0.0) {
										signedDistanceField[i][j][k] = -d0;
									}
								}
							} else if (d1 < d2) {
								if (d1 < Math.abs(signedDistanceField[i][j][k])) {
									signedDistanceField[i][j][k] = d1;
									if (Math.signum(sep1
											.dot(vertexNormals[t.v1])) < 0.0) {
										signedDistanceField[i][j][k] = -d1;
									}
								}
							} else {
								if (d2 < Math.abs(signedDistanceField[i][j][k])) {
									signedDistanceField[i][j][k] = d2;
									if (Math.signum(sep2
											.dot(vertexNormals[t.v2])) < 0.0) {
										signedDistanceField[i][j][k] = -d2;
									}
								}
							}
						}
					}
				}
			}
		}
	}