Determine point spread function of a single pixel image

Attached is an image of a single pixel lit-up on a micro-display. I also attached an image of the intensity line scan of the pixel, from which I can determine the FWHM. How do I determine the PSF from these images? Is there a MATLAB function that will allow me to do this. Thank you.

Answer 1

[Bjorn Gustavsson]

Your question is a bit confusing since you have two different pixel-domains where a PSF could be estimated - the micro-display and the image. For the micro-display it surely is just that single pixel (whose intensity might vary across its area but that seems pretty uninteresting with respect to PSF considerations - that would give you a PSF with sub-pixel width that still cannot be used to anything since the resolution of what can be displayed is limited to the pixel-size). For the image you would need to take the image when the size of the micro-display pixel would be less than a pixel wide (simple enough to estimate from the micro-display pixel-separation and the image field-of-view and the image resolution), here it appears to me that the micro-display pixel covers more than one pixel. That makes the problem challenging to solve - you effectively have 2 unknown spatial distributions first the PSF of your optical system and then the psf of the micro-display pixel but only one image with information about their convolution to work with.

Then when you have an image of the illuminating pixel that should be less than an image-pixel wide you can do this task by simple 2-D Gaussian fitting:

twoD_Gaussian = @(pars,x,y) pars(1)*exp(-(pars(2)-x).^2/pars(4)^2 - (pars(3)-y).^2/pars(5)^2)+pars(6);
errfcn = @(pars,x,y,img_reg,fitfcn) sum(sum((img_reg-fitfcn(pars,x,y)).^2));
img_reg = double(your_img(y_centre+[-30:30],x_centre+[-30:30]);
pars0 = [max(your_img(:)),x_centre,y_centre,5,5,0];
[x,y] = meshgrid(x_centre+[-30:30],y_centre+[-30:30]);
parsG = fminsearch(@(pars) errfcn(pars,x,y,img_reg,twoD_Gaussian),pars0);
sublpot(2,2,1),imagesc(img_reg-twoD_Gaussian(parsG,x,y))
subplot(2,2,2),imagesc(img_reg),hold on, contour(twoD_Gaussian(parsG,x,y),8,'k')

The twoD_Gaussian is just a function for a 2-D Gaussian with an additional background, the errfcn is just a general sum-of-squares function that takes an image-region a 2-D model-function and the inputs to that model-function as agruments. Then we cut out your region-of-interest, making sure that we get a double-array out, it should be a grayscale image first. Then set some reasonable start-guess parameters and generate a pair of coordinate matrices. Finaly we search for the best fitting parameters and check that the resulting fit is good. This would then give you an estimate of a PSF as a 2-D Gaussian. One can make more intricate PSF-model-functions (based on Bessel-functions, of Voigt-line-profiles etc), sometimes this is necessary especially if your imaging system has abberration-limited focus, but this is typically a good enough first start stab at this.

HTH

[LM]

Thank you for your answer. Ultimately I would like to deblur an image using the PSF. An example of a blurry Sieman's star image is attached (TOFA-003...jpg). I've used the script below that uses deconvlucy to sharpen the image. It's clearly not correct since I generate a black image. Could I please get your help with understanding how to use the deconvlucy function for this application? Thank you.

I = imread('TOFA-003 UV Cured - Pinwheel - Corrected.bmp');
FWHM = 10; % 10 pixels wide
sigma = (FWHM/2.35482);  %convert FWHM to sigma
B = imgaussfilt(I,sigma);
luc1 = uint8(B);
luc2 = deconvlucy(I,luc1,10);
imshow(luc2)

[Bjorn Gustavsson]

The second argument to deconvlucy should be the PSF (you send in a blurred/filtered version of the orignial image). Something like this:

Iz = zeros(31);            % This generates a
Iz(16,16) = 1;             % 31-by-31 
B = imgaussfilt(Iz,sigma); % Gaussian filter kernel,
sum(B(:))                  % check that it adds up to 1
luc2 = deconvlucy(I,B,1);  % First a single iteration...

Also have a look at deconvblind - it also generates updated estimates of the PSF, which can be nice if the initial psf-estimate is off.

[LM]

Thank you again for your help. I generated two images. The first is the original (TOFA-003 UV Cured...Corrected). The second is after applying deconvlucy (TOFA-003_deconvlucy). They look somewhat different but it is difficult to tell.

Do I use the MATLAB function Y = fft(X) to determine whether the deconvlucy image is sharper?

[Bjorn Gustavsson]

I'd typically do something like this (freehand-coded):

stopstopstop = 0;
sph = subplot(2,2,1);
imagesc(I)
sph(2) = subplot(2,2,2);
imagesc(luc2)
sph(3) = subplot(2,2,3);
imagesc(luc2(:,:,2)-I(:,:,2)) % just the green channel
sph(4) = subplot(2,2,4);
imagesc(luc2(:,:,2)./I(:,:,2))
linkaxes(sph)
zoom on
disp('Zoom to ROI, then push any key to continue...')
pause
disp('push left mouse-button to quit')
while stopstopstop ~= 1
  [x,y,stopstopstop] = ginput(1);
  if stopstopstop ~= 1
    luc2 = deconvlucy(luc2,PSF,1);
    axes(sph(2))
    imagesc(luc2)
    axes(sph(3))
    imagesc(luc2(:,:,2)-I(:,:,2))
    axes(sph(4))
    imagesc(luc2(:,:,2)./I(:,:,2))
    disp('Zoom around, then push any key to continue')
    pause
    disp('push left mouse-button to end this procedure')
  end
end

This should be some kind of start for a doodling-around with deconvolution on your images - but keep in mind that deconvolutions are inherentley high-frequency amplification operations and there is always a tendency of noise-amplification and ringing that comes with this. In my experience it has been possible to reduce the convolution-effect of about "half the PSF width", but this is something that obviously will vary from scene to scene and image to image.

[LM]

Thank you again for your comments.

Answer 2

[Matt J]

If you intend to model the PSF as a 2D Gaussian blob and you have the Optimization Toolbox then you can use this ready-made function for Gaussian fitting

https://www.mathworks.com/matlabcentral/fileexchange/69116-gaussfitn

[LM]

Thank you.