How to calculate image resolution (full-widt​h-at-half-​maximum) of a point source?

Dear Matlab Expert,

I have images obtained from Monte Carlo simulations of point source. I want to calculate the images resolution at full-width-at-maximum.

The image (picture and image.mat file) and code used is given below:

 V=max(max(img));
[M,N]=find(img==V);


[fwhmx, fwhmy, fwhm] =FWHM_calculation(img,N, M,pixelSize);


function [fwhm_x,fwhm_y,fwhmT]=FWHM_calculation(Image,cx,cy,pixelSize)

% this function is used to calculate the full width half maximum of image
% image: the input image, including a point source or square source in the
% center
% fwhm : the output value for fwhm
CC=ceil(cx);
CR=ceil(cy);
 N=25;
Column1=CC-N:CC+N ;
 Row1=CR-N:CR+N;
  ROIimage=Image(Row1,Column1); 
  profile_x=sum(ROIimage,1);
  X=1:1:2*N+1;

  V=max(profile_x);
	plot(X,profile_x,'r-x','LineWidth',0.7,'MarkerSize',5) 
   	hold on
  
  [fitresult1, gof1] =createFitSPECTfwhm2(X,profile_x,V,N);
  rmse=gof1.rmse;

  fwhm_x=2*sqrt(log(2))*pixelSize*fitresult1.c1
  
  profile_y=sum(ROIimage,2);
  plot(X,profile_y,'b-o','LineWidth',0.5,'MarkerSize',3)
  V1=max(profile_y);
  
  [fitresult2, gof2] =createFitSPECTfwhm2(X,profile_y,V1,N);

  fwhm_y=2*sqrt(log(2))*pixelSize*fitresult2.c1
  
  profile_yt=profile_y';
  profile=(profile_x+ profile_yt)/2;

   plot(X,profile,'k-*','LineWidth',0.7,'MarkerSize',5)
  V=max(profile);
  [fitresult2, gof2] =createFitSPECTfwhm2(X,profile,V,N);

  fwhmT=2*sqrt(log(2))*pixelSize*fitresult2.c1
  xlabel('Pixel Index')
 ylabel('Counts')
 title('Profile along Yaxis','fontSize',9)
 legend('Profile along X-axis','Profile along Y-axis','Average Profile')
  


function [fitresult, gof] = createFitSPECTfwhm2(X, profile,maxValue,Xcenter)
%CREATEFIT(X,PROFILE_X)
%  Create a fit.
%
%  Data for 'untitled fit 1' fit:
%      X Input : X
%      Y Output: profile_x
%  Output:
%      fitresult : a fit object representing the fit.
%      gof : structure with goodness-of fit info.
%
%  See also FIT, CFIT, SFIT.

%  Auto-generated by MATLAB on 26-Mar-2014 11:10:32

%% Fit: 'untitled fit 1'.

[xData, yData] = prepareCurveData( X, profile );

% Set up fittype and options.
ft = fittype( 'gauss1' );
% ft = fittype( 'gauss2' );
opts = fitoptions( ft );
opts.Display = 'Off';
% opts.Lower = [-Inf -Inf 0];
opts.Lower = [-Inf -Inf 0];
opts.StartPoint = [maxValue Xcenter 1];
opts.Upper = [Inf Inf Inf];
% opts.Upper = [Inf Inf Inf];

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
% Plot fit with data.
figure( 'Name', 'Gaussian fit 1' );
h = plot( fitresult, xData, yData );
legend( h, 'profile_x vs. X', 'Gaussian fit 1', 'Location', 'NorthEast' );
% Label axes
xlabel( 'X' );
ylabel( 'profile_x' );

Answer 1

[Bjorn Gustavsson]

I'd start with some simple models for your image response (2-D Gaussian, sinc^2(x).*sinc(y)^2) and fitted the parameters of those to fit your psf-image. Something like this:

fsinc2 = @(pars,x,y) pars(1)*sinc((x-pars(2))/pars(3)).^2.*sinc((y-pars(4))/pars(5)).^2;
fG2 = @(pars,x,y) pars(1)*exp(-(x-pars(2)).^2/pars(3)^2-(y-pars(4)).^2/pars(5));
err_fcn = @(pars,I,x,y,fcn) sum(sum((I-(fcn(pars,x,y)+pars(6))).^2));
pars0 = [2.9250e5 182 2 182 2 min(t(:))];                               
[x,y] = meshgrid(1:362);
parsG2 = fminsearch(@(pars) err_fcn(pars,t,x,y,fG2),pars0);
parsS2 = fminsearch(@(pars) err_fcn(pars,t,x,y,fsinc2),pars0);
disp(parsG2(2:end))
disp(parsS2(2:end))
subplot(2,3,1)
imagesc(t)
subplot(2,3,2)
imagesc(fG2(parsG2,x,y))
subplot(2,3,4)
subplot(2,3,5)
imagesc(t-fG2(parsG2,x,y))
subplot(2,3,3)
imagesc(fsinc2(parsG2,x,y)) 
imagesc(fsinc2(parsG2,x,y))
subplot(2,3,6)
imagesc(t-fsinc2(parsS2,x,y)) 
for i1 = 1:6,
  phs(i1) = subplot(2,3,i1);
end
linkaxes(phs)
% Zoom to your hearts content

After that you might want to spice up your model-function - the sinc^2^2 didn't do as well as I'd hope so something more is going on than a simple square apperture...

HTH

[Khalid Hussain]

Thank you @Bjorn Gustavsson for the comment. I try your piece of code but does not work.

May I request to @Image Analyst or others users to please guide me in this probelm

[Bjorn Gustavsson]

Which step does not work? It works perfectly fine on my computer.

[Khalid Hussain]

thank you dear @Bjorn Gustavsson,

I got this error of undefined t,

[Bjorn Gustavsson]

That's the name of the image in the file you attached.

whos -file img.mat
  Name        Size               Bytes  Class     Attributes

  t         362x362            1048352  double              

[Khalid Hussain]

Dear @Bjorn Gustavsson thank you for your help. I got output the below output. What will be value of fwhm (image resolution) from below values in the image? And secondly, what about the pixel size, as we have pixel size for each simulated images which may differ for others and probably it should be used while calculating the fwhm values.

Please see the 4th plot which is empty and Please check which value we can report as the fwhm of the given image.

Thank you very much

[Bjorn Gustavsson]

Look at the function specifications.

fsinc2 = @(pars,x,y) pars(1)*sinc((x-pars(2))/pars(3)).^2.*sinc((y-pars(4))/pars(5)).^2;
fG2 = @(pars,x,y) pars(1)*exp(-(x-pars(2)).^2/pars(3)^2-(y-pars(4)).^2/pars(5));

fsinc2 is the outer product of a sinc-square in x and one in y directions. For a sinc-square the FWHM is 2z for z such that sinc(z)^2-1/2 == 0. I get z to be about 0.4429. That should give you enough to figure out FWHM in the x and y directions, just use pars(3) and pars(5) to scale z. You will see that that gives you different values in these directions. For the 2-D Gaussian you should be able to figure it out once you look at the functions.

HTH

[Khalid Hussain]

Dear @Bjorn Gustavsson , sorry for too many questions,

Can you please eloborate the functions with comments or break the main functions in sub part so that I can understand? it will be your great favor,as I feel difficult to understand these compact functions. Thank you so much.

pars0 = [2.9250e5 182 2 182 2], second and fourth entry seems to be near half the image size while 3rd and 5th entries looks to be the dimensions of the image. But how can we get the first value of 2.9250e5? Please comment

How can I calculate the pars(1) ... pars5() or how function return these values

Please try some easy way so that my problem will be solved.

[Bjorn Gustavsson]

OK, I'll break it down for you.

fsinc2 = @(pars,x,y) pars(1)*sinc((x-pars(2))/pars(3)).^2.*sinc((y-pars(4))/pars(5)).^2;

In mathematical notation:

that is a 2-D sinc^2 centred at x_0 and y_0, where is pars(1), is pars(2), is pars(3), is pars(4), and is pars(5). The initial guess for pars(1) is obtained from the maximum of t, while the centre-point is rather close to the mid-point of the image and the default guess for the widths can be set to something suitably close to the spot-size.

fG2 = @(pars,x,y) pars(1)*exp(-(x-pars(2)).^2/pars(3)^2-(y-pars(4)).^2/pars(5));

Same but for a 2-D Gaussian, centred at x_0 and y_0

err_fcn = @(pars,I,x,y,fcn) sum(sum((I-fcn(pars,x,y)).^2));

err_fcn is just the squared sum of the differences between an image I and the values of the function fcn at the same locations. This can be minimized to find the best function parameters (I_max, x_0, dx, y_0 and dy)

pars0 = [2.9250e5 182 2 182 2];

Initial guess for the parameters mentioned above, same for both fsinc2 and fG2 (with different meanings with respect to FWHM since the functions are different...)

[x,y] = meshgrid(1:362);

2-D grid coordinates

parsG2 = fminsearch(@(pars) err_fcn(pars,t,x,y,fG2),pars0);
parsS2 = fminsearch(@(pars) err_fcn(pars,t,x,y,fsinc2),pars0);

Parameter optimization returning the best function-parameters in the parsG2 and parsS2 variables. Some care has to be taken to make sure that the search actually have reached a suitable target.

disp(parsG2(2:end))
disp(parsS2(2:end))
subplot(2,3,1)
imagesc(t)
subplot(2,3,2)
imagesc(fG2(parsG2,x,y))
subplot(2,3,4)
subplot(2,3,5)
imagesc(t-fG2(parsG2,x,y))
subplot(2,3,3)
imagesc(fsinc2(parsG2,x,y)) 
imagesc(fsinc2(parsG2,x,y))
subplot(2,3,6)
imagesc(t-fsinc2(parsS2,x,y)) 
for i1 = 1:6,
  phs(i1) = subplot(2,3,i1);
end
linkaxes(phs)
% Zoom to your hearts content

[Khalid Hussain]

Dear @Bjorn Gustavsson thank you so much for the detailed explation.

With due respect a little more help, I am still confused about the pixelSize, and fwhm.

Can you please calculate the fwhm for pixel size of 1.11mm and 1.00mm the image attached in this question? It will help me to calculate fwhm for the rest of the simulated images.

Secondly, Can we also get a graph of guassian or sinc fit function, something similar to this below graph.

Thank you much, really appreiate your cooperation.

[Bjorn Gustavsson]

Well the FWHM should be aroud 4 pixels, for the Gaussian function as defined above It should be something like this:

exp(-(dxFWHM/2)^2/parsG2(3)^2) = 1/2;

That should give you a simple enough equation to solve for dxFWHM, then you can scale it any which way you want as you see fit.

Just plot the function the same way as you plot the Gaussian. Something like this:

plot(1:size(t,2),fsinc2(parsS2,x(182,:),y(182,:)),'r')

also have now added a parameter for the back-ground intensity to the fitting.

[Khalid Hussain]

I solved the equation for fwhm and it returns a complex value, something like in below picture,. I am still unable to calculate the FWHM value. Please also check 1D gaussian, somthing like I shared at the start of question. Thank you

[Bjorn Gustavsson]

The parameters you get out of the fitting is:

parsG2(1) - peak intensity

parsG2(2) - centre-point in x-direction, in pixels

parsG2(3) - width-parameter in the x-direction for the Gaussian as specified above, in pixels

parsG2(4) - centre-point in y-direction, in pixels

parsG2(5) - width-parameter in y-direction of the Gaussian as specified above in pixels

This should give you the Gaussian FWHM in the x-direction as:

FWHM = 2*parsG2(3)*sqrt(log(2)); % In pixels.

You made a sign error when deriving your solution, the '-' should be inside the sqrt-function-call.

Then if you want to scale that to physical length you can do so by multiplying that FWHM-length in the unit of [pixels] with the something that has the dimension of [length/pixel] to arrive at something with dimension of [length]. When you look at the dimension of your FWHM you will see that you have [pixels]^2/[length/pixel] -> [pixels^3]/[length] - which is not right.

[Khalid Hussain]

Thank you for your detailed explaination.

I just multipled with pixelSize like FWHM = pixelSize*2*parsG2(3)*sqrt(log(2)), now it seems to be fine but for the unit is [pixel * Length] which may be not appropriate as it should be [Length]. So, Please check this, if these units can be okay.

For my other simulated images, I replaced 2.9250e5 with (max(max(img)) pars0 = [2.9250e5 182 2 182 2]; It works fine but for some images it gives error of Maximum number of function evaluations has been exceeded. If I increase the max funciton value to the suggested current function value then it also works fine.

Can you please guide how we can calculate the maximum number of evaluation, so that it works for all images and generate peoper Fit plot?

Thank you so much