Skip to main content

probability or statistics - FWHM of an image


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.



Answer



This is a fast discrete approach


i = Import["http://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["http://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


Comments

Post a Comment

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more