# Detecting and correcting gradational brightness variations in a grayscale image

I am trying to quantify the proportions of various components in a greyscale image (backscattered electron image of a polished rock sample). Here is the original image:

bse = Import@"http://i.stack.imgur.com/eQece.png";


The image is comprised of components (minerals) that should have the same intensity response. For example, the abundant 'mid-grey' component (plagioclase) should have a fixed grey-level across the entire image. It doesn't.

By masking the other components of the image and adjusting the contrast, the gradations in brightness become more obvious:

bsek = ImageSubtract[bse, ColorNegate@Binarize[bse, {0.16, 0.24}]];


There seems to be two effects: a first order decrease in brightness towards the bottom corners and a second-order periodic vertical striping (comprised of short wavelength gradients).

Why bother? - Correcting for variations in brightness is an important step prior to segmentation analysis of the image

By approximating the variation in grey-values for the plagioclase, the difference in birghtness between the top and bottom of the image can be modelled (albeit crudely).

ksim[im_] := Module[{kdf, kdl},
kdf = SmoothKernelDistribution[First@im];
kdl = SmoothKernelDistribution[Last@im];
Plot[{PDF[kdf, x], PDF[kdl, x]}, {x, 0, 0.3}, Filling -> Axis,
Exclusions -> None]]
ksim[ImageData[bsek]]


The following creates a gradational filter that causes the grey-level peaks to overlap (by adjusting the value of shift).

bsedim = ImageDimensions[bse]
shift = 0.2;
Transpose@
Array[Table[x, {x, -shift, 0., shift/bsedim[[1]]}] &, bsedim[[1]]];

ksim[ImageData[bsekm]]


This filter can then be applied to the original image:

bsem = ImageMultiply[bse, gradi]


This clumsy method does an OK job of levelling the first-order brightness defect. However, this effect looks to have a spherical gradient (dimmer in the bottom corners than the middle) which is not accounted for.

As yet, I have not found a way to correct the second order periodic effect.

Can anyone suggest a more straightforward way to detect and correct these kinds of brightness variation defects? Extra bonus gratitude to anyone who can help to correct the second-order periodic aberration.

Update

Here is a contrast stretched image showing the vertical aberrations:

 ImageAdjust[bse, {4.5, 3.85}]


-

Does this help? I've treated the error as additive, though you could easily change it to a multiplicative correction.

Starting with your bsek image, find the mean offset for rows and for columns:

id = ImageData[bsek];
horiz = (# - Mean[#]) &[Total[id]/Total[Unitize[id]]];
vert = (# - Mean[#]) &[Total[id, {2}]/Total[Unitize[id], {2}]];

ListLinePlot[horiz]
ListLinePlot[vert]


The horizontal plot shows an overall darkening at the left and right sides of the image, plus the periodic component. The vertical plot shows an overall darkening towards the bottom of the image, with a periodic component also.

Combining the horizontal and vertical offsets into an overall offset:

a = Outer[Plus, vert, horiz];


Now subtract the offset from the original data:

bse2 = Image[ImageData[bse] - a]


Using s.s.o.'s ClusteringComponents visualisation appears to show an improvement:

ClusteringComponents[bse, 10] // Colorize
ClusteringComponents[bse2, 10] // Colorize


-
This is a really great approach! Many Thanks! There seems to be a 'noise' issue with respect to the vertical and horizontal lines. This is visible if you adjust the contrast, e.g. ImageAdjust[bse2, {4.5, 3.85}]. Perhaps there is some way of smoothing the curves?.. –  geordie Aug 21 '13 at 12:01
I think the problem is that a correction determined at one brightness level does not necessarily apply at another. Ideally you would determine the detector response at multiple signal levels and apply a compensation curve rather than the single-point correction I have shown. Even then, you will still have residual error because of quantisation noise. Of course, there is also an assumption in my approach that the detector non-uniformity is separable into horizontal and vertical components, which won't in general be true. –  Simon Woods Aug 21 '13 at 13:50
@SimonWoods,I have a Question. Not that I know very well. If I look dat = ImageData@bse2; dat1 = ImageData@img; EuclideanDistance[dat, dat1] I get 2.50524 and only dat2 = ColorQuantize[img, 10, Dithering -> False]; EuclideanDistance[dat, ImageData@dat2] I get 2.35325. So, If the value is bigger does that mean you lose more information? Or how can we measure any ideas. Just out of my curiosity. –  s.s.o Aug 21 '13 at 14:23
@s.s.o, I'm not sure what you mean. Surely EuclideanDistance is just measuring how different the images are? –  Simon Woods Aug 21 '13 at 15:14
Simon, you know a lot more about this stuff than I. For my purposes your solution works well. I reprocessed the image using a mask generated from the darkest mineral (quartz). The v and h curves are a lot smoother and the correction seems to work better. Thanks again! –  geordie Aug 21 '13 at 23:49

I now know it is now what OP is looking for exactly, but till there is no better approach I will leave this answer.

dat = ImageData @ pic;
Histogram[dat // Flatten]


Image@Round[dat, .05]


-
Unfortunately not. I am trying to avoid losing 'fidelity' in the grey-levels. There are some very dark components that have very close grey-level ranges. You can see them in my raw and filtered images but they are lost from your filtered image. Thanks for trying. –  geordie Aug 21 '13 at 6:04
This is better for sure, but it is still downsampling the grey-levels(which i would prefer not to do). It might be good enough. I will test it one some other images and let you know. –  geordie Aug 21 '13 at 6:42
@geordie I do not expect it to be a the general approach, I wanted to try if this simple approach can be useful for you. About aberrations, are they on each image? Can you make a plain image with only those variations? If yes, normalize it and divide your measurement images by this so-called "flat". –  Kuba Aug 21 '13 at 6:45
Alas, I cannot obtain a plain image. That would be a great solution as it would probably remove both effects at once. –  geordie Aug 21 '13 at 6:59
I had a good look at the filtered image and, unfortunately, too much information is lost (even with the 0.5 rounding). –  geordie Aug 21 '13 at 7:09

Is this sort of what you need as a result? It's simple approach but I'll delete if it's not you asked for.

img = Import@"http://i.stack.imgur.com/eQece.png";

cl = ClusteringComponents[img];
im = Colorize[cl, ColorFunction -> GrayLevel]


you get

Edit 1

if you know roughly component numbers you can use

cl = ClusteringComponents[img, 10];
im = Colorize[cl]


Edit 2

Based on the solution by cormullion in his answer. (Sightly modified) You can quantify colors.

colorquantized = SortBy[Tally[Flatten[
ImageData[ColorQuantize[img, 12,
Dithering -> False]], 1]], Last]


{{0.0823529, 482}, {0.72549, 1306}, {0.333333, 3627}, {0.223529, 5561}, {0., 13718}, {0.145098, 21902}, {0.101961, 33696}, {0.172549, 37135}, {0.207843, 54467}, {0.196078, 55245}, {0.129412, 60045}, {0.113725, 103816}}

BarChart[colorquantized[[All, 2]],
ChartStyle -> GrayLevel /@ colorquantized[[All, 1]]]


Than you can use clustering to find colors to group to best fit to your purpose by:

FindClusters[colorquantized[[All, 2]], 7]


{{482, 1306}, {3627}, {5561}, {13718}, {21902}, {33696, 37135}, {54467, 55245, 60045, 103816}}

and if you like the number of clusters are ok. you use it

cl = ClusteringComponents[img,7];
im = Colorize[cl]


One disadvantage of this approach is choosing the number of clusters in a nice way (analytic way).

-
No, this approach results in a huge loss in 'fidelity'. The grey-levels are reduced to three levels. All the data for the other components are lost. Thanks anyway. –  geordie Aug 21 '13 at 9:26
Ok. By the way you can adjust the component numbers and get those... You can also try MeanShiftFilter... –  s.s.o Aug 21 '13 at 9:35
you may also interested in this answer about dominant colors as an alternative. mathematica.stackexchange.com/questions/24898/… –  s.s.o Aug 21 '13 at 9:41
Thanks for the suggestions. I have played around with MMA's morphological tools and filters. I will try to use these once the image is 'balanced'. As your second image clearly shows there are first and second order aberrations that interfere with the textural information in the image (the minerals). Removing these aberrations is a necessary first step before applying segmentation tools. –  geordie Aug 21 '13 at 10:10
Yes. The question is now clearer. Good luck with it. –  s.s.o Aug 21 '13 at 10:17