Quantify sharpness using frequencies

Question

I have several images of the same objects that I would like to rank based on sharpness. Say, for instance, to rank the quality of focus.

From theory I know that when you make an image less sharp you loose high frequencies. To test this I converted three images to the frequency domain:

It is clear the blurred images have more uniforme higher frequencies due to the noise. But what would be a good way to quantify this? I tried to get the mean of pixels at a certain distance away from the brightest pixel, but this does not show much difference

How would I use this method to find the sharpest picture.

Here is the code I used to make the above images:

%make gaussian filter
S=10;P=1.5;
x = [-S/2:1:S/2];
x = exp(-x.^2/(2*P^2));
filt = x'*x;filt = filt/sum(filt(:));

%laplacian kernel
LapK = [0 1 0; 1 -4 1 ; 0 1 0];

%get demo image
Iorg = imread('peppers.png');
Iorg = mean(Iorg,3);
Iorg = Iorg./max(Iorg(:));

%mesh grid
[X,Y]=meshgrid(1:size(Iorg,2),1:size(Iorg,1));
R = sqrt(X.^2+Y.^2);
Rmax = min(size(Iorg)./2);

%Add noise and fourier transform
I = Iorg+randn(size(Iorg))*0.01;I = I./max(I(:));
L = conv2(I,LapK); L1 = max(L(:));
F = fft2(Iorg);
%Radial mean
RadM1=nan([1,Rmax]);
for ct = 1:Rmax
    RadM1(ct) = mean(abs(F(R<=ct&R>(ct-1))));
end
figure(1);clf;colormap jet
subplot(3,3,1);imagesc(repmat(Iorg,[1,1,3]));axis image off;title('image');
subplot(3,3,2);imagesc(log(abs(fftshift(F))),[-4 10]);axis image off;title('fourier transform');
subplot(3,3,3);plot(1:Rmax,log(RadM1),'.-');title('radial mean');
ylabel('log(|F|)');ylim([2 12])

%filter
I=conv2(Iorg, filt, 'same');
%Add noise and fourier transform
I = I+randn(size(Iorg))*0.01;I = I./max(I(:));
L = conv2(I,LapK); L2 = max(L(:));
F = fft2(I);
%Radial mean
RadM2=nan([1,Rmax]);
for ct = 1:Rmax
    RadM2(ct) = mean(abs(F(R<=ct&R>(ct-1))));
end
subplot(3,3,4);imagesc(repmat(I,[1,1,3]));axis image off;title('filtered once');
subplot(3,3,5);imagesc(log(abs(fftshift(F))),[-4 10]);axis image off;title('fourier transform');
subplot(3,3,6);plot(1:Rmax,log(RadM2),'.-');title('radial mean');
ylabel('log(|F|)');ylim([2 12])

%filter twice
I=conv2(Iorg, filt, 'same');I=conv2(I, filt, 'same');
%Add noise and fourier transform
I = I+randn(size(Iorg))*0.01; I = I./max(I(:));
L = conv2(I,LapK); L3 = max(L(:));
F = fft2(I);
%Radial mean
RadM3=nan([1,Rmax]);
for ct = 1:Rmax
    RadM3(ct) = mean(abs(F(R<=ct&R>(ct-1))));
end
subplot(3,3,7);imagesc(repmat(I,[1,1,3]));axis image off;title('filtered twice');
subplot(3,3,8);imagesc(log(abs(fftshift(F))),[-4 10]);axis image off;title('fourier transform');
subplot(3,3,9);plot(1:Rmax,log(RadM3),'.-');title('radial mean');
ylabel('log(|F|)');ylim([2 12])
xlabel('radius(pix)');
fprintf(1,'1: %.2f \n2: %.2f \n3: %.2f\n',L1,L2,L3)

result

1: 1.04

2: 0.35

3: 0.27

[edit] This is indeed a duplicate. Convolution with a laplacian kernel and looking for the maximum works well in this case. [/edit]