5
$\begingroup$

Preamble:

I designed some code to help me take spectral images from a spectrometer and convert pixel number into wavelength. It does reasonably well at matching a calibration spectrum of a neon lamp to known spectral lines of neon, but only if I make sure to include all actual neon peaks from the calibration spectrum. Sometimes it's hard to tell if a peak in the spectrum is just noise/contamination or if it's real. If I exclude real peaks, my algorithm fails miserably. It's not a huge problem, but it'd be neat to automate this even further.

The Problem:

The essence of the problem is this: I have a series of known spectral line of neon. I take a spectrum of a neon lamp, and then want to match the lines I see in the image up to the known lines. However, I may not have exactly the same number of calibration lines as I do reference lines.

My reference data is this:

refdat = {{585.25, 27}, {588.19, 11}, {594.48, 29}, {597.55, 
7}, {603.00, 6}, {607.43, 20}, {609.62, 35}, {614.31, 
49.5}, {616.36, 16}, {621.73, 10.5}, {626.65, 24}, {630.48, 
11}, {633.44, 29}, {638.30, 43}, {640.22, 78.5}, {650.65, 
60}, {653.29, 19}, {659.90, 25.5}, {667.83, 39.5}, {671.70, 
30}, {692.95, 35.5}, {703.24, 100}, {717.39, 4.5}, {724.52, 50}, {743.89, 6}};

The x-values (wavelengths in nm) are the most reliable, so I would prefer to use these only. The intensity response (y-values) are pretty non-linear so I'm not sure they're of much use.

A good set of data looks like this:

good = {{11, 0.1032706710176553}, {44, 0.1146552044370179}, {115, 
  0.2113968543791952}, {150, 0.05051515444593649}, 
   {213, 0.052753899062279425}, {264, 0.13868347448613108}, {290, 
  0.22283763808371357}, 
   {344, 0.4252736853859361}, {368, 0.12755375457964316}, {430, 
  0.10566687181148425}, 
   {487, 0.2319607538664892}, {531, 0.08104101568883736}, {565, 
  0.2742127219013794}, 
   {622, 0.4568371231579498}, {644, 0.7877989366604489}, {765, 
  0.36850130722650787}, 
   {795, 0.13327028393091367}, {872, 0.1449608891795199}, {964, 
  0.20420376201460919}, 
   {1009, 0.13586396031680412}, {1256, 0.18521400598076912}, {1376, 
  1.}}

Good calibration data.

and a "bad" set of data might look like this:

{{11, 0.1032706710176553}, {44, 0.1146552044370179}, {115, 
  0.2113968543791952}, {264, 0.13868347448613108}, 
   {290, 0.22283763808371357}, {344, 0.4252736853859361}, {368, 
  0.12755375457964316}, 
   {430, 0.10566687181148425}, {487, 0.2319607538664892}, {565, 
  0.2742127219013794}, 
   {622, 0.4568371231579498}, {644, 0.7877989366604489}, {765, 
  0.36850130722650787}, 
   {795, 0.13327028393091367}, {872, 0.1449608891795199}, {964, 
  0.20420376201460919}, 
   {1009, 0.13586396031680412}, {1256, 0.18521400598076912}, {1376, 
  1.}}

Bad calibration data.

Basically, the threshold on my FindPeaks call just missed out on 3 peaks that were real. My algorithm is just using FindMinimum to do a linear fit of the pixel values to the reference data as in scale * pixel_values + offset:

results = Table[
   FindMinimum[
     Mean[(scale*pks[[All, 1]] - offset - 
        refdat[[i ;; i + Length@pks - 1, 1]])^2], 
      {{scale, 22.8}, {offset, 617.4}}],
    {i, Length@refdat - Length@pks + 1}];

Where pks is either good or bad. If I only have 22 peaks, but refdat has 25 points, I naively try refdat[[1;;22]], refdat[[2;;23]], refdat[[3;;24]], and refdat[[4;;25]]. This is fine if I happen to have selected every real peak but my spectrometer only collected the first or last $n$ peaks; it fails if I miss one in the middle somewhere.

Here is what happens if I have "good" data (it's pretty clear in this data what is a peak and what isn't, but that's not always the case):

Calibration of good data.

But if I miss those 3 points:

enter image description here

My Question:

Is there a reasonably efficient algorithm for linearly scaling one vector to match another with (a few; maybe 2-4) possible missed values?

Simply trying every combination of values seems prohibitive. My current code isn't too bad as I can usually pick out all the real values after a couple of tries with moving the FindPeaks threshold, though not always. As a last resort I can always pick out values manually, but it would be neat to know if there's a better solution. This seems like a problem that would probably be well-studied from a mathematical perspective, but I'm not well-versed enough to find it.

$\endgroup$
2
$\begingroup$

If most points in refdat do have corresponding points in good then you should be able to find a stretch of $n$ points that exist in both sets. Let's say that $n=5$. You can then find the scaling like this:

ref = Partition[First /@ refdat, 5, 1];
data = Partition[First /@ good, 5, 1];

findLinearModel[ref_, data_] := LinearModelFit[Transpose[{data, ref}], x, x]
findLinearModelError[ref_, data_] := Module[{lm},
  lm = findLinearModel[ref, data];
  Norm[ref - lm /@ data]
  ]

errors = Outer[findLinearModelError, ref, data, 1];
minPos = Position[errors, Min[errors]];

lms = Outer[findLinearModel, ref, data, 1];
lm = First@Extract[lms, minPos];

lines = InfiniteLine[{#, 0}, {0, 1}] & /@ (lm@*First /@ good);
intersections = Flatten[{x, y} /. Quiet@Solve[{x, y} ∈ #, {x, y} ∈ Line[refdat]] & /@ lines, 1];

ListLinePlot[
 refdat,
 Epilog -> {
   Red, PointSize[Large], Point[intersections]
   }]

Mathematica graphics

With this estimate based on five points, you can now exclude points you don't have correspondences for and find a new transform based on all of the ones you have, if you want.

x1 = refdat[[All, 1]];
x2 = lm /@ good[[All, 1]];
filteredRefDat = First /@ Nearest[x1, x2];
lm = LinearModelFit[Transpose[{First /@ good, filteredRefDat}], x, x];

lines = InfiniteLine[{#, 0}, {0, 1}] & /@ (lm@*First /@ good);
intersections = Flatten[{x, y} /. Quiet@Solve[{x, y} ∈ #, {x, y} ∈ Line[refdat]] & /@ lines, 1];

ListLinePlot[
 refdat,
 Epilog -> {
   Red, PointSize[Large], Point[intersections]
   }]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Thanks, this works pretty well! It's usually (though not always) possible to get 4 or 5 peaks in a row. I changed x1 and x2 in your filteredRefDat to refdat[[All, 1]] and intersection[[All, 1]] - I assume that was what was supposed to go there. $\endgroup$ – MassDefect Aug 17 at 19:20
  • $\begingroup$ @MassDefect Thanks, updated to fix that. $\endgroup$ – C. E. Aug 17 at 19:31

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.