Hi Yona,
I understand that you want to plot intensity along “x”, “y”, and “z” from a 3D Point Spread Function (PSF) stored in a .mat file and fit a Gaussian to the “z-axis” profile but the issue that you are encountering is that “max(data, [], 3)” loses depth information, making it difficult to extract and analyze the z-intensity profile correctly.
The following approach can be taken to solve this:
- Use “squeeze(data(y_center, :, z_slice))” for “x”, “squeeze(data(:, x_center, z_slice))” for “y”, and “squeeze(data(y_center, x_center, :))” for “z”.
- Then, fit a Gaussian function to the z-intensity profile using “lsqcurvefit” function, making sure that “z_range” and “intensity_z” are column vectors to match dimensions.
- Finally, plot the fitted curve to analyze the spread along “z”.
You may refer to the following corrected code:
% Load your data
% Get intensity projections
im_xy = max(data, [], 3); % Max intensity projection along Z
imshow(im_xy, []); % Display XY intensity projection
title('Max Projection (XY plane)');
% Define range
x_range = 1:size(data,2);
y_range = 1:size(data,1);
z_range = (1:size(data,3))'; % Ensure it's a column vector
% Select center for intensity profiles
x_center = round(size(data,2) / 2);
y_center = round(size(data,1) / 2);
% Extract intensity profiles
intensity_x = squeeze(data(y_center, :, round(size(data,3)/2))); % Along X
intensity_y = squeeze(data(:, x_center, round(size(data,3)/2))); % Along Y
intensity_z = squeeze(data(y_center, x_center, :)); % Along Z
% Plot intensity profiles
figure;
subplot(3,1,1);
plot(x_range, intensity_x, 'r', 'LineWidth', 1.5);
title('Intensity Profile along X');
xlabel('X'); ylabel('Intensity');
subplot(3,1,2);
plot(y_range, intensity_y, 'g', 'LineWidth', 1.5);
title('Intensity Profile along Y');
xlabel('Y'); ylabel('Intensity');
subplot(3,1,3);
plot(z_range, intensity_z, 'b', 'LineWidth', 1.5);
title('Intensity Profile along Z');
xlabel('Z'); ylabel('Intensity');
% Gaussian Fit for Z-axis intensity
gaussFit = @(p, x) p(1) * exp(-((x - p(2)).^2) / (2 * p(3)^2)) + p(4);
p0 = [max(intensity_z), mean(z_range), std(z_range), min(intensity_z)]; % Initial params
% Ensure input vectors are column vectors
z_range = z_range(:);
intensity_z = intensity_z(:);
p_opt = lsqcurvefit(gaussFit, p0, z_range, intensity_z);
% Plot fitted Gaussian
hold on;
plot(z_range, gaussFit(p_opt, z_range), 'k--', 'LineWidth', 2);
legend('Data', 'Gaussian Fit');
disp('Gaussian Fit Parameters (Z-axis):');
disp(p_opt);
Kindly refer to the following documentation for more details on the "lsqcurvefit" function in MATLAB:
Cheers & Happy Coding!
