Merge pull request #42 from Forceflow/dev

Dev

README.md

@@ -1,6 +1,6 @@
 [![Build Status](https://travis-ci.org/Forceflow/cuda_voxelizer.svg?branch=master)](https://travis-ci.org/Forceflow/cuda_voxelizer) ![](https://img.shields.io/github/license/Forceflow/cuda_voxelizer.svg) [![Donate](https://img.shields.io/badge/Donate-PayPal-green.svg)](https://www.paypal.me/Forceflow)
 
-# cuda_voxelizer v0.4.8
+# cuda_voxelizer v0.4.9
 A command-line tool to convert polygon meshes to (annotated) voxel grids.
  * Supported input formats: .ply, .off, .obj, .3DS, .SM and RAY
  * Supported output formats: .binvox, .obj, morton ordered grid
@@ -18,7 +18,7 @@ Program options:
    * `morton`: a binary file containing a Morton-ordered grid. This is a format I personally use for other tools.
  * `-cpu`: Force voxelization on the CPU instead of GPU. For when a CUDA device is not detected/compatible, or for very small models where GPU call overhead is not worth it. This is done multi-threaded, but will be slower for large models / grid sizes.
  * `-thrust` : Use Thrust library for copying the model data to the GPU, for a possible speed / throughput improvement. I found this to be very system-dependent. Default: disabled.
- * `-solid` : (Experimental) Use solid voxelization instead of voxelizing the mesh faces. Incompatible with `-cpu`. Needs a watertight input mesh.
+ * `-solid` : (Experimental) Use solid voxelization instead of voxelizing the mesh faces. Needs a watertight input mesh.
 
 
 ## Examples

msvc/vs2019/cuda_voxelizer.vcxproj

@@ -86,6 +86,8 @@ copy /y "$(SolutionDir)$(Platform)\$(Configuration)\$(TargetName).exe" "$(BINARY
       <FavorSizeOrSpeed>Speed</FavorSizeOrSpeed>
       <InlineFunctionExpansion>AnySuitable</InlineFunctionExpansion>
       <OpenMPSupport>true</OpenMPSupport>
+      <WholeProgramOptimization>false</WholeProgramOptimization>
+      <FloatingPointModel>Strict</FloatingPointModel>
     </ClCompile>
     <Link>
       <GenerateDebugInformation>true</GenerateDebugInformation>

src/cpu_voxelizer.cpp

@@ -1,6 +1,8 @@
 #include "cpu_voxelizer.h"
 #include <omp.h>
 
+#define float_error 0.000001
+
 namespace cpu_voxelizer {
 
 	// Set specific bit in voxel table
@@ -14,6 +16,8 @@ namespace cpu_voxelizer {
 		}
 	}
 
+
+
 	// Encode morton code using LUT table
 	uint64_t mortonEncode_LUT(unsigned int x, unsigned int y, unsigned int z) {
 		uint64_t answer = 0;
@@ -187,4 +191,137 @@ namespace cpu_voxelizer {
 		printf("[Debug] Marked %llu voxels as filled (includes duplicates!) on CPU \n", debug_n_voxels_marked);
 #endif
 	}
+
+	// use Xor for voxels whose corresponding bits have to flipped
+	void setBitXor(unsigned int* voxel_table, size_t index) {
+		size_t int_location = index / size_t(32);
+		unsigned int bit_pos = size_t(31) - (index % size_t(32)); // we count bit positions RtL, but array indices LtR
+		unsigned int mask = 1 << bit_pos;
+		//#pragma omp critical 
+		{
+			voxel_table[int_location] = (voxel_table[int_location] ^ mask);
+		}
+	}
+
+	bool TopLeftEdge(glm::vec2 v0, glm::vec2 v1) {
+		return ((v1.y < v0.y) || (v1.y == v0.y && v0.x > v1.x));
+	}
+
+	//check the triangle is counterclockwise or not
+	bool checkCCW(glm::vec2 v0, glm::vec2 v1, glm::vec2 v2) {
+		glm::vec2 e0 = v1 - v0;
+		glm::vec2 e1 = v2 - v0;
+		float result = e0.x * e1.y - e1.x * e0.y;
+		if (result > 0)
+			return true;
+		else
+			return false;
+	}
+
+	//find the x coordinate of the voxel
+	float get_x_coordinate(glm::vec3 n, glm::vec3 v0, glm::vec2 point) {
+		return (-(n.y * (point.x - v0.y) + n.z * (point.y - v0.z)) / n.x + v0.x);
+	}
+
+
+	//check the location with point and triangle
+	int check_point_triangle(glm::vec2 v0, glm::vec2 v1, glm::vec2 v2, glm::vec2 point) {
+		glm::vec2 PA = point - v0;
+		glm::vec2 PB = point - v1;
+		glm::vec2 PC = point - v2;
+
+		float t1 = PA.x * PB.y - PA.y * PB.x;
+		if (std::fabs(t1) < float_error && PA.x * PB.x <= 0 && PA.y * PB.y <= 0)
+			return 1;
+
+		float t2 = PB.x * PC.y - PB.y * PC.x;
+		if (std::fabs(t2) < float_error && PB.x * PC.x <= 0 && PB.y * PC.y <= 0)
+			return 2;
+
+		float t3 = PC.x * PA.y - PC.y * PA.x;
+		if (std::fabs(t3) < float_error && PC.x * PA.x <= 0 && PC.y * PA.y <= 0)
+			return 3;
+
+		if (t1 * t2 > 0 && t1 * t3 > 0)
+			return 0;
+		else
+			return -1;
+	}
+
+	// Mesh voxelization method
+	void cpu_voxelize_mesh_solid(voxinfo info, trimesh::TriMesh* themesh, unsigned int* voxel_table, bool morton_order) {
+		Timer cpu_voxelization_timer; cpu_voxelization_timer.start();
+
+		// PREPASS
+		// Move all vertices to origin (can be done in parallel)
+		trimesh::vec3 move_min = glm_to_trimesh<trimesh::vec3>(info.bbox.min);
+		for (int64_t i = 0; i < themesh->vertices.size(); i++) {
+			if (i == 0) { printf("[Info] Using %d threads \n", omp_get_num_threads()); }
+			themesh->vertices[i] = themesh->vertices[i] - move_min;
+		}
+
+		for (int64_t i = 0; i < info.n_triangles; i++) {
+			glm::vec3 v0 = trimesh_to_glm<trimesh::point>(themesh->vertices[themesh->faces[i][0]]);
+			glm::vec3 v1 = trimesh_to_glm<trimesh::point>(themesh->vertices[themesh->faces[i][1]]);
+			glm::vec3 v2 = trimesh_to_glm<trimesh::point>(themesh->vertices[themesh->faces[i][2]]);
+
+			// Edge vectors
+			glm::vec3 e0 = v1 - v0;
+			glm::vec3 e1 = v2 - v1;
+			glm::vec3 e2 = v0 - v2;
+			// Normal vector pointing up from the triangle
+			glm::vec3 n = glm::normalize(glm::cross(e0, e1));
+			if (std::fabs(n.x) < float_error) {
+				continue;
+			}
+
+			//Calculate the projection of three point into yoz plane
+			glm::vec2 v0_yz = glm::vec2(v0.y, v0.z);
+			glm::vec2 v1_yz = glm::vec2(v1.y, v1.z);
+			glm::vec2 v2_yz = glm::vec2(v2.y, v2.z);
+
+			//set the triangle counterclockwise
+			if (!checkCCW(v0_yz, v1_yz, v2_yz))
+			{
+				glm::vec2 v3 = v1_yz;
+				v1_yz = v2_yz;
+				v2_yz = v3;
+			}
+
+			// COMPUTE TRIANGLE BBOX IN GRID
+			// Triangle bounding box in world coordinates is min(v0,v1,v2) and max(v0,v1,v2)
+			glm::vec2 bbox_max = glm::max(v0_yz, glm::max(v1_yz, v2_yz));
+			glm::vec2 bbox_min = glm::min(v0_yz, glm::min(v1_yz, v2_yz));
+
+			glm::vec2 bbox_max_grid = glm::vec2(floor(bbox_max.x / info.unit.y - 0.5), floor(bbox_max.y / info.unit.z - 0.5));
+			glm::vec2 bbox_min_grid = glm::vec2(ceil(bbox_min.x / info.unit.y - 0.5), ceil(bbox_min.y / info.unit.z - 0.5));
+
+			for (int y = bbox_min_grid.x; y <= bbox_max_grid.x; y++)
+			{
+				for (int z = bbox_min_grid.y; z <= bbox_max_grid.y; z++)
+				{
+					glm::vec2 point = glm::vec2((y + 0.5) * info.unit.y, (z + 0.5) * info.unit.z);
+					int checknum = check_point_triangle(v0_yz, v1_yz, v2_yz, point);
+					if ((checknum == 1 && TopLeftEdge(v0_yz, v1_yz)) || (checknum == 2 && TopLeftEdge(v1_yz, v2_yz)) || (checknum == 3 && TopLeftEdge(v2_yz, v0_yz)) || (checknum == 0))
+					{
+						unsigned int xmax = int(get_x_coordinate(n, v0, point) / info.unit.x - 0.5);
+						for (int x = 0; x <= xmax; x++)
+						{
+							if (morton_order) {
+								size_t location = mortonEncode_LUT(x, y, z);
+								setBitXor(voxel_table, location);
+							}
+							else {
+								size_t location = static_cast<size_t>(x) + (static_cast<size_t>(y) * static_cast<size_t>(info.gridsize.y)) + (static_cast<size_t>(z) * static_cast<size_t>(info.gridsize.y) * static_cast<size_t>(info.gridsize.z));
+								setBitXor(voxel_table, location);
+							}
+							continue;
+						}
+					}
+				}
+			}
+
+		}
+		cpu_voxelization_timer.stop(); fprintf(stdout, "[Perf] CPU voxelization time: %.1f ms \n", cpu_voxelization_timer.elapsed_time_milliseconds);
+	}
 }

src/cpu_voxelizer.h

@@ -4,11 +4,13 @@
 #include <TriMesh.h>
 #include <glm/glm.hpp>
 #include <cstdio>
+#include <cmath> 
 #include "util.h"
 #include "timer.h"
 #include "morton_LUTs.h"
 
 
 namespace cpu_voxelizer {
 	void cpu_voxelize_mesh(voxinfo info, trimesh::TriMesh* themesh, unsigned int* voxel_table, bool morton_order);
+	void cpu_voxelize_mesh_solid(voxinfo info, trimesh::TriMesh* themesh, unsigned int* voxel_table, bool morton_order);
 }

src/main.cpp

@@ -20,7 +20,7 @@
 #include "cpu_voxelizer.h"
 
 using namespace std;
-string version_number = "v0.4.7";
+string version_number = "v0.4.9";
 
 // Forward declaration of CUDA functions
 float* meshToGPU_thrust(const trimesh::TriMesh *mesh); // METHOD 3 to transfer triangles can be found in thrust_operations.cu(h)
@@ -59,7 +59,7 @@ void printHelp(){
 	cout << " -o <output format: binvox, obj, obj_points or morton (default: binvox)>" << endl;
 	cout << " -thrust : Force using CUDA Thrust Library (possible speedup / throughput improvement)" << endl;
 	cout << " -cpu : Force CPU-based voxelization (slow, but works if no compatible GPU can be found)" << endl;
-	cout << " -solid : Force solid voxelization (experimental, needs watertight model, incompatible with -cpu)" << endl << endl;
+	cout << " -solid : Force solid voxelization (experimental, needs watertight model)" << endl << endl;
 	printExample();
 	cout << endl;
 }
@@ -244,17 +244,22 @@ int main(int argc, char* argv[]) {
 		if (solidVoxelization){
 			voxelize_solid(voxelization_info, device_triangles, vtable, useThrustPath, (outputformat == OutputFormat::output_morton));
 		}
-		else
-		{
+		else{
 			voxelize(voxelization_info, device_triangles, vtable, useThrustPath, (outputformat == OutputFormat::output_morton));
 		}
 	} else { 
 		// CPU VOXELIZATION FALLBACK
 		fprintf(stdout, "\n## CPU VOXELISATION \n");
 		if (!forceCPU) { fprintf(stdout, "[Info] No suitable CUDA GPU was found: Falling back to CPU voxelization\n"); }
 		else { fprintf(stdout, "[Info] Doing CPU voxelization (forced using command-line switch -cpu)\n"); }
+		// allocate zero-filled array
 		vtable = (unsigned int*) calloc(1, vtable_size);
-		cpu_voxelizer::cpu_voxelize_mesh(voxelization_info, themesh, vtable, (outputformat == OutputFormat::output_morton));
+		if (!solidVoxelization) {
+			cpu_voxelizer::cpu_voxelize_mesh(voxelization_info, themesh, vtable, (outputformat == OutputFormat::output_morton));
+		}
+		else {
+			cpu_voxelizer::cpu_voxelize_mesh_solid(voxelization_info, themesh, vtable, (outputformat == OutputFormat::output_morton));
+		}
 	}
 
 	//// DEBUG: print vtable

src/morton_LUTs.h

@@ -1,3 +1,4 @@
+#pragma once
 #include <stdint.h>
 
 // LUT tables to copy to GPU memory for quick morton decode / encode