17 views (last 30 days)

Show older comments

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' );

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

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

Start Hunting!