Estimate the "Blurry" distribution of an image

Question

I'm trying to analyze an image obtained by scanning electron microscope(SEM) where different parts of the same image has has different blurring radius. My goal is to determine how blurry different part of this image is.

Clearly, the central part of this image is clear and reveal a lot of details and the upper and down part is rather blurry.

The result of blurring radius needn't be accurate, but needs to show a trend that the central part is not blurry and the further from the central part, the blurrier the image will be.

Is this possible?

Thanks!

---

Some of my thoughts:

At the beginning I thought about MaxDetect method given here but later I realized this method may need a image with similar blurring radius everywhere, so this method probably won't work here.

---

A further goal

What if I would like to know its "approximate blurring radius" instead of just knowing how "blurry" it is?

Goal here is to approximately deternmine the blurring radius r in central part of this image (25%~60% counting from the bottom to the top) where the the object is only "slightly blurred" and still clearly visible.

So how to obtain the blurring radius r within an uncertainty of 30%?

---

Some other test images:

Answer 1

I think what you're looking is related to "scale selection" (Wikipedia link) in scale space theory. Simply put, the idea is: If you have an edge in your image that's been blurred with some filter size sigma, and you apply Gaussian derivative or Laplacian of Gaussian filters with varying sigmas to that image, you get the highest impulse response if the filter sizes match. (Think of it as template matching, although the theoretical reasoning is somewhat different.)

I'll illustrate my idea with a simple sample image, for example, a disk:

sample = N[DiskMatrix[32, 128]];

Now let's blur this disk with different (Gaussian) filters:

Table[
  (
   blurrySample = 
    GaussianFilter[sample, {3 s1, s1}, Method -> "Gaussian"]; (* continued *)

and apply a range of Laplacian of Gaussian filters to it:

   scaleSpaceStep = 1.1;
   scaleSpace = scaleSpaceStep^Range[Log[scaleSpaceStep, 40]];
   filter = 
    Table[s^2*
      Total[LaplacianGaussianFilter[blurrySample, {3 s, s}, 
         Method -> "Gaussian"]^2, ∞], {s, scaleSpace}]; (* continued *)

Note that I have to multiply the LoG filter with s^2 to make this work. This is a "normalization factor", and it depends on the image content, i.e. it's different for point-like, line-like or area-like features. We'll have to estimate this for your images. Let's look at the results:

   ListLinePlot[{scaleSpace, filter}\[Transpose], PlotRange -> All, 
    GridLines -> {{s1}, {}},
    Prolog -> {Inset[Image[blurrySample, ImageSize -> 40], 
       Scaled[{1, 1}], Scaled[{1, 1}]]}]
   ), {s1, 0, 29, 2.5}]

As you can see, the LoG filters have the strongest squared impulse response if the LoG filter's size (roughly) matches the size of the blurring filter.

Now let's try this on your images. First, I choose the LoG filter sizes I'll use:

scaleSpaceStep = 1.2;
scaleSpace = 
  scaleSpaceStep^
   Range[Round[Log[scaleSpaceStep, 1]], Log[scaleSpaceStep, 100]];

Then the scale selection is quite simple:

estimateScale[img_] := (
  (* apply LoG filters with different sizes *)
  log = Table[
    GaussianFilter[
     s^4 LaplacianGaussianFilter[ImageData[img], {3 s, s}, 
        Method -> "Gaussian"]^2, 50, Method -> "Gaussian"], {s, 
     scaleSpace}];
  (* and estimate the "best scale" from a weighted average *)
  perPixelMaxScale = (scaleSpace.log)/Total[log];
  (* fancy display stuff *)
  {
   Image[img, ImageSize -> 400],
   ArrayPlot[perPixelMaxScale, PlotLegends -> Automatic, 
    ColorFunction -> "ThermometerColors", ImageSize -> 400]
   }
  )

Applied to your images:

imgs = ColorConvert[Import[#], 
    "Grayscale"] & /@ {"https://i.stack.imgur.com/dNxZE.png", 
   "https://i.stack.imgur.com/CZ2sU.png", 
   "https://i.stack.imgur.com/x44HL.png", 
   "https://i.stack.imgur.com/ugMEg.png"};    
Grid[estimateScale /@ imgs]

Answer 2

You can define acutance to measure sharpness of an image part:

acutance[img_] := Mean@Flatten@ImageData@GradientFilter[img, 1]

You can then divide your image into blocks and calculate acutance for each of them. Blocks with higher acutance are sharper.

For your example:

blockSize = 50;
img = Import["https://i.stack.imgur.com/dNxZE.png"];
acutanceMap = Map[acutance, ImagePartition[img, blockSize], {2}];
MatrixPlot[acutanceMap, ColorFunction -> "TemperatureMap"]

you will get the following acutance map:

where darker red regions correspond to the sharpest parts.

You can also calculate acutance distribution along the vertical axis:

acutanceDistribution = Transpose[
    {blockSize*Range@Length@acutanceMap, Mean[Transpose@acutanceMap]}]

ListPlot[acutanceDistribution, PlotRange -> All, Joined -> True, Frame -> True]

The sharpest part of the image is located between the lines 350-500:

Results for the test images

Answer 3

We can use Blur or Sharpen to show the blurry model roughly.They will give a similar result.I will use Blur here.

img = Import["https://i.stack.imgur.com/dNxZE.png"];
model = PeronaMalikFilter[img - Blur[img], 10]

We need a border to locate the upper, cnenter and down as your demand.

border = Subdivide[Last[ImageDimensions[img]], 3]

{0,649/3,1298/3,649}

We can use the ListLinePlot to show the trend of blur.

data = MeanFilter[Mean@Transpose[
     ImageData[ImageAdjust[ColorConvert[model, "Grayscale"]]]], 10];
ListLinePlot[data, 
 Epilog -> {Red, Line[{{border[[2]], 0}, {border[[2]], 0.04}}], 
   Line[{{border[[3]], 0}, {border[[3]], 0.04}}]}, 
 Ticks -> {{{Mean[border[[;; 2]]], "Upper"}, {Mean[border[[2 ;; 3]]], 
     "Center"}, {Mean[border[[3 ;;]]], "down"}}}]

As we see,the clear part is not very near to the center,but close to the bottom,and its position is $62\%$ in the vertical direction.

N[First[Ordering[data, -1]]/Length[data]]

0.619414

Answer 4

imgs = Map[
  Import, {
   "https://i.stack.imgur.com/dNxZE.png", 
   "https://i.stack.imgur.com/CZ2sU.png", 
   "https://i.stack.imgur.com/x44HL.png"
    }]

With[{im = ImageAdjust[RangeFilter[#, 6], 3]},
   Show[#,
    Graphics[{Red, Line@First@ImageLines[im, MaxFeatures -> 1]}], 
    ImageSize -> 200]] & /@ imgs

Now we can rotate the images so the focal plane is horizontal

angle[{{x1_, y1_}, {x2_, y2_}}] := ArcTan[(y2 - y1)/(x2 - x1)]

rimg = With[{im = ImageAdjust[RangeFilter[#, 6], 3]},
   ImageRotate[#, -angle@First@ImageLines[im, MaxFeatures -> 1], 
    "MaxAreaCropping"]] & /@ imgs

We can separate strips with equal focus

Flatten@ImagePartition[#, {First@ImageDimensions[#], 41}] &@rimg[[2]]

Now we compare each sub image strip blurred over different pixel radius. You can see that focused strips differ more from Blurred for longer radii.

Table[
 Show[
    Plot[
     Evaluate@{
       LinearModelFit[
         Take[#, -10]
         , x, x][x],
       LinearModelFit[
         Prepend[Take[#, 2], {0, 0}]
         , x, x][x]
       }, {x, 0, 20}, PlotRange -> {0, 14000}],
    ListPlot[Prepend[#, {0, 0}], PlotRange -> {0, 10000}]
    ] &@Table[{n, 
    Total@Abs@
      Flatten[Subtract @@ (ImageData /@ {i, Blur[i, n]})]}, {n, 1, 20}]
 , {i, Flatten@ImagePartition[rimg[[1]], {866, 50}]}]

You can calculate the intersection of the two lines, i.e the radius at which the Blur error becomes shallow, which will be many pixels for sharp sub-image strips, and few for blurry sub-image strips.

Table[
 (x /.
     First@Solve[LinearModelFit[
          Take[#, -10]
          , x, x][x] ==
        LinearModelFit[
          Prepend[Take[#, 2], {0, 0}]
          , x, x][x], x]) &@
  Table[{n, 
    Total@Abs@
      Flatten[Subtract @@ (ImageData /@ {i, Blur[i, n]})]}, {n, 1, 20}]
 , {i, Flatten@ImagePartition[rimg[[1]], {866, 15}]}]