diff --git a/CMakeLists.txt b/CMakeLists.txt index a31b39c..9d2cc21 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,7 +11,9 @@ target_compile_features(pfe PRIVATE cxx_std_17) target_sources(pfe PRIVATE src/main.cc src/angles_filter.cc - src/angles_filter.h) + src/angles_filter.h + src/aspect_ratio_filter.cc + src/aspect_ratio_filter.h) target_link_libraries(pfe PRIVATE VTK::CommonCore diff --git a/src/aspect_ratio_filter.cc b/src/aspect_ratio_filter.cc new file mode 100644 index 0000000..b32019e --- /dev/null +++ b/src/aspect_ratio_filter.cc @@ -0,0 +1,69 @@ +#include "aspect_ratio_filter.h" + +#include +#include +#include +#include +#include + + +vtkStandardNewMacro(AspectRatioFilter); + + +int AspectRatioFilter::FillInputPortInformation(int port, vtkInformation *info) { + vtkUnstructuredGridAlgorithm::FillInputPortInformation(port, info); + // Ensure there is an "angles" array in the DataCell. + return 1; +} + + +vtkTypeBool AspectRatioFilter::RequestData( + vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) { + (void) request; + vtkUnstructuredGrid* input = + vtkUnstructuredGrid::GetData(inputVector[0]); + vtkUnstructuredGrid* output = + vtkUnstructuredGrid::GetData(outputVector); + + /* Pass the input structure and point data unmodified. */ + output->CopyStructure(input); + output->GetPointData()->PassData(input->GetPointData()); + + /* Pass the existing cell data unmodified. */ + vtkCellData *cellData = output->GetCellData(); + cellData->PassData(input->GetCellData()); + + /* Then add an array to it. */ + vtkNew aspectRatios; + aspectRatios->SetName("aspect_ratio"); + aspectRatios->SetNumberOfComponents(1); + aspectRatios->SetNumberOfTuples(input->GetNumberOfCells()); + + auto *angles = vtkDoubleArray::SafeDownCast(input->GetCellData()->GetArray("angles")); + + double cellAngles[12]; + double triangleRatios[4]; + for (vtkIdType cellId = 0; cellId < input->GetNumberOfCells(); cellId++) { + if (input->GetCellType(cellId) != VTK_TETRA) continue; + angles->GetTuple(cellId, cellAngles); + for (vtkIdType angleId = 0; + angleId < angles->GetNumberOfComponents(); + angleId += 3) { + double minAngle = std::min( + std::min(cellAngles[angleId], cellAngles[angleId + 1]), + cellAngles[angleId + 2]); + double maxAngle = std::max( + std::max(cellAngles[angleId], cellAngles[angleId + 1]), + cellAngles[angleId + 2]); + triangleRatios[angleId / 3] = minAngle / maxAngle; + } + double worst = std::min( + std::min(triangleRatios[0], triangleRatios[1]), + std::min(triangleRatios[2], triangleRatios[3])); + aspectRatios->SetTuple1(cellId, worst); + } + cellData->AddArray((vtkAbstractArray *) aspectRatios); + return true; +} diff --git a/src/aspect_ratio_filter.h b/src/aspect_ratio_filter.h new file mode 100644 index 0000000..e0d2f41 --- /dev/null +++ b/src/aspect_ratio_filter.h @@ -0,0 +1,23 @@ +#ifndef ASPECT_RATIO_FILTER_H +#define ASPECT_RATIO_FILTER_H + + +#include +#include + + +class AspectRatioFilter : public vtkUnstructuredGridAlgorithm { +public: + static AspectRatioFilter *New(); + vtkTypeMacro(AspectRatioFilter, vtkUnstructuredGridAlgorithm); + int FillInputPortInformation(int, vtkInformation *info) override; + vtkTypeBool RequestData(vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) override; + +protected: + ~AspectRatioFilter() override = default; +}; + + +#endif diff --git a/src/main.cc b/src/main.cc index a472af0..6548301 100644 --- a/src/main.cc +++ b/src/main.cc @@ -1,4 +1,5 @@ #include "angles_filter.h" +#include "aspect_ratio_filter.h" #include #include @@ -47,17 +48,16 @@ int main(int argc, char **argv) { return EXIT_FAILURE; } - vtkNew colors; - - std::array bkg{{26, 51, 102, 255}}; - colors->SetColor("BkgColor", bkg.data()); - - // vtkNew reader; + /* Reader */ vtkNew reader; reader->SetFileName(argv[1]); - // reader->Update(); + + + /* Angle filter */ vtkNew anglesFilter; anglesFilter->SetInputConnection(reader->GetOutputPort()); + + /* Display angles in console. */ anglesFilter->Update(); vtkUnstructuredGrid *grid = anglesFilter->GetOutput(); auto *angles = vtkDoubleArray::SafeDownCast(grid->GetCellData()->GetArray("angles")); @@ -68,16 +68,23 @@ int main(int argc, char **argv) { } std::cout << "\b\b \n"; } - // std::cout << "avg: " << average(angles) - // << ", stddev: " << standard_deviation(angles) - // << ", min: " << *std::min_element(angles.begin(), angles.end()) - // << ", max: " << *std::max_element(angles.begin(), angles.end()) << std::endl; + /* Aspect ratio */ + vtkNew aspectRatioFilter; + aspectRatioFilter->SetInputConnection(anglesFilter->GetOutputPort()); + + /* Display aspect ratios in console. */ + aspectRatioFilter->Update(); + grid = aspectRatioFilter->GetOutput(); + auto *aspectRatios = vtkDoubleArray::SafeDownCast(grid->GetCellData()->GetArray("aspect_ratio")); + for (ssize_t i = 0; i < aspectRatios->GetNumberOfTuples(); i++) { + std::cout << "cell " << i << ": " + << aspectRatios->GetTuple1(i) << "\n"; + } + + /* Volume rendering properties */ vtkNew volumeMapper; - volumeMapper->SetInputConnection(anglesFilter->GetOutputPort()); - // vtkNew mapper; - // mapper->SetInputConnection(reader->GetOutputPort()); - + volumeMapper->SetInputConnection(aspectRatioFilter->GetOutputPort()); vtkNew volume; volume->SetMapper(volumeMapper); vtkNew transferFunction; @@ -85,16 +92,12 @@ int main(int argc, char **argv) { transferFunction->AddPoint(1, 1); volume->GetProperty()->SetScalarOpacity(transferFunction); volume->GetProperty()->SetColor(transferFunction); - // vtkNew actor; - // actor->SetMapper(mapper); - - // actor->GetProperty()->SetColor( - // colors->GetColor4d("Tomato").GetData()); - // actor->RotateX(30.0); - // actor->RotateY(-45.0); + /* Renderer */ + vtkNew colors; + std::array bkg{{26, 51, 102, 255}}; + colors->SetColor("BkgColor", bkg.data()); vtkNew renderer; - // renderer->AddActor(actor); renderer->AddVolume(volume); renderer->SetBackground(colors->GetColor3d("BkgColor").GetData()); renderer->ResetCamera();