10

What is a resourceful method to calculate FWHM (full width at half the max) of every column in an image. Assuming the image is having a very high resolution for example 1944 x 2592.

Logic so far (applied to every column): sample image

gBW = ImageApply[Mean, zeroS];
data = ImageData[gBW];
datay = data[[All, 1]];

(* For column 1 through 1992 *)
datax = Table[i, {i, 1, Length[data]}];
ListLinePlot[data[[All, 1]], PlotStyle -> Thick,  PlotTheme -> "Detailed"];
(*Colum 1 Peak*)
datay = datay/Max[datay];
L = Length[datay];
Mag = 4;
PP = 2.2;
(*Find the centerindex of maximum*)
CenterIndex = Position[datay, Max[datay]][[1]]; (*Find index of max, first occurence*)
(*start searching lead trail*)
i = 2;
While[Sign[datay[[i]] - 0.5] == Sign[datay[[i - 1]] - 0.5],
      i = i + 1];
Interp = ((0.5 - datay[[i - 1]])/(datay[[i]] - datay[[i - 1]]))
Tlead = datax[[i - 1]] + Interp*(datax[[i]] - datax[[i - 1]])
i = CenterIndex[[1]] + 1

(* Start Searching for the next crossing at center *)
While[Sign[datay[[i]] - 0.5] == Sign[datay[[i - 1]] - 0.5] && (i <= L - 1),
       i = i + 1];
If[i != L, 
    Interp = (0.5 - datay[[i - 1]])/(datay[[i]] - datay[[i - 1]]);
    Ttrail = datax[[i - 1]] + Interp*(datax[[i]] - datax[[i - 1]]);
    FWHM = ((Ttrail - Tlead)/Mag)*PP]

Is there a better approach to do this or a logic to find FWHM along all the columns on the image and plot the same? ( assuming it is going to be a huge number like 2592 columns)

Application: Spectroscopy.

Dr. belisarius
  • 115,881
  • 13
  • 203
  • 453
Rene Duchamp
  • 1,419
  • 10
  • 20

2 Answers2

4

This is a fast discrete approach

i = Import["https://i.stack.imgur.com/D2GVe.png"];
tid = Transpose[ ImageData[ColorSeparate@ColorConvert[i, "HSB"] // Last]];
nm = Max /@ tid;
Histogram[2.2/4 MapThread[Count[Sign[#1 - #2/2], 1] &, {tid, nm}]]

Mathematica graphics

Addressing @BlacKow comments below, the following shows that the method works well for this data because all the included segments are contiguous to the maximum:

f[n_] := TakeWhile[Ordering[-tid[[n]]], tid[[n, #]] > Max@tid[[n]]/2 &] // 
                                                       Sort // Differences
fr = f /@ Range@Length@tid;
Times @@ (Flatten@fr)
(* 1 *)

So, the following is "safe" for almost any reasonable dataset:

Histogram[2.2/4 ((f[#] // Split // First // Length) + 1) & /@ Range@Length@tid]

Edit

A more careful treatment gives almost the same result:

i = Import["https://i.stack.imgur.com/D2GVe.png"];
tid = Transpose[ImageData[ColorSeparate@ColorConvert[i, "HSB"] // Last]];
nm = Max /@ tid;
rr = Unitize@MapThread[Sign[#1 - #2/2] + 1 &, {tid, nm}];
mask = {{-1, -1, -1, 1, 1, 1}};
dilated = Dilation[#, 3] & /@ ImageData@
                 HitMissTransform[Image@rr, #, Padding -> 0] & /@ {mask, -mask};
nums = (tid #) & /@ dilated;
ranges = ConstantArray[List /@ Range@Length@First@First@dilated,
                                                     {2, Length@First@dilated}];
lists = MapThread[List, {ranges, nums}, 3];
onlynums = lists /. {{_}, 0.} :> Sequence[];
ints = Transpose@Map[Interpolation, onlynums, {2}];
FWHM = -2.2/4 Subtract @@@ 
    MapThread[
     x /. FindRoot[#1[x] == #2, Evaluate@{x, Sequence @@ #1["Domain"][[1]]}] &, 
              {ints, Transpose@ConstantArray[nm/2, 2]}, 2];

Histogram@FWHM

Mathematica graphics

Dr. belisarius
  • 115,881
  • 13
  • 203
  • 453
  • So you are basically counting total number of points that are higher than half of the maximum, correct? So you are also counting points that don't belong to the peak, but just happen to be higher than half maximum. – BlacKow Jun 24 '15 at 20:51
  • @BlacKow Yes, but those points you mention don't seem to exist Manipulate[ ListPlot[tid[[i]], Epilog -> {Red, Line[{{0, #}, {82, #}}] &@Max[tid[[i]]/2]}], {i, 1, 1900}] – Dr. belisarius Jun 24 '15 at 20:55
  • Sure… For this particular dataset it can be ok, not very good approach for general solution. We somewhat discussed it here: http://mathematica.stackexchange.com/questions/85774/cutting-away-points-that-are-not-resonances/85793 – BlacKow Jun 24 '15 at 21:17
  • 2
    @BlacKow Speed and generality are opposing drivers only the OP knows how to balance – Dr. belisarius Jun 24 '15 at 21:41
  • You are exactly right! +1 – BlacKow Jun 24 '15 at 21:42
  • Speed is not an issue, I was looking if there was a more general way. Accuracy is more important. The reason why I want to treat every column as an image by itself. This is my need. My approach is correct, but I am looking for a method to do across the whole image and plot the variation of FWHM of each column to see if it lies anywhere between +/- 15 microns. – Rene Duchamp Jun 24 '15 at 22:12
  • @abhilashsukumari We're measuring pixels here. How does a pixel relates to a micron? Do you have the exact scale factor? – Dr. belisarius Jun 24 '15 at 22:46
  • 1
    as shown in the code, each pixel corresponds to [2.2/4]microns. As the magnification factor is 4x and 2.2 is the pixel pitch. – Rene Duchamp Jun 24 '15 at 23:19
  • 1
    @BlacKow See edits, please – Dr. belisarius Jun 25 '15 at 13:17
2

(Just an extension of belisarius' excellent answer which I have upvoted.)

The physically meaningful approach to the determination of FWHM requires working with actual physical intensity values, not gamma-corrected and artificially deformed intensity values which one gets with the usual photographic cameras generating images in the sRGB colorspace. If it is known that the photo is in the sRGB colorspace, the following conversion function originally published by Jari Paljakka can be used for conversion it into the linear RGB colorspace where channel values represent the actual physical intensity values:

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

When having actual physical intensity values, it is reasonable to define the overall intensity as a sum of linear intensities for all channels, so belisarius' code can be modifies as follows:

i = Import["https://i.stack.imgur.com/D2GVe.png"];
tid = Transpose@Total[srgb2linear[ImageData[i, Interleaving -> False]], 1];
nm = Max /@ tid;
Histogram[2.2/4 MapThread[Count[Sign[#1 - #2/2], 1] &, {tid, nm}]]

histogram

As one can see, working in the linear color space is crucially important in order to obtain correct estimates of FWHM.

Alexey Popkov
  • 61,809
  • 7
  • 149
  • 368