13
$\begingroup$

EDIT: Solutions by @Alexey Popkov and @Vitaliy Kaurov are very intuitive and both can be used to find a solution for the task.

I have a slit of 15-18 microns by 1 mm projected with magnified 4X onto a cmos sensor of 1944 x 2592 pixels. However, the width is not uniform and I want to find a method which shows how much variation occurs in the width along the whole length of the slit, which is 4 mm when magnified.

Example image of the slit:

enter image description here

$\endgroup$
  • $\begingroup$ How is your question related to the Mathematica software program ? $\endgroup$ – image_doctor Jun 19 '15 at 19:54
  • 1
    $\begingroup$ use of image correlation technique. I have done it before. I want to know if there is a way of doing a rectangular fit on the slit image. $\endgroup$ – Rene Duchamp Jun 19 '15 at 19:56
  • $\begingroup$ Calculate max value of every column, max/2 is threshold. Calculate number of pixels for all values above threshold times pixel pitch gives width of that col $\endgroup$ – Rene Duchamp Jun 19 '15 at 20:03
  • $\begingroup$ How would you define the width ? $\endgroup$ – image_doctor Jun 19 '15 at 20:04
  • 3
    $\begingroup$ ListPlot@DeleteCases[Total@ImageData[Binarize@img], 0] where img is the image. Scale the data as required to get desired units - as written this will be pixels. $\endgroup$ – ciao Jun 19 '15 at 23:23
30
$\begingroup$

I think the essence of the problem here is that width needs to be counted orthogonally to some best fit line going through the elongated shape. Even naked eye would estimate some non-zero slope. We need to make line completely horizontal on average. We could use ImageLines (see this compact example) but I suggest optimization approach. Import image:

i = Import["http://i.stack.imgur.com/BGbTa.jpg"];

See slope with this:

ImageAdd[#, ImageReflect[#]] &@ImageCrop[i]

enter image description here

Use this to devise a function to optimize:

f[x_Real] := Total[ImageData[
   ImageMultiply[#, ImageReflect[#]] &@ImageRotate[ImageCrop[i], x]], 3]

Realize where you are looking for a maximum:

Table[{x, f[x]}, {x, -.05, .05, .001}] // ListPlot

enter image description here

Find more precise maximum:

max = FindMaximum[f[x], {x, .02}]

{19073.462745098062, {x -> 0.02615131131124671}}

Use it to zero the slope

zeroS = ImageCrop[ImageRotate[i, max[[2, 1, 2]]]]

enter image description here

ListPlot3D[ImageData[ColorConvert[zeroS, "Grayscale"]], 
 BoxRatios -> {5, 1, 1}, Mesh -> False]

enter image description here

and get the width data (you can use different Binarize threshold or function):

data = Total[ImageData[Binarize[zeroS]]];
ListLinePlot[data, PlotTheme -> "Detailed", Filling -> Bottom, AspectRatio -> 1/4]

enter image description here

Get stats on your data:

N[#[data]] & /@ {Mean, StandardDeviation}

{14.28099910793934, 1.7445029175852613}

Remove narrowing end points outliers and find that your data are approximately under BinomialDistribution:

dis = FindDistribution[data[[5 ;; -5]]]
BinomialDistribution[19, 0.753676644441292`]

Show[Histogram[data[[5 ;; -5]], {8.5, 20, 1}, "PDF", PlotTheme -> "Detailed"],
 DiscretePlot[PDF[dis, k], {k, 7, 25}, PlotRange -> All, PlotMarkers -> Automatic]]

enter image description here

$\endgroup$
  • 4
    $\begingroup$ That was fun indeed! $\endgroup$ – chris Jun 20 '15 at 8:05
  • $\begingroup$ Finding more precise maximun and using it to find the slope, you have used min instead of max ? and my system gives {27089.4, {x -> 0.0241678}} for the same max = FindMaximum[f[x], {x, .02}]. Why?? $\endgroup$ – Rene Duchamp Jun 22 '15 at 23:05
  • $\begingroup$ @abhilashsukumari min was just a typo - I fixed it. Of course we are looking for max. $\endgroup$ – Vitaliy Kaurov Jun 23 '15 at 7:47
  • $\begingroup$ Got it. Is there a method to say how many of them lie in the range-{14,15} and exact 15 and range-{15,16} from data = Total[ImageData[Binarize[zeroS]]]. Depending on this can we display pass/fail criteria in mathematica ? I am trying to use bspec of histogram range. (should I add this as a new question?) $\endgroup$ – Rene Duchamp Jun 23 '15 at 22:20
  • $\begingroup$ @abhilashsukumari you should try experimenting with code. documentation is very detailed. just search for keywords and look through examples. something like this Length[Select[data, 14 <= # <= 15 &]] $\endgroup$ – Vitaliy Kaurov Jun 23 '15 at 23:14
18
$\begingroup$

The answer by Vitaliy is great but his approach has one drawback: ImageRotate introduces artifacts depending on the Resampling method which affects the final estimates for the slit width. A better solution would process the original image data without distorting it.

The following approach does not include any artificial manipulations with the original data and does not distort it in any way. The only crucial and arbitrary value it depends upon is the Binarize threshold (I use the default Otsu's algorithm, see some explanations here). One should keep in mind that selection of another threshold will shift the estimates for the slit width and may give more appropriate results depending on the actual task.


Convert the image data into the convenient form:

i = Import["http://i.stack.imgur.com/BGbTa.jpg"];

data = PixelValuePositions[Binarize[i], 1];

Now we can easily find the envelopes of the slit:

envelopes = 
  Sort /@ Transpose[
    Thread[{#[[1, 1]], MinMax[#[[;; , 2]]]}] & /@ 
     GatherBy[data, First]];

ListLinePlot[envelopes, PlotRange -> All, AspectRatio -> .1, 
 ImageSize -> 1000, Filling -> {1 -> {2}}, PlotTheme -> "Detailed"]

plot

Find the rotation angle of the slit:

Clear[a, b];

fit[theta_] := FindFit[data.RotationMatrix[theta], a + b x, {a, b}, x];

angle = FindRoot[(b /. fit[theta]) == 0, {theta, -.02}, Evaluated -> False]
{theta -> -0.0227125}

Rotate the envelopes and subtract the constant term:

envelopesCorr = 
  With[{a = a /. fit[angle[[1, 2]]]}, 
   TranslationTransform[{0, -a}] /@ envelopes.RotationMatrix[angle[[1, 2]]]];

ListLinePlot[envelopesCorr, PlotRange -> All, AspectRatio -> .1, 
 ImageSize -> 1000, Filling -> {1 -> {2}}, PlotTheme -> "Detailed"]

plot

Construct InterpolatingFunctions of the corrected envelopes and use them for upsampling the envelopes, then subtract from one another for obtaining the list of widths of the slit:

envelopesCorrF = 
  Interpolation[#, InterpolationOrder -> 1, 
     "ExtrapolationHandler" -> {Indeterminate &, 
       "WarningMessage" -> False}] & /@ envelopesCorr;

widthList = 
  Transpose[{#, Abs[Subtract @@ Through[envelopesCorrF@#]]}] &@
   Union[Flatten@envelopesCorr[[;; , ;; , 1]]];

ListLinePlot[widthList, PlotRange -> All, AspectRatio -> .1, 
 ImageSize -> 1000, Filling -> Bottom, PlotTheme -> "Detailed"]

Histogram[widthList[[;; , 2]], {11, 24, 1}, "PDF", PlotTheme -> "Detailed"]

plot

histogram


UPDATE: Selection of the Method option for Binarize

As I stated above, the final estimates of the slit width crucially depend on the selection of the Binarize threshold. Let us investigate it.

For the purposes of this analysis I prefer to work with the original lossless image of the slit instead of the lossy JPEG image shown in the question. The original image is in BMP format and of size 14.4 Mb. After converting to PNG we obtain exactly identical image of size 361 Kb which can be uploaded to imgur.com:

i=Import["http://i.stack.imgur.com/ZwVkc.png"];

We can reasonably assume that this image is in the sRGB colorspace. Let us examine the slit profile somewhere in the middle of the slit:

(* https://mathematica.stackexchange.com/a/20495/280 *)
iC = ImagePad[i, -BorderDimensions[i, 0.01]];

sRGBData = ImageData[iC, Interleaving -> False];
ListLinePlot[Transpose[sRGBData, {1, 3, 2}][[;; , 1200]], PlotRange -> All, 
 PlotStyle -> {Red, Green, Blue}]

sRGB plot

As we see, the intensity curves in the sRGB color space are not sharp and clear. Now let us convert the original image into the linear RGB colorspace where the channel values correspond to the actual physical intensity values. We will use the implementation of the corresponding formula from the Specification of sRGB published by Jari Paljakka:

srgb2linear = 
  Compile[{{Csrgb, _Real, 1}}, 
   With[{\[Alpha] = 0.055}, 
    Table[Piecewise[{{C/12.92, C <= 0.04045}, 
                     {((C + \[Alpha])/(1 + \[Alpha]))^2.4, C > 0.04045}}], 
          {C, Csrgb}]], RuntimeAttributes -> {Listable}];

linearRGBData = srgb2linear[sRGBData];

ListLinePlot[Transpose[linearRGBData, {1, 3, 2}][[;; , 1200]], PlotRange -> All, 
 PlotStyle -> {Red, Green, Blue}]

linear RGB curves

Now we see quite decent intensity curves which allow us to find the edges of the slit with good precision. The most informative is the Green channel which has almost straight slopes. Experimentation with available choices for the Method option of Binarize reveals that the most correct results for cropped image can be obtained with the method "Mean":

Show[ListLinePlot[Transpose[linearRGBData, {1, 3, 2}][[;; , 1200]], PlotRange -> All, 
  PlotStyle -> {Red, Green, Blue}], 
 ListLinePlot[ImageData[Binarize[iC, Method -> "Mean"]][[;; , 1200]], PlotStyle -> Gray]]

plot

$\endgroup$
  • $\begingroup$ Good for considering Linear RGB data for analysis. $\endgroup$ – Rene Duchamp Aug 27 '15 at 16:43
9
$\begingroup$

If you are looking for a quick estimate, you can erode the line until it disappears. The point at which it is gone gives an estimate of the thickness of the original line. In this implementation, adjust the slider until the line is just barely visible. The t parameter is the number of pixels eroded, and so corresponds to the thickness of the line. When I do it, I see it separating into multiple components by the 12th or 13th erosion and disappearing after the 14th or 15th.

Manipulate[Erosion[img, t], {t, 0, 30, 1}]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.