Hi,
While perfoming triangulation on XYZ data with 120 thousands of points I noticed the application crashes.
But it works fine with 100k of points.
It doesn’t fill much RAM so I don’t think the problem is connected with it.
Also I tested on python I can perform griddedinterpolation without problem.
#include <vtkActor.h>
#include <vtkDelaunay2D.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkPoints.h>
#include <vtkPolyData.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkVertexGlyphFilter.h>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>
// STD includes
#include <iostream>
#include <cassert>
#include <optional>
#include <filesystem>
namespace fs = std::filesystem;
// h5gt includes
#include <h5gt/H5File.hpp>
#include <h5gt/H5Group.hpp>
#include <h5gt/H5DataSet.hpp>
//Eigen includes
#include <Eigen/Dense>
int main(int, char*[])
{
vtkNew<vtkNamedColors> colors;
std::string h5FileName("C:/Users/my_name/dev/notebooks/predicted_data.h5");
std::string dsetPath("/data/block0_values");
size_t xcol = 0;
size_t ycol = 1;
size_t zcol = 2;
if (!fs::exists(h5FileName) || H5Fis_hdf5(h5FileName.c_str()) < 1){
std::cout << "Load: Input file is not an HDF5";
return -1;
}
std::optional<h5gt::File> h5File;
try {
h5File.emplace(h5FileName, h5gt::File::ReadOnly);
} catch (h5gt::Exception& err) {
std::cout << "Load: Unable to open HDF5 file: " << h5FileName;
return -1;
}
if (!h5File.has_value() ||
!h5File->hasObject(dsetPath, h5gt::ObjectType::Dataset)){
std::cout << "Load: Unable to find dataset: " << dsetPath;
return -1;
}
h5gt::DataSet dset = h5File->getDataSet(dsetPath);
std::vector dims = dset.getDimensions();
Eigen::VectorXd x(dims[0]);
Eigen::VectorXd y(dims[0]);
Eigen::VectorXd z(dims[0], 1);
try {
dset.select({0, xcol}, {dims[0], 1}).
read(x.data());
dset.select({0, ycol}, {dims[0], 1}).
read(y.data());
dset.select({0, zcol}, {dims[0], 1}).
read(z.data());
} catch (h5gt::Exception e) {
std::cout << "Load: Unable to read XYZ data from dataset: " << dsetPath;
return -1;
}
vtkIdType npts = x.size();
vtkNew<vtkPoints> points;
for (ptrdiff_t i = 0; i < 120000; i++){
points->InsertNextPoint(x(i), y(i), z(i)/3.28084);
if (i < 10){
std::cout << "x: " << x(i) << std::endl;
std::cout << "y: " << y(i) << std::endl;
std::cout << "z: " << z(i) << std::endl;
}
}
vtkNew<vtkDoubleArray> scalars;
scalars->SetNumberOfValues(z.size());
scalars->SetName("scalars");
for (size_t i = 0; i < scalars->GetNumberOfValues(); i++){
scalars->SetValue(i, z(i));
}
std::cout << "Scalars set" << std::endl;
// Add the grid points to a polydata object.
vtkNew<vtkPolyData> polydata;
polydata->SetPoints(points);
polydata->GetPointData()->SetScalars(scalars);
std::cout << "Polydata prepared" << std::endl;
// Triangulate the grid points
vtkNew<vtkDelaunay2D> delaunay;
// delaunay->SetAlpha(2.0);
delaunay->SetInputData(polydata);
delaunay->Update();
std::cout << "Delauney done" << std::endl;
// Visualize
vtkNew<vtkPolyDataMapper> meshMapper;
meshMapper->SetInputConnection(delaunay->GetOutputPort());
vtkNew<vtkActor> meshActor;
meshActor->SetMapper(meshMapper);
meshActor->GetProperty()->SetColor(colors->GetColor3d("Banana").GetData());
meshActor->GetProperty()->EdgeVisibilityOn();
vtkNew<vtkRenderer> renderer;
renderer->SetBackground(colors->GetColor3d("Mint").GetData());
vtkNew<vtkRenderWindow> renderWindow;
renderWindow->AddRenderer(renderer);
vtkNew<vtkRenderWindowInteractor> renderWindowInteractor;
renderWindowInteractor->SetRenderWindow(renderWindow);
renderer->AddActor(meshActor);
renderWindow->SetWindowName("Delaunay2D");
renderWindow->Render();
renderWindowInteractor->Start();
return EXIT_SUCCESS;
}
VTK: 9.1.20220125
Windows 11
MSVC 2019
