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

조회 수: 13(최근 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일
편집: Bjorn Gustavsson 2021년 6월 16일
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