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