elrnv / vtkio

Visualization ToolKit (VTK) file parser and writer
Apache License 2.0
57 stars 13 forks source link

Exported data for point cloud is not readable by Paraview 5.11 #38

Closed misterjoa closed 1 year ago

misterjoa commented 1 year ago

Hello ! Thank you for this nice piece of software which enabled exporting vtk directly from rust.

I have a problem to read point clouds created with vtkio (PolyData containing only points, no cells).

Files are opened only if those contain non-empty cells.

Minimal example (from doc)

readable by paraview

use vtkio::model::*; // import model definition of a VTK file

let mut vtk_bytes = Vec::<u8>::new();

Vtk {
    version: Version::new((2,0)),
    byte_order: ByteOrder::BigEndian,
    title: String::from("Triangle example"),
    file_path: None,
    data: DataSet::inline(PolyDataPiece {
        points: vec![0.0f32, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, -1.0].into(),
        polys: Some(VertexNumbers::Legacy {
            num_cells: 1,
            vertices: vec![3, 0, 1, 2]
        }),
        data: Attributes::new(),
        ..Default::default()
    })
}.write_xml(&mut vtk_bytes);

not readable by paraview (does load but no data is on screen)

use vtkio::model::*; // import model definition of a VTK file

let mut vtk_bytes = Vec::<u8>::new();

Vtk {
    version: Version::new((2,0)),
    byte_order: ByteOrder::BigEndian,
    title: String::from("Triangle example"),
    file_path: None,
    data: DataSet::inline(PolyDataPiece {
        points: vec![0.0f32, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, -1.0].into(),
        data: Attributes::new(),
        ..Default::default()
    })
}.write_xml(&mut vtk_bytes);

no data either for this example

use vtkio::model::*; // import model definition of a VTK file

let mut vtk_bytes = Vec::<u8>::new();

Vtk {
    version: Version::new((2,0)),
    byte_order: ByteOrder::BigEndian,
    title: String::from("Triangle example"),
    file_path: None,
    data: DataSet::inline(PolyDataPiece {
        points: vec![0.0f32, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, -1.0].into(),
        polys: Some(VertexNumbers::Legacy {
            num_cells: 0,
            vertices: vec![]
        }),
        data: Attributes::new(),
        ..Default::default()
    })
}.write_xml(&mut vtk_bytes);
elrnv commented 1 year ago

You should be able to export a point cloud using the verts field as follows:

use vtkio::model::*; // import model definition of a VTK file

let mut vtk_bytes = Vec::<u8>::new();

Vtk {
    version: Version::new((2,0)),
    byte_order: ByteOrder::BigEndian,
    title: String::from("Triangle example"),
    file_path: None,
    data: DataSet::inline(PolyDataPiece {
        points: vec![0.0f32, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, -1.0].into(),
        verts: Some(VertexNumbers::Legacy {
            num_cells: 3,
            vertices: vec![0, 1, 2]
        }),
        data: Attributes::new(),
        ..Default::default()
    })
}.write_xml(&mut vtk_bytes);

Does this work for you?

misterjoa commented 1 year ago

Hello, thank you for your answer!

It did not work with the solution that you proposed. I solved the problem, but there may be a need for a patch to manage cloud points properly without manual specification of connectivity and offset

In order to start debugging I wrote a file with no compression


    point_cloud_vtk = Vtk {...}

    let vtk_file = point_cloud_vtk
        .try_into_xml_format(vtkio::xml::Compressor::None, 0)
        .unwrap();
    let file = File::create("point_cloud.vtp").unwrap();
    xml::write(&vtk_file, file).unwrap();

Regarding debugging, I have a suggestion to make base64 encoding of arrays optional, as debugging is far easier with ascii-encoded arrays.

VertexNumbers::Legacy has an algorithm that generates (connectivity, offset) that seems not compatible with point clouds (the code tripped on some asserts and when commenting those, paraview could not read the file).

Using this resource https://discourse.paraview.org/t/xml-format-for-a-file-that-only-contains-points/7308/4

I managed to create a minimal cloud point file (.vtp otherwise paraview does not recognize it as such) as

<?xml version="1.0"?>
<VTKFile type='PolyData' version='0.1' byte_order='LittleEndian' header_type='UInt64'>
  <PolyData>
    <Piece NumberOfPoints='8' NumberOfVerts='8'>
      <Points>
        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="ascii">
                    0 0 0
            1 0 0
            1 1 0
            0 1 0
            0 0 1
            1 0 1
            1 1 1
            0 1 1
        </DataArray>
      </Points>
      <Verts>
        <DataArray type="Int32" Name="offsets">
        1 2 3 4 5 6 7 8
        </DataArray>
        <DataArray type="Int32" Name="connectivity">
        0 1 2 3 4 5 6 7
        </DataArray>
      </Verts>
    </Piece>
  </PolyData>
</VTKFile>

So in the end the following, "manual" solution worked

pub struct Node {
    pub id: i32,
    pub coordinates: Vec<f64>,
    pub dim: usize,
}

pub fn build_vtk_point_cloud(vert_list: &Vec<Node>) -> Vtk {
    let mut points: Vec<f32> = Vec::new();
    for vert in vert_list.iter() {
        points.extend(vert.coordinates.iter().map(|f| *f as f32));
    }

    Vtk {
        version: Version::new((2, 0)),
        byte_order: ByteOrder::BigEndian,
        title: "test".to_string(),
        file_path: None,
        data: DataSet::inline(PolyDataPiece {
            verts: Some(VertexNumbers::XML {
                connectivity: (0..vert_list.len() as u64).collect(),
                offsets: (1..(vert_list.len() + 1) as u64).collect(),
            }),
            points: vtkio::IOBuffer::F32(points),
            data: Attributes::new(),
            ..Default::default()
        }),
    }
}

In the generated file, I observed some anchors that are not related to anything (<PointData/> and <CellData/>), but the file was readable still :+1:

<?xml version="1.0"?>
<VTKFile type="PolyData" version="2.0" byte_order="BigEndian" header_type="UInt64">
  <PolyData>
    <Piece NumberOfPoints="3" NumberOfLines="0" NumberOfStrips="0" NumberOfPolys="0" NumberOfVerts="3">
      <PointData/>
      <CellData/>
      <Points>
        <DataArray type="Float32" format="binary" NumberOfComponents="3">AAAAAAAAACTERB3jwCLef0Od4ABFP795wBzg9EOqgABCcf+0w/JgAENoAAA=</DataArray>
      </Points>
      <Verts>
        <DataArray type="UInt64" Name="connectivity" format="binary" NumberOfComponents="1">AAAAAAAAABgAAAAAAAAAAAAAAAAAAAABAAAAAAAAAAI=</DataArray>
        <DataArray type="UInt64" Name="offsets" format="binary" NumberOfComponents="1">AAAAAAAAABgAAAAAAAAAAQAAAAAAAAACAAAAAAAAAAM=</DataArray>
      </Verts>
    </Piece>
  </PolyData>
</VTKFile>

I hope that this helps!

elrnv commented 1 year ago

Great! I am happy you got the issue resolved. I realized that the format for vertices should have a 1 prepended for every vertex index to be consistent with all other polygon types.

...
            vertices: vec![1, 0, 1, 1, 1, 2]
...

This is definitely unintuitive for point clouds, and should be documented somewhere, since it even tripped me up!

Also good point on base64 encoding! The functionality is there but the API is very thin for that, you would have to change update the desired DataArrayFormat manually for the time being.

Thanks for bringing this issue up, this helps a lot!