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

14 ビュー (過去 30 日間)
Khalid Hussain 2021 年 6 月 14 日
コメント済み: Khalid Hussain 2021 年 6 月 17 日
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.
%
% 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' );

サインインしてコメントする。

採用された回答

Bjorn Gustavsson 2021 年 6 月 14 日

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);
disp(parsG2(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)
for i1 = 1:6,
phs(i1) = subplot(2,3,i1);
end
% 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 件のコメント表示非表示 12 件の古いコメント
Khalid Hussain 2021 年 6 月 17 日
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

サインインしてコメントする。

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by