ccolin/pfe
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Compare commits

base: ccolin:36194f715e66bfe6371a4949c2cb024485605dc4
Branches Tags
ccolin:master
..
compare: ccolin:8ee09d686e6cbc3098428cd37f516675f497d2a5
Branches Tags
ccolin:master

No commits in common. "36194f715e66bfe6371a4949c2cb024485605dc4" and "8ee09d686e6cbc3098428cd37f516675f497d2a5" have entirely different histories.
36194f715e ... 8ee09d686e

12 changed files with 67 additions and 254 deletions
Show Stats Expand all files Collapse all files

7
CMakeLists.txt
View File

@ -12,7 +12,6 @@ set(VTK_COMPONENTS
VTK::IOGeometry
VTK::IOXML
VTK::FiltersModeling
VTK::FiltersGeometry
VTK::vtksys)
set(ENABLE_VIEWER OFF CACHE BOOL "Enable the 3D viewer, depends on Qt.")
if(ENABLE_VIEWER)
@ -63,11 +62,7 @@ target_sources(pfe PRIVATE
src/remove_external_cells_filter.cc
src/remove_external_cells_filter.h
src/relaxation_filter.cc
src/relaxation_filter.h
src/max_distance_filter.cc
src/max_distance_filter.h
src/closest_polymesh_point.cc
src/closest_polymesh_point.h)
src/relaxation_filter.h)
target_link_libraries(pfe PRIVATE ${VTK_COMPONENTS})

12
data/tetrahedron.obj
View File

@ -1,12 +0,0 @@
# Blender v2.79 (sub 0) OBJ File: ''
# www.blender.org
o Solid
v 0.000000 0.769688 0.000000
v 0.000000 0.000000 -0.769688
v 0.769688 0.000000 0.000000
v 0.000000 0.000000 0.000000
s off
f 1 2 3
f 1 3 4
f 1 4 2
f 2 4 3

60
src/closest_polymesh_point.cc
View File

@ -1,60 +0,0 @@
#include "closest_polymesh_point.h"
#include "point_tris_dist.h"
#include <limits>
vtkIdType findClosestPoint(const double *point, vtkPointSet *pointSet) {
vtkIdType minPoint = 0;
double minDist2 = std::numeric_limits<double>::infinity();
for (vtkIdType i = 0; i < pointSet->GetNumberOfPoints(); i++) {
double point2[3];
pointSet->GetPoint(i, point2);
double dist2 = vtkMath::Distance2BetweenPoints(point, point2);
if (dist2 < minDist2) {
minDist2 = dist2;
minPoint = i;
}
}
return minPoint;
}
void closestPolyMeshPoint(vtkPolyData *polyMesh,
const double *point,
vtkKdTree *kdTree,
vtkStaticCellLinks *links,
double *toClosestPoint,
double *distanceToClosestPoint) {
vtkIdType closestPolyMeshPoint;
{
double dummy;
// closestPolyMeshPoint = kdTree->FindClosestPoint((double *) point, dummy);
closestPolyMeshPoint = findClosestPoint(point, polyMesh);
}
vtkIdType *cellIds = links->GetCells(closestPolyMeshPoint);
vtkIdType nCells = links->GetNcells(closestPolyMeshPoint);
if (nCells == 0) {
*distanceToClosestPoint = -1;
return;
}
double minDist = std::numeric_limits<double>::infinity();
/* For each tri the closest point belongs to. */
for (vtkIdType i = 0; i < nCells; i++) {
vtkIdType cellId = cellIds[i];
double direction[3];
double dist = pointTriangleDistance(
(double *) point, polyMesh->GetCell(cellId), direction);
/* Find the closest one. */
if (dist < minDist) {
minDist = dist;
toClosestPoint[0] = direction[0];
toClosestPoint[1] = direction[1];
toClosestPoint[2] = direction[2];
}
}
vtkMath::MultiplyScalar(toClosestPoint, minDist);
*distanceToClosestPoint = minDist;
}

17
src/closest_polymesh_point.h
View File

@ -1,17 +0,0 @@
#ifndef CLOSEST_POLYMESH_POINT_H
#define CLOSEST_POLYMESH_POINT_H
#include <vtkPolyData.h>
#include <vtkKdTree.h>
#include <vtkStaticCellLinks.h>
void closestPolyMeshPoint(vtkPolyData *polyMesh,
const double *point,
vtkKdTree *kdTree,
vtkStaticCellLinks *links,
double *toClosestPoint,
double *distanceToClosestPoint);
#endif

10
src/dihedral_angles_filter.cc
View File

@ -1,7 +1,5 @@
#include "dihedral_angles_filter.h"
#include <algorithm>
#include <limits>
#include <vtkUnstructuredGrid.h>
#include <vtkPointData.h>
#include <vtkCellData.h>
@ -38,9 +36,6 @@ vtkTypeBool DihedralAnglesFilter::RequestData(
minDihedralAngleArray->SetNumberOfTuples(input->GetNumberOfCells());
double *minDihedralAngleBase = minDihedralAngleArray->GetPointer(0);
AverageMinDegrees = 0;
MinMinDegrees = std::numeric_limits<double>::infinity();
size_t i = 0;
vtkCellIterator *it = input->NewCellIterator();
for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
@ -67,14 +62,9 @@ vtkTypeBool DihedralAnglesFilter::RequestData(
minDihedralAngleBase[i / 6] =
*std::min_element(dihedralAnglesBase + i,
dihedralAnglesBase + i + 6);
AverageMinDegrees += minDihedralAngleBase[i / 6];
MinMinDegrees = std::min(MinMinDegrees, minDihedralAngleBase[i / 6]);
i += 6;
}
AverageMinDegrees /= input->GetNumberOfCells();
AverageMinDegrees = AverageMinDegrees * 180. / 3.141592653589793;
MinMinDegrees = MinMinDegrees * 180. / 3.141592653589793;
cellData->AddArray((vtkAbstractArray *) dihedralAnglesArray);
cellData->AddArray((vtkAbstractArray *) minDihedralAngleArray);
return true;

4
src/dihedral_angles_filter.h
View File

@ -12,13 +12,9 @@ public:
vtkTypeBool RequestData(vtkInformation *request,
vtkInformationVector **inputVector,
vtkInformationVector *outputVector) override;
vtkGetMacro(AverageMinDegrees, double);
vtkGetMacro(MinMinDegrees, double);
protected:
~DihedralAnglesFilter() override = default;
double AverageMinDegrees;
double MinMinDegrees;
};

24
src/main.cc
View File

@ -3,16 +3,13 @@
#include "surface_points_filter.h"
#include "project_surface_points_on_poly.h"
#include "relaxation_filter.h"
#include "max_distance_filter.h"
#include <vtkGeometryFilter.h>
#include <vtkCellArrayIterator.h>
#include <vtkCellData.h>
#include <vtkDataSetReader.h>
#include <vtkOBJReader.h>
#include <vtkPolyData.h>
#include <vtkPolyDataReader.h>
#include <vtkSetGet.h>
#include <vtkUnstructuredGrid.h>
#include <vtkUnstructuredGridReader.h>
#include <vtkXMLPolyDataReader.h>
@ -62,16 +59,14 @@ vtkSmartPointer<vtkAlgorithm> readerFromFileName(const char *fileName) {
int main(int argc, char **argv) {
if (!(argc == 3 || argc == 5)) {
std::cerr << "Usage: " << argv[0] << " tetmesh polydata [radiusScale relaxationIterCount]" << std::endl;
if (argc < 3 || argc > 5) {
std::cerr << "Usage: " << argv[0] << " tetmesh polydata [radiusScale]" << std::endl;
return EXIT_FAILURE;
}
double radiusScale = 2.;
int iterCount = 1;
if(argc > 3) {
if(argc == 4) {
radiusScale = std::stod(argv[3]);
iterCount = std::stoi(argv[4]);
}
auto tetMeshReader = readerFromFileName(argv[1]);
@ -93,31 +88,18 @@ int main(int argc, char **argv) {
project->SetInputConnection(1, polyMeshReader->GetOutputPort());
vtkNew<RelaxationFilter> relaxationFilter;
relaxationFilter->SetIterCount(iterCount);
relaxationFilter->SetInputConnection(project->GetOutputPort());
vtkNew<DihedralAnglesFilter> dihedralAnglesFilter;
dihedralAnglesFilter->SetInputConnection(
relaxationFilter->GetOutputPort());
vtkNew<vtkGeometryFilter> geometryFilter;
geometryFilter->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
vtkNew<MaxDistanceFilter> maxDistanceFilter;
maxDistanceFilter->SetInputConnection(0, geometryFilter->GetOutputPort());
maxDistanceFilter->SetInputConnection(1, polyMeshReader->GetOutputPort());
vtkNew<vtkXMLUnstructuredGridWriter> writer;
writer->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
writer->SetDataModeToAscii();
writer->SetFileName("out.vtu");
writer->Write();
maxDistanceFilter->Update();
std::cerr << "Max distance: " << maxDistanceFilter->GetMaxDist() << "\n"
<< "Average min angle: " << dihedralAnglesFilter->GetAverageMinDegrees() << "\n"
<< "Min min angle: " << dihedralAnglesFilter->GetMinMinDegrees() << "\n";
#ifdef USE_VIEWER
/* Volume rendering properties */
vtkNew<vtkOpenGLProjectedTetrahedraMapper> volumeMapper;

63
src/max_distance_filter.cc
View File

@ -1,63 +0,0 @@
#include "max_distance_filter.h"
#include "closest_polymesh_point.h"
#include <limits>
#include <vtkPolyData.h>
#include <vtkInformation.h>
#include <vtkKdTree.h>
vtkStandardNewMacro(MaxDistanceFilter);
MaxDistanceFilter::MaxDistanceFilter() {
SetNumberOfInputPorts(2);
SetNumberOfOutputPorts(0);
}
int MaxDistanceFilter::FillInputPortInformation(int port,
vtkInformation *info) {
if (port == 0) {
vtkPolyDataAlgorithm::FillInputPortInformation(port, info);
} else if (port == 1) {
info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData");
}
return 1;
}
vtkTypeBool MaxDistanceFilter::RequestData(vtkInformation *request,
vtkInformationVector **inputVector,
vtkInformationVector *outputVector) {
(void) request;
(void) outputVector;
vtkPolyData* input1 = vtkPolyData::GetData(inputVector[0]);
vtkPolyData* input2 = vtkPolyData::GetData(inputVector[1]);
vtkNew<vtkKdTree> tree1;
tree1->BuildLocatorFromPoints(input1->GetPoints());
vtkNew<vtkKdTree> tree2;
tree2->BuildLocatorFromPoints(input2->GetPoints());
vtkNew<vtkStaticCellLinks> links1;
links1->BuildLinks(input1);
vtkNew<vtkStaticCellLinks> links2;
links2->BuildLinks(input2);
MaxDist = 0;
for (vtkIdType i = 0; i < input1->GetNumberOfPoints(); i++) {
double point[3];
input1->GetPoint(i, point);
double vec[3];
double dist;
closestPolyMeshPoint(input2, point, tree2, links2, vec, &dist);
MaxDist = std::max(dist, MaxDist);
}
for (vtkIdType i = 0; i < input2->GetNumberOfPoints(); i++) {
double point[3];
input2->GetPoint(i, point);
double vec[3];
double dist;
closestPolyMeshPoint(input1, point, tree1, links1, vec, &dist);
MaxDist = std::max(dist, MaxDist);
}
return true;
}

23
src/max_distance_filter.h
View File

@ -1,23 +0,0 @@
#ifndef MAX_DISTANCE_FILTER_H
#define MAX_DISTANCE_FILTER_H
#include <vtkPolyDataAlgorithm.h>
class MaxDistanceFilter : public vtkPolyDataAlgorithm {
public:
static MaxDistanceFilter *New();
vtkTypeMacro(MaxDistanceFilter, vtkPolyDataAlgorithm);
MaxDistanceFilter();
int FillInputPortInformation(int, vtkInformation *info) override;
vtkTypeBool RequestData(vtkInformation *request,
vtkInformationVector **inputVector,
vtkInformationVector *outputVector) override;
vtkGetMacro(MaxDist, double);
protected:
double MaxDist;
};
#endif

37
src/project_surface_points_on_poly.cc
View File

@ -1,5 +1,6 @@
#include "point_tris_dist.h"
#include "project_surface_points_on_poly.h"
#include "closest_polymesh_point.h"
#include <vtkUnstructuredGrid.h>
#include <vtkPointData.h>
@ -35,6 +36,40 @@ int ProjectSurfacePointsOnPoly::FillInputPortInformation(int port, vtkInformatio
}
static void closestPolyMeshPoint(vtkPolyData *polyMesh,
const double *point,
vtkKdTree *kdTree,
vtkStaticCellLinks *links,
double *toClosestPoint,
double *distanceToClosestPoint) {
vtkIdType closestPolyMeshPoint;
{
double dummy;
closestPolyMeshPoint = kdTree->FindClosestPoint((double *) point, dummy);
}
vtkIdType *cellIds = links->GetCells(closestPolyMeshPoint);
vtkIdType nCells = links->GetNcells(closestPolyMeshPoint);
double minDist = std::numeric_limits<double>::infinity();
/* For each tri the closest point belongs to. */
for (vtkIdType i = 0; i < nCells; i++) {
vtkIdType cellId = cellIds[i];
double direction[3];
double dist = pointTriangleDistance(
(double *) point, polyMesh->GetCell(cellId), direction);
/* Find the closest one. */
if (dist < minDist) {
minDist = dist;
toClosestPoint[0] = direction[0];
toClosestPoint[1] = direction[1];
toClosestPoint[2] = direction[2];
}
}
vtkMath::MultiplyScalar(toClosestPoint, minDist);
*distanceToClosestPoint = minDist;
}
void findPointsWithinRadius(double radius, vtkPointSet *pointSet,
const double *point, vtkIdList *points) {
for (vtkIdType i = 0; i < pointSet->GetNumberOfPoints(); i++) {

60
src/relaxation_filter.cc
View File

@ -10,10 +10,6 @@
vtkStandardNewMacro(RelaxationFilter);
RelaxationFilter::RelaxationFilter()
:IterCount(1) {
}
vtkTypeBool RelaxationFilter::RequestData(
vtkInformation *request,
vtkInformationVector **inputVector,
@ -37,45 +33,43 @@ vtkTypeBool RelaxationFilter::RequestData(
output->BuildLinks();
double avg[3];
for (int iter = 0; iter < IterCount; iter++) {
for (vtkIdType i = 0; i < surfacePoints->GetNumberOfValues(); i++) {
vtkIdType id = surfacePoints->GetValue(i);
for (vtkIdType i = 0; i < surfacePoints->GetNumberOfValues(); i++) {
vtkIdType id = surfacePoints->GetValue(i);
output->GetPointCells(id, neighborCells);
output->GetPointCells(id, neighborCells);
if (neighborCells->GetNumberOfIds() != 0) {
if (neighborCells->GetNumberOfIds() != 0) {
for (vtkIdType j = 0; j < neighborCells->GetNumberOfIds(); j++) {
for (vtkIdType j = 0; j < neighborCells->GetNumberOfIds(); j++) {
vtkIdType cellId = neighborCells->GetId(j);
output->GetCellPoints(cellId, cellPoints);
vtkIdType cellId = neighborCells->GetId(j);
output->GetCellPoints(cellId, cellPoints);
for (vtkIdType k = 0; k < cellPoints->GetNumberOfIds(); k++) {
for (vtkIdType k = 0; k < cellPoints->GetNumberOfIds(); k++) {
vtkIdType neighborId = cellPoints->GetId(k);
if (isSurface->GetValue(neighborId)) {
neighbors.insert(neighborId);
}
vtkIdType neighborId = cellPoints->GetId(k);
if (isSurface->GetValue(neighborId)) {
neighbors.insert(neighborId);
}
cellPoints->Reset();
}
neighbors.erase(neighbors.find(id));
if (neighbors.size() != 0) {
avg[0] = 0; avg[1] = 0; avg[2] = 0;
for (const vtkIdType &j : neighbors) {
double point[3];
output->GetPoints()->GetPoint(j, point);
vtkMath::Add(point, avg, avg);
}
vtkMath::MultiplyScalar(avg, 1. / neighbors.size());
newPoints->SetPoint(id, avg);
}
cellPoints->Reset();
}
neighbors.erase(neighbors.find(id));
if (neighbors.size() != 0) {
avg[0] = 0; avg[1] = 0; avg[2] = 0;
for (const vtkIdType &j : neighbors) {
double point[3];
output->GetPoints()->GetPoint(j, point);
vtkMath::Add(point, avg, avg);
}
vtkMath::MultiplyScalar(avg, 1. / neighbors.size());
newPoints->SetPoint(id, avg);
}
neighborCells->Reset();
neighbors.clear();
}
output->SetPoints(newPoints);
neighborCells->Reset();
neighbors.clear();
}
output->SetPoints(newPoints);
return true;
}

4
src/relaxation_filter.h
View File

@ -9,16 +9,12 @@ class RelaxationFilter : public vtkUnstructuredGridAlgorithm {
public:
static RelaxationFilter *New();
vtkTypeMacro(RelaxationFilter, vtkUnstructuredGridAlgorithm);
RelaxationFilter();
vtkTypeBool RequestData(vtkInformation *request,
vtkInformationVector **inputVector,
vtkInformationVector *outputVector) override;
vtkSetMacro(IterCount, int);
vtkGetMacro(IterCount, int);
protected:
~RelaxationFilter() override = default;
int IterCount;
};