How can I plot an intensity profile along z axis in 3D data.

I have a 3D image data of mat. file. Those are pread point function of my optical setup. What I want to do is to plot the intensity over x, y and z and then fit them with a Gaussian fit. I have had a lot of problem so far.

What I do is :

load('C:\Users\Documents\data1015.mat');
im=max(data, [], 3);
imshow(im, [])

here I have the intensity information along x and y when I use the in the figure I can see the range where I need to plot my profile here I have the intensity information along x and y when I use the in the figure I can see the range where I need to plot my profile:

so I use the following code to plot instensity along x and y using in a specific range.

plot(im(:, :))
plot(im(: :))

How can I plot the intensity profeli along z? Well I know that the max operator looses the 3D information. Any help would be appreciated!

Answer 1

[sanidhyak]

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:

1. 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”.
2. 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.
3. 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:

• https://www.mathworks.com/help/optim/ug/lsqcurvefit.html

Cheers & Happy Coding!