Principal Component Analysis for Point Spread Function Measurement and Implementation of different scaling in each dimension

I made 3D measurements of beads (small spheres) with a microscope for getting the point spread function (psf). This can be approximated with a 3D-Gauß. I need the 3 main axes for this ellipsoid. So principal component analysis (pca) seems right for this task. The goal is to fit three 1D Gaussians through this axes.

What I have right now: Produce Sample Data:

%%produce sample data with 3D Gauß

data=zeros(200,200,200);
a=5;
b=2;
x0=100;
sigmax=1;
y0=100;
sigmay=5;
z0=100;
sigmaz=8;
alpha=1;
beta=0;
gamma=0;

for z=1:200
    for x=1:200
        for y=1:200

              data(x,y,z)= a*exp(-...
                  (...
                  ((-(z-z0)*sin(beta)+cos(beta)*((x-x0)*cos(gamma)-(y-y0)*sin(gamma))))^2/(2*sigmax^2)+...
                  ((cos(alpha)*((y-y0)*cos(gamma)+(x-x0)*sin(gamma))-sin(alpha)*((z-z0)*cos(beta)+sin(beta)*((x-x0)*cos(gamma)-(y-y0)*sin(gamma)))))^2/(2*sigmay^2)+...
                  ((sin(alpha)*((y-y0)*cos(gamma)+(x-x0)*sin(gamma))+cos(alpha)*((z-z0)*cos(beta)+sin(beta)*((x-x0)*cos(gamma)-(y-y0)*sin(gamma)))))^2/(2*sigmaz^2)...
                  ))+b;

          end
      end
  end
meanBead=data;

This produces a 3D-Gaussian which is rotated a bit. The pixel size in each direction can be different so we have to take the scaling into account.

metadata.scale.x=1;
metadata.scale.y=1;
metadata.scale.z=1;

Prepare the data for pca:

%get threshold of meanBead
[out, meanThr]=threshold(meanBead,'triangle');

%preallocate memory
PCAData=zeros(size(meanBead,1)*size(meanBead,2)*size(meanBead,3),4);

[xind, yind, zind]= ind2sub(size(data),1:numel(data));
PCAData(:,1:3) = cat(1,xind,yind,zind)';
PCAData(:,4) = meanBead(:);

%Delete rows with NaN in it
PCAData(isnan(PCAData(:,4)), :) = [];

%weigth the xyz values by intensity and scaling
PCADataWeights=PCAData(:,4);

PCAData=[PCAData(:,1)*metadata.scale.x, PCAData(:,2)*metadata.scale.y, PCAData(:,3)*metadata.scale.z];

%substract meanvalue of coordinates
PCAData = bsxfun(@minus,PCAData,mean(PCAData,1));

%calculate covariance matrix
[U,~,PCAVar]=pca(PCAData,'Weights',PCADataWeights);

The three columns of U are the directions with the most variance. They are orthogonal and normed. This script works totally fine, until the scaling changes slightly. If we change one value from 1 to 0.99 the vectors of U will get the normal axis (0,0,1) (0,1,0).. and so on.

The question now is, how do I need to take the scaling into account and implement it the right way here?

Thank you very much

[Image Analyst]

Can you post your image? The profile of your beads will be the PSF of the system convolved with a sombrero function, unless your beads were truly just a single pixel.

[Florian Vollrath]

yes, sure they look like a sombrero function, but the wiggelings' intensity is so small that it should fit a gaussian pretty well. Here are line treces through the measurement with 1D Gaussians. My problem is just the scaling through the traces. My code is the following:

%%linetrace through the pca lines
%number of interpolated points in the trace (whole image)
steps=100;
intVal=zeros(steps,6);
for j=1:3
t_min=-[cog(1)/U(1,j), cog(2)/U(2,j), cog(3)/U(3,j)];
t_max=[(size(meanBead,1)-cog(1)/U(1,j)), (size(meanBead,2)-cog(2)/U(2,j)), (size(meanBead,3)-cog(3)/U(3,j))];
%[range, minIdx]=min(abs(t_max-t_min));
[~, minIdx]=max(U(:,j));
range=abs(min([t_min(minIdx),t_max(minIdx)])-max([t_min(minIdx),t_max(minIdx)]));
stepsz=range/steps;
%adjust length and save it in intVal(:,->)
length=sqrt((scales(1)*U(1,j))^2+(scales(2)*U(2,j))^2+(scales(3)*U(3,j))^2);
% length=sqrt((U(1,j))^2+(U(2,j))^2+(U(3,j))^2);
k=1;
for i=min([t_min(minIdx),t_max(minIdx)]):stepsz:max([t_min(minIdx),t_max(minIdx)])
    interpolPos = U(:,j)*i+cog';
    intVal(k,2*j)=interp3(meanBead,interpolPos(2),interpolPos(1),interpolPos(3),'cubic');
    intVal(k,1+2*(j-1))=length*i;
    k=k+1;
end
end
%%fit 1D Gaussian
sigmasTH=[res_lat,res_lat,res_ax]*1.2/1000;
sigmas=zeros(1,3);
figure( 'Name', 'untitled fit 1' );
for i=1:3
% Set up fittype and options.
[xData, yData] = prepareCurveData( intVal(:,1+2*(i-1)), intVal(:,2*i) );
ft = fittype( 'a*exp(-(x-mu)^2/2/sigma^2)+c', 'independent', 'x', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.StartPoint = [nanmax(yData) nanmean(meanBead(:)) 0 sigmasTH(i)];
[fitresult, gof] = fit( xData, yData, ft, opts );

% Plot fit with data.

subplot(2,2,i);
h = plot( fitresult, xData, yData );
legend( h, 'y1 vs. x1', 'untitled fit 1', 'Location', 'NorthEast' );
% Label axes
xlabel x1
ylabel y1
grid on

%save sigma
sigmas(i)=abs(fitresult.sigma);
end
FWHMs=sigmas*2*sqrt(2*log(2));

[Florian Vollrath]

here is the mean measurement of one bead (scaling is (0.176,0.176,1.303)µm and three gaussian fits through lines