PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Proc Am Control Conf. Author manuscript; available in PMC 2010 May 7.
Published in final edited form as:
Proc Am Control Conf. 2009 June 10; 2009: 3848–3853.
doi:  10.1109/ACC.2009.5159832
PMCID: PMC2865653
NIHMSID: NIHMS133145

Comparison of Control Algorithms for a MEMS-based Adaptive Optics Scanning Laser Ophthalmoscope

Abstract

We compared four algorithms for controlling a MEMS deformable mirror of an adaptive optics (AO) scanning laser ophthalmoscope. Interferometer measurements of the static nonlinear response of the deformable mirror were used to form an equivalent linear model of the AO system so that the classic integrator plus wavefront reconstructor type controller can be implemented. The algorithms differ only in the design of the wavefront reconstructor. The comparisons were made for two eyes (two individuals) via a series of imaging sessions. All four controllers performed similarly according to estimated residual wavefront error not reflecting the actual image quality observed. A metric based on mean image intensity did consistently reflect the qualitative observations of retinal image quality. Based on this metric, the controller most effective for suppressing the least significant modes of the deformable mirror performed the best.

I. INTRODUCTION

Adaptive optics has received considerable attention for vision science applications since it was first applied to the eye in 1997 and shown the first images of single cone photoreceptors in a living human eye [1]. Similar to how the earth’s atmosphere degrades the image quality of ground-based telescopes, aberrations due to the eye’s optical imperfections degrade image quality making it difficult for clinicians and scientists to observe microscopic structures of the retina. Adaptive optics (AO) aims to remove most of these degradations by measuring the aberrations with a wavefront sensor (WFS), and through a feedback policy, adjust the surface profile of a deformable mirror (DM) to minimize the residual wavefront error.

The first AO retinal imager was a standard flood illuminated fundus camera [1, 2] and since then, AO has been successfully combined with other technologies such as the scanning laser ophthalmoscope (AOSLO) [3, 4] and optical coherence tomography [5]. But in comparison to astronomical AO systems, the refinement of system performance, particularly at the control system level, has not been rigorously addressed. As the number of applications increase, we can expect increases in the diversity of users and in the number of patients with more challenging optics (i.e. post-surgery, dry eyes, etc.). Improvements in system performance and robustness can significantly increase the clinical and scientific throughput (better quality images from a larger pool of patients) of an AO retinal imaging system.

Earlier control systems used in vision science AO systems often exhibited immediate clipping (saturation in one direction) and excessively long convergence times which made them impractical for clinical deployment [1, 4]. Many of these problems were patient dependent (fixation stability, tear quality, retinal reflectivity, etc.), so there is a large variability in retinal image quality among subjects which further complicates clinical investigations. The standard basis for AO control loop design, which involves modeling the plant (DM and WFS) with a static interaction matrix, was adopted somewhat recently in vision science [2]. The control law is an integrator in series with some type of inverse of the plant called the wavefront reconstructor. The interaction matrix is generated experimentally through a series of open loop measurements of each actuator’s spatial response, and either the standard pseudo-inverse or regularized inverse is used to compute the wavefront reconstructor.

We expanded on the standard AO controller design by: 1) incorporating the static nonlinear actuator response into an input linearization step and 2) implementing four different control algorithms on a AOSLO that use a MEMS DM [4]. Each algorithm optimizes a particular quadratic cost function in the design of the wavefront reconstructor and uses the standard integrator update law. We imaged two patients using each algorithm and quantified our findings using two different image quality metrics.

II. DESCRIPTION OF THE AOSLO

The control loop operates over the optical path of the AOSLO shown in Fig. (1). The infrared beam is provide by an 840 nm superluminescent diode (SLD) (Superlum Ltd., Russia) and a photomultiplier tube (PMT) (Hamamatsu, Japan) is used for light detection. The MEMS DM (Boston Micromachines Corporation, USA) has a 12 by 12 actuator array minus the corner pixels making a total of 140 inputs. Based on measurements made using a phase shifting interferometer (PSI), a single actuator has a stroke range of about 1.2 µm. The WFS is a Shack-Hartmann type with a subaperture diameter of 400 µm and a maximum frame rate of about 25 Hz. For a 6 mm diameter pupil, the wavefront is sampled at 213 locations.

Fig. 1
Schematic diagram of the AOSLO (PMT – photomultiplier tube, CP – confocal pinhole, BS – beam splitter, WFS – wavefront sensor, SLD – superluminescent diode, DM – deformable mirror, HS – horizontal ...

The infrared beam is coupled into the imaging path by a beam splitter and passes through a series of relay optics, DM, and scanning mirrors before being focused onto the retina. The reflected light from the retina returns along the same path reaching the same beam splitter. Most of the light passes through the beam splitter and is focused on the plane of the confocal pinhole. Light reaching the PMT is converted into a voltage signal that is digitized by a framegrabber to one 8-bit pixel value in the final image (512 by 512 pixels). Only a small area of the retina is illuminated at any point in time, so the pixel value depends on the amount of reflection from that area and the quality of the optics the light passes through before reaching the confocal pinhole. The PMT is synchronized with the scanning mirrors’ timing mechanism (HS and VS in Fig. (1)), so images are constructed pixel by pixel from raster scanning the retina.

The intensity point spread function (PSF) is the autocorrelation of the single-pass PSF (two-dimensional impulse response of the imaging system).

Hdp(x,y)=H(x',y')H(x'+x,y'+y)dx'dy'
(1)

Where Hdp and H are the double-pass and single-pass PSFs respectively. The double-pass PSF is imaged onto the plane of the circular confocal pinhole modeled here with a rectangle function:

c(x,y)=rect(x2+y2/D)
(2)

where D is pinhole diameter (50 to 80 µm). The PMT integrates light transmitted by the pinhole, so the value of a pixel is always proportional to the integrated intensity at that point in the raster scan [6]:

Ipixelc(x,y)Hdp(x,y)dxdy
(3)

Minimizing the residual wavefront error condenses the spatial distribution of the PSF so more light passes through the pinhole increasing the pixel value.

III. PROBLEM FORMULATION

A. AO System Loop

The wavefront from the eye is combined with the wavefront modulated by the DM surface profile to produce the residual wavefront seen at the WFS plane. According to the block diagram in Fig. (2), the error vector, which is the output of the WFS, is given by:

e(k)=H(z1)Gu(k)+H(z1)w
(4)

where u(k) and w are the input vector and the eye’s wave aberrations respectively. The WFS does not measure the wavefront directly but acquires a digital image, using a CCD camera, where wavefront gradients are estimated via an image processing algorithm. For the purpose of this study, we assume this algorithm to be sufficiently accurate. The best correction is achieved when the residual wavefront is flat or the gradient is zero. Since the DM is nearly instantaneous [7], the only plant dynamics are due to the CCD integration time of the WFS.

Fig. 2
Block diagram of the closed loop AO system. I is the 140 by 140 identity matrix and L is the wavefront reconstructor

The DM and WFS in eq. (4) are modeled by the interaction matrix T mapping the input vector to the wavefront gradients corresponding to the DM surface:

yT(k)=Tu(k1)
(5)

The one step delay is due to the integration time of the WFS. When combined with the gradients of the incoming wavefront from the eye, eq. (4) simplifies to:

e(k)=Tu(k1)+yw
(6)

where yw is the vector of sampled gradients of the eye’s aberrated wavefront. The standard integrator law employed by most AO systems is:

u(k)=u(k1)+L[e(k)+v(k)]
(7)

where L is the wavefront reconstructor matrix and v(k) is photon and CCD readout noise. From eq. (6) and eq. (7), we find the error update equation to be:

e(k+1)=[I+TL]e(k)+TLv(k)
(8)

When noise appears to be dominating the measurement signal, the CCD camera’s integration time is heuristically adjusted in real-time by the operator to detect more light from the retina increasing signal to noise ratio. The tradeoff is a reduction in temporal bandwidth which appears to be less critical than the accuracy in estimating the error vector.

B. Input linearization

The system described by eq. (6) through eq. (8) assumes the mapping between the input and measurement vectors to be linear. When the input is voltage, this model only applies to linear DMs such as the piezoelectric DMs employed in several other systems [1, 3]. For MEMS, a linearization step is required to approximate an equivalent linear system.

Deflection is achieved via electrostatic actuation, so a bias voltage needs be applied to the entire DM in order to actuate in both directions. The biased position should accommodate a maximum positive and negative single actuator deflection of equal magnitude, and the bias voltage was found using a PSI to be about 190 V. From this position, a single actuator is first released to 0 V and then driven incrementally to 265 V while measuring the deflection every 20 V (except from 260 V to 265 V). Deflections for six different actuators, averaged and normalized to range from −1 to 1, are shown in Fig. (3) with a second order polynomial fit. Defining the input vector as a set of normalized deflections better approximates the desired linear system model. The corresponding actuator voltages are found by solving for the roots of the fitted polynomial.

Fig. 3
Normalized deflection versus voltage curve. Error bars denote one standard deviation for six different actuators.

C. Interaction Matrix

Identification of the interaction matrix T is done by introducing a nearly flat wavefront into the system with a model eye (lens and a diffuse scatterer positioned at the focal point) and measuring the static response of all the actuators [8]. The DM is set at the biased state and each actuator is pushed and pulled while measurements are made with the WFS. Letting ti+ denote the gradients measured when an input of 1 (0 V) is applied to the ith actuator and ti the gradient measured when −1 (265 V) is applied to the same actuator, the ith column of the interaction matrix is:

ti=(ti+ti)/2
(9)

and the interaction matrix is defined as:

T=[t1  t2    t140]
(10)

Since there are 213 subapertures and each subaperture estimates both x and y derivatives, the dimensions of T are 426 by 140. The final model used is based the average of several generated interaction matrices.

IV. WAVEFRONT RECONSTRUCTION

A. Standard Regularization

A naive solution to the AO control problem seeks only to minimize the 2-norm of the error vector which results in a reconstructor that is the pseudo-inverse of the interaction matrix. However, the interaction matrix is ill-conditioned and the resultant reconstructor has been verified to be unstable in practice. A widely used technique for the inversion is the truncated singular value decomposition (SVD) where the smallest singular values of T are dropped. A more practical alternative to the SVD method is standard regularization where instead of dropping the smallest singular values, a constant regularization factor α is added to all the singular values [9]. This is equivalent to adding an input penalty to the standard cost function:

J(k)=e(k)22+α2u(k)22
(11)

Minimizing J with respect to u(k) obtains the reconstructor:

L=(TTT+α2I)1TT
(12)

Note that the design parameter α is the noise to signal ratio of the system.

B. Local Waffle Penalty

Waffle modes are created by driving adjacent actuators in opposite directions producing a voltage map resembling a checkerboard pattern. Patches of this pattern are often observed when SVD or standard regularization methods are used. Since they are not well sensed by the WFS, they can slowly build up in the control loop degrading retinal image quality in the process.

The suppression of waffle modes is a spatial frequency shaping problem. The following finite impulse response (FIR) filter is used to model local waffle structure [9]:

1111

which is simply a first derivative operation. By implementing this convolution operation as a matrix multiplication applied to the input vector u(k), the cost function and its corresponding reconstructor can be derived as:

J(k)=e(k)22+uT(k)[α2FTF+η2VVT]u(k)
(13)

L=(TTT+α2FTF+η2VVT)1TT
(14)

where F is the convolution matrix form of the FIR filter for local waffle. Matrix V is designed to span the nullspace of F (required for α2FTF + η2VVT [succeeds] 0). In principle, piston is the only mode that needs to be included in V but empirical observations of tip and tilt mode buildup during closed loop operation lead to their inclusion as well. The structures of piston, tip and tilt modes are shown in Fig. (4).

Fig. 4
From left to right: piston, tip and tilt input modes.

C. Kolmogorov Penalty

In astronomy, statistical models for noise and atmospheric turbulence are routinely used to optimize the design of the reconstructor matrix [9, 10]. Atmospheric turbulence is often modeled to follow Kolmogorov theory, but there is evidence that the spatial power spectrum of the eye’s wave aberrations also follow the classical Kolmogorov model [11]. Assuming the wavefront is proportional to the actuator commands, a sparse approximation for the inverse wavefront covariance matrix for the Kolmogorov model exists [12]:

Λww1CTC
(15)

where C is the convolution matrix form of the FIR filter for the discrete Laplacian operator:

010141010

Denoting the noise covariance matrix by Λvv, which we assume to be both white and constant in variance across the pupil, the cost function and its corresponding reconstructor become very similar to those for the previous design:

J(k)=eT(k)Λvv1e(k)+uT(k)[Λww1+η2VVT]u(k)=e(k)22+uT(k)[α2CTC+η2VVT]u(k)
(16)

L=(TTT+α2CTC+η2VVT)1TT
(17)

The nullspace of C is spanned by piston, tip and tilt so these modes make up the columns of V. It is worth noting that this reconstructor under integral control is the minimum variance controller for nearly ideal imaging conditions [13].

D. Zernike Polynomials

Zernike polynomials, a set of two dimensional polynomials that are orthogonal over the unit circle, are almost exclusively used to quantify the eye’s wave aberrations [11, 14, 15]. If a finite set of Zernike modes can accurately represent the eye’s aberrations, projecting the input vector onto a Zernike spanned subspace should improve system robustness during less-than-ideal experimental conditions. Furthermore, this method can allow for the correction of certain Zernike modes (i.e. defocus, astigmatism, etc.) while leaving all others intact which may be useful for certain applications such as vision performance testing.

An input vector defined by the first N Zernike modes can be described by:

u=i=0N1ciZi=Zc
(18)

Where ci and Zi are the ith Zernike coefficient and mode vector (evaluated at each actuator point) respectively, and the matrices are defined as:

Z=[Z3  Z4    ZN1]140×N
(19)

c=[c3  c4    cN1]T140
(20)

Z0 through Z2 correspond to piston, tip and tilt and are therefore not included in Z. N is 66 for the 10th order correction used in this study. Substituting eq. (18) into eq. (6) obtains the relationship between the error vector and the Zernike coefficient vector:

e(k)=TZc(k1)+yw
(21)

The Zernike polynomial reconstructor minimizes the cost function with respect to the Zernike coefficient vector:

J(k)=e(k)22+α2c(k)22
(22)

L=Z[(TZ)T(TZ)+α2I]1(TZ)T
(23)

It is worth noting from an implementation standpoint that the Zernike polynomials lose their orthogonality when discretized and extrapolated over a nearly square actuator array. For this reason, matrix Z is ill-conditioned, so the Gram-Schmidt procedure is applied to the columns of Z before evaluating eq. (23).

V. STABILITY

Stability for the closed loop system described by eqs. (7) and (8) can be addressed by analyzing the behavior of the Lyapunov function:

V(e(k+1))=eT(k+1)e(k+1)=eT(k)[I+TL]T[I+TL]e(k)
(24)

Letting W be the weighting matrix on u(k) (c(k) for the Zernike polynomial case), it can be shown that:

V(e(k+1))=eT(k)e(k)e(k)LT(TTT+2W)Le(k)ΔV(e(k))=eT(k)LT(TTT+2W)Le(k)

Matrix W must be positive definite because the lack of weighting on input u(t) would create unbounded input magnitudes leading to actuator saturation. For the standard regularization and Zernike polynomial cases, W was proportional to identity so positive definiteness was trivial. For the other two designs, we manually identified specific modes that needed to be explicitly penalized in order to establish positive definiteness. It follows immediately that:

TTT+2W0
(25)

But since L is a weighted inverse of T, it cannot have full column rank (L [set membership] R140×426). Therefore,

LT(TTT+2W)L0ΔV(e(k))0
(26)

so the system is stable in the sense of Lyapunov regardless of which reconstructor is used. However, this is enough to guarantee that the control signal and wavefront error signals do not go unbounded.

A more thorough stability analysis would require accurate modeling of electrostatic actuation coupled with the membrane deformation properties of the DM. A mathematical model of the type of MEMS device used in this study has been assessed [16]. It was not adopted for this study because the model’s predicted membrane response did not match the actual membrane response of our MEMS device at high voltages.

VI. RESULTS

The imaging sessions were kept short (~20 seconds) and administered minutes apart to minimize subject fatigue which may bias the comparisons. The center of the raster scan was placed approximately 0.4 degrees outside of the subjects preferred fixation point. The scanning was performed over approximately 0.9 degrees retinal eccentricity which corresponds to about 0.265 mm for subject 1 and 0.279 mm for subject 2 since eye sizes differ among individuals. The two image quality metrics used to quantify system performance are the root-mean-square (RMS) wavefront error and the mean pixel value of the retinal image. The RMS wavefront error was computed in real-time and logged for each experiment. Mean pixel values were computed offline.

The sampled wavefront gradients from output vector e(k) were integrated to estimate the two-dimensional wavefront maps used to compute the RMS wavefront error [17]. According to Fig. (5), the RMS error converged in all four cases for both subjects and remained near the best corrected state over the entire imaging session. However, performance traces based on the mean pixel value of each acquire frame, shown in Fig. (6), displayed a much greater level of diversity. The retinal image pixel values are more direct indicators of image quality since certain aberration profiles from either the DM or the eye can lie beyond the sampling capabilities of the WFS.

Fig. 5
Performance based on the estimated RMS of the residual wavefront error for subject 1 (left) and subject 2 (right). SR – standard regularization; LWP – local waffle penalty; KP – Kolmogorov penalty; ZP – Zernike polynomials ...
Fig. 6
Performance based on the mean pixel value of the retinal images

The temporal mean and standard deviation for both image quality metrics during closed loop operation are tabulated in Fig. (7). The results suggest that the Kolmogorov penalty reconstructor modeled the power spectrum of the eye’s wave aberrations most accurately out of the four tested reconstructors. The Zernike polynomial reconstructor performed the worst overall for the two subjects in this study, but we have experienced imaging conditions where measurements were poor, and only the Zernike polynomial reconstructor was robust enough to make an effective correction.

Fig. 7
Temporal mean and standard deviation of the RMS wavefront error and the mean pixel value beginning from about two seconds after closing the control loop

A set of acquired cone mosaic images are shown in Fig. (8). They correspond to the data in Fig. (5) and Fig. (6) for subject 2, stabilized and frame averaged (100 manually selected frames) to improve contrast. The round spots with varying intensities packed in a nearly hexagonal array are cone photoreceptors. For this particular subject, these features are noticeably more blurry in images acquired using the Zernike polynomials and the local waffle penalty reconstructors, which is consistent with the observed mean pixel value traces. Cones become much smaller and more tightly packed at the fovea center [18], and further system refinements are needed to reliably resolve them.

Fig. 8
Examples of AOSLO images from subject 2 acquired with each of the four control algorithms (top left: SR, top right: LWP, bottom left: KP and bottom right: ZP). Each image subtends from approximately 0 (fovea center) to 0.8 degrees (0.25 mm) eccentricity. ...

VII. CONCLUSION

The work presented in this paper marks the first step to improving the resolution of a vision science AO imaging systems by using more advanced controllers; an important design component that is often overlooked based on relevant literature. We have addressed the AO control problem focusing primarily on implementation and testing and demonstrated that sharper images of the human cone mosaic can be obtained by improving the control system. For countering the static nonlinear properties of the MEMS device, we have added an input linearization step based on deflection measurements made using a PSI. A critical design component of the integrator type controller for AO systems is the wavefront reconstructor. Four reconstructors based on the optimization of quadratic cost functions were described and tested on two real eyes. In addition to the standard regularization method, we considered two designs using spatial frequency shaped DM modes and also a modal reconstructor based on a finite set of Zernike polynomials [9, 12, 19].

Two quantitative image quality metrics were used to evaluate the performance of the control algorithms: 1) RMS residual wavefront error and the 2) mean pixel value of the acquired retinal image. Even though the four controllers performed similarly according to the computed RMS wavefront error, they did not all produce retinal images of similar quality. The mean pixel value is a more sensitive indicator of retinal image quality because it is directly related to the system’s PSF. The Kolmogorov penalty reconstructor performed the best according to the mean pixel value suggesting that it makes a reasonable approximation of the statistics of the eye’s wave aberrations. Even though the Zernike polynomial reconstructor did not outperform the other methods in most cases, it can be a practical alternative during less than ideal imaging conditions.

ACKNOWLEDGMENT

The authors thank Yuhua Zhang for the development of the MEMS-based AOSLO, Curtis Vogel for answering many mathematics questions, and the two individual who volunteered to be imaged.

This work was supported by the NIH Bioengineering Research Partnership Grant EY014375 and by the NSF Science and Technology Center for Adaptive Optics, managed by the University of California at Santa Cruz under cooperative agreement AST-9876783

Contributor Information

Kaccie Y. Li, vision science group in the School of Optometry at the University of California at Berkeley, USA.

Sandipan Mishra, department of Mechanical and Industrial Engineering at the University of Illinois at Urbana-Champaign.

Pavan Tiruveedhula, vision science group in the School of Optometry at the University of California at Berkeley, USA.

Austin Roorda, vision science group in the School of Optometry at the University of California at Berkeley, USA.

REFERENCES

1. Liang JZ, Williams DR, Miller DT. Supernormal vision and high-resolution retinal imaging through adaptive optics. JOSA A. 1997 Nov;vol. 14:2884–2892. [PubMed]
2. Hofer H, Chen L, Yoon GY, Singer B, Yamauchi Y, Williams DR. Improvement in retinal image quality with dynamic correction of the eye's aberrations. Optics Express. 2001 May 21;vol. 8:631–643. [PubMed]
3. Roorda A, Romero-Borja F, Donnelly WJ, Queener H, Hebert TJ, Campbell MCW. Adaptive optics scanning laser ophthalmoscopy. Optics Express. 2002 May 6;vol. 10:405–412. [PubMed]
4. Zhang YH, Poonja S, Roorda A. MEMS-based adaptive optics scanning laser ophthalmoscopy. Optics Letters. 2006 May 1;vol. 31:1268–1270. [PubMed]
5. Hermann B, Fernandez EJ, Unterhuber A, Sattmann H, Fercher AF, Drexler W, Prieto PM, Artal P. Adaptive-optics ultrahigh-resolution optical coherence tomography. Optics Letters. 2004 Sep 15;vol. 29:2142–2144. [PubMed]
6. Venkateswaran K, Roorda A, Romero-Borja F. Theoretical modeling and evaluation of the axial resolution of the adaptive optics scanning laser ophthalmoscope. JBO. 2004 Jan–Feb;vol. 9:132–138. [PubMed]
7. Bifano T, Bierden P, Perreault J. Micromachined deformable mirrors for dynamic wavefront control. Proceedings of SPIE. 2004;vol. 5553:1–16.
8. Boyer C, Michau V, Rousset G. Adaptive optics: interaction matrix measurement and real-time control algorithms for the Come-On project. Proceedings of SPIE. 1990;vol. 1542:46–61.
9. Gavel D. Suppressing anomalous localized waffle behavior in least squares wavefront reconstructors. Proceedings of the SPIE. 2002;vol. 4839:972–980.
10. van Dam MA, Le Mignant D, Macintosh BA. Performance of the Keck Observatory adaptive-optics system. Applied Optics. 2004 Oct 10;vol. 43:5458–5467. [PubMed]
11. Cagigal MP, Canales VF, Castejon-Mochon JF, Prieto PM, Lopez-Gil N, Artal P. Statistical description of wave-front aberration in the human eye. Optics Letters. 2002 Jan 1;vol. 27:37–39. [PubMed]
12. Ellerbroek BL. Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques. JOSA A. 2002 Sep;vol. 19:1803–1816. [PubMed]
13. Looze D. Minimum variance control structure for adaptive optics systems. Proceedings of the 2005 American Control Conference; June 8–10 2005; pp. 1466–1471.
14. Porter J, Guirao A, Cox IG, Williams DR. Monochromatic aberrations of the human eye in a large population. JOSA A. 2001 Aug;vol. 18:1793–1803. [PubMed]
15. Castejon-Mochon JF, Lopez-Gil N, Benito A, Artal P. Ocular wave-front aberration statistics in a normal young population. Vision Research. 2002 Jun;vol. 42:1611–1617. [PubMed]
16. Vogel CR, Yang Q. Modeling, simulation, and open-loop control of a continuous facesheet MEMS deformable mirror. JOSA A. 2006 May;vol. 23:1074–1081. [PubMed]
17. Herrmann J. Least-Squares Wave-Front Errors of Minimum Norm. JOSA. 1980;vol. 70:28–35.
18. Curcio CA, Sloan KR, Kalina RE, Hendrickson AE. Human Photoreceptor Topography. Journal of Comparative Neurology. 1990 Feb 22;vol. 292:497–523. [PubMed]
19. van Dam M. In: The reconstruction matrix. Li KY, editor. Kamuela: 2003.