Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Multiregion IDs in .ele files following tetrahedralization #12

Open
bobbydawson opened this issue Oct 24, 2023 · 0 comments
Open

Multiregion IDs in .ele files following tetrahedralization #12

bobbydawson opened this issue Oct 24, 2023 · 0 comments

Comments

@bobbydawson
Copy link

bobbydawson commented Oct 24, 2023

I have been trying to implement TetGen in libigl and am struggling to produce an .ele file similar to what TetGen produces as a standalone program. Specifically, I would like the tetrahedra to have accurate region IDs. When I run the following code, the .ele file that is produced has exactly one tetrahedron with region ID = 2, with all others equal to 1. I am using the seed point to initialize the two regions as I understand it, but this might be the issue. The "Trabeculae_shell.obj" and "Endosteum_shell_inner.obj" are watertight, non-intersecting nested meshes, with the latter enveloping the former (both attached). Any help would be much appreciated.

#include <igl/read_triangle_mesh.h>
#include <igl/copyleft/cgal/CSGTree.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/jet.h>
#include <Eigen/Core>

#include <igl/copyleft/tetgen/tetrahedralize.h>

int main(int argc, char * argv[])
{
  using namespace Eigen;
  using namespace igl::copyleft::cgal;
  using namespace std;
  using namespace igl;

  // Define Mesh structure
  struct Mesh {
      Eigen::MatrixXd V;
      Eigen::MatrixXi F;
  };


  // Load meshes using read_triangle_mesh
  std::vector<std::string> mesh_files = {"Trabeculae_shell.obj", "Endosteum_shell_inner.obj"};
  std::vector<Mesh> meshes;

  for (const auto& file : mesh_files) {
      Mesh mesh;
      igl::read_triangle_mesh(file, mesh.V, mesh.F);
      meshes.push_back(mesh);
  }

  Eigen::MatrixXd V_all;
  Eigen::MatrixXi F_all;

  // Merge all meshes into V_all and F_all
  for (auto &mesh : meshes) {
      int offset = V_all.rows();
      V_all.conservativeResize(V_all.rows() + mesh.V.rows(), 3);
      F_all.conservativeResize(F_all.rows() + mesh.F.rows(), 3);

      V_all.block(offset, 0, mesh.V.rows(), 3) = mesh.V;
      F_all.block(offset, 0, mesh.F.rows(), 3) = mesh.F.array() + offset;
  }

  Eigen::MatrixXd V_tet;
  Eigen::MatrixXi T;

  Eigen::VectorXi VM,FM;
  Eigen::MatrixXd H;
  Eigen::MatrixXd R(2, 5);
  Eigen::VectorXi TM,TR,PT;
  Eigen::MatrixXi FT,TN;
  int num_regions;

  Eigen::RowVector3d trabec_point(0.146617, 0.09489, 0.204463);
  Eigen::RowVector3d endo_point(0.00708, 0.065277, 0.058339);

  R.row(0) << trabec_point(0), trabec_point(1), trabec_point(2), 1, 0;
  R.row(1) << endo_point(0), endo_point(1), endo_point(2), 2, 0;

  // Tetrahedralize
  igl::copyleft::tetgen::tetrahedralize(V_all, F_all, H, VM, FM, R, "pYA", V_tet, T, F_all, TM, TR, TN, PT, FT, num_regions);
  std::cout << "number of regions: " << num_regions << std::endl;
  int x = 0;
  int num_reg_1 = 0;
  int num_reg_2 = 0;
  for (int i = 0; i < T.rows(); ++i) {
    if (TR(i) != x) {
      std::cout << "Tetrahedron " << i << " is in Region " << TR(i) << std::endl;
      x = TR(i);
    }
    if (TR(i) == 1) {
      num_reg_1++;
    }
    else if (TR(i) == 2) {
      num_reg_2++;
    }
  }
  std::cout << "number of tetrahedra in region 1: " << num_reg_1 << std::endl;
  std::cout << "number of tetrahedra in region 2: " << num_reg_2 << std::endl;

  // Prepare the TetGenIO structure for saving
  tetgenio outmesh;

  // Fill with vertices
  outmesh.numberofpoints = V_tet.rows();
  outmesh.pointlist = new REAL[outmesh.numberofpoints * 3];
  for (int i = 0; i < V_tet.rows(); ++i) {
      for (int j = 0; j < 3; ++j) {
          outmesh.pointlist[i * 3 + j] = V_tet(i, j);
      }
  }

  // Fill with tetrahedra
  outmesh.numberoftetrahedra = T.rows();
  outmesh.tetrahedronlist = new int[outmesh.numberoftetrahedra * 4];
  for (int i = 0; i < T.rows(); ++i) {
      for (int j = 0; j < 4; ++j) {
          outmesh.tetrahedronlist[i * 4 + j] = T(i, j) + 1;  // TetGen is 1-indexed
      }
  }

  outmesh.numberofregions = num_regions;
  outmesh.regionlist = new REAL[outmesh.numberofregions * 5];

  for (int i = 0; i < outmesh.numberofregions; ++i) {
      for (int j = 0; j < 5; ++j) {
          outmesh.regionlist[i * 5 + j] = R(i, j);
      }
  }

  outmesh.numberoftetrahedronattributes = 1;
  outmesh.tetrahedronattributelist = new REAL[outmesh.numberoftetrahedra];
  for (int i = 0; i < outmesh.numberoftetrahedra; ++i) {
      outmesh.tetrahedronattributelist[i] = TR(i);
  }


  // Save to .node and .ele files
  outmesh.save_nodes("output_file");
  outmesh.save_elements("output_file");

  // Cleanup
  delete[] outmesh.pointlist;
  delete[] outmesh.tetrahedronlist;

  delete[] outmesh.regionlist;
  delete[] outmesh.tetrahedronattributelist;



  igl::opengl::glfw::Viewer viewer;
  viewer.data().set_mesh(V_tet, F_all);
  viewer.data().add_points(trabec_point, Eigen::RowVector3d(1,0,0)); // red for trabec_point
  viewer.data().add_points(endo_point, Eigen::RowVector3d(0,1,0));   // green for endo_point
  viewer.launch();

}

Endosteum_shell_inner.zip
Trabeculae_shell.zip

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant