Skip to content

Commit

Permalink
#734 #735 Implement orientFacets method for meshes loaded from a .vtp…
Browse files Browse the repository at this point in the history 
… file
  • Loading branch information
@mobernabeu
mobernabeu committed Mar 27, 2020
1 parent 3d46232 commit a4764ed
Show file tree
Hide file tree
Showing 5 changed files with 119 additions and 40 deletions.
75 changes: 60 additions & 15 deletions Code/redblood/Mesh.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,8 @@
#include "constants.h"
#include <vtkXMLPolyDataReader.h>
#include <vtkSmartPointer.h>
#include <vtkPolyDataNormals.h>
#include <vtkCellData.h>

namespace hemelb
{
Expand Down Expand Up @@ -141,27 +143,36 @@ namespace hemelb
return result;
}

std::shared_ptr<MeshData> readVTKMesh(std::string const &filename)
std::shared_ptr<MeshData> readVTKMesh(std::string const &filename, bool fixFacetOrientation)
{
log::Logger::Log<log::Debug, log::Singleton>("Reading red blood cell from %s",
filename.c_str());

if (!util::file_exists(filename.c_str()))
throw Exception() << "Red-blood-cell mesh file '" << filename.c_str() << "' does not exist";

std::shared_ptr<MeshData> meshData;
vtkSmartPointer<vtkPolyData> polyData;
std::tie(meshData, polyData) = readMeshDataFromVTKPolyData(filename);

if (fixFacetOrientation)
{
orientFacets(*meshData, *polyData);
}

return meshData;
}

std::tuple<std::shared_ptr<MeshData>, vtkSmartPointer<vtkPolyData> > readMeshDataFromVTKPolyData(std::string const &filename)
{
log::Logger::Log<log::Debug, log::Singleton>("Reading red blood cell from VTK polydata file");

// Read in VTK polydata object
vtkSmartPointer<vtkXMLPolyDataReader> reader = vtkSmartPointer<vtkXMLPolyDataReader>::New();
reader->SetFileName(filename.c_str());
reader->Update();
// TODO: check that the file read succeeded, e.g. wrong format
vtkPolyData* polydata = reader->GetOutput();

return readVTKMesh(polydata);
}

std::shared_ptr<MeshData> readVTKMesh(vtkPolyData* polydata, bool fixFacetOrientation)
{
log::Logger::Log<log::Debug, log::Singleton>("Reading red blood cell from polydata object");
vtkSmartPointer<vtkPolyData> polydata(reader->GetOutput());

// Number of vertices
unsigned int num_vertices = polydata->GetNumberOfPoints();
Expand Down Expand Up @@ -207,12 +218,7 @@ namespace hemelb
num_vertices,
num_facets);

if (fixFacetOrientation)
{
// TODO: Write new implementation of orientFacets that uses the VTK algorithm to work out what's the outward pointing normal
orientFacets(*result);
}
return result;
return std::make_tuple(result, polydata);
}


Expand Down Expand Up @@ -757,6 +763,8 @@ namespace hemelb
}
void orientFacets(MeshData &mesh, bool outward)
{
log::Logger::Log<log::Warning, log::Singleton>("orientFacets method for meshes constructed from a .msh file has been deprecated. Consider using VTK input files instead.");

// create a new mesh at slighly smaller scale
MeshData smaller(mesh);
auto const scale = 0.99;
Expand Down Expand Up @@ -784,7 +792,44 @@ namespace hemelb
std::swap(facet[0], facet[2]);
}
}
}
unsigned orientFacets(MeshData &mesh, vtkPolyData &polydata, bool outward)
{
vtkSmartPointer<vtkPolyDataNormals> normalCalculator = vtkSmartPointer<vtkPolyDataNormals>::New();
normalCalculator->SetInputData(&polydata);
normalCalculator->ComputePointNormalsOff();
normalCalculator->ComputeCellNormalsOn();
normalCalculator->SetAutoOrientNormals(1);
normalCalculator->Update();
auto normals = normalCalculator->GetOutput()->GetCellData()->GetNormals();

assert(normals->GetNumberOfComponents() == 3);
assert(normals->GetNumberOfTuples() == mesh.facets.size());

// Loop over each facet, checks orientation and modify as appropriate
unsigned normalId = 0;
unsigned numSwapped = 0;
for (auto &facet : mesh.facets)
{
auto const &v0 = mesh.vertices[facet[0]];
auto const &v1 = mesh.vertices[facet[1]];
auto const &v2 = mesh.vertices[facet[2]];

double vtkNormal[3];
normals->GetTuple(normalId, vtkNormal);
LatticePosition direction(vtkNormal[0], vtkNormal[1], vtkNormal[2]);

if ( ( (v0 - v1).Cross(v2 - v1).Dot(direction) > 0e0) xor outward)
{
std::swap(facet[0], facet[2]);
++numSwapped;
}

++normalId;
}

return numSwapped;
}

}
} // hemelb::redblood
8 changes: 5 additions & 3 deletions Code/redblood/Mesh.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,8 +56,10 @@ namespace hemelb
LatticeVolume volume(MeshData::Vertices const &vertices, MeshData::Facets const &facets);
LatticeArea area(MeshData const &mesh);
LatticeArea area(MeshData::Vertices const &vertices, MeshData::Facets const &facets);
//! Orients facet inward, or inward
//! DEPRECATED. Orients facets outward, or inward. Algorithm cannot handle case of facet being coplanar with mesh barycenter.
void orientFacets(MeshData &mesh, bool outward = true);
//! Orients facets inwards/outwards using VTK algorithm to determining outward facing direction. MeshData object should have been constructed from vtkPolyData object. See readMeshDataFromVTKPolyData.
unsigned orientFacets(MeshData &mesh, vtkPolyData &polydata, bool outward = true);

//! Holds raw topology data
class MeshTopology
Expand Down Expand Up @@ -240,8 +242,8 @@ namespace hemelb
std::shared_ptr<MeshData> readMesh(std::istream &stream, util::UnitConverter const &);

//! Read VTK mesh from file
std::shared_ptr<MeshData> readVTKMesh(std::string const &filename);
std::shared_ptr<MeshData> readVTKMesh(vtkPolyData* polydata, bool fixFacetOrientation=true);
std::shared_ptr<MeshData> readVTKMesh(std::string const &filename, bool fixFacetOrientation=true);
std::tuple<std::shared_ptr<MeshData>, vtkSmartPointer<vtkPolyData> > readMeshDataFromVTKPolyData(std::string const &filename);

//! Write mesh from file
//! Format is from T. Krueger's thesis
Expand Down
1 change: 1 addition & 0 deletions Code/unittests/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,4 +91,5 @@ add_to_resources(resources/four_cube.gmy resources/four_cube.xml
resources/rbc_ico_720.msh
resources/rbc_ico_720.vtp
resources/992Particles_rank3_26_t992.msh
resources/992Particles_rank3_26_t992.vtp
)
55 changes: 33 additions & 22 deletions Code/unittests/redblood/RedBloodMeshVTKDataIOTests.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
#include "resources/Resource.h"
#include "util/UnitConverter.h"
#include "unittests/redblood/Fixtures.h"
#include <vtkPolyData.h>

namespace hemelb
{
Expand All @@ -31,14 +32,13 @@ namespace hemelb
{
CPPUNIT_TEST_SUITE (RedBloodMeshVTKDataIOTests);
CPPUNIT_TEST (testVTPReadMesh);
//CPPUNIT_TEST (testVTPWriteMesh);
CPPUNIT_TEST (testOrientOriginalRBCMesh);
CPPUNIT_TEST (testOrientTimmSimRBCMesh);
CPPUNIT_TEST_SUITE_END();

public:
void setUp()
{
std::string filename = resources::Resource("rbc_ico_720.vtp").Path();
mesh = readVTKMesh(filename);
}

void tearDown()
Expand All @@ -47,6 +47,9 @@ namespace hemelb

void testVTPReadMesh()
{
std::string filename = resources::Resource("rbc_ico_720.vtp").Path();
std::shared_ptr<MeshData> mesh = readVTKMesh(filename);

CPPUNIT_ASSERT(mesh);
CPPUNIT_ASSERT_EQUAL(mesh->vertices.size(), 362ul);
CPPUNIT_ASSERT_EQUAL(mesh->facets.size(), 720ul);
Expand All @@ -65,23 +68,33 @@ namespace hemelb
CPPUNIT_ASSERT(any(mesh->facets.back(), 157));
}

// void testWriteMesh()
// {
// std::ostringstream output;
// writeMesh(output, *mesh, util::UnitConverter(1, 1, LatticePosition(0, 0, 0)));
// std::istringstream input(output.str());
// std::shared_ptr<MeshData> other = readMesh(input);
// CPPUNIT_ASSERT(other->vertices.size() == mesh->vertices.size());
// CPPUNIT_ASSERT(other->facets.size() == mesh->facets.size());
// CPPUNIT_ASSERT(compare(mesh->vertices.front() - other->vertices.front()));
// CPPUNIT_ASSERT(compare(mesh->vertices.back() - other->vertices.back()));
//
// for (size_t i(0); i < 3; ++i)
// {
// CPPUNIT_ASSERT(mesh->facets.front()[i] == other->facets.front()[i]);
// CPPUNIT_ASSERT(mesh->facets.back()[i] == other->facets.back()[i]);
// }
// }
void testOrientOriginalRBCMesh()
{
// This file has 684 out of 720 faces oriented inwards
std::string filename = resources::Resource("rbc_ico_720.vtp").Path();

std::shared_ptr<MeshData> meshData;
vtkSmartPointer<vtkPolyData> polyData;
std::tie(meshData, polyData) = readMeshDataFromVTKPolyData(filename);

auto numSwaps = orientFacets(*meshData, *polyData);

CPPUNIT_ASSERT_EQUAL(numSwaps, 684u);
}

void testOrientTimmSimRBCMesh()
{
// This file has all 720 faces oriented inwards
std::string filename = resources::Resource("992Particles_rank3_26_t992.vtp").Path();

std::shared_ptr<MeshData> meshData;
vtkSmartPointer<vtkPolyData> polyData;
std::tie(meshData, polyData) = readMeshDataFromVTKPolyData(filename);

auto numSwaps = orientFacets(*meshData, *polyData);

CPPUNIT_ASSERT_EQUAL(numSwaps, 720u);
}

static bool compare(util::Vector3D<double> const &in)
{
Expand All @@ -92,8 +105,6 @@ namespace hemelb
return vec[0] == value or vec[1] == value or vec[2] == value;
}

private:
std::shared_ptr<MeshData> mesh;
};


Expand Down
Loading

0 comments on commit a4764ed

Please sign in to comment.