Using phase contrast MRI, I recorded images of harmonic propagating shear waves in in vivo brain tissue in a 2D axial slice (126x126 pixel at 1.6mm x 1.6mm resolution). I recorded 8 images over one vibration cycle, from which I extracted the harmonic motion (complex valued wave field) using the temporal Fourier transform.

Antomic image of 2D axial slice Real part of complex wave field

Here you see the animated wave in a different slice for illustration:

Animated wave field

Next, I would like to decompose this wave field into different propagation directions or plane waves (locally plane) in order to calculate the local phase gradient on those plane waves. This does not work for superimposed waves. Therefore I calculated the spatial Fourier transform of the wave field and multiplied with directional filters (gaussian filter) and performed the inverse Fourier transform to end up with decomposed waves in for example 8 directions. However I am afraid that the filter in Fourier space, as it is shown in the image, introduces artifacts to the filtered wave fields. As a final step, I use the decomposed wave fields to reconstruct a 2D map of local wave numbers and wave speed based on the gradient of the phase of the complex waves. The result is illustrated here:

enter image description here

How would you decompose a complex valued shear wave field into it's main propagation direction such that the plane wave assumption fo the phase gradient calculation is fullfilled?

Directional Filter The Matlab code to decompose the wavefield into 8 propagation direction is given below. Moreover example data for Matlab (waveField_single.mat) and a plot script can be downloaded here: https://drive.google.com/drive/folders/1g3_so9WTx5N4YjeSGl94kvCrqQD3QS76?usp=sharing

% Filter waveField into 8 propagation directions

n = 126;
reso = 1.6;
nTheta = 8;% number of directions

fourierWF = fftshift( fft2(waveField) );    % in-plane fourier transformation
k1 = -( (0:n-1)-fix(n/2) ) / (n*reso) * 2 * pi;%[rad/m] wavenumber in 1st direction
k2 =  ( (0:n-1)-fix(n/2) ) / (n*reso) * 2 * pi;%[rad/m] wavenumber in 2nd direction
[theta, rho] = cart2pol( repmat(k2,[n,1]), repmat(k1',[1, n]) );% transform to polar coordinates

% calculate the angles of directions
sigmaTheta = 2 * pi / nTheta;
thetaValue= 2 * pi * linspace( 0, 1-1/nTheta, nTheta);%[rad]
filter = zeros( n ,n , nTheta );

for iTheta = 1 : nTheta % loop over directions
    currentThetaValue = thetaValue( iTheta );%[rad] selected theta value
    currentTheta = angle( exp( 1i * (theta-currentThetaValue) ) ); %[rad]  rotated theta
    filter(:,:,iTheta) = 0.4 * exp( -1/2 * ( currentTheta/sigmaTheta ).^2 );% gaussian function with selected theta (0.4=theta function(1/sqrt(e))
    filterdFourierWF = fourierWF .* filter(:,:,iTheta);
    waveField_filt(:, :, 1, iTheta) = ifft2(ifftshift(filterdFourierWF));% inverse fft
  • $\begingroup$ What's wrong with your field decomposition? You seem to have picked a rather high temporal frequency component (compared to the frequency content of your original signal) for the images shown, but they look nice! $\endgroup$ Oct 9, 2021 at 5:06
  • $\begingroup$ For clarification, do you mean that in the bottom pictures, the wavelength of the decomposed images (right side) appears shorter than the wavelength of the input image on the left? $\endgroup$
    – Helge
    Oct 11, 2021 at 18:56
  • $\begingroup$ Yes! The frequency content seems much higher in the decomposed field than in the original complexed value field. Were all the images computed with the same selected temporal frequency ( of 30 Hz, as you mention)? $\endgroup$ Oct 11, 2021 at 20:19
  • $\begingroup$ Yes, it's all the same frequency of 30 Hz. That's why I wonder if the Fourier decomposition using these radial cones alters the wavelength and therefore introduces an error. $\endgroup$
    – Helge
    Oct 12, 2021 at 8:33
  • $\begingroup$ I would assume somethnig's wrong with your algorithm, then... $\endgroup$ Oct 12, 2021 at 19:10

1 Answer 1


The spatial Fourier transform precisely does a decomposition of your field into plane waves. A local peak in the Fourier decomposition with $k_{x}$, $k_y$ coordinates indicate that the field has a plane wave component propagating at angle $\theta=atan\frac{k_{y}}{k_{x}}$ and wavenumber $k^{2}=k_{x}^2+k_{y}^{2}$. The gradient of the phase of a plane wave points along the propagation direction (given by the wavenumber vector) and with norm, the wavenumber. So I suggest you work directly in the Fourier domain to obtain the wavenumbers associated with the relevant peaks.

EDIT : this is the output of your code on my computer, the colormap is the standard "jet" from matlab

enter image description here

  • $\begingroup$ Thank you for your reply. Indeed that would give me the dominant frequencies and wave numbers present in my image. However I would like to obtain a 2D map of local wavenumbers, from which I can compute the local wave speed. Therefore I need to go back to image space if I am right. I added a final map of shear wave speed to my original post as illustration. $\endgroup$
    – Helge
    Oct 12, 2021 at 9:04
  • $\begingroup$ This is avery interesting problem, I know this is what they do in MRE to obtain maps of the shear moduli in elastic tissues. I never quite understood how the algorithm used for this application works to be honest, and I'm a bit dubious about the obtention of such a high resolution mapping in $V_{s}$! $\endgroup$ Oct 12, 2021 at 19:33
  • $\begingroup$ Suggestion: could you not detect isophase contours directly from your complex-valued harmonic field and take the gradient along the local normal to these contours? $\endgroup$ Oct 12, 2021 at 19:39
  • $\begingroup$ Thank you for your suggestion, I will give it a try. Indeed I am working in the field of MRE. $\endgroup$
    – Helge
    Oct 13, 2021 at 9:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.