k-Wave

1. DominikWaibel
   

   Member
   Joined: Apr '17
   Posts: 1

   Hi,
   is there any way to create the needed HDF5 file with a C++ code? I am writing a plugin in C++ which should execute the "kspaceFirstOrder3D-OMP.exe" and therefore needs the correct HDF5-File.
   Optimal for me would be a code which creates the HDF5 File for a given .nrrd-file containing the initial pressure distribution.
   Thank you very much for your help!
   Best,
   Dominik

   Posted 7 years ago #

2. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Dominik,

   We don't currently have any C++ code to write the input file (although you could potentially modify the i/o component of kspaceFirstOrder3D-OMP used to write the output file). The complete list of fields is given in the k-Wave manual. It shouldn't be too difficult if you familiarise yourself with the HDF5 library.

   Brad.

   Posted 7 years ago #

3. Jiri Jaros
   

   Developer
   Joined: Feb '12
   Posts: 115

   Hi Dominik,
   we've been working on a Python code to generate input files for k-Wave. We've been also working on voxelization of 3D models stored in the OBJ format. I could potentially send you our codes for inspiration. Any way, python has a very simple interface for the HDF5 so it's not so difficult to prepare the input file.

   Jiri

   Posted 7 years ago #

4. lbruju
   

   Member
   Joined: Mar '18
   Posts: 19

   Hi Jiri Jaros,

   I am also working on a Python code to generate input files for k-wave. As you propose, could you please post your Python code which generates the input files for inspiration?

   Thanks a lot,

   lbruju

   Posted 6 years ago #

5. JakobD
   

   Member
   Joined: Mar '18
   Posts: 1

   Hi everyone

   I don't have matlab, so I also want to generate a hdf5 file with python (for after execute kSpaceFirstOrder3D-OMP). But I have problems when generating this file.
   (1) Is it necessary to add one by one all the attributes for each element? It seems to be quite long and not very optimized...
   (2) I found a hdf5 generated with matlab (on the forum) and I don't understand what représente the variables "ddx_k_shift_neg_r" (and same with neg and for the other axis y and z)
   (3) With the matlab code, lots of functions permit to create different shape of source or sensor (makeArc, makeBowl...). How it could be possible to mimic these functions in hdf5 file?
   (4) and last question, did I well understand that there is no equivalent function of kWaveDiffusion in C++?

   Finally, I will be very grateful if someone can share an example of python script generating a hdf5 file !

   Jakob

   Posted 6 years ago #

6. Jiri Jaros
   

   Developer
   Joined: Feb '12
   Posts: 115

   Hi Jakob, Ibruju,
   I will upload a sample matlab script later today.

   But first, let's answer your questions.
   (1) Yes, it is necessary to add HDF5 attributes to all datasets. The reason is to double check the data has been generated correctly.
   (2) ddx_k_shift_neg_r - these are the spatial shifts in the x, y and z direction for staggered grids. The _r suffix says, the size of the domain is reduced in this direction. Since the pressure/velocity matrices are real, the FFTs are symmetric, and we can only store the first half. Consequently, the ddx_k_shift_neg is truncated from N to N/2 + 1 complex numbers. The reason why this is generated by Matlab not C++ is that we wanted some freedom in the staggered grid definition. - the code to generate these values will be part of the script I'm gonna upload.
   (3) The C++ code only supports an index sensor mask. You will have to identify the grid points on the arc manually and fill the index sensor mask by yourself.
   (4) There's no diffusion code in C++ at the moment. We thought it was so fast it wasn't worth to reimplemented it in C++ at that moment.

   Best
   Jiri

   Posted 6 years ago #

7. lbruju
   

   Member
   Joined: Mar '18
   Posts: 19

   Hi,

   Thanks a lot for your answer.

   I think I manage to write a script which generates a hdf5 file. But I obtain surprising results.
   In 'p_source_index' are the index of half of a sphere.
   In 'p_source_input', I defined a list 'source' by: 'source=amp*amp*np.sin(2*np.pi*freq*t_array)' with 'amp=10.0 #Pa'.
   I defined all other dataset. At the end, I obtained a pressure field around 1E8, that seems to be very high...
   Do you know what could be the origin of the error?

   Thanks

   lbruju

   Posted 6 years ago #

8. lbruju
   

   Member
   Joined: Mar '18
   Posts: 19

   Hi again,
   I have another question. I don't get the difference between defining a transducer and defining a source mask and a time varying pressure?

   In my study case, I want to mimic the evolution of pressure with a focused transducer. Do I need to define a transducer (with transducer_source_flag = 1) or to define a source (with 'p_source_index' and 'p_source_input')?

   If I need to define a transducer, can you confirm that "u_source_index" is equivalent to "p_source_index" and "transducer_source_input" is equivalent to "p_source_input". But what is the role of "delay_mask"? Is that for defining the shape of the transducer?

   Thanks a lot,

   lbruju

   Posted 6 years ago #

9. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi lbruju,

   The transducer source allows a single time series to be used for multiple grid points within the simulation. It always uses a velocity source in the x-direction (u_source_index is for a velocity source, and p_source_index is for a mass or pressure source). For details, see the Ultrasound Imaging examples. I wouldn't recommend using this, as we will likely implement this functionality in a slightly different (and more general) way in a future release. So in short, define your source using p_source_index and p_source_input.

   Brad.

   Posted 6 years ago #

10. Jiri Jaros
    

    Developer
    Joined: Feb '12
    Posts: 115

    Hi Ibruju,
    My PhD student, Filip Vaverka, has created a prototype of python generator of k-wave input files. You can download it here.

    http://www.fit.vutbr.cz/~jarosjir/download/k-wave-generator-example.zip

    Jiri

    Posted 6 years ago #

11. lbruju
    

    Member
    Joined: Mar '18
    Posts: 19

    Hi,

    Thanks for the files.
    I am lost regarding indexing and shape. It's written in the k-wave_user_manual in the chapter "Format of the C++ HDF5 Files" that the dimension sizes are given following the MATLAB indexing convention (which is row-major ordering) and for creating files outside
    MATLAB, dataset dimensions should be given as (Nz, Ny, Nx) .
    But in the python script you shared, the shape are given and for example, for ddx_k_shift_neg the shape (1,1,Nz) is indicated. I don't understand why the shape is not (Nz,1,1) ?

    After executing kspaceFirstOrder3D-OMP.exe, the resulting acoustic pressure field is in a matrix of size (np.size(p_source_input), Nx*Ny*Nz) . I want to reshape the acoustic pressure field to have in final_press[t][x,y,z] the resulting pressure at time t, at point (x,y,z) . For each time, do I need to reshape among (Nz,Ny,Nx) (i.e p.reshape(pressure[time,0:Nx*Ny*Nz],(Nz,Ny,Nx))) or among (Nx,Ny,Nz) (i.e p.reshape(pressure[time,0:Nx*Ny*Nz],(Nz,Ny,Nx)))?

    Thanks a lot,

    lbruju

    Posted 6 years ago #

12. Jiri Jaros
    

    Developer
    Joined: Feb '12
    Posts: 115

    Hi lbruju,
    HDF5 is transparent to C++ (row-major ordering) / Matlab (column-major ordering) and respects the convention of the host language. Consequently, when creating a matrix in Matlab/Fortran, you provide the dimensions in the Matlab ordering (Nx, Ny, Nz). However, when doing the same thing in C++ you work with (Nz, Ny, Nz). HDF5 is responsible for data rotation when reading from different languages. The k-Wave manual is meant for Matlab users and thus respects the Matlab conventions.
    I'd recommend to use h5ls and h5dump commandline utilities to see the organization.

    Now the question is what about Python. And the answer is, it depends :) If you use NumPy, you can choose whether you want to work with C-like of Fortran-like ordering by the parameter order.
    When using the C-ordering, which is default, you provide dimension size in this order [Nz, Ny, Nx].

    We have now realized that there are unnecessary transpositions in the script I shared. Simply said, the X-Z are swapped and the data rotated, which is actually the same thing as before. So this is not necessary (you don't need to transpose data, rotate indices, etc.)

    Best
    Jiri

    Posted 6 years ago #

13. lbruju
    

    Member
    Joined: Mar '18
    Posts: 19

    Hi Jiri Jaros,

    Thanks a lot for your answer, I removed all the transpositions and swapped x and z axis, and my simulation seems to be correct now.

    Best regards,
    lbruju

    Posted 6 years ago #

14. Tristan2
    

    Member
    Joined: Jan '19
    Posts: 3

    Hi Jiri Jaros,

    thank you very much for the python input generator. I get some errors if i try to use the set_u_source method.

    I'm not sure if i completely understand what the script is doing but is it possible that something like

    "self.input_data_sets['velocity_source_terms']['u_source_many'].value = 1 if ux_source_series.shape[0] > 1 or ux_source_series.shape[0] > 1 or ux_source_series.shape[0] > 1 else 0"

    is missing from the set_u_source method?

    Best regards,
    Tristan

    Posted 5 years ago #

15. Filip Vaverka
    

    Developer
    Joined: Mar '19
    Posts: 4

    Yes, you are right it is missing from that method. Also "self.input_data_sets['simulation_flags']['ux_source_flag'].value" (and other "X_source_flag") are wrongly set to length of the signal series instead to index of number of timesteps before the source series begins to be introduced in the domain.

    Posted 5 years ago #

16. Patrick
    

    Member
    Joined: Jun '19
    Posts: 14

    Hi,

    I have attempted to use your Python script to generate the input file, however, running the exact same configuration in Matlab gives drastically different results. (p_max nearly 100% different). Moreover, with identical timesteps, the input file generated by the Python script will produce NaN pressure values, while the Matlab code will run without problem. What could be the cause of this? And are there any plans for full support of this Python script? I know many of us would highly appreciate it.

    Thanks,

    PC

    Posted 5 years ago #

17. Jiri Jaros
    

    Developer
    Joined: Feb '12
    Posts: 115

    Hi pcook,
    And are there any plans for full support of this Python script? - Yes, I've reduced the complexity of the input file so that the derivative operators, pml and shift variables are generated inside the c++ code. This will be in the next release.

    The more complicated question is about wrong results. There must be nonsence in the data. NaN means some variables are ill-defined leading to divisions by zero, etc. I'd suspect some PML setting or derivative operators ddx_shift....

    I'd recommend cross-checking the python input file with the Matlab input file. Simply generate those two files, then open matlab and read dataset by dataset using h5read and compare the contents.

    Best
    Jiri

    Posted 5 years ago #

18. Patrick
    

    Member
    Joined: Jun '19
    Posts: 14

    Dr. Jaros,

    Fantastic, that's great news. I look forward to the next release.

    I will read through both generated files and see where the differences were and try to track down my error.

    Thanks for your help!
    PC

    Posted 5 years ago #

19. djps
    

    Member
    Joined: Sep '15
    Posts: 5

    I am trying to run some diagnostic simulations without matlab, and can't write hdf5 files in octave, so have to do pretty much everything in python. It makes cross-checking between things generated pretty difficult!

    For clarity: from the files provided, where is it necessary to remove the transpositions - in the example_generate_input.py or kwave_input_file.py? If this is removed, shouldn't arrays be defined as np.zeros((Nz,Ny,Nx)) though?

    Thanks, david

    Posted 4 years ago #

20. Filip Vaverka
    

    Developer
    Joined: Mar '19
    Posts: 4

    Hi djps,

    The issue with Python/numpy is that it defaults to C ordered arrays, whereas MATLAB is F ordered. It seems that there are two options to deal with this:
    1) Reorder the data that are stored in the file (this leads to transpositions as "array.transpose(2, 1, 0)).
    2) Flip whole simulation instead (this is done by swapping x/z axis and keeping the data in C order).

    Either way, we prepared new version (its complete rewrite - hence the delay) of the Python script with some examples, which should work with latest k-Wave release.
    The script should support both "Flip" and "Reorder" (by default and bit more tested) modes. The scripts require Python 3, h5py, numpy, scipy and matplotlib for "driver example" and can be downloaded from https://www.fit.vut.cz/~ivaverka/k-wave-gen-lite.zip.
    Once again this is not thoroughly tested and therefore provided "as is". Any comments and questions can be posted here or you can email me directly.

    Posted 4 years ago #

RSS feed for this topic

Reply »

You must log in to post.