How to read a binary structured_points data file generated in Fortran by C++ and VTK?

Question

I am new to VTK, and I am trying to read a binary structured_points (image data) file in vtk but getting the "Error reading binary data" error with a blank window. I am using a simple fortran program to create the data file (gaussian field data) as follows,

 program gaussian

            implicit none
            integer i, j, k
            character(1) c

    open(unit=100, file='energyDensity.vtk',
 1        form="unformatted",access="stream")
    write(100) '# vtk DataFile Version 3.0', new_line(c)
    write(100) 'First time trying vtk import \n', new_line(c)
    write(100) 'BINARY', new_line(c)
    write(100) 'DATASET STRUCTURED_POINTS', new_line(c)
    write(100) 'DIMENSIONS 101 101 101', new_line(c)
    write(100) 'ORIGIN 0 0 0', new_line(c)
    write(100) 'SPACING 1 1 1', new_line(c)
    write(100) 'POINT_DATA 1030301', new_line(c)
    write(100) 'SCALARS volume_scalars double 1', new_line(c)
    write(100) 'LOOKUP_TABLE default', new_line(c)


            do k = -50,50
            do j = -50,50
            do i = -50,50

            write(100) 50.*exp(float((-(i*i+j*j+k*k))/25)) 

            enddo
            enddo
            enddo
    close(100)

endprogram

VTK reads and plots the data nicely if the data is in ASCII format (see the picture below)

I am using the following code C++ in VTK to read the data (without any setting for endian),

vtkNew<vtkStructuredPointsReader> reader;
  reader->SetFileName (argv[1]);
  reader->Update();

I have tried searching a lot on internet but couldn't find the correct way of reading structured_points data in binary. It seems there is no way of setting the endian thing for structured points reader as well. I am not sure what do here. Any help would be really appreciated.

Are you sure the data is correct? Does Paraview or Visit open your binary data correctly? — Vladimir F
Vladimir, Paraview managed to open the file but the plot looked nothing like the picture above. Paraview showed four weird spheres (regularly spaced), which doesn't make sense. So, no, paraview does not read the binary file correctly. I have no idea how to resolve this issue. The way I produce binary files in Fortran (as in the code above) has always worked fine for me when importing in softwares like Mathematica. — singularity

Vladimir F Vladimir F · Accepted Answer · 2019-02-22T13:15:02

There are several issues in your writing procedure.

First, you are writing single precision data (by using float()) but you claim that the data is double in the VTK header.

Second, the binary data in legacy VTK files should be big-endian.

You are also doing integer division in (-(i*i+j*j+k*k))/25) but that might be intentional (or not).

This works OK

 program gaussian

   use iso_fortran_env

            implicit none
            integer i, j, k
            character(1) c

    open(unit=100, file='energyDensity.vtk', &
         form="unformatted",access="stream")
    write(100) '# vtk DataFile Version 3.0', new_line(c)
    write(100) 'First time trying vtk import \n', new_line(c)
    write(100) 'BINARY', new_line(c)
    write(100) 'DATASET STRUCTURED_POINTS', new_line(c)
    write(100) 'DIMENSIONS 101 101 101', new_line(c)
    write(100) 'ORIGIN 0 0 0', new_line(c)
    write(100) 'SPACING 1 1 1', new_line(c)
    write(100) 'POINT_DATA 1030301', new_line(c)
    write(100) 'SCALARS volume_scalars double 1', new_line(c)
    write(100) 'LOOKUP_TABLE default', new_line(c)


            do k = -50,50
            do j = -50,50
            do i = -50,50

            write(100) SwapB64(exp(real((-(i*i+j*j+k*k)), real64) / 25))

            enddo
            enddo
            enddo
    close(100)
contains

    elemental function SwapB64(x) result(res)
      real(real64) :: res
      real(real64),intent(in) :: x
      character(8) :: bytes
      integer(int64) :: t
      real(real64) :: rbytes, rt
      equivalence (rbytes, bytes)
      equivalence (t, rt)

      rbytes = x

      t = ichar(bytes(8:8),int64)

      t = ior( ishftc(ichar(bytes(7:7),int64),8),  t )

      t = ior( ishftc(ichar(bytes(6:6),int64),16), t )

      t = ior( ishftc(ichar(bytes(5:5),int64),24), t )

      t = ior( ishftc(ichar(bytes(4:4),int64),32), t )

      t = ior( ishftc(ichar(bytes(3:3),int64),40), t )

      t = ior( ishftc(ichar(bytes(2:2),int64),48), t )

      t = ior( ishftc(ichar(bytes(1:1),int64),56), t )

      res = rt

    end function
endprogram

How to read a binary structured_points data file generated in Fortran by C++ and VTK?

2 Answers