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

32 Ansichten (letzte 30 Tage)
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' );

Akzeptierte Antwort

Bjorn Gustavsson
Bjorn Gustavsson am 14 Jun. 2021
Bearbeitet: Bjorn Gustavsson am 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 Kommentare
Bjorn Gustavsson
Bjorn Gustavsson am 17 Jun. 2021
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
Khalid Hussain am 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

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Produkte


Version

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by