Error with fft-based-hauptmann in run_all.py

[alexpattyn]

Description

Looks like there is a divide by zero error in fft_based_hauptmann_2018.py when the user runs run_all.py. This should just work out of the box. I can take a closer look.

Error

/var/home/apattyn/Projects/standardised-image-reconstruction/image_reconstruction/reconstruction_utils/beamforming/fft_based_hauptmann_2018.py:248: RuntimeWarning: divide by zero encountered in true_divide
  sf = np.divide(sf, 2*w)
/var/home/apattyn/Projects/standardised-image-reconstruction/image_reconstruction/reconstruction_utils/beamforming/fft_based_hauptmann_2018.py:248: RuntimeWarning: invalid value encountered in true_divide
  sf = np.divide(sf, 2*w)

[alexpattyn]

Seems like this is expected? In fft_based_hauptmann_2018.py we have the following around line 241:

        w = c*kgrid_back_c.kx
        w_new = kgrid_back.c*kgrid_back.k

        # scalig factor
        sf = np.square(w/c) - np.square(kgrid_back_c.ky)
        sf = c*c*np.sqrt(sf.astype(np.complex))

        sf = np.divide(sf, 2*w)

Where w is causing the divide by zero issue, and it seems to be sound speed multiplied by the x-axis which goes from -X:X. Which will cause the error above. Numpy returns inf on divide by zero, so perhaps this should be replaced by a 1?

[bencox]

Shouldn't be a 1 but should be the limit as w goes to zero. The matlab version has the following line to correct the divide by zero: sf(w == 0 & kgrid.ky == 0 & kgrid.kz == 0) = c ./ 2; so probably the inf should be replaced by c/2. We should check with Jenni Poimala, who I think wrote this version.

[alexpattyn]

Once we get clarification I can make a PR for the fix.

[alexpattyn]

The above code is a little confusing to me, since it doesn't fix the divide by zero error. it just overwrites the regions where w == 0 & kgrid.ky == 0 & kgrid.kz == 0 within sf and may leave some Inf values. I.e. there may be some indexes where w==0 but kgrid.ky!=0.

It would be simpler if we had something like w(w==0)=foo, but I am not sure how the original algorithm is written.

I'll wait until Jenni comments.

[jpoimala]

My code is based on the matlab version, so Inf values are handled in the same way as in the Ben's comment.

[alexpattyn]

@jpoimala I can now see that below the code I took a snip of, thanks! :+1:

idx1 = np.where((w == 0) & (kgrid_back_c.ky == 0))
sf[idx1] = c/2

idx2 = np.where(np.abs(w) < np.abs(c*kgrid_back_c.ky))
sf[idx2] = 0

idx = np.where(np.isnan(sf))
sf[idx] = 0

Still would be nice to fix the error at the source (using w(w==0)=foo). Does anyone know the paper the code was based off of? Wanted to take a peak and see if they say anything about this.

[jgroehl]

We collect the references in the reconstruction_algorithm classes. For this case, the cited reference is the following:

Hauptmann, Andreas, et al.
"Approximate k-space models and deep learning for fast photoacoustic reconstruction."
International Workshop on Machine Learning for Medical Image Reconstruction. Springer, Cham, 2018.

Since these references are pretty hidden in the code, should we maybe add a list of all implemented algorithms and the respective citations in the README?

[bencox]

The original paper from 2001 is here: https://doi.org/10.1088/0031-9155/46/7/309. The same algorithm was also published in 2002 here: https://doi.org/10.1109/TMI.2002.801172

[alexpattyn]

Thanks all :+1: I'll review that paper to see if there is a better way to handle that divide by zero issue.