Czy istnieje sposób na wykrycie, czy obraz jest nieostry?

Zastanawiałem się, czy istnieje sposób, aby ustalić, czy obraz jest rozmyty, czy nie, analizując dane obrazu.

Tak to jest. Oblicz szybką transformatę Fouriera i przeanalizuj wynik. Transformacja Fouriera informuje, jakie częstotliwości występują na obrazie. Jeśli występuje duża liczba wysokich częstotliwości, wówczas obraz jest rozmazany.

Zdefiniowanie terminów „niski” i „wysoki” zależy od Ciebie.

Edytuj :

Jak stwierdzono w komentarzach, jeśli chcesz pojedynczego elementu zmiennoprzecinkowego reprezentującego rozmycie danego obrazu, musisz opracować odpowiednią metrykę.

Odpowiedź nikie dostarcza takiej miary. Przekształć obraz w jądro Laplaciana:

   1
1 -4  1
   1

I użyj solidnej maksymalnej metryki na wyjściu, aby uzyskać liczbę, której możesz użyć do progowania. Staraj się unikać nadmiernego wygładzania obrazów przed obliczeniem Laplaciana, ponieważ dowiesz się tylko, że wygładzony obraz jest rzeczywiście rozmazany :-).

Innym bardzo prostym sposobem oszacowania ostrości obrazu jest użycie filtra Laplace'a (lub LoG) i po prostu wybranie maksymalnej wartości. Zastosowanie solidnej miary, takiej jak kwantyl 99,9%, jest prawdopodobnie lepsze, jeśli spodziewasz się szumu (tj. Wybranie N-tego najwyższego kontrastu zamiast najwyższego kontrastu.) Jeśli spodziewasz się różnej jasności obrazu, powinieneś również dołączyć etap wstępnego przetwarzania w celu normalizacji jasności obrazu / kontrast (np. wyrównanie histogramu).

Zaimplementowałem sugestię Simona i tę w Mathematica i wypróbowałem ją na kilku obrazach testowych:

Pierwszy test rozmywa obrazy testowe przy użyciu filtra Gaussa o różnej wielkości jądra, a następnie oblicza FFT zamazanego obrazu i przyjmuje średnią z 90% najwyższych częstotliwości:

testFft[img_] := Table[
  (
   blurred = GaussianFilter[img, r];
   fft = Fourier[ImageData[blurred]];
   {w, h} = Dimensions[fft];
   windowSize = Round[w/2.1];
   Mean[Flatten[(Abs[
       fft[[w/2 - windowSize ;; w/2 + windowSize, 
         h/2 - windowSize ;; h/2 + windowSize]]])]]
   ), {r, 0, 10, 0.5}]

Wynik w logarytmicznym wykresie:

5 linii reprezentuje 5 obrazów testowych, oś X reprezentuje promień filtra Gaussa. Wykresy maleją, więc FFT jest dobrym miernikiem ostrości.

Oto kod estymatora rozmycia „najwyższego LoG”: Po prostu stosuje filtr LoG i zwraca najjaśniejszy piksel w wyniku filtru:

testLaplacian[img_] := Table[
  (
   blurred = GaussianFilter[img, r];
   Max[Flatten[ImageData[LaplacianGaussianFilter[blurred, 1]]]];
   ), {r, 0, 10, 0.5}]

Wynik w logarytmicznym wykresie:

Rozpiętość dla niewyraźnych obrazów jest tutaj nieco lepsza (2,5 vs 3,3), głównie dlatego, że ta metoda wykorzystuje tylko najsilniejszy kontrast na obrazie, podczas gdy FFT jest zasadniczo wartością średnią dla całego obrazu. Funkcje również maleją szybciej, dlatego łatwiej ustawić próg „rozmycia”.

Podczas pracy z obiektywem z automatyczną regulacją ostrości natknąłem się na bardzo przydatny zestaw algorytmów do wykrywania ostrości obrazu . Jest zaimplementowany w MATLAB, ale większość funkcji jest dość łatwa do przeniesienia do OpenCV za pomocą filter2D .

Jest to w zasadzie implementacja ankietowa wielu algorytmów pomiaru ostrości. Jeśli chcesz przeczytać oryginalne artykuły, w kodzie znajdują się odniesienia do autorów algorytmów. Artykuł z 2012 r. Autorstwa Pertuza i in. Analiza operatorów miary ostrości dla kształtu z ogniska (SFF) daje świetny przegląd wszystkich tych miar, a także ich wydajności (zarówno pod względem szybkości, jak i dokładności w odniesieniu do SFF).

EDYCJA: Dodano kod MATLAB na wypadek śmierci łącza.

function FM = fmeasure(Image, Measure, ROI)
%This function measures the relative degree of focus of 
%an image. It may be invoked as:
%
%   FM = fmeasure(Image, Method, ROI)
%
%Where 
%   Image,  is a grayscale image and FM is the computed
%           focus value.
%   Method, is the focus measure algorithm as a string.
%           see 'operators.txt' for a list of focus 
%           measure methods. 
%   ROI,    Image ROI as a rectangle [xo yo width heigth].
%           if an empty argument is passed, the whole
%           image is processed.
%
%  Said Pertuz
%  Abr/2010


if ~isempty(ROI)
    Image = imcrop(Image, ROI);
end

WSize = 15; % Size of local window (only some operators)

switch upper(Measure)
    case 'ACMO' % Absolute Central Moment (Shirvaikar2004)
        if ~isinteger(Image), Image = im2uint8(Image);
        end
        FM = AcMomentum(Image);

    case 'BREN' % Brenner's (Santos97)
        [M N] = size(Image);
        DH = Image;
        DV = Image;
        DH(1:M-2,:) = diff(Image,2,1);
        DV(:,1:N-2) = diff(Image,2,2);
        FM = max(DH, DV);        
        FM = FM.^2;
        FM = mean2(FM);

    case 'CONT' % Image contrast (Nanda2001)
        ImContrast = inline('sum(abs(x(:)-x(5)))');
        FM = nlfilter(Image, [3 3], ImContrast);
        FM = mean2(FM);

    case 'CURV' % Image Curvature (Helmli2001)
        if ~isinteger(Image), Image = im2uint8(Image);
        end
        M1 = [-1 0 1;-1 0 1;-1 0 1];
        M2 = [1 0 1;1 0 1;1 0 1];
        P0 = imfilter(Image, M1, 'replicate', 'conv')/6;
        P1 = imfilter(Image, M1', 'replicate', 'conv')/6;
        P2 = 3*imfilter(Image, M2, 'replicate', 'conv')/10 ...
            -imfilter(Image, M2', 'replicate', 'conv')/5;
        P3 = -imfilter(Image, M2, 'replicate', 'conv')/5 ...
            +3*imfilter(Image, M2, 'replicate', 'conv')/10;
        FM = abs(P0) + abs(P1) + abs(P2) + abs(P3);
        FM = mean2(FM);

    case 'DCTE' % DCT energy ratio (Shen2006)
        FM = nlfilter(Image, [8 8], @DctRatio);
        FM = mean2(FM);

    case 'DCTR' % DCT reduced energy ratio (Lee2009)
        FM = nlfilter(Image, [8 8], @ReRatio);
        FM = mean2(FM);

    case 'GDER' % Gaussian derivative (Geusebroek2000)        
        N = floor(WSize/2);
        sig = N/2.5;
        [x,y] = meshgrid(-N:N, -N:N);
        G = exp(-(x.^2+y.^2)/(2*sig^2))/(2*pi*sig);
        Gx = -x.*G/(sig^2);Gx = Gx/sum(Gx(:));
        Gy = -y.*G/(sig^2);Gy = Gy/sum(Gy(:));
        Rx = imfilter(double(Image), Gx, 'conv', 'replicate');
        Ry = imfilter(double(Image), Gy, 'conv', 'replicate');
        FM = Rx.^2+Ry.^2;
        FM = mean2(FM);

    case 'GLVA' % Graylevel variance (Krotkov86)
        FM = std2(Image);

    case 'GLLV' %Graylevel local variance (Pech2000)        
        LVar = stdfilt(Image, ones(WSize,WSize)).^2;
        FM = std2(LVar)^2;

    case 'GLVN' % Normalized GLV (Santos97)
        FM = std2(Image)^2/mean2(Image);

    case 'GRAE' % Energy of gradient (Subbarao92a)
        Ix = Image;
        Iy = Image;
        Iy(1:end-1,:) = diff(Image, 1, 1);
        Ix(:,1:end-1) = diff(Image, 1, 2);
        FM = Ix.^2 + Iy.^2;
        FM = mean2(FM);

    case 'GRAT' % Thresholded gradient (Snatos97)
        Th = 0; %Threshold
        Ix = Image;
        Iy = Image;
        Iy(1:end-1,:) = diff(Image, 1, 1);
        Ix(:,1:end-1) = diff(Image, 1, 2);
        FM = max(abs(Ix), abs(Iy));
        FM(FM<Th)=0;
        FM = sum(FM(:))/sum(sum(FM~=0));

    case 'GRAS' % Squared gradient (Eskicioglu95)
        Ix = diff(Image, 1, 2);
        FM = Ix.^2;
        FM = mean2(FM);

    case 'HELM' %Helmli's mean method (Helmli2001)        
        MEANF = fspecial('average',[WSize WSize]);
        U = imfilter(Image, MEANF, 'replicate');
        R1 = U./Image;
        R1(Image==0)=1;
        index = (U>Image);
        FM = 1./R1;
        FM(index) = R1(index);
        FM = mean2(FM);

    case 'HISE' % Histogram entropy (Krotkov86)
        FM = entropy(Image);

    case 'HISR' % Histogram range (Firestone91)
        FM = max(Image(:))-min(Image(:));


    case 'LAPE' % Energy of laplacian (Subbarao92a)
        LAP = fspecial('laplacian');
        FM = imfilter(Image, LAP, 'replicate', 'conv');
        FM = mean2(FM.^2);

    case 'LAPM' % Modified Laplacian (Nayar89)
        M = [-1 2 -1];        
        Lx = imfilter(Image, M, 'replicate', 'conv');
        Ly = imfilter(Image, M', 'replicate', 'conv');
        FM = abs(Lx) + abs(Ly);
        FM = mean2(FM);

    case 'LAPV' % Variance of laplacian (Pech2000)
        LAP = fspecial('laplacian');
        ILAP = imfilter(Image, LAP, 'replicate', 'conv');
        FM = std2(ILAP)^2;

    case 'LAPD' % Diagonal laplacian (Thelen2009)
        M1 = [-1 2 -1];
        M2 = [0 0 -1;0 2 0;-1 0 0]/sqrt(2);
        M3 = [-1 0 0;0 2 0;0 0 -1]/sqrt(2);
        F1 = imfilter(Image, M1, 'replicate', 'conv');
        F2 = imfilter(Image, M2, 'replicate', 'conv');
        F3 = imfilter(Image, M3, 'replicate', 'conv');
        F4 = imfilter(Image, M1', 'replicate', 'conv');
        FM = abs(F1) + abs(F2) + abs(F3) + abs(F4);
        FM = mean2(FM);

    case 'SFIL' %Steerable filters (Minhas2009)
        % Angles = [0 45 90 135 180 225 270 315];
        N = floor(WSize/2);
        sig = N/2.5;
        [x,y] = meshgrid(-N:N, -N:N);
        G = exp(-(x.^2+y.^2)/(2*sig^2))/(2*pi*sig);
        Gx = -x.*G/(sig^2);Gx = Gx/sum(Gx(:));
        Gy = -y.*G/(sig^2);Gy = Gy/sum(Gy(:));
        R(:,:,1) = imfilter(double(Image), Gx, 'conv', 'replicate');
        R(:,:,2) = imfilter(double(Image), Gy, 'conv', 'replicate');
        R(:,:,3) = cosd(45)*R(:,:,1)+sind(45)*R(:,:,2);
        R(:,:,4) = cosd(135)*R(:,:,1)+sind(135)*R(:,:,2);
        R(:,:,5) = cosd(180)*R(:,:,1)+sind(180)*R(:,:,2);
        R(:,:,6) = cosd(225)*R(:,:,1)+sind(225)*R(:,:,2);
        R(:,:,7) = cosd(270)*R(:,:,1)+sind(270)*R(:,:,2);
        R(:,:,7) = cosd(315)*R(:,:,1)+sind(315)*R(:,:,2);
        FM = max(R,[],3);
        FM = mean2(FM);

    case 'SFRQ' % Spatial frequency (Eskicioglu95)
        Ix = Image;
        Iy = Image;
        Ix(:,1:end-1) = diff(Image, 1, 2);
        Iy(1:end-1,:) = diff(Image, 1, 1);
        FM = mean2(sqrt(double(Iy.^2+Ix.^2)));

    case 'TENG'% Tenengrad (Krotkov86)
        Sx = fspecial('sobel');
        Gx = imfilter(double(Image), Sx, 'replicate', 'conv');
        Gy = imfilter(double(Image), Sx', 'replicate', 'conv');
        FM = Gx.^2 + Gy.^2;
        FM = mean2(FM);

    case 'TENV' % Tenengrad variance (Pech2000)
        Sx = fspecial('sobel');
        Gx = imfilter(double(Image), Sx, 'replicate', 'conv');
        Gy = imfilter(double(Image), Sx', 'replicate', 'conv');
        G = Gx.^2 + Gy.^2;
        FM = std2(G)^2;

    case 'VOLA' % Vollath's correlation (Santos97)
        Image = double(Image);
        I1 = Image; I1(1:end-1,:) = Image(2:end,:);
        I2 = Image; I2(1:end-2,:) = Image(3:end,:);
        Image = Image.*(I1-I2);
        FM = mean2(Image);

    case 'WAVS' %Sum of Wavelet coeffs (Yang2003)
        [C,S] = wavedec2(Image, 1, 'db6');
        H = wrcoef2('h', C, S, 'db6', 1);   
        V = wrcoef2('v', C, S, 'db6', 1);   
        D = wrcoef2('d', C, S, 'db6', 1);   
        FM = abs(H) + abs(V) + abs(D);
        FM = mean2(FM);

    case 'WAVV' %Variance of  Wav...(Yang2003)
        [C,S] = wavedec2(Image, 1, 'db6');
        H = abs(wrcoef2('h', C, S, 'db6', 1));
        V = abs(wrcoef2('v', C, S, 'db6', 1));
        D = abs(wrcoef2('d', C, S, 'db6', 1));
        FM = std2(H)^2+std2(V)+std2(D);

    case 'WAVR'
        [C,S] = wavedec2(Image, 3, 'db6');
        H = abs(wrcoef2('h', C, S, 'db6', 1));   
        V = abs(wrcoef2('v', C, S, 'db6', 1));   
        D = abs(wrcoef2('d', C, S, 'db6', 1)); 
        A1 = abs(wrcoef2('a', C, S, 'db6', 1));
        A2 = abs(wrcoef2('a', C, S, 'db6', 2));
        A3 = abs(wrcoef2('a', C, S, 'db6', 3));
        A = A1 + A2 + A3;
        WH = H.^2 + V.^2 + D.^2;
        WH = mean2(WH);
        WL = mean2(A);
        FM = WH/WL;
    otherwise
        error('Unknown measure %s',upper(Measure))
end
 end
%************************************************************************
function fm = AcMomentum(Image)
[M N] = size(Image);
Hist = imhist(Image)/(M*N);
Hist = abs((0:255)-255*mean2(Image))'.*Hist;
fm = sum(Hist);
end

%******************************************************************
function fm = DctRatio(M)
MT = dct2(M).^2;
fm = (sum(MT(:))-MT(1,1))/MT(1,1);
end

%************************************************************************
function fm = ReRatio(M)
M = dct2(M);
fm = (M(1,2)^2+M(1,3)^2+M(2,1)^2+M(2,2)^2+M(3,1)^2)/(M(1,1)^2);
end
%******************************************************************

Kilka przykładów wersji OpenCV:

// OpenCV port of 'LAPM' algorithm (Nayar89)
double modifiedLaplacian(const cv::Mat& src)
{
    cv::Mat M = (Mat_<double>(3, 1) << -1, 2, -1);
    cv::Mat G = cv::getGaussianKernel(3, -1, CV_64F);

    cv::Mat Lx;
    cv::sepFilter2D(src, Lx, CV_64F, M, G);

    cv::Mat Ly;
    cv::sepFilter2D(src, Ly, CV_64F, G, M);

    cv::Mat FM = cv::abs(Lx) + cv::abs(Ly);

    double focusMeasure = cv::mean(FM).val[0];
    return focusMeasure;
}

// OpenCV port of 'LAPV' algorithm (Pech2000)
double varianceOfLaplacian(const cv::Mat& src)
{
    cv::Mat lap;
    cv::Laplacian(src, lap, CV_64F);

    cv::Scalar mu, sigma;
    cv::meanStdDev(lap, mu, sigma);

    double focusMeasure = sigma.val[0]*sigma.val[0];
    return focusMeasure;
}

// OpenCV port of 'TENG' algorithm (Krotkov86)
double tenengrad(const cv::Mat& src, int ksize)
{
    cv::Mat Gx, Gy;
    cv::Sobel(src, Gx, CV_64F, 1, 0, ksize);
    cv::Sobel(src, Gy, CV_64F, 0, 1, ksize);

    cv::Mat FM = Gx.mul(Gx) + Gy.mul(Gy);

    double focusMeasure = cv::mean(FM).val[0];
    return focusMeasure;
}

// OpenCV port of 'GLVN' algorithm (Santos97)
double normalizedGraylevelVariance(const cv::Mat& src)
{
    cv::Scalar mu, sigma;
    cv::meanStdDev(src, mu, sigma);

    double focusMeasure = (sigma.val[0]*sigma.val[0]) / mu.val[0];
    return focusMeasure;
}

Nie ma gwarancji, czy te środki są najlepszym wyborem dla twojego problemu, ale jeśli wyśledzisz dokumenty związane z tymi środkami, mogą dać ci więcej wglądu. Mam nadzieję, że kod Ci się przyda! Wiem, że tak.

Opierając się na odpowiedzi Nike. Zaimplementowanie metody opartej na laplacianach przy pomocy opencv jest proste:

short GetSharpness(char* data, unsigned int width, unsigned int height)
{
    // assumes that your image is already in planner yuv or 8 bit greyscale
    IplImage* in = cvCreateImage(cvSize(width,height),IPL_DEPTH_8U,1);
    IplImage* out = cvCreateImage(cvSize(width,height),IPL_DEPTH_16S,1);
    memcpy(in->imageData,data,width*height);

    // aperture size of 1 corresponds to the correct matrix
    cvLaplace(in, out, 1);

    short maxLap = -32767;
    short* imgData = (short*)out->imageData;
    for(int i =0;i<(out->imageSize/2);i++)
    {
        if(imgData[i] > maxLap) maxLap = imgData[i];
    }

    cvReleaseImage(&in);
    cvReleaseImage(&out);
    return maxLap;
}

Zwróci krótki wskazujący maksymalną wykrytą ostrość, która w oparciu o moje testy na próbkach ze świata rzeczywistego jest całkiem dobrym wskaźnikiem tego, czy aparat jest ostry, czy nie. Nic dziwnego, że normalne wartości są zależne od sceny, ale znacznie mniej niż metoda FFT, która musi być zbyt wysoka na poziomie fałszywie dodatnim, aby była użyteczna w mojej aplikacji.

Wymyśliłem zupełnie inne rozwiązanie. Musiałem przeanalizować nieruchome klatki wideo, aby znaleźć najostrzejszą z każdej (X) klatki. W ten sposób wykryłbym rozmycie ruchu i / lub nieostre obrazy.

Skończyło się na wykrywaniu Canny Edge i uzyskałem BARDZO dobre wyniki z prawie każdym rodzajem wideo (metodą Nikie miałem problemy ze zdigitalizowanymi filmami VHS i ciężkimi filmami z przeplotem).

Zoptymalizowałem wydajność, ustawiając obszar zainteresowania (ROI) na oryginalnym obrazie.

Za pomocą EmguCV:

//Convert image using Canny
using (Image<Gray, byte> imgCanny = imgOrig.Canny(225, 175))
{
    //Count the number of pixel representing an edge
    int nCountCanny = imgCanny.CountNonzero()[0];

    //Compute a sharpness grade:
    //< 1.5 = blurred, in movement
    //de 1.5 à 6 = acceptable
    //> 6 =stable, sharp
    double dSharpness = (nCountCanny * 1000.0 / (imgCanny.Cols * imgCanny.Rows));
}

Dzięki nikie za wspaniałą sugestię Laplace'a. Dokumenty OpenCV skierowały mnie w tym samym kierunku: używając Pythona, CV2 (OpenCv 2.4.10) i Numpy ...

gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) numpy.max(cv2.convertScaleAbs(cv2.Laplacian(gray_image,3)))

wynik jest między 0-255. Odkryłem, że wszystko powyżej 200ish jest bardzo ostre, a po 100 jest wyraźnie rozmyte. maksimum nigdy nie osiąga dużo poniżej 20, nawet jeśli jest całkowicie rozmyte.

Jeden ze sposobów, którego obecnie używam, mierzy rozkład krawędzi na obrazie. Poszukaj tego papieru:

@ARTICLE{Marziliano04perceptualblur,
    author = {Pina Marziliano and Frederic Dufaux and Stefan Winkler and Touradj Ebrahimi},
    title = {Perceptual blur and ringing metrics: Application to JPEG2000,” Signal Process},
    journal = {Image Commun},
    year = {2004},
    pages = {163--172} }

Zwykle jest za zaporą, ale widziałem kilka darmowych kopii. Zasadniczo lokalizują pionowe krawędzie obrazu, a następnie mierzą szerokość tych krawędzi. Uśrednienie szerokości daje końcowy wynik oceny rozmycia obrazu. Szersze krawędzie odpowiadają rozmytym obrazom i odwrotnie.

Ten problem należy do dziedziny oceny jakości obrazu bez odniesienia . Jeśli spojrzysz na to w Google Scholar, otrzymasz mnóstwo przydatnych referencji.

EDYTOWAĆ

Oto wykres oszacowań rozmycia uzyskanych dla 5 zdjęć w poście nikie. Wyższe wartości odpowiadają większemu rozmyciu. Użyłem filtru Gaussa 11 x 11 o stałym rozmiarze i zmieniłem standardowe odchylenie (używając convertpolecenia imagemagick do uzyskania rozmytych obrazów).

Jeśli porównujesz obrazy o różnych rozmiarach, nie zapomnij o normalizacji według szerokości obrazu, ponieważ większe obrazy będą miały szersze krawędzie.

Wreszcie znaczącym problemem jest rozróżnienie między rozmyciem artystycznym a niepożądanym rozmyciem (spowodowanym brakiem ostrości, kompresją, względnym ruchem obiektu do aparatu), ale wykracza to poza proste podejścia, takie jak ten. Na przykład artystyczne rozmycie, spójrz na obraz Lenny: odbicie Lenny w lustrze jest rozmyte, ale jej twarz jest idealnie skupiona. Przyczynia się to do większego oszacowania rozmycia obrazu Lenna.

Wypróbowałem rozwiązanie oparte na filtrze Laplaciana z tego postu. To mi nie pomogło. Więc wypróbowałem rozwiązanie z tego postu i było dobre dla mojej sprawy (ale jest powolne):

import cv2

image = cv2.imread("test.jpeg")
height, width = image.shape[:2]
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

def px(x, y):
    return int(gray[y, x])

sum = 0
for x in range(width-1):
    for y in range(height):
        sum += abs(px(x, y) - px(x+1, y))

Mniej zamazany obraz ma maksymalną sumwartość!

Możesz także dostroić prędkość i dokładność, zmieniając krok, np

ta część

for x in range(width - 1):

możesz zastąpić tym

for x in range(0, width - 1, 10):

Powyższe odpowiedzi wyjaśniły wiele rzeczy, ale myślę, że użyteczne jest rozróżnienie pojęciowe.

Co zrobić, jeśli zrobisz idealnie rozogniskowane zdjęcie zamazanego obrazu?

Problem wykrywania rozmycia jest dobrze postawiony tylko wtedy, gdy masz referencję . Jeśli potrzebujesz zaprojektować np. System automatycznego ustawiania ostrości, porównujesz sekwencję zdjęć wykonanych przy różnym stopniu rozmycia lub wygładzenia i próbujesz znaleźć punkt minimalnego rozmycia w tym zestawie. Innymi słowy, musisz odnieść się do różnych obrazów przy użyciu jednej z technik zilustrowanych powyżej (w zasadzie - z różnymi możliwymi poziomami udoskonalenia w podejściu - szukając jednego obrazu o najwyższej zawartości wysokiej częstotliwości).

Kod Matlaba dwóch metod opublikowanych w bardzo cenionych czasopismach (Transakcje IEEE dotyczące przetwarzania obrazu) jest dostępny tutaj: https://ivulab.asu.edu/software

sprawdź algorytmy CPBDM i JNBM. Jeśli sprawdzisz kod, przeniesienie go nie będzie trudne, a przy okazji bazuje na metodzie Marzialiano jako podstawowej funkcji.

zaimplementowałem to użyj fft w Matlabie i sprawdź histogram średniej obliczeniowej FFT i STD, ale można również wykonać funkcję dopasowania

fa =  abs(fftshift(fft(sharp_img)));
fb = abs(fftshift(fft(blured_img)));

f1=20*log10(0.001+fa);
f2=20*log10(0.001+fb);

figure,imagesc(f1);title('org')
figure,imagesc(f2);title('blur')

figure,hist(f1(:),100);title('org')
figure,hist(f2(:),100);title('blur')

mf1=mean(f1(:));
mf2=mean(f2(:));

mfd1=median(f1(:));
mfd2=median(f2(:));

sf1=std(f1(:));
sf2=std(f2(:));

Tak robię w Opencv, aby wykryć jakość ogniskowania w regionie:

Mat grad;
int scale = 1;
int delta = 0;
int ddepth = CV_8U;
Mat grad_x, grad_y;
Mat abs_grad_x, abs_grad_y;
/// Gradient X
Sobel(matFromSensor, grad_x, ddepth, 1, 0, 3, scale, delta, BORDER_DEFAULT);
/// Gradient Y
Sobel(matFromSensor, grad_y, ddepth, 0, 1, 3, scale, delta, BORDER_DEFAULT);
convertScaleAbs(grad_x, abs_grad_x);
convertScaleAbs(grad_y, abs_grad_y);
addWeighted(abs_grad_x, 0.5, abs_grad_y, 0.5, 0, grad);
cv::Scalar mu, sigma;
cv::meanStdDev(grad, /* mean */ mu, /*stdev*/ sigma);
focusMeasure = mu.val[0] * mu.val[0];