- Email: [email protected]
Shift and intensity modeling in spectroscopy—general concept and applications
Chemometrics and Intelligent Laboratory Systems 45 Ž1999. 361–370
Shift and intensity modeling in spectroscopy—general concept and applications Frank Westad b
a,)
, Harald Martens
b
a Camo ASA, P.O. Box 1405 Vika, N-0115 Oslo, Norway Norwegian UniÕersity of Science and Technology, Institute of Physical Chemistry, N-7034 Trondheim, Norway
Received 7 January 1998; revised 11 June 1998; accepted 11 June 1998
Abstract Multivariate analysis and calibration in spectroscopy have been widely used in many fields of applications in recent years. The number of components to use in such applications for sufficiently describing the spectra is dependent on the stability of positions of the peaks of interest. Unwanted shifts give rise to a more complex model, i.e., more components. Various kinds of shift correction can be applied as preprocessing tools to cope with these effects. In this paper we describe a general concept of how to model intensity Žpeak height. and position as separate phenomena. We described methods and theory from image modeling, and apply them on one-dimensional signals such as traditional Raman spectroscopy. q 1999 Elsevier Science B.V. All rights reserved. Keywords: Shift and intensity modeling; Spectroscopy
Contents 1. Introduction .
...................................................
2. Theory . . . . . . . . . . . . . . . 2.1. Shift and intensity modeling 2.2. Motion estimation . . . . . .
........................................ ........................................ ........................................
362 362 364
............................................
368
...................................................
370
3. Results from experiments . 4. Conclusions .
Acknowledgements . References
)
362
.................................................
370
......................................................
370
Corresponding author.
0169-7439r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 9 - 7 4 3 9 Ž 9 8 . 0 0 1 4 4 - 0
362
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
1. Introduction In image modeling and analysis, there has been a strong focus on methods to find the motion fields between a sequence of images Žframes.. This enables the possibility of defining objects as being pixels that moves coherently for a sequence of frames. When such motions are found for a number of objects, the intensities of the objects can be modeled separately from their motions. This concept of separating motion and intensities as two different phenomena is a concept that has a very general nature. As motion is related to time, we will rather use the term shift to describe objects that change their positions, e.g., peaks shifting on a wavelength axis. Some examples in spectroscopy are shifts in NMR due to chemical andror physical effects, and broadening or narrowing effects in Raman and IR spectroscopy as a function of temperature. In the case of NMR spectroscopy, Brekke et al. w1x, Spraul et al. w2x and Vogels et al. w3x have addressed shift as a problem in multivariate analysis. Vogels et al. solve the problem by a technique called partial linear fit ŽPLF. where a region of interest is shifted w.r.t. a reference until a minimum RMS Žroot mean square. error is reached. Brown and Stoyanova w4x determine frequency and phase shifts for single resonant peaks in NMR which have the same lineshape. In this paper we introduce a more general concept for all types of shifts as a tool for better interpretation and more parsimonious models compared to straight forward intensity modeling.
2. Theory To have a compact notation we introduce the following notation is used in this paper: Ir Reference spectrum In Spectrum Žin position of In. ˆIn Reconstructed spectrum Žby moving Ir to In position. InAtr Spectrum reconstructed by moving In back to reference position DI Difference spectrum ŽIn y Ir. DIAtn Difference spectrum in position of In DIAtr Difference spectrum in reference position dx Shift vector
Here, e.g., DIAtr should be read as ‘Difference spectrum after moving In with shift vector d x Žs InAtr. and subtract it from reference spectrum, Ir’. 2.1. Shift and intensity modeling Assume that a number of spectra have been collected with the purpose of establishing a multivariate model ŽPCA or regression type model. to describe the changes in intensity for the peaks of interest, and that there is a shift effect present in regions of the spectra. Fig. 1 shows an example of such data. The original spectrum used to generate the data was the Raman spectrum of cyclohexanol in the region 758–662 cmy1 , a total of 801 points. The regions 758–752, 707–702 and 667–662 cmy1 were set to 0, and the first peak was shifted in step of 8 units, five steps in both directions, yielding a total of 11 spectra. The second peak was scaled in equal steps to give pure intensity changes ranging from 0.1–1.1 A.U. A smoothing average with 11 points as smoothing basis was performed on the original spectrum. In Figs. 1–3 the abscissa shows the point numbers only, as we are focusing on the shift effects. The first four loadings from a PCA on mean centered spectra are shown in Fig. 2a; Fig. 2b shows the loadings scaled with the singular values. It is clear that the shifts give rise to more components Žthe loadings for the non-shifted data would show only one component.. In fact, the 10th loadings will have a complex structure instead of showing noisy behavior, although the eigenvalue is rather small. Fig. 3 shows the eigenvalues and percent residual variance for the two PCA models on shifted and non-shifted spectra. From these figures we conclude that it is much to gain if we can model the shifts as a separate phenomenon. Let one spectrum in a series of spectra be denoted as reference spectrum, Ir, and let any other spectrum be denoted as In. Fig. 4a shows Ir, In and the difference spectrum DI s In y Ir. Here the left peak in In is uniformly shifted eight points to the right compared to the reference. Thus the spectrum In can be expressed as a function of Ir, the motion vector d x, and error, E: In s move Ž Ž Ir q DIAtr . ,d x . q E s move Ž Ir q Ž InAtr y Ir . ,d x . q E,
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
Fig. 1. Simulated shifts for spectral data.
Fig. 2. Ža. Loadings for first four PCs of data in Fig. 1. Žb. Loadings scaled by singular values.
363
364
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
tion in Fig. 1, these spectra will now have only intensity differences, which can be modeled by one component if the motion vector d x is sufficiently well estimated. 2.2. Motion estimation
Fig. 3. Eigenvalues from PCA on data with Ž1. and without shift Ž2..
where DIAtr means the difference spectrum in reference position. The function moÕe effectuates interpolation to ensure continuity in the resulting spectrum. InAtr is In moved back to reference with motion vector d x. We can also reconstruct In by moving Ir with the motion vector, d x, yielding ˆIn. The difference ˆIn y In is denoted DIAtn. Assume that we have N spectra, and that motion vectors have been estimated for N y 1 of these, with spectrum one as reference, Ir. Then move all spectra back in reference position to give InAtr for all spectra. Compared to the situa-
In image analysis there exist several methods for estimating the motion of pixels for a given number of frames. The main two types are Ži. Pixel matching based methods and Žii. Gradient methods. As image modeling and analysis is a typical ad-hoc science with little deterministic models involved, the implemented algorithms for this purpose vary both in complexity and the degree of sound mathematical foundation. For that reason a large number of parameters are required, or optional, as input to avoid divergence problems and to estimate true, consistent motion. In this paper we will focus on one method by Horn and Schunck w5,6x that requires few input parameters and has a simple algorithm. 2.2.1. The Horn and Schunck algorithm This algorithm is a gradient based method which assumes that intensities are preserved for a number of
Fig. 4. Ir, In and DI before and after shift-correction.
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
365
Fig. 5. Ža. In and Ir. Žb. True and found shifts. Shifts are given as ordinate values.
adjacent points that moves between samples. At first we will describe the method in an image context, i.e., for 2D data, and then show the 1D implementation that is used on our data. We advocate for an initial 2D description because spectroscopic imaging and electrophoresis are other relevant fields of applications. Let the intensity of one point Ž x, y . for a given sample Ži.e., frame. be denoted as I Ž x, y, t .. Assuming that the intensity is constant, dI dt
s0
and using the chain rule we get
EI d x E x dt
EI d y q
E y dt
EI q
Et
s0
and let
us
dx dt
and Õ s
dy dt
we have a linear equation with two unknowns, u and Õ, I x u q I y Õ q It s 0 where we have introduced abbreviations for the partial derivatives, I x , I y , and It . As each point has a number of possible independent motions, we can not solve the above equation without adding some constraint. With the assumption that neighboring points have similar motion, we can express this as to mini-
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
366
Fig. 6. Ža. Blow up of In, Ir and ˆIn. Žb. Difference spectrum before and after moving Ir to In position ŽˆIn..
mize the square of the magnitude of the gradient of the motion:
Eu
2
Eu
2
EÕ
2
EÕ
2
ž / ž / ž / ž / Ex
q
and
Ey
Ex
q
Ey
,
which will be used as a smoothness measure. Laplacians of u and Õ are defined as 2
= 2us
2
E u Ex
2
E u q
Ey
2
2
and = 2 Õ s
Ex
The
2
E Õ 2
1
E Õ q
E y2
To estimate the partial derivatives, it is important that I x , I y and It are consistent, i.e., they should re-
1 For 1D signals this is not a requirement for solving the equations, but we will still impose a smoothness constraint because peaks in spectroscopy are smooth by nature.
fer to the same point in the image at the same time. This is achieved by averaging the differences for adjacent points in a cube consisting of eight points: 1 I x f Ii , jq1, k y Ii , j, k q Iiq1, jq1, k y Iiq1, j, k 4 qIi , jq1, kq1 y Ii , j, kq1 q Iiq1, jq1, kq1 yIiq1, j, kq1 4 1 I y f Iiq1, j, k y Ii , j, k q Iiq1, jq1, k y Ii , jq1, k 4 qIiq1, j, kq1 y Ii , j, kq1 q Iiq1, jq1, kq1 yIi , jq1, kq1 4 1 It f Ii , j, kq1 y Ii , j, k q Iiq1, j, kq1 y Iiq1, j, k 4 qIi , jq1, kq1 y Iiq1, jq1, k q Iiq1, jq1, kq1 yIi , jq1, k 4
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
367
Fig. 7. Ža. Original data. Žb. Shift vectors. Žc. Shift corrected data in reference position ŽInAtr.. Abscissa in cmy1 . Ordinate in Raman intensity.
where i, j and k are indices for x, y and t, respectively. In the equations below we will also need the Laplacians of u and Õ, approximated by = 2 u s Ž u i, j, k y u i, j, k . and = 2 Õ s Ž Õi, j, k y Õi, j, k .. The Laplacian for one point is estimated by a weighted average of neighboring points. For points in a 2D image, the weights for the local averages are 1r12 1r6 1r12
1r6 y1 1r6
1r12 1r6 1r12
The problem then is to minimize both
´ b s E x u q E y Õ q Et and
´c s
Eu
2
Eu
2
EÕ
2
EÕ
2
ž / ž / ž / ž / Ex
q
Ey
q
Ex
q
Ey
with a suitable weighing factor, a , for the ‘true’ motion values Ž u, Õ .. Since the intensity Žbrightness. is
368
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
Fig. 8. Eigenvalues for Ž1. original data Ž2. InAtr.
assumed to have an uncertainty or quantization error, the total error is given by
´ 2 s HH Ž a 2´ c2 q ´ b2 . d xd y where a 2 is proportional to the noise in the measurement. Using the calculus of variation we obtain I x2 u q I x I y Õ s a 2= 2 u y I x It I x I y u q I y2 u s a 2= 2n y I y It Introducing the Laplacians and solving for u and Õ we get
Ž a 2 q Ix2 q I y2 . Ž u y u . s yIx Ž Ix u q I y Õ q It . Ž a 2 q Ix2 q I y2 . Ž Õ y Õ . s yI y Ž Ix u q I y Õ q It . These two equations can be solved by an iterative method like Gauss–Seidel in the form u nq 1 s u n y I x Ž I x u n q I y Õ n q It . r Ž a 2 q I x2 q I y2 . Õ nq 1 s Õ n y I y Ž I x u n q I y Õ n q It . r Ž a 2 q I x2 q I y2 . In the case of 1D data this reduces to u nq 1 s u n y I x Ž I x u n q It . r Ž a 2 q I x2 .
which is the formula that can be solved using Matlab 4.0. We notice that the a value will stabilize regions were the partial derivatives are close to zero.
3. Results from experiments To evaluate the concept of shift correction, we will show two examples. The first one concerns the simulated shift-data shown in Fig. 1. Let Ir and In be the spectra shown in Fig. 5a, and apply H & S motion estimation on these data. H & S was developed in image modeling, with intensities typical in the range 0–255. For that reason, we will scale the data to have equivalent intensities because the alpha value is related to this range of intensity. This will, however, primarily affect the convergence properties, and not the final shift vector. The found shift vector is shown in Fig. 5b. Although the found and true shifts differ for large parts of the spectrum, it does not mean that the reconstruction is wrong for all parts. In regions where the spectrum is constant, e.g., the leftmost points, an incorrect shift is not critical, whereas in the peak regions the reconstruction will be severely affected by
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
incorrect shift estimates. Fig. 6a shows In, Ir and ˆIn for a range of the first peak, and Fig. 6b shows DI and DIAtn. As a measure of error, the RMSE Žroot mean square error. values were computed; yielding 15.3 for DI and 0.8 for DIAtn. As mentioned above, peaks in spectra can also change their shape Žnarrowrwide. due to, e.g., temperature effects, as well as having asymmetric shapes. In such situations, a more flexible model than mere peak shifting is required. To illustrate situations where shift and broadening both are present, we performed motion estimation on a selected region of 12 Raman spectra of chlorocyclohexane collected at temperatures 233–288 K. The mean spectrum was used as the reference spectrum, and the data were scaled to range 0–255 with an alpha value of five as input to the H & S algorithm. The objective of this
369
study was to compare the complexity of the data moved back to reference position ŽInAtr., and the original data. As H & S assumes smoothness in intensities, there will be artifacts around peak maxima, especially when the intensity differences increase. For that reason, the motion vectors were median filtered by 13 points before the moveback operation was applied. Fig. 7 shows two of the peaks in the original data Ža., the unfiltered shift vectors Žb., and InAtr Žc., for four of the spectra. The shift correction effect from the filtered motion vectors dx has reduced the complexity, and thus the sizes of eigenvalues of the data, which are shown in Fig. 8 for the whole spectral range. It is worth to mention that there are artifacts in the moved back spectra for large motions, and intensity differences can also have an impact on the estimated motion vectors. As seen from Fig. 9, which
Fig. 9. Loadings for first principal component. Ža. Original data. Žb. InAtr.
370
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
shows the first loadings from mean centered PCA, the loadings for InAtr Žb., are jagged compared to the original data, but the typical shift is no longer dominating.
Acknowledgements The authors are grateful to Dr. Egil Nodland at the University of Bergen, Norway, for providing the Raman spectra analyzed in this paper.
4. Conclusions References The image-related concept of separate shift and intensity modeling has been proven useful also for one-dimensional spectroscopic data. Methods for motion estimation like, e.g., Horn and Schunck give estimates that can move the spectra with respect to some selected reference spectrum back to reference position. This applies to shifted, asymmetric and broadrnarrow peaks. The intensity models are of lesser complexity compared to the original data. There is also an important aspect of parsimony in separating position and intensity information.
w1x T. Brekke, O.M. Kvalheim, E. Sletten, Anal. Chim. Acta 223 Ž1989. 12–134. w2x M. Spraul, P. Neidig, U. Klauck, P. Kessler, E. Holmes, J.K. Nickolson, B.C. Sweatman, S.R. Salman, R.D. Farrant, E. Rahr, C.R. Beddel, J.C. Lindon, J. Pharm. Biomed. Anal. 12 Ž1994. 1215–1225. w3x J.T.W.E. Vogels, A.C. Tas, J. Venekamp, J. Van Der Greef, Journal of Chemometrics 10 Ž1996. 425–438. w4x T.R. Brown, R. Stoyanova, Journal of Magnetic Resonance, Series B 112 Ž1996. 32–43. w5x B.K.P. Horn, B.G. Schunck, Artif. Intell. 17 Ž1981. 185–203. w6x B.K.P. Horn, B.G. Schunck, Artif. Intell. 59 Ž1993. 81–87.
Shift and intensity modeling in spectroscopy—general concept and applications Frank Westad b
a,)
, Harald Martens
b
a Camo ASA, P.O. Box 1405 Vika, N-0115 Oslo, Norway Norwegian UniÕersity of Science and Technology, Institute of Physical Chemistry, N-7034 Trondheim, Norway
Received 7 January 1998; revised 11 June 1998; accepted 11 June 1998
Abstract Multivariate analysis and calibration in spectroscopy have been widely used in many fields of applications in recent years. The number of components to use in such applications for sufficiently describing the spectra is dependent on the stability of positions of the peaks of interest. Unwanted shifts give rise to a more complex model, i.e., more components. Various kinds of shift correction can be applied as preprocessing tools to cope with these effects. In this paper we describe a general concept of how to model intensity Žpeak height. and position as separate phenomena. We described methods and theory from image modeling, and apply them on one-dimensional signals such as traditional Raman spectroscopy. q 1999 Elsevier Science B.V. All rights reserved. Keywords: Shift and intensity modeling; Spectroscopy
Contents 1. Introduction .
...................................................
2. Theory . . . . . . . . . . . . . . . 2.1. Shift and intensity modeling 2.2. Motion estimation . . . . . .
........................................ ........................................ ........................................
362 362 364
............................................
368
...................................................
370
3. Results from experiments . 4. Conclusions .
Acknowledgements . References
)
362
.................................................
370
......................................................
370
Corresponding author.
0169-7439r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 9 - 7 4 3 9 Ž 9 8 . 0 0 1 4 4 - 0
362
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
1. Introduction In image modeling and analysis, there has been a strong focus on methods to find the motion fields between a sequence of images Žframes.. This enables the possibility of defining objects as being pixels that moves coherently for a sequence of frames. When such motions are found for a number of objects, the intensities of the objects can be modeled separately from their motions. This concept of separating motion and intensities as two different phenomena is a concept that has a very general nature. As motion is related to time, we will rather use the term shift to describe objects that change their positions, e.g., peaks shifting on a wavelength axis. Some examples in spectroscopy are shifts in NMR due to chemical andror physical effects, and broadening or narrowing effects in Raman and IR spectroscopy as a function of temperature. In the case of NMR spectroscopy, Brekke et al. w1x, Spraul et al. w2x and Vogels et al. w3x have addressed shift as a problem in multivariate analysis. Vogels et al. solve the problem by a technique called partial linear fit ŽPLF. where a region of interest is shifted w.r.t. a reference until a minimum RMS Žroot mean square. error is reached. Brown and Stoyanova w4x determine frequency and phase shifts for single resonant peaks in NMR which have the same lineshape. In this paper we introduce a more general concept for all types of shifts as a tool for better interpretation and more parsimonious models compared to straight forward intensity modeling.
2. Theory To have a compact notation we introduce the following notation is used in this paper: Ir Reference spectrum In Spectrum Žin position of In. ˆIn Reconstructed spectrum Žby moving Ir to In position. InAtr Spectrum reconstructed by moving In back to reference position DI Difference spectrum ŽIn y Ir. DIAtn Difference spectrum in position of In DIAtr Difference spectrum in reference position dx Shift vector
Here, e.g., DIAtr should be read as ‘Difference spectrum after moving In with shift vector d x Žs InAtr. and subtract it from reference spectrum, Ir’. 2.1. Shift and intensity modeling Assume that a number of spectra have been collected with the purpose of establishing a multivariate model ŽPCA or regression type model. to describe the changes in intensity for the peaks of interest, and that there is a shift effect present in regions of the spectra. Fig. 1 shows an example of such data. The original spectrum used to generate the data was the Raman spectrum of cyclohexanol in the region 758–662 cmy1 , a total of 801 points. The regions 758–752, 707–702 and 667–662 cmy1 were set to 0, and the first peak was shifted in step of 8 units, five steps in both directions, yielding a total of 11 spectra. The second peak was scaled in equal steps to give pure intensity changes ranging from 0.1–1.1 A.U. A smoothing average with 11 points as smoothing basis was performed on the original spectrum. In Figs. 1–3 the abscissa shows the point numbers only, as we are focusing on the shift effects. The first four loadings from a PCA on mean centered spectra are shown in Fig. 2a; Fig. 2b shows the loadings scaled with the singular values. It is clear that the shifts give rise to more components Žthe loadings for the non-shifted data would show only one component.. In fact, the 10th loadings will have a complex structure instead of showing noisy behavior, although the eigenvalue is rather small. Fig. 3 shows the eigenvalues and percent residual variance for the two PCA models on shifted and non-shifted spectra. From these figures we conclude that it is much to gain if we can model the shifts as a separate phenomenon. Let one spectrum in a series of spectra be denoted as reference spectrum, Ir, and let any other spectrum be denoted as In. Fig. 4a shows Ir, In and the difference spectrum DI s In y Ir. Here the left peak in In is uniformly shifted eight points to the right compared to the reference. Thus the spectrum In can be expressed as a function of Ir, the motion vector d x, and error, E: In s move Ž Ž Ir q DIAtr . ,d x . q E s move Ž Ir q Ž InAtr y Ir . ,d x . q E,
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
Fig. 1. Simulated shifts for spectral data.
Fig. 2. Ža. Loadings for first four PCs of data in Fig. 1. Žb. Loadings scaled by singular values.
363
364
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
tion in Fig. 1, these spectra will now have only intensity differences, which can be modeled by one component if the motion vector d x is sufficiently well estimated. 2.2. Motion estimation
Fig. 3. Eigenvalues from PCA on data with Ž1. and without shift Ž2..
where DIAtr means the difference spectrum in reference position. The function moÕe effectuates interpolation to ensure continuity in the resulting spectrum. InAtr is In moved back to reference with motion vector d x. We can also reconstruct In by moving Ir with the motion vector, d x, yielding ˆIn. The difference ˆIn y In is denoted DIAtn. Assume that we have N spectra, and that motion vectors have been estimated for N y 1 of these, with spectrum one as reference, Ir. Then move all spectra back in reference position to give InAtr for all spectra. Compared to the situa-
In image analysis there exist several methods for estimating the motion of pixels for a given number of frames. The main two types are Ži. Pixel matching based methods and Žii. Gradient methods. As image modeling and analysis is a typical ad-hoc science with little deterministic models involved, the implemented algorithms for this purpose vary both in complexity and the degree of sound mathematical foundation. For that reason a large number of parameters are required, or optional, as input to avoid divergence problems and to estimate true, consistent motion. In this paper we will focus on one method by Horn and Schunck w5,6x that requires few input parameters and has a simple algorithm. 2.2.1. The Horn and Schunck algorithm This algorithm is a gradient based method which assumes that intensities are preserved for a number of
Fig. 4. Ir, In and DI before and after shift-correction.
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
365
Fig. 5. Ža. In and Ir. Žb. True and found shifts. Shifts are given as ordinate values.
adjacent points that moves between samples. At first we will describe the method in an image context, i.e., for 2D data, and then show the 1D implementation that is used on our data. We advocate for an initial 2D description because spectroscopic imaging and electrophoresis are other relevant fields of applications. Let the intensity of one point Ž x, y . for a given sample Ži.e., frame. be denoted as I Ž x, y, t .. Assuming that the intensity is constant, dI dt
s0
and using the chain rule we get
EI d x E x dt
EI d y q
E y dt
EI q
Et
s0
and let
us
dx dt
and Õ s
dy dt
we have a linear equation with two unknowns, u and Õ, I x u q I y Õ q It s 0 where we have introduced abbreviations for the partial derivatives, I x , I y , and It . As each point has a number of possible independent motions, we can not solve the above equation without adding some constraint. With the assumption that neighboring points have similar motion, we can express this as to mini-
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
366
Fig. 6. Ža. Blow up of In, Ir and ˆIn. Žb. Difference spectrum before and after moving Ir to In position ŽˆIn..
mize the square of the magnitude of the gradient of the motion:
Eu
2
Eu
2
EÕ
2
EÕ
2
ž / ž / ž / ž / Ex
q
and
Ey
Ex
q
Ey
,
which will be used as a smoothness measure. Laplacians of u and Õ are defined as 2
= 2us
2
E u Ex
2
E u q
Ey
2
2
and = 2 Õ s
Ex
The
2
E Õ 2
1
E Õ q
E y2
To estimate the partial derivatives, it is important that I x , I y and It are consistent, i.e., they should re-
1 For 1D signals this is not a requirement for solving the equations, but we will still impose a smoothness constraint because peaks in spectroscopy are smooth by nature.
fer to the same point in the image at the same time. This is achieved by averaging the differences for adjacent points in a cube consisting of eight points: 1 I x f Ii , jq1, k y Ii , j, k q Iiq1, jq1, k y Iiq1, j, k 4 qIi , jq1, kq1 y Ii , j, kq1 q Iiq1, jq1, kq1 yIiq1, j, kq1 4 1 I y f Iiq1, j, k y Ii , j, k q Iiq1, jq1, k y Ii , jq1, k 4 qIiq1, j, kq1 y Ii , j, kq1 q Iiq1, jq1, kq1 yIi , jq1, kq1 4 1 It f Ii , j, kq1 y Ii , j, k q Iiq1, j, kq1 y Iiq1, j, k 4 qIi , jq1, kq1 y Iiq1, jq1, k q Iiq1, jq1, kq1 yIi , jq1, k 4
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
367
Fig. 7. Ža. Original data. Žb. Shift vectors. Žc. Shift corrected data in reference position ŽInAtr.. Abscissa in cmy1 . Ordinate in Raman intensity.
where i, j and k are indices for x, y and t, respectively. In the equations below we will also need the Laplacians of u and Õ, approximated by = 2 u s Ž u i, j, k y u i, j, k . and = 2 Õ s Ž Õi, j, k y Õi, j, k .. The Laplacian for one point is estimated by a weighted average of neighboring points. For points in a 2D image, the weights for the local averages are 1r12 1r6 1r12
1r6 y1 1r6
1r12 1r6 1r12
The problem then is to minimize both
´ b s E x u q E y Õ q Et and
´c s
Eu
2
Eu
2
EÕ
2
EÕ
2
ž / ž / ž / ž / Ex
q
Ey
q
Ex
q
Ey
with a suitable weighing factor, a , for the ‘true’ motion values Ž u, Õ .. Since the intensity Žbrightness. is
368
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
Fig. 8. Eigenvalues for Ž1. original data Ž2. InAtr.
assumed to have an uncertainty or quantization error, the total error is given by
´ 2 s HH Ž a 2´ c2 q ´ b2 . d xd y where a 2 is proportional to the noise in the measurement. Using the calculus of variation we obtain I x2 u q I x I y Õ s a 2= 2 u y I x It I x I y u q I y2 u s a 2= 2n y I y It Introducing the Laplacians and solving for u and Õ we get
Ž a 2 q Ix2 q I y2 . Ž u y u . s yIx Ž Ix u q I y Õ q It . Ž a 2 q Ix2 q I y2 . Ž Õ y Õ . s yI y Ž Ix u q I y Õ q It . These two equations can be solved by an iterative method like Gauss–Seidel in the form u nq 1 s u n y I x Ž I x u n q I y Õ n q It . r Ž a 2 q I x2 q I y2 . Õ nq 1 s Õ n y I y Ž I x u n q I y Õ n q It . r Ž a 2 q I x2 q I y2 . In the case of 1D data this reduces to u nq 1 s u n y I x Ž I x u n q It . r Ž a 2 q I x2 .
which is the formula that can be solved using Matlab 4.0. We notice that the a value will stabilize regions were the partial derivatives are close to zero.
3. Results from experiments To evaluate the concept of shift correction, we will show two examples. The first one concerns the simulated shift-data shown in Fig. 1. Let Ir and In be the spectra shown in Fig. 5a, and apply H & S motion estimation on these data. H & S was developed in image modeling, with intensities typical in the range 0–255. For that reason, we will scale the data to have equivalent intensities because the alpha value is related to this range of intensity. This will, however, primarily affect the convergence properties, and not the final shift vector. The found shift vector is shown in Fig. 5b. Although the found and true shifts differ for large parts of the spectrum, it does not mean that the reconstruction is wrong for all parts. In regions where the spectrum is constant, e.g., the leftmost points, an incorrect shift is not critical, whereas in the peak regions the reconstruction will be severely affected by
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
incorrect shift estimates. Fig. 6a shows In, Ir and ˆIn for a range of the first peak, and Fig. 6b shows DI and DIAtn. As a measure of error, the RMSE Žroot mean square error. values were computed; yielding 15.3 for DI and 0.8 for DIAtn. As mentioned above, peaks in spectra can also change their shape Žnarrowrwide. due to, e.g., temperature effects, as well as having asymmetric shapes. In such situations, a more flexible model than mere peak shifting is required. To illustrate situations where shift and broadening both are present, we performed motion estimation on a selected region of 12 Raman spectra of chlorocyclohexane collected at temperatures 233–288 K. The mean spectrum was used as the reference spectrum, and the data were scaled to range 0–255 with an alpha value of five as input to the H & S algorithm. The objective of this
369
study was to compare the complexity of the data moved back to reference position ŽInAtr., and the original data. As H & S assumes smoothness in intensities, there will be artifacts around peak maxima, especially when the intensity differences increase. For that reason, the motion vectors were median filtered by 13 points before the moveback operation was applied. Fig. 7 shows two of the peaks in the original data Ža., the unfiltered shift vectors Žb., and InAtr Žc., for four of the spectra. The shift correction effect from the filtered motion vectors dx has reduced the complexity, and thus the sizes of eigenvalues of the data, which are shown in Fig. 8 for the whole spectral range. It is worth to mention that there are artifacts in the moved back spectra for large motions, and intensity differences can also have an impact on the estimated motion vectors. As seen from Fig. 9, which
Fig. 9. Loadings for first principal component. Ža. Original data. Žb. InAtr.
370
F. Westad, H. Martensr Chemometrics and Intelligent Laboratory Systems 45 (1999) 361–370
shows the first loadings from mean centered PCA, the loadings for InAtr Žb., are jagged compared to the original data, but the typical shift is no longer dominating.
Acknowledgements The authors are grateful to Dr. Egil Nodland at the University of Bergen, Norway, for providing the Raman spectra analyzed in this paper.
4. Conclusions References The image-related concept of separate shift and intensity modeling has been proven useful also for one-dimensional spectroscopic data. Methods for motion estimation like, e.g., Horn and Schunck give estimates that can move the spectra with respect to some selected reference spectrum back to reference position. This applies to shifted, asymmetric and broadrnarrow peaks. The intensity models are of lesser complexity compared to the original data. There is also an important aspect of parsimony in separating position and intensity information.
w1x T. Brekke, O.M. Kvalheim, E. Sletten, Anal. Chim. Acta 223 Ž1989. 12–134. w2x M. Spraul, P. Neidig, U. Klauck, P. Kessler, E. Holmes, J.K. Nickolson, B.C. Sweatman, S.R. Salman, R.D. Farrant, E. Rahr, C.R. Beddel, J.C. Lindon, J. Pharm. Biomed. Anal. 12 Ž1994. 1215–1225. w3x J.T.W.E. Vogels, A.C. Tas, J. Venekamp, J. Van Der Greef, Journal of Chemometrics 10 Ž1996. 425–438. w4x T.R. Brown, R. Stoyanova, Journal of Magnetic Resonance, Series B 112 Ž1996. 32–43. w5x B.K.P. Horn, B.G. Schunck, Artif. Intell. 17 Ž1981. 185–203. w6x B.K.P. Horn, B.G. Schunck, Artif. Intell. 59 Ž1993. 81–87.