Decompose complex shear wave field into propagation directions

Question

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
load('waveField_single.mat')

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

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