MATLAB Answers

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

17 views (last 30 days)
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' );

Accepted Answer

Bjorn Gustavsson
Bjorn Gustavsson on 14 Jun 2021
Edited: Bjorn Gustavsson on 16 Jun 2021
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
  13 Comments
Khalid Hussain
Khalid Hussain on 17 Jun 2021
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

Sign in to comment.

More Answers (0)

Products


Release

R2018a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by