# Estimate the "Blurry" distribution of an image

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:

• What is the definition of "approximate blurring radius"? May 27, 2017 at 10:34
• If the image is given by blurimg=Blur[sharpimg, radius] I would like the to have a matrix telling me this blurimg has a blurring radius of radius in all parts of this image. and in this case, I would like to have a matrix telling me something like: the clear part in the center is with a blurring radius of 1px, while the part next to it seems like some image been blurred by 5px. Is this explanation satisfying?
– Wjx
May 27, 2017 at 10:40
• Note that at present, both answers cannot cope with the last sample image properly, in the last image, the lighter parts, especially the right bottom one, is blurry, while the central, darker parts are clear. Now both algorithm failed to acknoledge this because both results are dominanted by those edges with extremely large gradient instead of shapness of small details, sometimes with smaller gradient, which determine our human perception of "sharp" while with .
– Wjx
May 27, 2017 at 10:53
• Some local variance by linear filtering should do it, but I am a happy beginner in the language so I should probably not be the one to demonstrate how to do it using the language. May 27, 2017 at 11:01
• @mathreadler maybe you can demonstrate your idea in the comment in a more mathematical form? I may complete the language part and post it as an community wiki~
– Wjx
May 27, 2017 at 11:43

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]
}
)


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]


• This is truly awesome!
– Wjx
May 28, 2017 at 1:31
• A small advice: maybe some preprocessing like img1=img/Blur[img, 50] would help to make the result of the last image more reasonable.
– Wjx
May 28, 2017 at 3:29
• Also, some form of caliberation can help to achieve better radius estimation.
– Wjx
May 29, 2017 at 14:55

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.

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

– Wjx
May 27, 2017 at 9:44
• A small suggestion: maybe change the second parameter of "GradientFilter" to something like 2 would improve the noise resilliance of this approach.
– Wjx
May 27, 2017 at 10:49
• Note that the last acutance map is contradictory to human perception, which is exactly the case I've mentioned in the comment of the question. The central part is in focus while the brighter parts are out of focus. It's just the opposite to your acutance map.
– Wjx
May 27, 2017 at 15:26
• @Wjx Compared to the previous three images, the last one has almost uniform sharpnes. This can be seen on the acutance map. You can try to adjust GradientFilter parameters and block sizes to get the result closer to what you want. May 27, 2017 at 16:04

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[
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

– Wjx
May 27, 2017 at 9:44
• @Wjx Be my guest and feel free,please.Just choose what you want to choose. :)
– yode
May 27, 2017 at 10:11
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}]}]


• Facinating~ so how to convert the value you've got to the blurry radius?
– Wjx
Jun 1, 2017 at 0:55
• @Wjx I don't think that can be done directly. You could deduce it by knowing the Depth of Filed and deducing the angle. If the sample is flat, then should be linear with the distance to focus. Jun 1, 2017 at 14:16