[Small Feature] Add VTK read/write support for Linear_cell_complex

This commit is contained in:
Yliess Bellargui 2025-07-28 07:37:54 +02:00
parent 09365799e9
commit c37745641a
8 changed files with 866 additions and 649 deletions

File diff suppressed because one or more lines are too long

View File

@ -1,29 +0,0 @@
# vtk DataFile Version 2.0
Example Tetrahedron with Scalars
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 4 float
0.0 0.0 0.0
1.0 0.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0
CELLS 1 5
4 0 1 2 3
CELL_TYPES 1
10
POINT_DATA 4
SCALARS vertex_scalars float 1
LOOKUP_TABLE default
1.0
2.0
3.0
4.0
CELL_DATA 1
SCALARS volume_scalars float 1
LOOKUP_TABLE default
42.0

View File

@ -1,52 +1,29 @@
/*!
\ingroup PkgLinearCellComplexExamples
\brief Minimal example: read a `.3map` file, compute per-volume vertex count, and write to VTK.
\brief Example showing how to read and write a 3D linear cell complex from/to a VTK file.
This example loads a volumetric mesh from a `.vtk` file into a `Linear_cell_complex_for_combinatorial_map<3>`,
optionally reads scalar fields (vertex/volume), and writes the modified structure back to a `.vtk` file.
\cgalFeature{Linear_cell_complex_vtk_io}
\cgalExampleOutput{
Loaded LCC from data/lcc_input.vtk
- 80 vertices
- 52 volumes
Wrote LCC to data/lcc_output.vtk
}
This example loads a 3D linear cell complex from a `.3map` file using
\cgal function `CGAL::load_combinatorial_map()`. It computes for each
3-cell (volume) the number of incident vertices (0-cells), stores these values
in a `std::vector<std::size_t>`, and writes the result to a `.vtk` file with
`CGAL::write_lcc_to_vtk()`, using the computed values as scalars for each volume.
*/
#include <CGAL/Linear_cell_complex_for_combinatorial_map.h>
#include <CGAL/Linear_cell_complex_traits.h>
#include <CGAL/Linear_cell_complex_vtk_io.h>
#include <CGAL/IO/io.h>
#include <iostream>
#include <CGAL/Combinatorial_map_save_load.h>
#include <vector>
#include <fstream>
#include <cstdlib>
typedef CGAL::Linear_cell_complex_for_combinatorial_map<3> LCC;
int main()
{
using LCC = CGAL::Linear_cell_complex_for_combinatorial_map<3, 3, CGAL::Linear_cell_complex_traits<3>>;
LCC lcc;
std::vector<float> vertex_scalars, volume_scalars;
const char* input_filename = "data/lcc_input.vtk";
if (!CGAL::read_vtk(lcc, input_filename, &vertex_scalars, &volume_scalars)) {
std::cerr << "Failed to read: " << input_filename << std::endl;
return EXIT_FAILURE;
}
std::cout << "Loaded LCC from " << input_filename << std::endl;
std::cout << " - " << lcc.number_of_vertex_attributes() << " vertices" << std::endl;
std::cout << " - " << lcc.template one_dart_per_cell<3>().size() << " volumes" << std::endl;
const char* output_filename = "data/lcc_output.vtk";
if (!CGAL::write_vtk(lcc, output_filename, &vertex_scalars, &volume_scalars)) {
std::cerr << "Failed to write: " << output_filename << std::endl;
return EXIT_FAILURE;
}
std::cout << "Wrote LCC to " << output_filename << std::endl;
return EXIT_SUCCESS;
LCC lcc; CGAL::load_combinatorial_map("data/beam-with-mixed-cells.3map", lcc);
std::vector<std::size_t> v;
for(auto it = lcc.template one_dart_per_cell<3>().begin(), itend = lcc.template one_dart_per_cell<3>().end(); it != itend; ++it)
v.push_back(std::distance(lcc.template one_dart_per_incident_cell<0,3>(lcc.dart_descriptor(*it)).begin(), lcc.template one_dart_per_incident_cell<0,3>(lcc.dart_descriptor(*it)).end()));
CGAL::write_lcc_to_vtk(lcc, "data/beam-with-mixed-cells.vtk", nullptr, &v);
return EXIT_SUCCESS;
}

View File

@ -261,8 +261,7 @@ public:
prev_dart =lcc.null_descriptor;
}
void add_vertex_to_facet(size_type i)
{
void add_vertex_to_facet(size_type i, std::vector<DH>* tabdarts = nullptr) {
CGAL_assertion(i<vertex_map.size());
// std::cout<<i<<" "<<std::flush;
DH cur_dart=Add_vertex_to_face<LCC>::run(lcc, vertex_map[i], prev_dart);
@ -289,6 +288,7 @@ public:
{ first_dart=cur_dart; min_vertex=max_vertex=i; min_dart=cur_dart; }
prev_dart=cur_dart;
if(tabdarts != nullptr) { tabdarts->push_back(cur_dart); }
}
// End of the facet. Return the first dart of this facet.
@ -325,13 +325,14 @@ public:
return first_dart;
}
DH add_facet(std::initializer_list<size_type> l)
{
begin_facet();
for (size_type i:l)
{ add_vertex_to_facet(i); }
return end_facet();
}
DH add_facet(std::initializer_list<size_type> l, std::vector<DH>* tabdarts = nullptr)
{
if(tabdarts != nullptr) { tabdarts->reserve(tabdarts->size() + l.size()); }
begin_facet();
for (size_type i:l)
{ add_vertex_to_facet(i, tabdarts); }
return end_facet();
}
void begin_surface()
{
@ -404,5 +405,197 @@ private:
} //namespace CGAL
///////////////////////////////////////////////////////////////////////////////
/* Create an hexahedron, given the indices of its vertices (in the following
* order), the vertex must already have been added in the incremental builder.
* 3
* /|\
* 0-|-2
* \|/
* 1
*/
template<typename IncrementalBuilder>
typename IncrementalBuilder::LCC::Dart_descriptor
make_tetrahedron_with_builder(IncrementalBuilder& ib,
std::size_t i0,
std::size_t i1,
std::size_t i2,
std::size_t i3,
std::vector<typename IncrementalBuilder::LCC::Dart_descriptor>*
tabdarts=nullptr)
{
ib.begin_surface();
ib.add_facet({i0,i1,i2}, tabdarts);
ib.add_facet({i1,i0,i3}, tabdarts);
ib.add_facet({i2,i1,i3}, tabdarts);
ib.add_facet({i0,i2,i3}, tabdarts);
return ib.end_surface();
}
///////////////////////////////////////////////////////////////////////////////
/* 4
* /|\
* 0-|-3
* | | |
* 1---2
*/
template<typename IncrementalBuilder>
typename IncrementalBuilder::LCC::Dart_descriptor
make_pyramid_with_builder(IncrementalBuilder& ib,
std::size_t i0,
std::size_t i1,
std::size_t i2,
std::size_t i3,
std::size_t i4,
std::vector<typename IncrementalBuilder::LCC::Dart_descriptor>*
tabdarts=nullptr)
{
ib.begin_surface();
ib.add_facet({i0,i1,i2,i3}, tabdarts);
ib.add_facet({i1,i0,i4}, tabdarts);
ib.add_facet({i2,i1,i4}, tabdarts);
ib.add_facet({i3,i2,i4}, tabdarts);
ib.add_facet({i0,i3,i4}, tabdarts);
return ib.end_surface();
}
///////////////////////////////////////////////////////////////////////////////
/* 3
* /|\
* 4---5
* | | |
* | 0 |
* |/ \|
* 1---2
*/
template<typename IncrementalBuilder>
typename IncrementalBuilder::LCC::Dart_descriptor
make_prism_with_builder(IncrementalBuilder& ib,
std::size_t i0,
std::size_t i1,
std::size_t i2,
std::size_t i3,
std::size_t i4,
std::size_t i5,
std::vector<typename IncrementalBuilder::LCC::Dart_descriptor>*
tabdarts=nullptr)
{
ib.begin_surface();
ib.add_facet({i0,i1,i2}, tabdarts);
ib.add_facet({i1,i0,i3,i4}, tabdarts);
ib.add_facet({i2,i1,i4,i5}, tabdarts);
ib.add_facet({i0,i2,i5,i3}, tabdarts);
ib.add_facet({i5,i4,i3}, tabdarts);
return ib.end_surface();
}
///////////////////////////////////////////////////////////////////////////////
/* 7----6
* /| /|
* 4----5 |
* | 3--|-2
* |/ |/
* 0----1
*/
template<typename IncrementalBuilder>
typename IncrementalBuilder::LCC::Dart_descriptor
make_hexahedron_with_builder(IncrementalBuilder& ib,
std::size_t i0,
std::size_t i1,
std::size_t i2,
std::size_t i3,
std::size_t i4,
std::size_t i5,
std::size_t i6,
std::size_t i7,
std::vector<typename IncrementalBuilder::LCC::Dart_descriptor>*
tabdarts=nullptr)
{
ib.begin_surface();
ib.add_facet({i0,i1,i2,i3}, tabdarts);
ib.add_facet({i1,i0,i4,i5}, tabdarts);
ib.add_facet({i2,i1,i5,i6}, tabdarts);
ib.add_facet({i3,i2,i6,i7}, tabdarts);
ib.add_facet({i0,i3,i7,i4}, tabdarts);
ib.add_facet({i7,i6,i5,i4}, tabdarts);
return ib.end_surface();
}
///////////////////////////////////////////////////////////////////////////////
template<typename IncrementalBuilder>
typename IncrementalBuilder::LCC::Dart_descriptor
make_pentagonal_prism_with_builder(IncrementalBuilder& ib,
std::size_t i0,
std::size_t i1,
std::size_t i2,
std::size_t i3,
std::size_t i4,
std::size_t i5,
std::size_t i6,
std::size_t i7,
std::size_t i8,
std::size_t i9,
std::vector<typename IncrementalBuilder::LCC::Dart_descriptor>*
tabdarts=nullptr)
{
ib.begin_surface();
ib.add_facet({i0,i1,i2,i3,i4}, tabdarts);
ib.add_facet({i1,i0,i5,i6}, tabdarts);
ib.add_facet({i2,i1,i6,i7}, tabdarts);
ib.add_facet({i3,i2,i7,i8}, tabdarts);
ib.add_facet({i4,i3,i8,i9}, tabdarts);
ib.add_facet({i0,i4,i9,i5}, tabdarts);
ib.add_facet({i9,i8,i7,i6,i5}, tabdarts);
return ib.end_surface();
}
///////////////////////////////////////////////////////////////////////////////
template<typename IncrementalBuilder>
typename IncrementalBuilder::LCC::Dart_descriptor
make_hexagonal_prism_with_builder(IncrementalBuilder& ib,
std::size_t i0,
std::size_t i1,
std::size_t i2,
std::size_t i3,
std::size_t i4,
std::size_t i5,
std::size_t i6,
std::size_t i7,
std::size_t i8,
std::size_t i9,
std::size_t i10,
std::size_t i11,
std::vector<typename IncrementalBuilder::LCC::Dart_descriptor>*
tabdarts=nullptr)
{
ib.begin_surface();
ib.add_facet({i0,i1,i2,i3,i4,i5}, tabdarts);
ib.add_facet({i1,i0,i6,i7}, tabdarts);
ib.add_facet({i2,i1,i7,i8}, tabdarts);
ib.add_facet({i3,i2,i8,i9}, tabdarts);
ib.add_facet({i4,i3,i9,i10}, tabdarts);
ib.add_facet({i5,i4,i10,i11}, tabdarts);
ib.add_facet({i0,i5,i11,i6}, tabdarts);
ib.add_facet({i11,i10,i9,i8,i7,i6}, tabdarts);
return ib.end_surface();
}
///////////////////////////////////////////////////////////////////////////////
template<typename IncrementalBuilder>
typename IncrementalBuilder::LCC::Dart_descriptor
make_generic_cell_with_builder(IncrementalBuilder& ib,
const std::vector<std::size_t>& faces,
std::vector<typename IncrementalBuilder::LCC::Dart_descriptor>*
tabdarts=nullptr)
{
ib.begin_surface();
std::size_t i=1, end; // Start to 1 because faces[0] is the number of faces
for(; i<faces.size(); )
{
end=i+1+faces[i]; // faces[i] is the number of vertices of the face; +i is the index of the end
++i; // I prefer to increment i after its use!
ib.begin_facet();
for(; i<end; ++i)
{ ib.add_vertex_to_facet(faces[i], tabdarts); }
ib.end_facet();
}
return ib.end_surface();
}
///////////////////////////////////////////////////////////////////////////////
#endif // CGAL_LINEAR_CELL_COMPLEX_INCREMENTAL_BUILDER_3_H //
// EOF //

View File

@ -12,18 +12,25 @@
#ifndef CGAL_LINEAR_CELL_COMPLEX_VTK_IO_H
#define CGAL_LINEAR_CELL_COMPLEX_VTK_IO_H 1
#include "Linear_cell_complex_incremental_builder_3.h"
#include <CGAL/assertions.h>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include <unordered_map>
#include <vector>
namespace CGAL {
/** \file Linear_cell_complex_vtk_io.h
* Functions to import/export 3D Linear_cell_complex from/to VTK legacy ASCII format.
*
*
* Only supports:
* - Linear_cell_complex_for_combinatorial_map<3,3>
* - Linear_cell_complex_for_combinatorial_map<3,3>
* - VTK legacy ASCII format (.vtk files)
* - Optional scalar fields for vertices and volumes
*
*
* Supported VTK cell types:
* - VTK_TETRA (10): Tetrahedron
* - VTK_VOXEL (11): Voxel (special hexahedron ordering)
@ -35,6 +42,10 @@ namespace CGAL {
* - VTK_POLYHEDRON (42): Generic polyhedron
*/
// ============================================================================
// Declarations
// ============================================================================
/**
* \brief Read a VTK legacy ASCII file and load it into a 3D Linear_cell_complex.
* \ingroup PkgLinearCellComplexExamples
@ -48,11 +59,11 @@ namespace CGAL {
* If provided, will be resized to match number of volumes.
* \return `true` if loading was successful, `false` otherwise
*/
template <typename LCC>
bool read_vtk(LCC& alcc,
template <typename LCC, typename ScalarType = float>
bool read_lcc_from_vtk(LCC& alcc,
const char* filename,
std::vector<float>* vertex_scalars = nullptr,
std::vector<float>* volume_scalars = nullptr);
std::vector<ScalarType>* vertex_scalars = nullptr,
std::vector<ScalarType>* volume_scalars = nullptr);
/**
* \brief Write a 3D Linear_cell_complex to a VTK legacy ASCII file.
@ -67,14 +78,593 @@ bool read_vtk(LCC& alcc,
* same size as number of 3-cells in the LCC.
* \return `true` if writing was successful, `false` otherwise
*/
template <typename LCC>
bool write_vtk(const LCC& alcc,
template <typename LCC, typename ScalarType = float>
bool write_lcc_to_vtk(const LCC& alcc,
const char* filename,
const std::vector<float>* vertex_scalars = nullptr,
const std::vector<float>* volume_scalars = nullptr);
const std::vector<ScalarType>* vertex_scalars = nullptr,
const std::vector<ScalarType>* volume_scalars = nullptr);
// Advanced versions with functors
template <typename LCC, typename VertexScalarReader, typename CellScalarReader>
bool read_lcc_from_vtk(LCC& alcc, const char* filename, VertexScalarReader vertex_reader, CellScalarReader cell_reader);
template <typename LCC, typename PointFunctor, typename CellFunctor>
bool write_lcc_to_vtk(const LCC& alcc, const char* filename, PointFunctor ptval, CellFunctor cellval);
// ============================================================================
// Implementation details
// ============================================================================
namespace internal {
// Helper: read a scalar of the right type from stream
template <typename Functor>
bool read_scalar_by_vtk_type(std::istream& is, const std::string& vtk_type, std::size_t n, Functor f) {
if(vtk_type == "float") {
for(std::size_t i = 0; i < n; ++i) {
float v;
if(!(is >> v))
return false;
f(i, v);
}
} else if(vtk_type == "double") {
for(std::size_t i = 0; i < n; ++i) {
double v;
if(!(is >> v))
return false;
f(i, v);
}
} else if(vtk_type == "int") {
for(std::size_t i = 0; i < n; ++i) {
int v;
if(!(is >> v))
return false;
f(i, v);
}
} else if(vtk_type == "unsigned_int") {
for(std::size_t i = 0; i < n; ++i) {
unsigned int v;
if(!(is >> v))
return false;
f(i, v);
}
} else if(vtk_type == "short") {
for(std::size_t i = 0; i < n; ++i) {
short int v;
if(!(is >> v))
return false;
f(i, v);
}
} else if(vtk_type == "unsigned_short") {
for(std::size_t i = 0; i < n; ++i) {
unsigned short int v;
if(!(is >> v))
return false;
f(i, v);
}
} else if(vtk_type == "char") {
for(std::size_t i = 0; i < n; ++i) {
char v;
int tmp;
if(!(is >> tmp))
return false;
v = static_cast<char>(tmp);
f(i, v);
}
} else if(vtk_type == "unsigned_char") {
for(std::size_t i = 0; i < n; ++i) {
unsigned char v;
int tmp;
if(!(is >> tmp))
return false;
v = static_cast<unsigned char>(tmp);
f(i, v);
}
} else if(vtk_type == "long") {
for(std::size_t i = 0; i < n; ++i) {
long int v;
if(!(is >> v))
return false;
f(i, v);
}
} else if(vtk_type == "unsigned_long") {
for(std::size_t i = 0; i < n; ++i) {
unsigned long int v;
if(!(is >> v))
return false;
f(i, v);
}
} else {
std::cerr << "[ERROR] read_lcc_from_vtk: unsupported scalar type: " << vtk_type << std::endl;
return false;
}
return true;
}
// VTK type name mapping
template <typename T> struct gettype
{
static std::string name() { return "unknown"; }
};
template <> struct gettype<bool>
{
static std::string name() { return "bit"; }
};
template <> struct gettype<unsigned char>
{
static std::string name() { return "unsigned_char"; }
};
template <> struct gettype<char>
{
static std::string name() { return "char"; }
};
template <> struct gettype<unsigned short int>
{
static std::string name() { return "unsigned_short"; }
};
template <> struct gettype<short int>
{
static std::string name() { return "short"; }
};
template <> struct gettype<unsigned int>
{
static std::string name() { return "unsigned_int"; }
};
template <> struct gettype<int>
{
static std::string name() { return "int"; }
};
template <> struct gettype<unsigned long int>
{
static std::string name() { return "unsigned_long"; }
};
template <> struct gettype<long int>
{
static std::string name() { return "long"; }
};
template <> struct gettype<float>
{
static std::string name() { return "float"; }
};
template <> struct gettype<double>
{
static std::string name() { return "double"; }
};
// VTK cell type constants
enum VTK_Cell_Type {
VTK_TETRA = 10,
VTK_VOXEL = 11,
VTK_HEXAHEDRON = 12,
VTK_WEDGE = 13,
VTK_PYRAMID = 14,
VTK_PENTAGONAL_PRISM = 15,
VTK_HEXAGONAL_PRISM = 16,
VTK_POLYHEDRON = 42
};
// Helper: detect VTK cell type from a 3-cell
template <typename LCC, typename Dart> inline VTK_Cell_Type get_vtk_cell_type(const LCC& lcc, Dart itvol) {
// Heuristic: count number of vertices and faces
std::size_t nbv = 0, nbf = 0;
for(auto itv = lcc.template one_dart_per_incident_cell<0, 3>(itvol).begin(),
itvend = lcc.template one_dart_per_incident_cell<0, 3>(itvol).end();
itv != itvend; ++itv)
++nbv;
for(auto itf = lcc.template one_dart_per_incident_cell<2, 3>(itvol).begin(),
itfend = lcc.template one_dart_per_incident_cell<2, 3>(itvol).end();
itf != itfend; ++itf)
++nbf;
if(nbv == 4 && nbf == 4)
return VTK_TETRA;
if(nbv == 5 && nbf == 5)
return VTK_PYRAMID;
if(nbv == 6 && nbf == 5)
return VTK_WEDGE;
if(nbv == 8 && nbf == 6)
return VTK_HEXAHEDRON;
if(nbv == 10 && nbf == 7)
return VTK_PENTAGONAL_PRISM;
if(nbv == 12 && nbf == 8)
return VTK_HEXAGONAL_PRISM;
return VTK_POLYHEDRON;
}
template <typename LCC, typename VertexScalarReader, typename CellScalarReader>
bool read_lcc_from_vtk_ascii(std::istream& is, LCC& alcc, VertexScalarReader vertex_reader, CellScalarReader cell_reader) {
static_assert(LCC::dimension == 3 && LCC::ambient_dimension == 3,
"read_lcc_from_vtk() only supports 3D Linear_cell_complexes (3,3)");
using Point = typename LCC::Point;
using FT = typename LCC::FT;
alcc.clear();
Linear_cell_complex_incremental_builder_3<LCC> ib(alcc);
std::string line, tmp;
std::size_t npoints, ncells;
// Skip to POINTS section
while(std::getline(is, line)) {
if(line.find("POINTS") != std::string::npos)
break;
}
if(is.eof()) {
std::cerr << "[ERROR] read_lcc_from_vtk: POINTS section not found" << std::endl;
return false;
}
std::stringstream ss(line);
std::getline(ss, tmp, ' '); // skip "POINTS"
ss >> npoints;
// Read points
std::vector<typename LCC::Vertex_attribute_descriptor> points(npoints);
for(std::size_t i = 0; i < npoints; ++i) {
FT x, y, z;
if(!(is >> x >> y >> z)) {
std::cerr << "[ERROR] read_lcc_from_vtk: failed to read point " << i << std::endl;
return false;
}
points[i] = ib.add_vertex(Point(x, y, z));
}
// Skip to CELLS section
while(std::getline(is, line)) {
if(line.find("CELLS") != std::string::npos)
break;
}
if(is.eof()) {
std::cerr << "[ERROR] read_lcc_from_vtk: CELLS section not found" << std::endl;
return false;
}
ss = std::stringstream(line);
std::getline(ss, tmp, ' '); // skip "CELLS"
ss >> ncells;
// Read connectivity
std::vector<std::vector<std::size_t>> faces(ncells);
std::size_t points_per_cell;
for(std::size_t i = 0; i < ncells; ++i) {
if(!(is >> points_per_cell)) {
std::cerr << "[ERROR] read_lcc_from_vtk: failed to read cell " << i << std::endl;
return false;
}
faces[i].resize(points_per_cell);
for(std::size_t j = 0; j < points_per_cell; ++j) {
if(!(is >> faces[i][j])) {
std::cerr << "[ERROR] read_lcc_from_vtk: failed to read cell " << i << " vertex " << j << std::endl;
return false;
}
}
}
// Skip to CELL_TYPES section
while(std::getline(is, line)) {
if(line.find("CELL_TYPES") != std::string::npos)
break;
}
if(is.eof()) {
std::cerr << "[ERROR] read_lcc_from_vtk: CELL_TYPES section not found" << std::endl;
return false;
}
// Create cells based on types
std::size_t cell_type;
for(std::size_t i = 0; i < ncells; ++i) {
if(!(is >> cell_type)) {
std::cerr << "[ERROR] read_lcc_from_vtk: failed to read cell type " << i << std::endl;
return false;
}
const auto& v = faces[i];
switch(cell_type) {
case 10: // TETRA
if(v.size() == 4)
make_tetrahedron_with_builder(ib, v[0], v[1], v[2], v[3]);
break;
case 11: // VOXEL
if(v.size() == 8)
make_hexahedron_with_builder(ib, v[0], v[1], v[3], v[2], v[4], v[5], v[7], v[6]);
break;
case 12: // HEXAHEDRON
if(v.size() == 8)
make_hexahedron_with_builder(ib, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7]);
break;
case 13: // PRISM (WEDGE)
if(v.size() == 6)
make_prism_with_builder(ib, v[0], v[1], v[2], v[3], v[4], v[5]);
break;
case 14: // PYRAMID
if(v.size() == 5)
make_pyramid_with_builder(ib, v[0], v[1], v[2], v[3], v[4]);
break;
case 15: // PENTAGONAL_PRISM
if(v.size() == 10)
make_pentagonal_prism_with_builder(ib, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8], v[9]);
break;
case 16: // HEXAGONAL_PRISM
if(v.size() == 12)
make_hexagonal_prism_with_builder(ib, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8], v[9], v[10], v[11]);
break;
case 42: // GENERIC CELL
make_generic_cell_with_builder(ib, v);
break;
default:
std::cerr << "[ERROR] read_lcc_from_vtk: type " << cell_type << " unknown." << std::endl;
}
}
// Clean up unused vertex attributes
for(auto itv = alcc.vertex_attributes().begin(); itv != alcc.vertex_attributes().end();) {
if(alcc.template dart_of_attribute<0>(itv) == alcc.null_descriptor) {
alcc.erase_vertex_attribute(itv);
} else {
++itv;
}
}
// Read POINT_DATA scalars if present
while(std::getline(is, line)) {
if(line.find("POINT_DATA") != std::string::npos) {
std::size_t ndata;
ss = std::stringstream(line);
std::getline(ss, tmp, ' '); // skip "POINT_DATA"
ss >> ndata;
std::getline(is, line); // SCALARS line
std::stringstream scalars_line(line);
std::string scalars, name, vtk_type;
scalars_line >> scalars >> name >> vtk_type;
std::getline(is, line); // LOOKUP_TABLE line
if(!read_scalar_by_vtk_type(is, vtk_type, ndata, vertex_reader)) {
std::cerr << "[ERROR] read_lcc_from_vtk: failed to read POINT_DATA" << std::endl;
}
break;
}
}
// Read CELL_DATA scalars if present
while(std::getline(is, line)) {
if(line.find("CELL_DATA") != std::string::npos) {
std::size_t ndata;
ss = std::stringstream(line);
std::getline(ss, tmp, ' '); // skip "CELL_DATA"
ss >> ndata;
std::getline(is, line); // SCALARS line
std::stringstream scalars_line(line);
std::string scalars, name, vtk_type;
scalars_line >> scalars >> name >> vtk_type;
std::getline(is, line); // LOOKUP_TABLE line
if(!read_scalar_by_vtk_type(is, vtk_type, ndata, cell_reader)) {
std::cerr << "[ERROR] read_lcc_from_vtk: failed to read CELL_DATA" << std::endl;
}
break;
}
}
return true;
}
template <typename LCC, typename PointFunctor, typename CellFunctor>
bool write_lcc_to_vtk_ascii(std::ostream& os, const LCC& alcc, PointFunctor ptval, CellFunctor cellval) {
static_assert(LCC::dimension == 3 && LCC::ambient_dimension == 3,
"write_lcc_to_vtk() only supports 3D Linear_cell_complexes (3,3)");
// Write VTK header
os << "# vtk DataFile Version 2.0\n";
os << "CGAL Linear_cell_complex\n";
os << "ASCII\n";
os << "DATASET UNSTRUCTURED_GRID\n\n";
// Build vertex index map
std::unordered_map<typename LCC::Vertex_attribute_const_descriptor, std::size_t> vertex_index;
std::size_t nbpts = 0;
for(auto itv = alcc.vertex_attributes().begin(), itvend = alcc.vertex_attributes().end(); itv != itvend; ++itv) {
vertex_index[itv] = nbpts++;
}
// Write points
os << "POINTS " << nbpts << " double\n";
for(auto itv = alcc.vertex_attributes().begin(), itvend = alcc.vertex_attributes().end(); itv != itvend; ++itv) {
const auto& p = itv->point();
os << p.x() << " " << p.y() << " " << p.z() << "\n";
}
os << "\n";
// Count cells and build connectivity
std::size_t nbcells = 0;
std::size_t total_size = 0;
std::ostringstream cell_stream;
std::ostringstream type_stream;
for(auto itvol = alcc.template one_dart_per_cell<3>().begin(), itvolend = alcc.template one_dart_per_cell<3>().end();
itvol != itvolend; ++itvol)
{
++nbcells;
VTK_Cell_Type cell_type = get_vtk_cell_type(alcc, itvol);
type_stream << static_cast<int>(cell_type) << "\n";
if(cell_type == VTK_POLYHEDRON) {
// Generic polyhedron format write as face-vertex connectivity
std::vector<std::vector<std::size_t>> faces;
std::size_t cell_size = 1; // Start with 1 for number of faces
for(auto itface = alcc.template one_dart_per_incident_cell<2, 3>(itvol).begin(),
itfaceend = alcc.template one_dart_per_incident_cell<2, 3>(itvol).end();
itface != itfaceend; ++itface)
{
faces.push_back(std::vector<std::size_t>());
auto& face = faces.back();
for(auto itvert = alcc.template darts_of_orbit<1>(itface).begin(),
itvertend = alcc.template darts_of_orbit<1>(itface).end();
itvert != itvertend; ++itvert)
{
face.push_back(vertex_index[alcc.vertex_attribute(itvert)]);
}
cell_size += face.size() + 1; // +1 for face size
}
cell_stream << cell_size << " " << faces.size();
for(const auto& face : faces) {
cell_stream << " " << face.size();
for(auto v : face) {
cell_stream << " " << v;
}
}
cell_stream << "\n";
total_size += cell_size + 1; // +1 for cell size
} else {
// Standard cell types write vertex connectivity directly
std::vector<std::size_t> vertices;
for(auto itvert = alcc.template one_dart_per_incident_cell<0, 3>(itvol).begin(),
itvertend = alcc.template one_dart_per_incident_cell<0, 3>(itvol).end();
itvert != itvertend; ++itvert)
{
vertices.push_back(vertex_index[alcc.vertex_attribute(itvert)]);
}
cell_stream << vertices.size();
for(auto v : vertices) {
cell_stream << " " << v;
}
cell_stream << "\n";
total_size += vertices.size() + 1;
}
}
// Write cells section
os << "CELLS " << nbcells << " " << total_size << "\n";
os << cell_stream.str();
os << "\n";
// Write cell types
os << "CELL_TYPES " << nbcells << "\n";
os << type_stream.str();
os << "\n";
// Write vertex scalars if ptval is not nullptr
if constexpr(!std::is_same_v<PointFunctor, std::nullptr_t>) {
os << "POINT_DATA " << nbpts << "\n";
os << "SCALARS vertex_scalars " << gettype<decltype(ptval(alcc, alcc.null_dart_descriptor))>::name() << " 1\n";
os << "LOOKUP_TABLE default\n";
for(auto itv = alcc.vertex_attributes().begin(), itvend = alcc.vertex_attributes().end(); itv != itvend; ++itv) {
auto dart = alcc.template dart_of_attribute<0>(itv);
CGAL_assertion(dart != alcc.null_dart_descriptor);
os << ptval(alcc, dart) << "\n";
}
}
// Write cell scalars if cellval is not nullptr/pointer
if constexpr(!std::is_same_v<CellFunctor, std::nullptr_t> && !std::is_pointer_v<CellFunctor>) {
os << "CELL_DATA " << nbcells << "\n";
os << "SCALARS volume_scalars " << gettype<decltype(cellval(alcc, alcc.null_dart_descriptor))>::name() << " 1\n";
os << "LOOKUP_TABLE default\n";
for(auto itvol = alcc.template one_dart_per_cell<3>().begin(),
itvolend = alcc.template one_dart_per_cell<3>().end();
itvol != itvolend; ++itvol)
{
os << cellval(alcc, itvol) << "\n";
}
}
return true;
}
} // namespace internal
// ============================================================================
// Public interface implementation
// ============================================================================
// Functor-based versions
template <typename LCC, typename VertexScalarReader, typename CellScalarReader>
inline bool read_lcc_from_vtk(LCC& alcc, const char* filename, VertexScalarReader vertex_reader, CellScalarReader cell_reader) {
CGAL_assertion(filename != nullptr);
std::ifstream file(filename);
if(!file.is_open()) {
std::cerr << "[ERROR] read_lcc_from_vtk: cannot open file " << filename << std::endl;
return false;
}
return internal::read_lcc_from_vtk_ascii(file, alcc, vertex_reader, cell_reader);
}
template <typename LCC, typename PointFunctor, typename CellFunctor>
inline bool write_lcc_to_vtk(const LCC& alcc, const char* filename, PointFunctor ptval, CellFunctor cellval) {
CGAL_assertion(filename != nullptr);
std::ofstream file(filename);
if(!file.is_open()) {
std::cerr << "[ERROR] write_lcc_to_vtk: cannot open file " << filename << std::endl;
return false;
}
return internal::write_lcc_to_vtk_ascii(file, alcc, ptval, cellval);
}
// Vector-based versions (convenience wrappers)
template <typename LCC, typename ScalarType>
inline bool read_lcc_from_vtk(LCC& alcc,
const char* filename,
std::vector<ScalarType>* vertex_scalars,
std::vector<ScalarType>* volume_scalars) {
auto v_writer = [&](std::size_t i, ScalarType val) {
if(vertex_scalars) {
if(vertex_scalars->size() <= i)
vertex_scalars->resize(i + 1);
(*vertex_scalars)[i] = val;
}
};
auto c_writer = [&](std::size_t i, ScalarType val) {
if(volume_scalars) {
if(volume_scalars->size() <= i)
volume_scalars->resize(i + 1);
(*volume_scalars)[i] = val;
}
};
return read_lcc_from_vtk(alcc, filename, v_writer, c_writer);
}
template <typename LCC, typename ScalarType>
inline bool write_lcc_to_vtk(const LCC& alcc,
const char* filename,
const std::vector<ScalarType>* vertex_scalars,
const std::vector<ScalarType>* volume_scalars) {
// Build index maps
std::unordered_map<typename LCC::Vertex_attribute_const_descriptor, std::size_t> vertex_indices;
std::size_t idx = 0;
for(auto itv = alcc.vertex_attributes().begin(), itvend = alcc.vertex_attributes().end(); itv != itvend; ++itv)
vertex_indices[itv] = idx++;
std::unordered_map<typename LCC::Dart_const_descriptor, std::size_t> volume_indices;
idx = 0;
for(auto itvol = alcc.template one_dart_per_cell<3>().begin(), itvolend = alcc.template one_dart_per_cell<3>().end();
itvol != itvolend; ++itvol)
volume_indices[itvol] = idx++;
return write_lcc_to_vtk(
alcc, filename,
[vertex_scalars, &vertex_indices](const LCC& lcc, typename LCC::Dart_const_descriptor d) -> ScalarType {
if(vertex_scalars)
return (*vertex_scalars)[vertex_indices.at(lcc.attribute<0>(d))];
return ScalarType();
},
[volume_scalars, &volume_indices](const LCC& lcc, typename LCC::Dart_const_descriptor d) -> ScalarType {
if(volume_scalars)
return (*volume_scalars)[volume_indices.at(d)];
return ScalarType();
});
}
} // namespace CGAL
#include <CGAL/Linear_cell_complex_vtk_io_impl.h>
#endif // CGAL_LINEAR_CELL_COMPLEX_VTK_IO_H

View File

@ -1,524 +0,0 @@
// Copyright (c) 2025 CNRS and LIRIS' Establishments (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Guillaume Damiand <guillaume.damiand@liris.cnrs.fr>
#ifndef CGAL_LINEAR_CELL_COMPLEX_VTK_IO_IMPL_H
#define CGAL_LINEAR_CELL_COMPLEX_VTK_IO_IMPL_H 1
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include <unordered_map>
#include <vector>
#include <CGAL/assertions.h>
namespace CGAL {
namespace internal {
// VTK cell type constants
enum VTK_Cell_Type {
VTK_TETRA = 10,
VTK_VOXEL = 11,
VTK_HEXAHEDRON = 12,
VTK_WEDGE = 13,
VTK_PYRAMID = 14,
VTK_PENTAGONAL_PRISM = 15,
VTK_HEXAGONAL_PRISM = 16,
VTK_POLYHEDRON = 42
};
// Helper function to create a wedge/prism using the incremental builder
template<typename LCC>
void create_wedge(LCC& lcc,
typename LCC::Vertex_attribute_descriptor p0,
typename LCC::Vertex_attribute_descriptor p1,
typename LCC::Vertex_attribute_descriptor p2,
typename LCC::Vertex_attribute_descriptor p3,
typename LCC::Vertex_attribute_descriptor p4,
typename LCC::Vertex_attribute_descriptor p5)
{
// Create wedge using incremental builder
// A wedge has 2 triangular faces and 3 quadrilateral faces
typename LCC::Dart_descriptor d1, d2, d3, d4, d5;
// Create the triangular faces
d1 = lcc.make_triangle(p0, p1, p2);
d2 = lcc.make_triangle(p3, p5, p4); // reversed order for proper orientation
// Create quadrilateral faces
d3 = lcc.make_quadrangle(p0, p3, p4, p1);
d4 = lcc.make_quadrangle(p1, p4, p5, p2);
d5 = lcc.make_quadrangle(p2, p5, p3, p0);
// Now we need to sew these faces together
// This is a simplified approach in practice, you'd need to carefully
// manage the dart orientations and sewing operations
// For now, we'll leave them as separate faces
}
// Helper function to create a pyramid
template<typename LCC>
void create_pyramid(LCC& lcc,
typename LCC::Vertex_attribute_descriptor p0,
typename LCC::Vertex_attribute_descriptor p1,
typename LCC::Vertex_attribute_descriptor p2,
typename LCC::Vertex_attribute_descriptor p3,
typename LCC::Vertex_attribute_descriptor p4)
{
// Create pyramid using incremental builder
// A pyramid has 1 quadrilateral base and 4 triangular faces
// Create the quadrilateral base
typename LCC::Dart_descriptor base = lcc.make_quadrangle(p0, p1, p2, p3);
// Create triangular faces
lcc.make_triangle(p0, p4, p1);
lcc.make_triangle(p1, p4, p2);
lcc.make_triangle(p2, p4, p3);
lcc.make_triangle(p3, p4, p0);
// This creates separate faces. For a proper pyramid, you'd need
// to sew these faces together using the sew operations of the LCC
}
// Helper to determine cell topology
template<typename LCC>
VTK_Cell_Type get_vtk_cell_type(const LCC& lcc,
typename LCC::Dart_const_descriptor dart)
{
// Count vertices and faces to determine cell type
std::size_t num_vertices = 0;
std::size_t num_faces = 0;
for (auto it = lcc.template one_dart_per_incident_cell<0,3>(dart).begin(),
itend = lcc.template one_dart_per_incident_cell<0,3>(dart).end();
it != itend; ++it) {
++num_vertices;
}
for (auto it = lcc.template one_dart_per_incident_cell<2,3>(dart).begin(),
itend = lcc.template one_dart_per_incident_cell<2,3>(dart).end();
it != itend; ++it) {
++num_faces;
}
// Simple heuristic based on vertex/face count
if (num_vertices == 4 && num_faces == 4) return VTK_TETRA;
if (num_vertices == 8 && num_faces == 6) return VTK_HEXAHEDRON;
if (num_vertices == 5 && num_faces == 5) return VTK_PYRAMID;
if (num_vertices == 6 && num_faces == 5) return VTK_WEDGE;
if (num_vertices == 10 && num_faces == 7) return VTK_PENTAGONAL_PRISM;
if (num_vertices == 12 && num_faces == 8) return VTK_HEXAGONAL_PRISM;
return VTK_POLYHEDRON; // Default to generic polyhedron
}
template<typename LCC>
bool read_vtk_ascii(std::istream& is,
LCC& alcc,
std::vector<float>* vertex_scalars,
std::vector<float>* volume_scalars)
{
static_assert(LCC::dimension == 3 && LCC::ambient_dimension == 3,
"read_vtk() only supports 3D Linear_cell_complexes (3,3)");
using Point = typename LCC::Point;
using FT = typename LCC::FT;
alcc.clear();
std::string line, tmp;
std::size_t npoints, ncells;
// Skip to POINTS section
while (std::getline(is, line)) {
if (line.find("POINTS") != std::string::npos) break;
}
if (is.eof()) {
std::cerr << "[ERROR] read_vtk: POINTS section not found" << std::endl;
return false;
}
std::stringstream ss(line);
std::getline(ss, tmp, ' '); // skip "POINTS"
ss >> npoints;
// Read points
std::vector<typename LCC::Vertex_attribute_descriptor> points(npoints);
FT x, y, z;
for (std::size_t i = 0; i < npoints; ++i) {
if (!(is >> x >> y >> z)) {
std::cerr << "[ERROR] read_vtk: failed to read point " << i << std::endl;
return false;
}
points[i] = alcc.create_vertex_attribute(Point(x, y, z));
}
// Skip to CELLS section
while (std::getline(is, line)) {
if (line.find("CELLS") != std::string::npos) break;
}
if (is.eof()) {
std::cerr << "[ERROR] read_vtk: CELLS section not found" << std::endl;
return false;
}
ss = std::stringstream(line);
std::getline(ss, tmp, ' '); // skip "CELLS"
ss >> ncells;
// Read connectivity
std::vector<std::vector<std::size_t>> faces(ncells);
std::size_t points_per_cell;
for (std::size_t i = 0; i < ncells; ++i) {
if (!(is >> points_per_cell)) {
std::cerr << "[ERROR] read_vtk: failed to read cell " << i << std::endl;
return false;
}
faces[i].resize(points_per_cell);
for (std::size_t j = 0; j < points_per_cell; ++j) {
if (!(is >> faces[i][j])) {
std::cerr << "[ERROR] read_vtk: failed to read cell " << i
<< " vertex " << j << std::endl;
return false;
}
}
}
// Skip to CELL_TYPES section
while (std::getline(is, line)) {
if (line.find("CELL_TYPES") != std::string::npos) break;
}
if (is.eof()) {
std::cerr << "[ERROR] read_vtk: CELL_TYPES section not found" << std::endl;
return false;
}
// Create cells based on types
std::size_t cell_type;
for (std::size_t i = 0; i < ncells; ++i) {
if (!(is >> cell_type)) {
std::cerr << "[ERROR] read_vtk: failed to read cell type " << i << std::endl;
return false;
}
const auto& cell_vertices = faces[i];
switch (cell_type) {
case VTK_TETRA:
if (cell_vertices.size() == 4) {
alcc.make_tetrahedron(points[cell_vertices[0]], points[cell_vertices[1]],
points[cell_vertices[2]], points[cell_vertices[3]]);
}
break;
case VTK_VOXEL:
if (cell_vertices.size() == 8) {
// VTK voxel has different vertex ordering than standard hexahedron
alcc.make_hexahedron(points[cell_vertices[0]], points[cell_vertices[1]],
points[cell_vertices[3]], points[cell_vertices[2]],
points[cell_vertices[4]], points[cell_vertices[5]],
points[cell_vertices[7]], points[cell_vertices[6]]);
}
break;
case VTK_HEXAHEDRON:
if (cell_vertices.size() == 8) {
alcc.make_hexahedron(points[cell_vertices[0]], points[cell_vertices[1]],
points[cell_vertices[2]], points[cell_vertices[3]],
points[cell_vertices[4]], points[cell_vertices[5]],
points[cell_vertices[6]], points[cell_vertices[7]]);
}
break;
case VTK_WEDGE:
if (cell_vertices.size() == 6) {
// Use helper function to create wedge
create_wedge(alcc, points[cell_vertices[0]], points[cell_vertices[1]],
points[cell_vertices[2]], points[cell_vertices[3]],
points[cell_vertices[4]], points[cell_vertices[5]]);
}
break;
case VTK_PYRAMID:
if (cell_vertices.size() == 5) {
// Use helper function to create pyramid
create_pyramid(alcc, points[cell_vertices[0]], points[cell_vertices[1]],
points[cell_vertices[2]], points[cell_vertices[3]],
points[cell_vertices[4]]);
}
break;
case VTK_POLYHEDRON:
// For generic polyhedron, we'd need more complex construction
// This is a simplified version in practice, VTK polyhedron format
// includes face connectivity information that we're not parsing here
std::cerr << "[WARNING] read_vtk: VTK_POLYHEDRON not fully supported yet" << std::endl;
break;
default:
std::cerr << "[ERROR] read_vtk: unsupported cell type " << cell_type << std::endl;
break;
}
}
// Clean up unused vertex attributes
for (auto itv = alcc.vertex_attributes().begin();
itv != alcc.vertex_attributes().end(); ) {
if (alcc.template dart_of_attribute<0>(itv) == alcc.null_descriptor) {
alcc.erase_vertex_attribute(itv);
} else {
++itv;
}
}
// Try to read scalar data if requested
if (vertex_scalars != nullptr) {
vertex_scalars->clear();
// Look for POINT_DATA section
while (std::getline(is, line)) {
if (line.find("POINT_DATA") != std::string::npos) {
std::size_t ndata;
ss = std::stringstream(line);
std::getline(ss, tmp, ' '); // skip "POINT_DATA"
ss >> ndata;
// Skip SCALARS and LOOKUP_TABLE lines
std::getline(is, line); // SCALARS line
std::getline(is, line); // LOOKUP_TABLE line
vertex_scalars->resize(ndata);
for (std::size_t i = 0; i < ndata; ++i) {
if (!(is >> (*vertex_scalars)[i])) {
std::cerr << "[WARNING] read_vtk: failed to read vertex scalar " << i << std::endl;
vertex_scalars->clear();
break;
}
}
break;
}
}
}
if (volume_scalars != nullptr) {
volume_scalars->clear();
// Reset stream or continue from current position
while (std::getline(is, line)) {
if (line.find("CELL_DATA") != std::string::npos) {
std::size_t ndata;
ss = std::stringstream(line);
std::getline(ss, tmp, ' '); // skip "CELL_DATA"
ss >> ndata;
// Skip SCALARS and LOOKUP_TABLE lines
std::getline(is, line); // SCALARS line
std::getline(is, line); // LOOKUP_TABLE line
volume_scalars->resize(ndata);
for (std::size_t i = 0; i < ndata; ++i) {
if (!(is >> (*volume_scalars)[i])) {
std::cerr << "[WARNING] read_vtk: failed to read volume scalar " << i << std::endl;
volume_scalars->clear();
break;
}
}
break;
}
}
}
return true;
}
template<typename LCC>
bool write_vtk_ascii(std::ostream& os,
const LCC& alcc,
const std::vector<float>* vertex_scalars,
const std::vector<float>* volume_scalars)
{
static_assert(LCC::dimension == 3 && LCC::ambient_dimension == 3,
"write_vtk() only supports 3D Linear_cell_complexes (3,3)");
// Write VTK header
os << "# vtk DataFile Version 2.0\n";
os << "CGAL Linear_cell_complex\n";
os << "ASCII\n";
os << "DATASET UNSTRUCTURED_GRID\n\n";
// Build vertex index map
std::unordered_map<typename LCC::Vertex_attribute_const_descriptor, std::size_t> vertex_index;
std::size_t nbpts = 0;
for (auto itv = alcc.vertex_attributes().begin(),
itvend = alcc.vertex_attributes().end(); itv != itvend; ++itv) {
vertex_index[itv] = nbpts++;
}
// Write points
os << "POINTS " << nbpts << " double\n";
for (auto itv = alcc.vertex_attributes().begin(),
itvend = alcc.vertex_attributes().end(); itv != itvend; ++itv) {
const auto& p = itv->point();
os << p.x() << " " << p.y() << " " << p.z() << "\n";
}
os << "\n";
// Count cells and build connectivity
std::size_t nbcells = 0;
std::size_t total_size = 0;
std::ostringstream cell_stream;
std::ostringstream type_stream;
for (auto itvol = alcc.template one_dart_per_cell<3>().begin(),
itvolend = alcc.template one_dart_per_cell<3>().end();
itvol != itvolend; ++itvol) {
++nbcells;
VTK_Cell_Type cell_type = get_vtk_cell_type(alcc, itvol);
type_stream << static_cast<int>(cell_type) << "\n";
if (cell_type == VTK_POLYHEDRON) {
// Generic polyhedron format write as face-vertex connectivity
std::vector<std::vector<std::size_t>> faces;
std::size_t cell_size = 1; // Start with 1 for number of faces
for (auto itface = alcc.template one_dart_per_incident_cell<2,3>(itvol).begin(),
itfaceend = alcc.template one_dart_per_incident_cell<2,3>(itvol).end();
itface != itfaceend; ++itface) {
faces.push_back(std::vector<std::size_t>());
auto& face = faces.back();
for (auto itvert = alcc.template darts_of_orbit<1>(itface).begin(),
itvertend = alcc.template darts_of_orbit<1>(itface).end();
itvert != itvertend; ++itvert) {
face.push_back(vertex_index[alcc.vertex_attribute(itvert)]);
}
cell_size += face.size() + 1; // +1 for face size
}
cell_stream << cell_size << " " << faces.size();
for (const auto& face : faces) {
cell_stream << " " << face.size();
for (auto v : face) {
cell_stream << " " << v;
}
}
cell_stream << "\n";
total_size += cell_size + 1; // +1 for cell size
} else {
// Standard cell types write vertex connectivity directly
std::vector<std::size_t> vertices;
for (auto itvert = alcc.template one_dart_per_incident_cell<0,3>(itvol).begin(),
itvertend = alcc.template one_dart_per_incident_cell<0,3>(itvol).end();
itvert != itvertend; ++itvert) {
vertices.push_back(vertex_index[alcc.vertex_attribute(itvert)]);
}
cell_stream << vertices.size();
for (auto v : vertices) {
cell_stream << " " << v;
}
cell_stream << "\n";
total_size += vertices.size() + 1;
}
}
// Write cells section
os << "CELLS " << nbcells << " " << total_size << "\n";
os << cell_stream.str();
os << "\n";
// Write cell types
os << "CELL_TYPES " << nbcells << "\n";
os << type_stream.str();
os << "\n";
// Write vertex scalars if provided
if (vertex_scalars != nullptr) {
if (vertex_scalars->size() != nbpts) {
std::cerr << "[ERROR] write_vtk: vertex_scalars size (" << vertex_scalars->size()
<< ") does not match number of vertices (" << nbpts << ")" << std::endl;
return false;
}
os << "POINT_DATA " << nbpts << "\n";
os << "SCALARS vertex_scalars float 1\n";
os << "LOOKUP_TABLE default\n";
for (float val : *vertex_scalars) {
os << val << "\n";
}
os << "\n";
}
// Write volume scalars if provided
if (volume_scalars != nullptr) {
if (volume_scalars->size() != nbcells) {
std::cerr << "[ERROR] write_vtk: volume_scalars size (" << volume_scalars->size()
<< ") does not match number of cells (" << nbcells << ")" << std::endl;
return false;
}
os << "CELL_DATA " << nbcells << "\n";
os << "SCALARS volume_scalars float 1\n";
os << "LOOKUP_TABLE default\n";
for (float val : *volume_scalars) {
os << val << "\n";
}
}
return true;
}
} // namespace internal
// Public interface implementation
template<typename LCC>
bool read_vtk(LCC& alcc,
const char* filename,
std::vector<float>* vertex_scalars,
std::vector<float>* volume_scalars)
{
CGAL_assertion(filename != nullptr);
std::ifstream file(filename);
if (!file.is_open()) {
std::cerr << "[ERROR] read_vtk: cannot open file " << filename << std::endl;
return false;
}
return internal::read_vtk_ascii(file, alcc, vertex_scalars, volume_scalars);
}
template<typename LCC>
bool write_vtk(const LCC& alcc,
const char* filename,
const std::vector<float>* vertex_scalars,
const std::vector<float>* volume_scalars)
{
CGAL_assertion(filename != nullptr);
std::ofstream file(filename);
if (!file.is_open()) {
std::cerr << "[ERROR] write_vtk: cannot open file " << filename << std::endl;
return false;
}
return internal::write_vtk_ascii(file, alcc, vertex_scalars, volume_scalars);
}
} // namespace CGAL
#endif // CGAL_LINEAR_CELL_COMPLEX_VTK_IO_IMPL_H

View File

@ -1,48 +1,45 @@
#include <CGAL/Linear_cell_complex_for_combinatorial_map.h>
#include <CGAL/Linear_cell_complex_vtk_io.h>
#include <cstdlib>
#include <fstream>
#include <cassert>
#include <cstdio>
#include <string>
#include <unordered_map>
#include <vector>
typedef CGAL::Linear_cell_complex_for_combinatorial_map<3, 3> LCC;
typedef LCC::Point Point;
int main()
{
LCC lcc_out;
int main() {
LCC lcc1, lcc2;
std::vector<float> v_scalars1, vol_scalars1;
std::vector<float> v_scalars2, vol_scalars2;
// Create a hexahedron
auto d = lcc_out.make_hexahedron(
lcc_out.create_vertex_attribute(Point(0, 0, 0)),
lcc_out.create_vertex_attribute(Point(1, 0, 0)),
lcc_out.create_vertex_attribute(Point(1, 1, 0)),
lcc_out.create_vertex_attribute(Point(0, 1, 0)),
lcc_out.create_vertex_attribute(Point(0, 0, 1)),
lcc_out.create_vertex_attribute(Point(1, 0, 1)),
lcc_out.create_vertex_attribute(Point(1, 1, 1)),
lcc_out.create_vertex_attribute(Point(0, 1, 1))
);
const std::string input_file = "data/beam-with-mixed-cells.vtk";
assert(CGAL::read_lcc_from_vtk(lcc1, input_file.c_str(), &v_scalars1, &vol_scalars1));
std::vector<float> v_scalars = {1,2,3,4,5,6,7,8};
std::vector<float> vol_scalars = {42.0f};
// Build index maps
std::unordered_map<LCC::Vertex_attribute_const_descriptor, std::size_t> vertex_indices;
std::size_t idx = 0;
for(auto itv = lcc1.vertex_attributes().begin(), itvend = lcc1.vertex_attributes().end(); itv != itvend; ++itv)
vertex_indices[itv] = idx++;
std::unordered_map<LCC::Dart_const_descriptor, std::size_t> volume_indices;
idx = 0;
for(auto itvol = lcc1.one_dart_per_cell<3>().begin(), itvolend = lcc1.one_dart_per_cell<3>().end(); itvol != itvolend;
++itvol)
volume_indices[itvol] = idx++;
const char* fname = "tmp_test_lcc_vtk.vtk";
const char* tmp_file = "tmp_test_lcc_vtk.vtk";
assert(CGAL::write_lcc_to_vtk(
lcc1, tmp_file,
[&v_scalars1, &vertex_indices](const LCC& lcc, LCC::Dart_const_descriptor d) {
return v_scalars1[vertex_indices.at(lcc.attribute<0>(d))];
},
[&vol_scalars1, &volume_indices](const LCC& lcc, LCC::Dart_const_descriptor d) {
return vol_scalars1[volume_indices.at(d)];
}));
assert(CGAL::read_lcc_from_vtk(lcc2, tmp_file, &v_scalars2, &vol_scalars2));
assert(CGAL::write_vtk(lcc_out, fname, &v_scalars, &vol_scalars));
assert(lcc1.is_isomorphic_to(lcc2, false, true, true));
LCC lcc_in;
std::vector<float> v_scalars_in, vol_scalars_in;
assert(CGAL::read_vtk(lcc_in, fname, &v_scalars_in, &vol_scalars_in));
assert(lcc_in.number_of_vertex_attributes() == lcc_out.number_of_vertex_attributes());
assert(lcc_in.template number_of_cells<3>() == lcc_out.template number_of_cells<3>());
assert(v_scalars_in.size() == v_scalars.size());
assert(vol_scalars_in.size() == vol_scalars.size());
std::remove(fname);
std::remove(tmp_file);
return EXIT_SUCCESS;
}

View File

@ -0,0 +1,11 @@
# vtk DataFile Version 2.0
CGAL Linear_cell_complex
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 0 double
CELLS 0 0
CELL_TYPES 0