Fitting 3D Gaussian to frequency signals in a 3D Power Spectrum

Hello,

I have a 3d power spectrum from calculating the FFT from a 3d volume by fftn(X). The power spectrum is loaded in a 3d array containing the intensity values, which are spherical symmetric to the array center. The aim is basically to measure the position (either in x,y,z or spherical coordinates) of signal intensities from the array center as exctly as possible. A parameter transformation can be simply done by cart2sph(x,y,z) since measuring distance and angle of the signals from the center is the most reasonable way.

I took out a small Volume (cube of 9 voxels per edge) from the whole power spectrum, containing two signals that I'm interested in (attached tiff-file "Test_volume_3D_Power_Spectrum"). In this way, I don't need to load the whole Power Spectrum, which is really a large file and the fit function should stay in this volume during the fit process. Since I have the coordinates of the box I just need to add the according vector from the center voxel.

My idea was to fit some kind of 3d gaussian density using mvnpdf(X,mu,Sigma), estimate the 3d vector "mu" and the covariance matrix "Sigma" as starting values from the intensity maxima in the dataset. To build the fit function for the start is not a problem either before or after the parameter transformation. but fitting it to the signal is not really straight forward, since this is 4d dataset and a visualization of the fit result is not simply done. This is how far I came.

%import tiff file to double array
input_name = 'Test_volume_3D_Power_Spectrum';
vol = importTIF(strcat(input_name,'.tif'),0); %uses the import function in the attached file

%start values estimated from raw data
mu1 = [5 6 6];
sigma1 = eye(3); 
A2 = 3.0e+10; %scaling factor for the gaussian
mu2 = [7 5 3];
sigma2 = eye(3);
A2 = 3.0e+10; %scaling factor for the gaussian

%setting up meash grid. Note: I think for a reasonable fit at subvoxel range, the steps must be at least 0.1 instead of 1. I put it to 1 
%in order to see, if the fit function is built correctly.
%also I want to have x, y and z to keep the option for non cubic volumes
x = 1:1:9;
y = 1:1:9;
z = 1:1:9;

[X,Y,Z] = meshgrid(x,y,z);

grid_values = [X(:) Y(:) Z(:)];

%building the fit function
PDF1 = mvnpdf(grid_values, mu1, sigma1);
PDF2 = mvnpdf(grid_values, mu2, sigma2);

Fit1 = A1*PDF1;
Fit2 = A2*PDF2;

%Prepare the built function for export for checking
Fit_data_1 = reshape(Fit1, length(x), length(y), length(z));
Fit_data_1 = reshape(Fit2, length(x), length(y), length(z));

Fit_combined = Fit_data_1 + Fit_data_2; 
Fit_combined = int32(Fit_combined);

%export to tif
save_name = strcat(ipnut_name,'_3D_Gauss.tif');
var = Tiff(save_name,'w');

tagstruct.ImageLength = size(Fit_combined,1);
tagstruct.ImageWidth = size(Fit_combined,2);
tagstruct.SampleFormat = 1;
tagstruct.Photometric = Tiff.Photometric.MinIsBlack;
tagstruct.BitsPerSample = 32;
tagstruct.SamplesPerPixel = 1;
tagstruct.Compression = Tiff.Compression.Deflate;
tagstruct.PlanarConfiguration = Tiff.PlanarConfiguration.Chunky;
var.setTag(tagstruct)

    if length(size(Fit_combined)) == 3
        for ii=1:size(data,3)
            var.setTag(tagstruct);
            var.write(data(:,:,ii));
            var.writeDirectory;
        end
    else
        var.write(data);
    end
    
close(var);

Now I want to fit the built function to my test volume. I studied fitgmdist(X,k) which should principally do, what I want, but I couldn't figure out, how to adapt it. Furthermore, I didn't find a way to extract the the parameters I want to have. But maybe I don't understand this method correctly.

Also some kind of surface plot fit should be possible, but the existing functions seem to restricted to 2D surfaces in 3D datasets. What I would need is a 3D surface fit within a 4D dataset. In general fitting a function on a smal cube of 8 voxels at the edge is quite challenging within a very large volume of several hundred voxels.

Maybe there are much easier ways to achieve my aim, but I got stuck at this point. I will keep on trying to find a way, but some advise would be very helpful. If something is missing, that you, just tell me and I will correct the code.

Thank you,