refactor angle filter

This commit is contained in:
ccolin 2022-02-17 14:02:20 +01:00
parent 41502f0c33
commit 8e0c231708
4 changed files with 141 additions and 121 deletions

View File

@ -8,7 +8,11 @@ add_executable(pfe)
target_compile_features(pfe PRIVATE cxx_std_17)
target_sources(pfe PRIVATE src/main.cc)
target_sources(pfe PRIVATE
src/main.cc
src/angles_filter.cc
src/angles_filter.h)
target_link_libraries(pfe PRIVATE
VTK::CommonCore
VTK::ViewsCore

87
src/angles_filter.cc Normal file
View File

@ -0,0 +1,87 @@
#include "angles_filter.h"
#include <vtkUnstructuredGrid.h>
#include <vtkPointData.h>
#include <vtkCellData.h>
#include <vtkDoubleArray.h>
#include <vtkCellIterator.h>
vtkStandardNewMacro(AnglesFilter);
vtkTypeBool AnglesFilter::RequestData(
vtkInformation *request,
vtkInformationVector **inputVector,
vtkInformationVector *outputVector) {
(void) request;
vtkUnstructuredGrid* input =
vtkUnstructuredGrid::GetData(inputVector[0]);
vtkUnstructuredGrid* output =
vtkUnstructuredGrid::GetData(outputVector);
output->CopyStructure(input);
output->GetPointData()->PassData(input->GetPointData());
vtkCellData *cellData = output->GetCellData();
cellData->PassData(input->GetCellData());
vtkNew<vtkDoubleArray> anglesArray;
anglesArray->SetName("angles");
anglesArray->SetNumberOfComponents(12);
anglesArray->SetNumberOfTuples(input->GetNumberOfCells());
double *anglesBase = anglesArray->GetPointer(0);
size_t i = 0;
auto it = input->NewCellIterator();
for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
if (it->GetCellType() != VTK_TETRA) continue;
cellAngles(input, it->GetPointIds(), anglesBase + i);
i += 12;
}
cellData->AddArray((vtkAbstractArray *) anglesArray);
return true;
}
void cellAngles(vtkDataSet *dataSet, vtkIdList *idList, double *angles) {
// std::cout << "nb points: " << idList->GetNumberOfIds() << std::endl;
double a[3], b[3], c[3], d[3];
dataSet->GetPoint(idList->GetId(0), a);
dataSet->GetPoint(idList->GetId(1), b);
dataSet->GetPoint(idList->GetId(2), c);
dataSet->GetPoint(idList->GetId(3), d);
// std::cout << "ids " << idList->GetId(0)
// << " " << idList->GetId(1)
// << " " << idList->GetId(2)
// << " " << idList->GetId(3) << std::endl;
// std::cout << "coords" << std::endl
// << a[0] << ", " << a[1] << ", " << a[2] << std::endl
// << b[0] << ", " << b[1] << ", " << b[2] << std::endl
// << c[0] << ", " << c[1] << ", " << c[2] << std::endl
// << d[0] << ", " << d[1] << ", " << d[2] << std::endl;
double ab[3], ac[3], ad[3],
ba[3], bc[3], bd[3],
ca[3], cb[3], cd[3],
da[3], db[3], dc[3];
vtkMath::Subtract(b, a, ab);
vtkMath::Subtract(c, a, ac);
vtkMath::Subtract(d, a, ad);
vtkMath::Subtract(a, b, ba);
vtkMath::Subtract(c, b, bc);
vtkMath::Subtract(d, b, bd);
vtkMath::Subtract(a, c, ca);
vtkMath::Subtract(b, c, cb);
vtkMath::Subtract(d, c, cd);
vtkMath::Subtract(a, d, da);
vtkMath::Subtract(b, d, db);
vtkMath::Subtract(c, d, dc);
angles[0] = vtkMath::AngleBetweenVectors(ab, ac);
angles[1] = vtkMath::AngleBetweenVectors(ac, ad);
angles[2] = vtkMath::AngleBetweenVectors(ad, ab);
angles[3] = vtkMath::AngleBetweenVectors(ba, bc);
angles[4] = vtkMath::AngleBetweenVectors(bc, bd);
angles[5] = vtkMath::AngleBetweenVectors(bd, ba);
angles[6] = vtkMath::AngleBetweenVectors(ca, cb);
angles[7] = vtkMath::AngleBetweenVectors(cb, cd);
angles[8] = vtkMath::AngleBetweenVectors(cd, ca);
angles[9] = vtkMath::AngleBetweenVectors(da, db);
angles[10] = vtkMath::AngleBetweenVectors(db, dc);
angles[11] = vtkMath::AngleBetweenVectors(dc, da);
}

25
src/angles_filter.h Normal file
View File

@ -0,0 +1,25 @@
#ifndef ANGLES_FILTER_H
#define ANGLES_FILTER_H
#include <vtkUnstructuredGridAlgorithm.h>
#include <vtkIdList.h>
void cellAngles(vtkDataSet *dataSet, vtkIdList *idList, double *angles);
class AnglesFilter : public vtkUnstructuredGridAlgorithm {
public:
static AnglesFilter *New();
vtkTypeMacro(AnglesFilter, vtkUnstructuredGridAlgorithm);
vtkTypeBool RequestData(vtkInformation *request,
vtkInformationVector **inputVector,
vtkInformationVector *outputVector) override;
protected:
~AnglesFilter() override = default;
};
#endif

View File

@ -1,15 +1,9 @@
#include <vtkCell.h>
#include <vtkPointData.h>
#include "angles_filter.h"
#include <vtkCellData.h>
#include <vtkUnstructuredGrid.h>
#include <algorithm>
#include <vtkActor.h>
#include <vtkCamera.h>
#include <vtkCylinderSource.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkVolumeProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
@ -18,11 +12,7 @@
#include <vtkVolume.h>
#include <vtkOpenGLProjectedTetrahedraMapper.h>
#include <vtkUnstructuredGridReader.h>
#include <vtkPolyDataReader.h>
#include <vtkXMLPolyDataReader.h>
#include <vtkPiecewiseFunction.h>
#include <vtkCellIterator.h>
#include <vtkUnstructuredGridAlgorithm.h>
#include <vtkDoubleArray.h>
#include <array>
@ -30,111 +20,25 @@
#include <algorithm>
void cellAngles(vtkDataSet *dataSet, vtkIdList *idList, double *angles) {
// std::cout << "nb points: " << idList->GetNumberOfIds() << std::endl;
double a[3], b[3], c[3], d[3];
dataSet->GetPoint(idList->GetId(0), a);
dataSet->GetPoint(idList->GetId(1), b);
dataSet->GetPoint(idList->GetId(2), c);
dataSet->GetPoint(idList->GetId(3), d);
// std::cout << "ids " << idList->GetId(0)
// << " " << idList->GetId(1)
// << " " << idList->GetId(2)
// << " " << idList->GetId(3) << std::endl;
// std::cout << "coords" << std::endl
// << a[0] << ", " << a[1] << ", " << a[2] << std::endl
// << b[0] << ", " << b[1] << ", " << b[2] << std::endl
// << c[0] << ", " << c[1] << ", " << c[2] << std::endl
// << d[0] << ", " << d[1] << ", " << d[2] << std::endl;
double ab[3], ac[3], ad[3],
ba[3], bc[3], bd[3],
ca[3], cb[3], cd[3],
da[3], db[3], dc[3];
vtkMath::Subtract(b, a, ab);
vtkMath::Subtract(c, a, ac);
vtkMath::Subtract(d, a, ad);
vtkMath::Subtract(a, b, ba);
vtkMath::Subtract(c, b, bc);
vtkMath::Subtract(d, b, bd);
vtkMath::Subtract(a, c, ca);
vtkMath::Subtract(b, c, cb);
vtkMath::Subtract(d, c, cd);
vtkMath::Subtract(a, d, da);
vtkMath::Subtract(b, d, db);
vtkMath::Subtract(c, d, dc);
angles[0] = vtkMath::AngleBetweenVectors(ab, ac);
angles[1] = vtkMath::AngleBetweenVectors(ac, ad);
angles[2] = vtkMath::AngleBetweenVectors(ad, ab);
angles[3] = vtkMath::AngleBetweenVectors(ba, bc);
angles[4] = vtkMath::AngleBetweenVectors(bc, bd);
angles[5] = vtkMath::AngleBetweenVectors(bd, ba);
angles[6] = vtkMath::AngleBetweenVectors(ca, cb);
angles[7] = vtkMath::AngleBetweenVectors(cb, cd);
angles[8] = vtkMath::AngleBetweenVectors(cd, ca);
angles[9] = vtkMath::AngleBetweenVectors(da, db);
angles[10] = vtkMath::AngleBetweenVectors(db, dc);
angles[11] = vtkMath::AngleBetweenVectors(dc, da);
}
// template <typename T>
// T average(const std::vector<T> &data) {
// T avg = 0;
// for (const T &t : data) {
// avg += t;
// }
// return avg / data.size();
// }
class AddAngleComputationAlgorithm : public vtkUnstructuredGridAlgorithm {
public:
static AddAngleComputationAlgorithm *New();
vtkTypeMacro(AddAngleComputationAlgorithm, vtkUnstructuredGridAlgorithm);
vtkTypeBool RequestData(vtkInformation *request,
vtkInformationVector **inputVector,
vtkInformationVector *outputVector) override {
(void) request;
vtkUnstructuredGrid* input =
vtkUnstructuredGrid::GetData(inputVector[0]);
vtkUnstructuredGrid* output =
vtkUnstructuredGrid::GetData(outputVector);
output->CopyStructure(input);
output->GetPointData()->PassData(input->GetPointData());
vtkCellData *cellData = output->GetCellData();
cellData->PassData(input->GetCellData());
vtkNew<vtkDoubleArray> anglesArray;
anglesArray->SetName("angles");
anglesArray->SetNumberOfComponents(12);
anglesArray->SetNumberOfTuples(input->GetNumberOfCells());
double *anglesBase = anglesArray->GetPointer(0);
size_t i = 0;
auto it = input->NewCellIterator();
for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
if (it->GetCellType() != VTK_TETRA) continue;
cellAngles(input, it->GetPointIds(), anglesBase + i);
i += 12;
}
cellData->AddArray((vtkAbstractArray *) anglesArray);
return true;
}
protected:
~AddAngleComputationAlgorithm() override = default;
};
vtkStandardNewMacro(AddAngleComputationAlgorithm);
template <typename T>
T average(const std::vector<T> &data) {
T avg = 0;
for (const T &t : data) {
avg += t;
}
return avg / data.size();
}
template <typename T>
T standard_deviation(const std::vector<T> &data) {
T avg = average(data);
T stddev = 0;
for (const T &t : data) {
stddev += (t - avg) * (t - avg);
}
return std::sqrt(stddev / data.size());
}
// template <typename T>
// T standard_deviation(const std::vector<T> &data) {
// T avg = average(data);
// T stddev = 0;
// for (const T &t : data) {
// stddev += (t - avg) * (t - avg);
// }
// return std::sqrt(stddev / data.size());
// }
int main(int argc, char **argv) {
@ -152,10 +56,10 @@ int main(int argc, char **argv) {
vtkNew<vtkUnstructuredGridReader> reader;
reader->SetFileName(argv[1]);
// reader->Update();
vtkNew<AddAngleComputationAlgorithm> addAngles;
addAngles->SetInputConnection(reader->GetOutputPort());
addAngles->Update();
vtkUnstructuredGrid *grid = addAngles->GetOutput();
vtkNew<AnglesFilter> anglesFilter;
anglesFilter->SetInputConnection(reader->GetOutputPort());
anglesFilter->Update();
vtkUnstructuredGrid *grid = anglesFilter->GetOutput();
auto *angles = vtkDoubleArray::SafeDownCast(grid->GetCellData()->GetArray("angles"));
for (ssize_t i = 0; i < angles->GetNumberOfTuples(); i++) {
std::cout << "cell " << i << ": ";
@ -170,7 +74,7 @@ int main(int argc, char **argv) {
// << ", max: " << *std::max_element(angles.begin(), angles.end()) << std::endl;
vtkNew<vtkOpenGLProjectedTetrahedraMapper> volumeMapper;
volumeMapper->SetInputConnection(addAngles->GetOutputPort());
volumeMapper->SetInputConnection(anglesFilter->GetOutputPort());
// vtkNew<vtkPolyDataMapper> mapper;
// mapper->SetInputConnection(reader->GetOutputPort());