# Decompose complex shear wave field into propagation directions

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.

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

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:

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?

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
clearvars

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
end

• 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! Oct 9, 2021 at 5:06
• 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? Oct 11, 2021 at 18:56
• 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)? Oct 11, 2021 at 20:19
• 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. Oct 12, 2021 at 8:33
• I would assume somethnig's wrong with your algorithm, then... Oct 12, 2021 at 19:10

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.
• 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}$! Oct 12, 2021 at 19:33