CMake and VTK support

CMakeLists.txt

@@ -0,0 +1,36 @@
+cmake_minimum_required(VERSION 3.8)
+
+project(quartet LANGUAGES CXX)
+
+option(QUARTET_ENABLE_VIEWER "Build the mesh viewer" ON)
+
+add_library(quartetcore
+    src/predicates.cpp
+    src/geometry_queries.cpp
+    src/sdf.cpp
+    src/trimesh.cpp
+    src/tet_mesh.cpp
+    src/feature.cpp
+    src/read_obj.cpp
+    src/tet_quality.cpp
+    src/match_features.cpp
+    src/optimize_tet_mesh.cpp
+    src/make_signed_distance.cpp
+    src/make_tet_mesh.cpp
+)
+
+add_executable(quartet src/main.cpp)
+target_link_libraries(quartet quartetcore)
+
+if (QUARTET_ENABLE_VIEWER)
+
+    set(Open_GL_PREFERENCE "GLVND")
+    find_package(OpenGL REQUIRED)
+    find_package(GLUT REQUIRED)
+
+    add_library(quartetview src/gluvi.cpp)
+    target_link_libraries(quartetview PUBLIC OpenGL::OpenGL OpenGL::GLU GLUT::GLUT)
+
+    add_executable(view-tet src/view_tet.cpp)
+    target_link_libraries(view-tet quartetcore quartetview)
+endif()

src/main.cpp

@@ -1,17 +1,16 @@
-#include "make_tet_mesh.h"
-#include "make_signed_distance.h"
-#include "read_obj.h"
-#include "feature.h"
-#include "trimesh.h"
-#include <sstream>
-#include <fstream>
 #include <cstdio>
 #include <cstring>
 #include <ctime>
+#include <fstream>
+#include <sstream>
 
-int
-main(int argc,
-     char **argv)
+#include "feature.h"
+#include "make_signed_distance.h"
+#include "make_tet_mesh.h"
+#include "read_obj.h"
+#include "trimesh.h"
+
+int main(int argc, char** argv)
 {
     //
     // Parse input arguments
@@ -20,7 +19,7 @@ main(int argc,
     bool usage = false, failure = false;
     bool feature = false, automatic = false, surface = false,
          intermediate = false, optimize = false, unsafe = false;
-    char* featureFile = NULL, *surfaceFile = NULL;
+    char *featureFile = NULL, *surfaceFile = NULL;
     float angleThreshold;
 
     if (argc < 4 || argc > 12) {
@@ -64,7 +63,7 @@ main(int argc,
                     automatic = true;
                 } else {
                     std::printf("Error interpreting [%s] as angle threshold\n",
-                                argv[arg]);
+                        argv[arg]);
                     usage = true;
                     --arg;
                 }
@@ -116,6 +115,29 @@ main(int argc,
         }
     }
 
+    //
+    // Check and prepare output filename before the actual computation. There
+    // is nothing more frustrating than doing a long computation just to find
+    // out that we cannot write the output because the output name is wrong.
+    //
+    std::string outfile(argv[3]);
+    size_t extpos = outfile.find_last_of(".");
+    std::string outformat;
+
+    if (extpos == std::string::npos) {
+        failure = true;
+        std::printf("Error: output file extension required\n");
+    } else {
+        outformat = outfile.substr(extpos);
+        if (!(outformat == ".vtk" || outformat == ".tet")) {
+            std::printf("Error: supported output formats are .tet or .vtk\n");
+            failure = true;
+        }
+    }
+
+    // Strip off file extension for creating other files.
+    std::string filenameStr = outfile.substr(0, outfile.find_last_of('.'));
+
     if (feature && automatic) {
         std::printf("Error: Feature file and automatic feature generation "
                     "selected.  Disabling automatic feature generation "
@@ -124,7 +146,7 @@ main(int argc,
     }
 
     if (usage) {
-        std::printf("usage: quartet <input.obj> <dx> <output.tet> "
+        std::printf("usage: quartet <input.obj> <dx> <output.tet|output.vtk> "
                     "[-f <feature.feat>] [-a <angle_threshold>] "
                     "[-s <surface.obj>] [-i] [-o] \n"
                     "  -f -- Match the features in the given "
@@ -143,7 +165,6 @@ main(int argc,
         return 1;
     }
 
-
     //
     // Run quartet mesh generation
     //
@@ -161,46 +182,42 @@ main(int argc,
     }
 
     // Get the grid spacing
-    float dx=1;
+    float dx = 1;
     if (sscanf(argv[2], "%f", &dx) != 1) {
-        std::printf("error interpreting [%s] as a float grid spacing\n",
-                    argv[2]);
+        std::printf("error interpreting [%s] as a float grid spacing\n", argv[2]);
         return 4;
     }
 
     // Find bounding box
-    Vec3f xmin=surf_x[0], xmax=surf_x[0];
-    for (size_t i=1; i<surf_x.size(); ++i)
+    Vec3f xmin = surf_x[0], xmax = surf_x[0];
+    for (size_t i = 1; i < surf_x.size(); ++i)
         update_minmax(surf_x[i], xmin, xmax);
-    
+
     // Build triangle mesh data structure
     TriMesh trimesh(surf_x, surf_tri);
 
-
     // Make the level set
-	
-	// Determining dimensions of voxel grid.
-	// Round up to ensure voxel grid completely contains bounding box.
-	// Also add padding of 2 grid points around the bounding box.
-	// NOTE: We add 5 here so as to add 4 grid points of padding, as well as
-	// 1 grid point at the maximal boundary of the bounding box 
-	// ie: (xmax-xmin)/dx + 1 grid points to cover one axis of the bounding box
-	Vec3f origin=xmin-Vec3f(2*dx);
-	int ni = (int)std::ceil((xmax[0]-xmin[0])/dx)+5,
-		nj = (int)std::ceil((xmax[1]-xmin[1])/dx)+5,
-		nk = (int)std::ceil((xmax[2]-xmin[2])/dx)+5;
-	
+
+    // Determining dimensions of voxel grid.
+    // Round up to ensure voxel grid completely contains bounding box.
+    // Also add padding of 2 grid points around the bounding box.
+    // NOTE: We add 5 here so as to add 4 grid points of padding, as well as
+    // 1 grid point at the maximal boundary of the bounding box
+    // ie: (xmax-xmin)/dx + 1 grid points to cover one axis of the bounding box
+    Vec3f origin = xmin - Vec3f(2 * dx);
+    int ni = (int)std::ceil((xmax[0] - xmin[0]) / dx) + 5,
+        nj = (int)std::ceil((xmax[1] - xmin[1]) / dx) + 5,
+        nk = (int)std::ceil((xmax[2] - xmin[2]) / dx) + 5;
+
     SDF sdf(origin, dx, ni, nj, nk); // Initialize signed distance field.
     std::printf("making %dx%dx%d level set\n", ni, nj, nk);
     make_signed_distance(surf_tri, surf_x, sdf);
-
-
+
     // Then the tet mesh
     TetMesh mesh;
 
     // Read in feature file, if given
     if (feature || automatic) {
-
         FeatureSet featureSet;
 
         if (feature) {
@@ -219,93 +236,81 @@ main(int argc,
         std::printf("making tet mesh\n");
         make_tet_mesh(mesh, sdf, featureSet, optimize, intermediate, unsafe);
         std::printf("done tet mesh\n");
-        
+
         time = clock() - time;
-        std::printf("mesh generation took %f seconds.\n", 
-                    ((float)time)/CLOCKS_PER_SEC);
+        std::printf("mesh generation took %f seconds.\n",
+            ((float)time) / CLOCKS_PER_SEC);
     } else {
         clock_t time = clock();
 
         // Make tet mesh without features
         std::printf("making tet mesh\n");
-		make_tet_mesh(mesh, sdf, optimize, intermediate, unsafe);
-		std::printf("done tet mesh\n");
-        
+        make_tet_mesh(mesh, sdf, optimize, intermediate, unsafe);
+        std::printf("done tet mesh\n");
+
         time = clock() - time;
-        std::printf("mesh generation took %f seconds.\n", 
-                    ((float)time)/CLOCKS_PER_SEC);
+        std::printf("mesh generation took %f seconds.\n",
+            ((float)time) / CLOCKS_PER_SEC);
     }
 
-    // Write it out
-    std::printf("Writing mesh to file: %s\n", argv[3]);
-    bool result = mesh.writeToFile(argv[3]);
-    if (!result)
-    {
-        std::printf("ERROR: Failed to write mesh to file.\n");
+    // Write resulting mesh
+    std::printf("Writing mesh to file: %s\n", outfile.c_str());
+    bool result = false;
+    if (outformat == ".vtk") {
+        result = mesh.writeToFileVTK(outfile.c_str());
+    } else if (outformat == ".tet") {
+        result = mesh.writeToFile(outfile.c_str());
     }
 
-    // Strip off file extension for creating other files.
-    std::string filenameStr(argv[3]);
-    if (filenameStr.find_last_of('.') != std::string::npos)
-    {
-        filenameStr = filenameStr.substr(0, filenameStr.find_last_of('.'));
+    if (!result) {
+        std::printf("ERROR: Failed to write mesh to file.\n");
     }
 
     // Write out an "info" file for the mesh as well.
     std::stringstream ss;
     ss << filenameStr << ".info";
     std::printf("Writing mesh info to file: %s\n", ss.str().c_str());
     result = mesh.writeInfoToFile(ss.str().c_str());
-    if (!result)
-    {
+    if (!result) {
         std::printf("ERROR: Failed to write mesh info to file.\n");
     }
-    
+
     // Write out a Matlab file containing the dihedral angles of the mesh.
     std::vector<float> dihedralAngles;
-    dihedralAngles.reserve(mesh.tSize()*6);
-    for (size_t tIdx = 0; tIdx < mesh.tSize(); ++tIdx)
-    {
+    dihedralAngles.reserve(mesh.tSize() * 6);
+    for (size_t tIdx = 0; tIdx < mesh.tSize(); ++tIdx) {
         Tet currTet = mesh.getTet(tIdx);
         std::vector<float> currAngles;
         currTet.dihedralAngles(currAngles);
-        for (size_t i = 0; i < currAngles.size(); ++i)
-        {
+        for (size_t i = 0; i < currAngles.size(); ++i) {
             dihedralAngles.push_back(currAngles[i] * 180.0 / M_PI);
         }
     }
     std::stringstream ss2;
     ss2 << filenameStr << "_hist.m";
     std::ofstream out(ss2.str().c_str());
-    if (out.good())
-    {   
+    if (out.good()) {
         std::printf("Writing histogram to file: %s\n", ss2.str().c_str());
         out << "histValues = [ ";
-        for (size_t i = 0; i < dihedralAngles.size(); ++i)
-        {
+        for (size_t i = 0; i < dihedralAngles.size(); ++i) {
             out << dihedralAngles[i] << " ";
         }
         out << "];" << std::endl;
-    }
-    else
-    {
+    } else {
         std::printf("Failed to write histogram file: %s\n", ss2.str().c_str());
     }
 
     // Write out the exterior as a trimesh, if desired.
-    if (surface)
-    {
+    if (surface) {
         std::vector<Vec3f> surf_verts;
         std::vector<Vec3i> surf_tris;
         mesh.getBoundary(surf_verts, surf_tris);
         std::printf("Writing surface mesh to file: %s\n", surfaceFile);
         result = write_objfile(surf_verts, surf_tris, surfaceFile);
-        if (!result)
-        {
+        if (!result) {
             std::printf("ERROR: Failed to write surface mesh file.\n");
         }
     }
 
     return 0;
 }
-

src/tet_mesh.cpp

@@ -436,6 +436,49 @@ bool TetMesh::writeToFile(const char* filename) const
     return true;
 }
 
+bool TetMesh::writeToFileVTK(const char *filename) const
+{
+  std::ofstream file(filename);
+  if (!file.is_open()) {
+      // Handle error: file could not be opened
+      return false;
+  }
+
+  // Write header
+  file << "# vtk DataFile Version 2.0\n";
+  file << "quartet\n";
+  file << "ASCII\n";
+  file << "DATASET UNSTRUCTURED_GRID\n";
+
+  // Write number of points
+  file << "POINTS " << v.size() << " double\n";
+
+  // Write point coordinates
+  for (const Vec3f& point : v) {
+    file << point[0] << " " << point[1] << " " << point[2] << "\n";
+  }
+
+  // Write number of cells
+  int numCells = t.size();
+  int numCellPoints = 1 + 4;
+  int totalNumCellPoints = numCells * numCellPoints;
+  file << "CELLS " << numCells << " " << totalNumCellPoints << "\n";
+
+  // Write cell connectivity
+  for (const Vec4i& tet : t) {
+    file << "4 " << tet[0] << " " << tet[1] << " " << tet[2] << " " << tet[3] << "\n";
+  }
+
+  // Write cell types (all tetrahedra)
+  int cellType = 10; // VTK cell type for tetrahedron
+  file << "CELL_TYPES " << numCells << "\n";
+  for (int i = 0; i < numCells; ++i) {
+    file << cellType << "\n";
+  }
+
+  file.close();
+  return true;
+}
 
 bool TetMesh::writeInfoToFile(const char* filename) const
 {

src/tet_mesh.h

@@ -106,6 +106,8 @@ struct TetMesh
     // etc.) and write it to the file with the given name.
     bool writeInfoToFile(const char* filename) const;
 
+
+    bool writeToFileVTK(const char *filename) const;
 private:
     std::vector<Vec3f> v;
     std::vector<Vec4i> t;