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 Героям слава, Feb 22 '19 at 09:11
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, Feb 22 '19 at 10:14

Vladimir F Героям слава · Answer 1 · 2019-02-24T08:23:14.640

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

Vladimir, thanks for the response. The VTK user guide is not very clear about legacy formats only accepting big_endian. By the way, there is a very simple trick for converting data to big_endian in fortran using the "convert" parameter in the "open" statement. I will post the code that worked for me below. Also, changing float to real(kind=8) helped. The float had worked fine with ASCII format and so I forgot to change it to double. And yes, I am deliberately dividing by 25 (gaussian decays fast this way). Thanks a lot. — singularity, Feb 23 '19 at 18:44
@phd-applicant I do know about that function and I do not like it. Various compilers offer various bodges, but nothing is really portable. — Vladimir F Героям слава, Feb 23 '19 at 18:52
@VladimirF, `FLOAT` is a standard specific name for `REAL` taking a default integer input. Table 16.3 in some editions of the standard lists it. — user5713492, Feb 24 '19 at 01:44

singularity · Accepted Answer · 2019-02-24T04:58:11.653

Following Vladimir's suggestions about legacy formats only accepting big_endian, a few simple changes to my fortran code above did the trick for me. Setting the "convert" parameter in the "open" file statement to "big_endian" solved the biggest problem. Also, float should be changed to "real(real64)" to be consistent with double data type in the vtk data file.

    program gaussian

    use iso_fortran_env

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

    open(unit=100, file='energyDensity.vtk',
 1        form="unformatted",access="stream",convert="big_endian")
    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

            x=50.*exp(real((-(i*i+j*j+k*k))/25,real64))

            write(100) x

            enddo
            enddo
            enddo
    close(100)

    endprogram

Do note that `real(kind=8)` is not portable and does NOT mean 8-byte. See https://stackoverflow.com/questions/838310/fortran-90-kind-parameter and https://stackoverflow.com/questions/3170239/fortran-integer4-vs-integer4-vs-integerkind-4?noredirect=1&lq=1 I used the `real64` parameter for a reason. — Vladimir F Героям слава, Feb 23 '19 at 18:53
Vladimir, thanks. I replaced kind=8 with real64. I will keep the convert parameter in the open statement for now as it is a bit compact. Also, I only ever use intel and gfortran compilers, which seem to support the convert keyword. — singularity, Feb 24 '19 at 04:54

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

2 Answers2