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 "kspaceFirstOrder3DOMP.exe" and therefore needs the correct HDF5File.
Optimal for me would be a code which creates the HDF5 File for a given .nrrdfile containing the initial pressure distribution.
Thank you very much for your help!
Best,
Dominik
kWave
A MATLAB toolbox for the timedomain
simulation of acoustic wave fields
Creating the HDF5 File in C++
(18 posts) (8 voices)
Posted 2 years ago #

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 kspaceFirstOrder3DOMP used to write the output file). The complete list of fields is given in the kWave manual. It shouldn't be too difficult if you familiarise yourself with the HDF5 library.
Brad.
Posted 2 years ago # 
Hi Dominik,
we've been working on a Python code to generate input files for kWave. 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 2 years ago # 
Hi Jiri Jaros,
I am also working on a Python code to generate input files for kwave. As you propose, could you please post your Python code which generates the input files for inspiration?
Thanks a lot,
lbruju
Posted 1 year ago # 
Hi everyone
I don't have matlab, so I also want to generate a hdf5 file with python (for after execute kSpaceFirstOrder3DOMP). 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 1 year ago # 
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
JiriPosted 1 year ago # 
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 1 year ago # 
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 1 year ago # 
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 xdirection (
u_source_index
is for a velocity source, andp_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 usingp_source_index
andp_source_input
.Brad.
Posted 1 year ago # 
Hi Ibruju,
My PhD student, Filip Vaverka, has created a prototype of python generator of kwave input files. You can download it here.http://www.fit.vutbr.cz/~jarosjir/download/kwavegeneratorexample.zip
Jiri
Posted 1 year ago # 
Hi,
Thanks for the files.
I am lost regarding indexing and shape. It's written in the kwave_user_manual in the chapter "Format of the C++ HDF5 Files" that the dimension sizes are given following the MATLAB indexing convention (which is rowmajor 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, forddx_k_shift_neg
the shape(1,1,Nz)
is indicated. I don't understand why the shape is not(Nz,1,1)
?After executing kspaceFirstOrder3DOMP.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 infinal_press[t][x,y,z]
the resulting pressure at time t, at point(x,y,z)
. For eachtime
, do I need to reshape among(Nz,Ny,Nx)
(i.ep.reshape(pressure[time,0:Nx*Ny*Nz],(Nz,Ny,Nx))
) or among(Nx,Ny,Nz)
(i.ep.reshape(pressure[time,0:Nx*Ny*Nz],(Nz,Ny,Nx))
)?Thanks a lot,
lbruju
Posted 1 year ago # 
Hi lbruju,
HDF5 is transparent to C++ (rowmajor ordering) / Matlab (columnmajor 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 kWave 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 Clike of Fortranlike ordering by the parameter order.
When using the Cordering, 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 XZ 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
JiriPosted 1 year ago # 
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,
lbrujuPosted 1 year ago # 
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,
TristanPosted 6 months ago # 
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 months ago # 
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 2 months ago # 
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 illdefined leading to divisions by zero, etc. I'd suspect some PML setting or derivative operators ddx_shift....
I'd recommend crosschecking 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
JiriPosted 1 month ago # 
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!
PCPosted 1 month ago #
Reply
You must log in to post.