HYPERSPECTRAL IMAGING-BASED SPATIALLY-RESOLVED TECHNIQUE FOR ACCURATE MEASUREMENT OF THE OPTICAL PROPERTIES OF HORTICULTRAL PRODUCTS By Haiyan Cen A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Biosystems Engineering 2011 ABSTRACT HYPERSPECTRAL IMAGING-BASED SPATIALLY-RESOLVED TECHNIQUE FOR ACCURATE MEASUREMENT OF THE OPTICAL PROPERTIES OF HORTICULTRAL PRODUCTS By Haiyan Cen Hyperspectral imaging-based spatially-resolved technique is promising for determining the optical properties and quality attributes of horticultural and food products. However, considerable challenges still exist for accurate determination of spectral absorption and scattering properties from intact horticultural products. The objective of this research was, therefore, to develop and optimize hyperspectral imaging-based spatially-resolved technique for accurate measurement of the optical properties of horticultural products. Monte Carlo (MC) simulations and experiments for model samples of known optical properties were performed to optimize the inverse algorithm of a single-layer diffusion model and the optical designs, for extracting the absorption (µa) and reduced scattering (µs') coefficients from spatially-resolved reflectance profiles. The logarithm and integral data transformation and the relative weighting methods were found to greatly improve the parameter estimation accuracy with the relative errors of 10.4%, 10.7%, and 11.4% for µa, and 6.6%, 7.0%, and 7.1% for µs', respectively. More accurate measurements of optical properties were obtained when the light beam was of Gaussian type with the diameter of less than 1 mm, and the minimum and maximum source-detector distances were 1.5 mm and 10-20 transport mean free paths, respectively. An optical property measuring prototype was built based on the optimization results. The instrument was used to measure the optical properties, and assess quality/maturity, of 500 ‘Redstar’ peaches and 1039 ‘Golden Delicious’ (GD) and 1040 ‘Delicious’ (RD) apples. A separate study was also conducted on confocal laser scanning and scanning electron microscopic image analysis and compression test of fruit tissue specimens to measure the structural and mechanical properties of GD and ‘Granny Smith’ (GS) apples under accelerated softening at high temperature (22 ºC)/high humidity (95%) for up to 30 days. The absorption spectra of peach and apple fruit were featured with the absorption peaks of major pigments (i.e., chlorophylls and anthocyanin) and water, while the reduced scattering coefficient generally decreased with the increase of wavelength. Partial least squares regression resulted in various levels of correlation of µa and µs' with the firmness, soluble solids content, and skin and flesh color parameters of peaches (r = 0.204-0.855) and apples (r = 0.460-0.885), and the combination of the two optical parameters generally gave higher correlations (up to 0.893). The mean value of µa and µs' for GD and GS apples for each storage date was positively correlated with acoustic/impact firmness, Young’s modulus, and cell parameters (r = 0.585-0.948 for GD and r = 0.292-0.993 for GS). A two-layer diffusion model for determining the optical properties of fruit skin and flesh was further investigated through solid model samples. The average errors of determining two and four optical parameters were 6.8% and 15.3%, respectively, for MC reflectance data. The errors of determining the first layer of model samples were approximately 23.0% for µa and 18.4% for µs', indicating the difficulty and also potential in applying the two-layer diffusion model for fruit. This research has demonstrated the usefulness of hyperspectral imaging-based spatiallyresolved technique for determining the optical properties and maturity/quality of fruits. However, further research is needed to reduce measurement variability or error caused by irregular or rough surface of fruit and the presence of fruit skin, and apply the technique to other foods and biological materials. This work is dedicated to my husband Dahai Gao and my son Nicholas Jiatong Gao iv    ACKNOWLEDGMENTS I could not have reached the height in my life without tremendous help and support from the others. I sincerely thank my advisor and major professor, Dr. Renfu Lu, for his guidance, encouragement, understanding, patience, and friendship throughout my Ph.D. study at Michigan State University, which was essential to the completion of this research. I would also like to warmly thank my committee members, Dr. Daniel E. Guyer, co-advisor, for his help on my course and research work and valuable comments on improvement of my presentation skills, and Dr. Randolph Beaudry and Dr. Lalita Udpa for their time and enormously helpful advice on my research project. My thanks also go to my talented colleagues and friends  Dr. Diwan P. Ariana, Dr. Fernando Mendoza, Dr. Akira Mizushima, Mr. Benjamin Bailey, Mr. Irwin Donis Gonzalez, and Mr. Ahmed Rady  for their help, support and advice in my Ph.D. project. In particular, I would like to acknowledge Dr. Diwan P. Ariana and Mr. Benjamin Bailey for their significant contributions to the development of the optical property measuring instrument, and Dr. Fernando Mendoza for his help on image analysis. I am also grateful to Dr. Kirk Dolan and Dr. James V. Beck for their contributions in solving inverse light transport problems. I especially would like to thank Dr. Fred W. Bakker-Arkema, Dr. James F. Steffe, and Dr. Ajit K. Srivastava for their help during my Ph.D. study at Michigan State University. I also thank the Department of Biosystems and Agricultural Engineering and the College of Agriculture and Natural Resources for providing me with a fellowship to complete this dissertation. Last but not least, I would like to thank my husband Dahai Gao for his intellectual help, moral support, and deep love. My deepest gratitude and love go to my parents, Qingmiao Cen, v    and Meidi Tu, for their dedication and the many years of support during my study. Thanks also go to my parents-in-law, Feng Gao and Yang Yang, and my sister, Jie Cen. vi    TABLE OF CONTENTS LIST OF TABLES .........................................................................................................................x  LIST OF FIGURES .................................................................................................................... xii  KEY TO SYMBOLS OR ABBREVIATIONS ...................................................................... xviii  CHAPTER 1  INTRODUCTION ............................................................................................1  1.1  Background .......................................................................................................................1  1.2  Objectives ..........................................................................................................................6  CHAPTER 2  LITERATURE REVIEW ................................................................................7  2.1  Principles of Light Interaction with Biological Materials .................................................7  2.1.1  Absorption Coefficient...............................................................................................7  2.1.2  Scattering Coefficient ................................................................................................9  2.1.3  Anisotropy Factor ......................................................................................................9  2.1.4  Reduced Scattering Coefficient ...............................................................................10  2.1.5  Refractive Index .......................................................................................................11  2.2  Modeling of Light Transport in Tissues ..........................................................................12  2.2.1  Radiation Transport Model ......................................................................................12  2.2.2  Diffusion Model .......................................................................................................14  2.2.3  Adding-doubling Model...........................................................................................15  2.2.4  Monte Carlo Simulation ...........................................................................................17  2.3  Techniques for Measuring Optical Properties of Biological Materials...........................19  2.3.1  Time-resolved Technique ........................................................................................19  2.3.2  Frequency-domain Technique .................................................................................22  2.3.3  Spatially-resolved Technique...................................................................................25  2.4  Approaches to Inverse Problems for Determining Optical Parameters ..........................31  2.4.1  Sensitivity Analysis .................................................................................................31  2.4.2  Inverse Algorithm ....................................................................................................32  2.4.3  Data Transformation and Weighting Methods ........................................................33  2.5  Nondestructive Techniques for Measurement of Maturity and Quality .........................36  2.5.1  Maturity and Quality of Fruits .................................................................................36  2.5.2  Mechanical Techniques ...........................................................................................38  2.5.3  Optical Techniques ..................................................................................................40  CHAPTER 3  OPTIMIZATION OF INVERSE ALGORITHM FOR ESTIMATING OPTICAL PROPERTIES OF BIOLOGICAL MATERIALS ................................................46  3.1  Introduction .....................................................................................................................46  3.2  Forward Problem .............................................................................................................48  3.2.1  One-layer Diffusion Model for the Steady-state Case .............................................48  3.2.2  Monte Carlo Simulation ...........................................................................................49  3.3  Inverse Problem...............................................................................................................50  3.3.1  Nonlinear Least Squares Inverse Algorithm ............................................................50  vii    3.3.2  Data Transformation and Weighting Methods ........................................................51  3.3.3  Sensitivity Analysis .................................................................................................51  3.3.4  Statistical Analysis and Model Assessment .............................................................52  3.4  Simulation Experiments ..................................................................................................53  3.5  Results and Discussion ....................................................................................................55  3.5.1  Comparison between the Diffusion Model and Monte Carlo Simulations ..............55  3.5.2  Sensitivity Coefficients ............................................................................................57  3.5.3  Estimation of Optical Parameters based on Monte Carlo Simulation Data .............60  3.5.4  Statistical Analysis ...................................................................................................63  3.6  Conclusions .....................................................................................................................68  CHAPTER 4  OPTIMIZATION OF THE HYPERSPECTRAL IMAGING-BASED SPATIALLY-RESOLVED SYSTEM ........................................................................................70  4.1  Introduction .....................................................................................................................70  4.2  Materials and Methods ....................................................................................................72  4.2.1  Hyperspectral Imaging-based Spatially-resolved Technique ..................................72  4.2.2  Model Samples Preparation .....................................................................................75  4.2.3  Reference Methods ..................................................................................................76  4.2.4  Optimization of the Light Beam and Source-detector Distance ..............................80  4.2.5  Assessment for Accuracy, Precision/Reproducibility, and Sensitivity ....................82  4.3  Results and Discussion ....................................................................................................83  4.3.1  Accuracy of the Reference Methods ........................................................................83  4.3.2  Light Beam Characteristics ......................................................................................85  4.3.3  Source-detector Distance .........................................................................................89  4.3.4  Accuracy, Precision/Reproducibility, and Sensitivity .............................................94  4.4  Conclusions ...................................................................................................................100  CHAPTER 5  DEVELOPMENT OF AN OPTICAL PROPERTY MEASUREMENT PROTOTYPE…………….........................................................................................................101  5.1  Introduction ...................................................................................................................101  5.2  System Setup and Software Development ....................................................................102  5.2.1  Optical Property Analyzer Setup ...........................................................................102  5.2.2  Optical Property Analyzer Software Development ...............................................105  5.3  System Calibration and Evaluation ...............................................................................110  CHAPTER 6  MATURITY/QUALITY ASSESSMENT FOR APPLE AND PEACH FRUIT USING OPTICAL PROPERTIES..............................................................................114  6.1  Introduction ...................................................................................................................114  6.2  Materials and Methods ..................................................................................................115  6.2.1  Fruit Samples .........................................................................................................115  6.2.2  Optical Properties Measurement ............................................................................116  6.2.3  Acoustic Measurement...........................................................................................118  6.2.4  Destructive Maturity/Quality Measurement ..........................................................119  6.2.5  Data Analysis .........................................................................................................120  6.3  Results and Discussion ..................................................................................................121  viii    6.3.1  Statistical Data of Fruit Maturity/Quality ..............................................................121  6.3.2  Peach Maturity/Quality Evaluation........................................................................123  6.3.3  Apple Internal Quality Evaluation .........................................................................131  6.4  Conclusions ...................................................................................................................135  CHAPTER 7  RELATIONSHIP BETWEEN THE OPTICAL AND MECHANICAL/STRUCTURAL PROPERTIES OF APPLE TISSUE ..............................137  7.1  Introduction ...................................................................................................................137  7.2  Materials and Methods ..................................................................................................138  7.2.1  Apple Samples .......................................................................................................138  7.2.2  Optical Properties Measurement ............................................................................139  7.2.3  Acoustic Firmness Measurement and Compression Test for Measuring Mechanical Properties .........................................................................................................................139  7.2.4  Confocal Laser Scanning Microscopy and Scanning Electron Microscopy for Microstructural Analysis ..................................................................................................141  7.2.5  Data Analysis .........................................................................................................142  7.3  Results and Discussion ..................................................................................................146  7.3.1  Changes in the Optical Properties of Apples during Storage ................................146  7.3.2  Changes in Acoustic Firmness and Mechanical Properties during Storage ..........149  7.3.3  Microstructural Changes in Apple Tissue during Storage .....................................152  7.3.4  Correlations between Optical Properties and Mechanical/Structural Properties ...159  7.4  Conclusions ...................................................................................................................162  CHAPTER 8  DETERMINATION OF THE OPTICAL PROPERTIES OF TWOLAYER TURBID MATEIRALS ..............................................................................................164  8.1  Introduction ...................................................................................................................164  8.2  Materials and Methods ..................................................................................................165  8.2.1  Two-layer Diffusion Model ...................................................................................165  8.2.2  Model Samples Preparation and Measurement......................................................167  8.3  Validation of the Two-layer Diffusion Model and Inverse Algorithm .........................168  8.4  Results and Discussion ..................................................................................................170  8.4.1  Monte Carlo Simulations of Spatially-resolved Diffuse Reflectance ....................170  8.4.2  Sensitivity Coefficient Analysis ............................................................................171  8.4.3  Extraction of Optical Properties from Monte Carlo Simulation Data ...................174  8.4.4  Optical Properties of Model Samples Measured with Integrating Sphere .............178  8.4.5  Optical Properties of Model Samples Determined from Hyperspectral Imaging Measurements ..................................................................................................................179  8.5  Conclusions ...................................................................................................................182  CONCLUSION AND FUTURE RESEARCH ........................................................................184  BIBLIOGRAPHY ......................................................................................................................189  ix    LIST OF TABLES Table 2.1 Summary of the performance of time-resolved, frequency-domain and spatiallyresolved measurements (TS: time-resolved, FD: frequency-domain, SR: spatiallyresolved, POM: polyoxymethylene). .......................................................................... 29  Table 3.1 Thirty-six combinations of the absorption (µa) and reduced scattering coefficients (µs') and their corresponding transport mean free path ( mfp' ) used in Monte Carlo -1 simulations (unit: cm for µa & µs', and mm for mfp'). ............................................ 54  Table 3.2 Statistical results for estimating the optical parameters using the LTDM method. ..... 67  Table 5.1 Components of the Optical Property Analyzer control software. ............................... 106  Table 6.1 Statistics of fruit maturity/quality parameters for 500 ‘Redstar’ peaches measured by standard methods. ..................................................................................................... 122  Table 6.2 Statistics of the firmness (i.e., maximum Magness-Taylor force or FM) and soluble solids content (SSC) of ‘Golden Delicious’ (GD) and ‘Delicious’ (RD) apples for the freshly harvested, after-storage, and combined groups. ........................................... 123  Table 6.3 Partial least squares (PLS) prediction results of ‘Redstar’ peach maturity/quality parameters using absorption coefficient (µa), reduced scattering coefficient (µs') and their three combinations (µa & µs', µa × µs', and µeff).* .......................................... 129  Table 6.4 Firmness prediction results for ‘Golden Delicious’ and ‘Delicious’ apples for the freshly harvested, after-storage and combined groups.* .......................................... 132  Table 6.5 Prediction results for the soluble solids content of ‘Golden Delicious’ and ‘Delicious’ apples for the freshly harvested, after-storage and combined groups.* .................... 133  Table 7.1 Mean and standard error of the cell size/shape parameters for ‘Golden Delicious’ and ‘Granny Smith’ apple tissues for different storage days.* ........................................ 154  Table 7.2 Morphological features extracted from the scanning electron microscopic (SEM) images acquired at the low resolution for ‘Golden Delicious’ and ‘Granny Smith’ apple tissues.* ........................................................................................................... 158  Table 7.3 Correlations of selected optical parameters with acoustic and impact firmness and Young’s modulus for ‘Golden Delicious’ and ‘Granny Smith’ apples.* ................. 160  Table 7.4 Correlations of optical parameters with the cell size parameters extracted from the confocal laser scanning microscopic images of ‘Golden Delicious’ and ‘Granny Smith’ apples. ........................................................................................................... 161  x    Table 7.5 Correlation of optical parameters with the fractal analysis parameters extracted from the scanning electron microscopic images of ‘Golden Delicious’ and ‘Granny Smith’ apples. ....................................................................................................................... 161  Table 8.1 Optical properties of two two-layer model samples at the wavelengths of 535 nm and 700 nm determined by the integrating sphere and the adding-doubling method (Prahl et al., 1993). .............................................................................................................. 179  xi    LIST OF FIGURES Figure 2.1 Light interaction with matter: (a) absorption, (b) scattering, (c) scattering event with anisotropy factor g, and (d) refraction (For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this dissertation). .................................................................................................................. 8  Figure 2.2 Light incident at an angle v that is reflected and transmitted by a slab at an angle v ' in the adding-doubling method (Welch and van Gemert, 1995)..................................... 17  Figure 2.3 Measurement principle of (a) time-resolved and (b) frequency-domain techniques. . 20  Figure 2.4 Measurement principle for spatially-resolved technique. ........................................... 26  Figure 2.5 Force/deformation (F/D) curves for cylindrical tissue specimens from firm and soft apples under constant strain-rate (0.42 mm/s) compression (Abbott, 1999). ............. 39  Figure 2.6 Diagram of the whole electromagnetic spectrum (Mohsenin, 1984). ......................... 41  Figure 2.7 Three different modes for NIR spectroscopic measurements: a) reflectance, b) transmittance and c) interactance. ............................................................................... 42  Figure 3.1 Comparison of the spatially-resolved diffuse reflectance obtained from the diffusion model (symbols) and Monte Carlo simulation (solid lines): (a) µs'/ µa = 5, and (b) µs'/ µa = 100. ..................................................................................................................... 56  -1 - Figure 3.2 Scaled sensitivity coefficients of the optical parameters [µa = 0.06 cm & µs' = 4 cm 1 -1 -1 for the left pane of plots (a1, b1, c1, d1), and µa = 0.57 cm & µs' = 40 cm for the right pane of plots (a2, b2, c2, d2)] as functions of the source-detector distance for ODM (a1, a2), LTDM (b1, b2), ITDM (c1, c2), and RWDM (d1, d2). Solid curves are for R, dash curves for µa and dot curves for µs' (a.u. = arbitrary unit)................. 58  Figure 3.3 Relative errors of estimating 29 groups of µa (a) and µs' (b) shown in Table 3.1 by the original model, and the data transformation and relative weighting methods: ODM (·), LTDM (◦), ITDM (*) and RWDM (∆). ...................................................................... 61  Figure 3.4 Absolute values of the relative errors of estimating the optical properties (µa = 0.40 -1 -1 cm , µs' = 20 cm ) for different noise levels (each noise level includes 10 replications). Solid circles represent the absorption coefficient, while open circles are for the reduced scattering coefficient. ......................................................................... 63  Figure 3.5 Residual plots for the reflectance data versus source-detector distances from (a) ODM, -1 -1 (b) LTDM, (c) ITDM, and (d) RWDM with µa = 0.06 cm & µs' = 4 cm (dot curve) -1 -1 and µa = 0.57 cm & µs' = 40 cm (asterisk curve). ................................................. 64  xii    Figure 3.6 Residual histograms for the reflectance data from ODM (a), LTDM (b), ITDM (c), -1 -1 and RWDM (d) with 1) µa = 0.06 cm & µs' = 4 cm for the left pane of plots (a1, -1 -1 b1, c1, d1), and 2) µa = 0.57 cm & µs' = 40 cm for the right pane of plots (a2, b2, c2, d2). ........................................................................................................................ 65  Figure 3.7 3-D plot of sum of squares for (a) group 25 and (b) group 30 using the LTDM method. ..................................................................................................................................... 67  Figure 4.1 Hyperspectral imaging-based spatially-resolved method: (a) schematic showing the major components of the system; (b) top view of multi-line scanning mode for acquiring spatially-resolved reflectance profiles. ....................................................... 73  Figure 4.2 Hyperspectral reflectance image of a liquid model sample: (a) 2-dimensional (2-D) display with intensities being indicated by pseudo colors, and (b) original spatiallyresolved reflectance profiles at 570 nm and 700 nm. ................................................. 75  Figure 4.3 Collimated transmittance measurement for absorption using a miniature fiber optic spectrometer. ............................................................................................................... 77  Figure 4.4 Single integrating sphere configurations for: (a) total diffuse reflectance measurement, and (b) total transmittance measurement. ................................................................... 78  Figure 4.5 Spectra of the absorption coefficient for the three standard solutions measured by the transmittance method [asterisks () for standard 0.1, triangles () for standard 0.5, and squares () for standard 0.8] and their true vlaues (solid lines). ........................ 84  Figure 4.6 Average differences of the reduced scattering coefficient for the model samples at different concentrations of Intralipid for 500-900 nm obtained from the integrating sphere measurement and the empirical equation. ....................................................... 85  Figure 4.7 Comparison of spatially-resolved reflectance produced by an infinitely small beam, a flat beam with diameter of d = 1 mm, and a Gaussian beam with d = 1 and 2 mm. .. 87  Figure 4.8 Error analysis for different beam sizes: (a) average errors of estimating six sets of -1 -1 optical properties with 0.060 ≤ µa ≤ 2.000 cm and 4.0 ≤ µs' ≤40.0 cm , and the ratios of µs'/µa = 20 and µs'/µa = 70; and (b) and (c) relative errors for three sets of optical properties with increased µa and µs' from A to C with µa = 0.006, 0.029, 0.057 -1 -1 cm and µs' = 0.4, 2.0, 4.0 cm . ............................................................................... 87  Figure 4.9 three-dimensional profiles (a and b) and 2-dimensional (c and d) contours of the incident light beam at wavelengths of 650 nm and 950 nm, where D1 is the direction along the scan line and D2 is perpendicular to the scan line. ..................................... 89  Figure 4.10 Relative errors of estimating µa (squares) and µs' (asterisks) from the spatiallyresolved reflectance data generated by Monte Carlo simulations when using different xiii    -1 minimum and maximum source-detector distances: (a1) and (b1) µa = 0.290 cm -1 -1 -1 & µs' = 20.0 cm , and (a2) and (b2) µa = 0.430 cm & µs' = 30.0 cm . .................... 90  Figure 4.11 Signal-to-noise ratio measurement of the hyperspectral imaging system: (a) average spatially-resolved reflectance profile of 10 measurements of a model sample at 650 nm, and (b) signal-to-noise ratio of the measurements within 10 mm source-detector distance. ...................................................................................................................... 92  -1 -1 Figure 4.12 Errors of estimating µa = 1.00 cm and µs' = 20.0 cm introduced by different spatial resolutions relative to the optical properties obtained for the resolution of 0.01 mm. ............................................................................................................................. 93  Figure 4.13 Spectra of absorption and reduced scattering coefficients of three model samples with (a) blue dye, (b) green dye, and (c) mixed dye as absorbers measured by the hyperspectral imaging and reference methods............................................................ 95  Figure 4.14 Confocal laser scanning microscopy images of Intralipid solutions with (a) blue dye, and (b) green dye. ....................................................................................................... 97  Figure 4.15 Coefficient of variation versus reordered ascending absorption coefficients of a model sample at different wavelengths....................................................................... 99  Figure 5.1 Optical Property Analyzer (or OPA) for measuring optical properties and acquiring hyperspectral images................................................................................................. 102  Figure 5.2 Schematic of the Optical Property Analyzer (or OPA) for measuring the optical properties of biological materials.............................................................................. 103  Figure 5.3 Main window of the Optical Property Analyzer software......................................... 107  Figure 5.4 Display window of the Optical Property Analyzer for the image acquisition setup. 108  Figure 5.5 Display window of the Optical Property Analyzer for optical properties computation. ................................................................................................................................... 109  Figure 5.6 Display windows of the Optical Property Analyzer for (a) calibration information and (b) sample size information. ..................................................................................... 110  Figure 5.7 Reproducibility of the Optical Property Analyzer for measuring absorption and reduced scattering coefficients of a model sample at 555 nm at each measurement day with respect to the average value calculated over 5 days. ........................................ 112  Figure 5.8 Coefficient of variation versus reordered ascending absorption coefficients of a model sample at different wavelengths................................................................................ 113  Figure 6.1 Hyperspectral reflectance image and optical property spectra of a peach sample: (a) 2D display of the original reflectance image, (b) a raw spectrum extracted for a xiv    specific scattering distance, (c) a spatially-resolved reflectance profile extracted for a selected wavelength, (d) pre-processed or averaged spatially-resolved reflectance profile at the selected wavelength, (e) the spectrum of absorption coefficient (µa), and (f) the spectrum of reduced scattering coefficient (µs'). ........................................... 117  Figure 6.2 Spectra of (a) absorption coefficient for three peaches with different skin colors (light red, red, dark red), and of (b) reduced scattering coefficient for three peaches at different firmness levels (hard, medium, soft). ......................................................... 124  Figure 6.3 Distributions of (a) absorption (µa) and (b) reduced scattering coefficient (µs') for 500 ‘Redstar’ peaches at four wavelengths (525 nm for anthocyanin, 620 nm for chlorophyll-b, 675 nm for chlorophyll-a, and 970 nm for water). ........................... 127  Figure 6.4 (a) Prediction of Magness-Taylor (MT) firmness (FM) of ‘Redstar’ peaches using effective attenuation coefficient (µeff) with partial least squares (PLS) and (b) correlation between acoustic data and MT firmness (FM)....................................... 130  Figure 6.5 Spectra of (a) absorption and (b) reduced scattering coefficients for four ‘Golden Delicious’ (GD) apples and four ‘Delicious’ (RD) apples. ...................................... 131  Figure 6.6 Prediction of fruit firmness (a, c) and soluble solids content or SSC (b, d) using the best combinations of µa and µs' for the combined group of ‘Golden Delicious’ (GD) and ‘Delicious’ (RD) apples. .................................................................................... 134  Figure 7.1 Experimental procedures for compression test of apple tissue specimens: (a) cylindrical specimen taken from the fruit flesh of the apple; (b) cork borer and double-blade knife used for cutting tissue specimens; (c) final tissue specimen with 13.9 × 13 mm (D × L); and (d) compression test setup with the Texture Analyzer. 140  Figure 7.2 Procedures of image processing for extracting cell parameters. ............................... 143  Figure 7.3 Mean spectra of 10 apple samples each for the five times of storage: (a1) absorption and (a2) reduced scattering coefficients of ‘Golden Delicious’ (GD) apples; and (b1) absorption and (b2) reduced scattering coefficients of ‘Granny Smith’ (GS) apples. ................................................................................................................................... 147  Figure 7.4 Changes in a) acoustic firmness (FI) and impact firmness (IF), and (b) fruit weight, for the same 10 fruit each of ‘Golden Delicious’ (GD) and ‘Granny Smith’ (GS) measured at five storage times (the vertical bars denote two standard deviations). . 150  Figure 7.5 Stress/strain curves obtained from the compression test for the flesh specimens of (a) five ‘Golden Delicious’ apples and (b) five ‘Granny Smith’ apples, one for each storage time. .............................................................................................................. 150  Figure 7.6 Mean values of Young’s modulus for 10 ‘Golden Delicious’ (GD) (columns with capital letters) and ‘Granny Smith’ (columns with small-cap letters) apples for each of xv    the five storage times (columns with different letters for the same cultivar are significantly different with p < 0.05). ....................................................................... 151  Figure 7.7 Confocal laser scanning microscopic images of a cross-section of apple flesh tissue 10 mm beneath the skin for ‘Golden Delicious’ apple at day 1 (a1) and day 30 (a2), and for ‘Granny Smith’ apple at day 1 (b1) and day 30 (b2) ( ‘ → ’ denotes the connection between cells, and ‘×’ denotes the broken cells). ............................ 152  Figure 7.8 Scanning electron microscopy images of one cross-section of apple flesh tissue 10 mm beneath the skin for ‘Golden Delicious’ at day 1 (a1) and day 30 (a2) acquired at low resolution for an overview of the cell structures, and at day 1 (b1) and day 30 (b2) acquired at high resolution for a more detailed view of the cell-to-cell connection [the red lines in (a1) denotes five cell connections counted manually as examples]. ...... 156  Figure 7.9 Scanning electron microscopy images of one cross-section of apple flesh tissue 10 mm beneath the skin for ‘Granny Smith’ at day 1 (a1) and day 30 (a2) acquired at low resolution for an overview of the cell structures, and at day 1 (b1) and day 30 (b2) acquired at high resolution for a more detailed view of the cell-to-cell connection [the red lines in (a1) denote five cell connections counted manually as examples]. ....... 157  Figure 8.1 Comparison of spatially-resolved diffuse reflectance obtained from the two-layer diffusion model (symbols) and Monte Carlo simulations (solid lines): (a) µa1/µa2 = 0.50 and µs1'/µs2' = 0.86 (0.05, 12, 0.10, 14, 1), (b) µa1/µa2 = 6.50 and µs1'/µs2' = 1.80(2.60, 36, 0.40, 20, 0.5), (c) µa1/µa2 = 0.80 and µs1'/µs2' = 1.58 (0.32, 19, 0.40, 12, 1), and (d) µa1/µa2 = 2 and µs1'/µs2' = 0.75 (1.00, 15, 0.50, 20, 1). Values in the -1 parentheses are µa1, µs1', µa2, µs2', d with the units of cm for optical properties and mm for the thickness of the first layer. ..................................................................... 171  Figure 8.2 Scaled sensitivity coefficients as functions of source-detector distances for four combinations of optical properties (a) µa1/µa2 = 0.50 and µs1'/µs2' = 0.86 (0.05, 12, 0.10, 14, 1), (b) µa1/µa2 = 6.50 and µs1'/µs2' = 1.80 (2.60, 36, 0.40, 20, 0.5), (c) µa1/µa2 = 0.80 and µs1'/µs2' = 1.58 (0.32, 19, 0.40, 12, 1), and (d) µa1/µa2 = 2 and µs1'/µs2' = 0.75 (1.00, 15, 0.50, 20, 1), and the values in the bracket are µa1, µs1', µa2, -1 µs2', d with the unit cm for optical properties and mm for the thickness of the first layer (‘-’, ‘·’, ‘+’, ‘ₒ’, and ‘∆’denote reflectance, and scaled sensitivity coefficients of µa1, µs1', µa2, and µs2', respectively)........................................................................ 172  Figure 8.3 Scaled sensitivity coefficients as functions of source-detector distances for µa1 = 0.05 -1 -1 -1 -1 cm , µs1' = 12 cm , µa2 = 0.10 cm , µs2' = 14 cm , and the thickness of the first layer: (a) d = 0.85 mm, (b) d = 2 mm, (c) d = 4 mm, and (d) d = 6 mm (‘-’, ‘·’, ‘+’, ‘ₒ’, and ‘∆’denote reflectance, and scaled sensitivity coefficients of µa1, µs1', µa2, and µs2', respectively). ..................................................................................................... 173  xvi    Figure 8.4 Estimated absorption and reduced scattering coefficients of the first layer from fitting the diffusion model to the Monte Carlo simulation data: (a) µa1 varies 0.20 - 1.20 cm 1 -1 -1 -1 with µs1' = 12 cm , µa2 = 0.50 cm , µs2' = 11 cm , and d = 2 mm, and (b) µs1' -1 -1 -1 varies 15-45 cm-1 with µa1= 0.10 cm , µa2 = 0.05 cm , µs2' = 11 cm , and d = 2 mm. ........................................................................................................................... 175  Figure 8.5 Relative errors of estimating the optical properties of two layers from the Monte Carlo diffuse reflectance for four combination sets of optical properties: µa1/µa2 = 0.50 and µs1'/µs2' = 0.86 (0.05, 12, 0.10, 14, 1), µa1/µa2 = 6.50 and µs1'/µs2' = 1.80 (2.60, 36, 0.40, 20, 0.5), µa1/µa2 = 0.80 and µs1'/µs2' = 1.58 (0.32, 19, 0.40, 12, 1), and µa1/µa2 = 2 and µs1'/µs2' = 0.75 (1.00, 15, 0.50, 20, 1). Values in the parentheses are µa1, µs1', -1 µa2, µs2', d with the units of cm for optical properties and mm for the thickness of the first layer. ............................................................................................................ 177  Figure 8.6 Absorption and reduced scattering spectra of the homogenous model disk samples measured by the integrating sphere........................................................................... 178  Figure 8.7 Comparison of diffuse reflectance from measurements (symbols) and the two-layer diffusion model (solid lines) for sample 1 at 535 nm (a) and 700 nm (b), and sample 2 at 535 nm (c) and at 700 nm (d) (See Table 8.1 for the two samples)...................... 180  Figure 8.8 Relative errors of the estimated optical properties for two model samples (denoted as S1 and S2) at wavelengths of 535 nm and 700 nm (see Table 8.1 for the optical property data for the two model samples)................................................................. 181                                      xvii    KEY TO SYMBOLS OR ABBREVIATIONS CA Controlled atmosphere CCD Charge-coupled device CI Confidence interval CLSM Confocal laser scanning microscopy CV Coefficient of variation EM Electromagnetic EMCCD Electron multiplying charge-coupled device FD Frequency-domain F/D Force/deformation GD Golden Delicious GPGPU General-purpose computing on graphics processing units GS Granny Smith HISR Hyperspectral imaging-based spatially-resolved IAD Inverse adding-doubling IF Impact firmness IMC Inverse Monte Carlo xviii    ITDM Integral-transformed diffusion model LTDM Logarithm-transformed diffusion model MAP Maximum a posteriori MC Monte Carlo MDI Multiple document interface ML Maximum likelihood MT Magness-Taylor NIRS Near-infrared spectroscopy OBF Order-blocking filter ODM Original diffusion model OLS Ordinary least squares OPA Optical Property Analyzer PCG Preconditioned conjugate gradients PLS Partial least squares PRESS Predicted residual error sum of squares RD Delicious RT Radiation transport xix    RWDM Diffusion model with relative weighting SDK Software developer kit SEM Scanning electron microscopy SEP Standard error of prediction SNR Signal-to-noise ratio SR Spatially-resolved SS Sum of squares SSC Soluble solids content TCSPC Time-correlated single-photon counting TR Time-resolved A Internal reflection coefficient, or beam area a' Transport albedo c Light velocity C * Chroma D Diffusion coefficient, or fractal dimension E Young’s modulus xx    f Frequency FI acoustic firmness index g Anisotropy factor h° Hue angle I Light intensity * * L , a , and b * CIELAB space LVs Latent variables J0 Zeroth-order Bessel function M Modulation m Mass mfp' Transport mean free path n Refractive index P Beam perimeter Q Source term in radiation transport equation R Reflectance r Spatial distance, or correlation coefficient xxi    Rd Roundness S Isotropic source T Transmittance θ Scattering angle, or phase angle Ω Solid angle Φ Fluence rate τ' Optical thickness σ Standard error, or failure stress λ Wavelength ε Relative error, or failure strain a Absorption coefficient s' Reduced scattering coefficient s Scattering coefficient µt' Total attenuation coefficient µeff Effective attenuation coefficient xxii    CHAPTER 1 INTRODUCTION 1.1 Background Quality of fruits and vegetables is characterized by appearance (color, size and/or shape), flavor (e.g., sugar content or soluble solids content), texture (e.g., firmness and crispness) and aroma. The relative importance of each quality index depends on the commodity and its final use. Many factors such as growing condition, maturity at harvest, and postharvest handling, storage, and transportation can directly or indirectly affect the quality of fresh produce. With the increasing demand for all year-round supply of high quality fruits and vegetables in both domestic and export markets, it is critical that producers supply fresh horticultural products with desired quality attributes in order to assure consumer satisfaction and confidence. As such, the horticultural industry is highly interested in adopting new, cost effective objective methods and techniques for analyzing and monitoring fruit and vegetable quality. Conventional methods for quality measurement of horticultural products include Magness-Taylor (MT) firmness test and Brix refractometry for sugar content (or soluble solids content). These techniques are destructive and only can test a few product pieces, and they are not suitable for sorting and grading each piece of fruits and vegetables. Over years, many nondestructive methods have been developed for quality evaluation of food and agricultural products. Among the many types of nondestructive sensing techniques (i.e., optical, mechanical, and electrical), optical techniques are particularly useful because they are rapid, nondestructive, cost effective, and generally safe to use. Abundant literature is available on optical methods and technologies for nondestructive quality measurement of horticultural products. However, few studies have been reported about 1   light propagation in plant materials, and of the determination of optical properties of horticultural products, which are critical for developing effective optical sensing techniques. Optical techniques, like imaging and near-infrared spectroscopy (NIRS), are widely used for nondestructive quality assessment of agricultural and food products (Abbott et al., 1997; Blasco et al., 2003; McGlone et al., 2003). Imaging or machine vision is mainly used to inspect the external features of product items. NIRS, on the other hand, has been applied for internal quality evaluation of agricultural and food products as well as nonfood materials. NIRS measurement provides approximate quantification of absorption properties of the sample (Dahm and Dahm, 2001), but it does not separate absorption properties from scattering properties, which could limit its capability in predicting such quality attributes as fruit firmness. Biological materials like fruits are optically inhomogeneous, and thus light will be attenuated in the fruit through absorption and scattering, which are affected by the structural and compositional characteristics of the tissue. Therefore, it would be advantageous and/or desirable to measure both absorption and scattering properties. Light propagation in turbid media is characterized by the absorption coefficient (a) and the reduced scattering coefficient (s') (Tuchin, 2000). Light absorption is mainly related to the chemical compositions of the medium, while scattering is influenced by the structural properties (density, particle size, and cellular structures). Changes in the chemical composition and physical structures are accompanied with the changes in the quality of fruit and vegetables (Lu et al., 2009; Qin and Lu, 2008). Hence, quantification of the absorption and scattering properties can greatly expand our understanding of light interaction with biological materials and enable us to assess the physiological state, properties or characteristics, and thus quality of agricultural and food products. 2   The propagation of light in a medium can be described by radiation transport theory, which is also known as the Boltzmann equation. When scattering is dominant over absorption (i.e, s'a), the radiation transport equation may be simplified by diffusion approximation (Durduran et al., 1997; Ishimaru, 1978). The diffusion equation has been widely used for modeling light propagation in single-layer or homogeneous and multi-layer turbid media. With appropriate boundary conditions and known optical parameters, direct solutions to the diffusion equation, referred to as forward problems, can provide quantitative description of light propagation in the turbid medium. Conversely, the diffusion equation can be used to estimate the absorption and reduced scattering coefficients of media for specific boundary conditions, which is called inverse problems. In solving inverse problems, we often resort to numerical methods (e.g., Monte Carlo, asymptotic approximation, and finite element) because they are able to deal with complex boundary conditions (Gonzalez-Rodriguez and Kim, 2008; Schweiger and Arridge, 1997; Seo et al., 2007). However, these methods are computationally intensive. For this reason, analytical solutions to the diffusion equation, coupled with a proper inverse algorithm, are preferred in determining the optical properties of turbid materials, which is also referred to as parameter estimation problems. Inverse light transport problems are much more complicated than forward problems, and sometimes they are ill-posed, especially for the two-layer diffusion model, in which four or five parameters (a and s' for each layer plus the thickness of the first layer) need to be estimated. Therefore, in order to have accurate estimation of optical properties, it is important to understand the intrinsic properties of the diffusion model, evaluate the feasibility of estimating all parameters and the uniqueness and stability of the solution, and consequently develop a reliable inverse algorithm. 3   Based on the radiation transport theory or diffusion approximation, several advanced optical techniques, which mainly include time-resolved (TR), frequency-domain (FD), and spatiallyresolved (SR), have been developed for measuring the optical properties of biological materials. These techniques have attracted considerable attention in the biomedical field for such applications as disease diagnosis, photodynamic therapy, and noninvasive monitoring (Bykov et al., 2006; Welzel et al., 2004; Wilson and Patterson, 2008). However, only limited research has been reported on measuring the optical properties of agricultural and food products using these techniques (Anderson et al., 2007; Cubeddu et al., 2001a; Lu et al., 2009; Qin and Lu, 2008). Time-resolved and frequency-domain techniques require sophisticated instrumentation, and they currently are still expensive and inconvenient or time consuming in measurement with limited wavelength range selections. In comparison, the spatially-resolved technique is less expensive, easier to use, and faster in measurement. It is, thus, more viable for agricultural and food applications. Spatially-resolved technique is commonly implemented in one of the two sensing configurations, i.e., multiple fiber arrays connected to spectrometers and non-contact reflectance imagery (Doornbos et al., 1999; Pilz et al., 2008). The former needs good contact between the detecting probes and the measured medium, which may not be desirable for food and agricultural products because of the safety and sanitation requirements. Most research using the non-contact sensing configuration can only provide optical property information at single or several wavelengths. In spite of great progress in instrumentation and theories over the past decade, we still lack a reliable and convenient method for accurate measurement of the optical properties of biological materials in general and food and agricultural products in particular. No commercial instruments based on the spatially-resolved spectroscopic principle are currently available. Over the past 4   decade, hyperspectral imaging has emerged as a powerful technique for quality and safety inspection of food and agricultural products (Lu and Chen, 1998; Lu and Peng, 2006; Nicolai et al., 2006; Park et al., 2002). Hyperspectral imaging combines conventional imaging and spectroscopy techniques to acquire spatial and spectral information simultaneously, and it is thus ideally suited for measuring spatially-resolved diffuse reflectance profiles for a broad spectral region. In light of those advantages, our laboratory recently developed a hyperspectral imagingbased spatially-resolved technique for measuring the spectra of optical properties of biological materials. Studies by Qin et al. (2007) showed that the technique is promising for optical characterization of food and agricultural products. However, limitations in application of this new technique were observed in our previous studies, and several critical issues need to be addressed before the technique is suitable for practical applications. These issues include the development of a robust light propagation model and inverse algorithm, optimization of the optical design of light beam and source-detector distance in the hyperspectral imaging-based spatially-resolved system, and the establishment of standard procedures for evaluating the instrument capability. Moreover, cost and convenience in use of this advanced optical technique in food and agriculture are also of great concern. Therefore, it is highly desirable to develop an instrument integrated with robust and user-friendly algorithms, which is capable of automatically acquiring and processing hyperspectral scattering image data to extract the optical absorption and scattering properties. In short, optical characterization of the fundamental absorption and scattering properties will provide opportunities in the development of new sensing technologies for nondestructive quality evaluation of horticultural and food products. 5   1.2 Objectives This research was aimed at developing a hyperspectral imaging-based spatially-resolved technique for accurate and reliable measurement of the optical properties of horticultural and food products and for maturity/quality assessment of peach and apple. The specific objectives were to: (1) Develop an improved inverse algorithm, based on a single-layer diffusion model, for extracting the absorption and reduced scattering coefficients from the hyperspectral reflectance profiles, and optimize the algorithm using Monte Carlo simulation and statistical analysis; (2) Evaluate and optimize the design of light beam and source-detector distance so as to achieve accurate acquisition of spatially-resolved spectral scattering profiles for determination of the optical properties; (3) Develop, test and evaluate an optical property measuring prototype, with integration of image acquisition and processing programs and the optimized inverse algorithm and optical designs, for automatic measurement of the optical properties of horticultural and food products; (4) Measure absorption and reduced scattering spectra of peaches and apples for assessing fruit maturity and/or quality; (5) Correlate the spectral absorption and scattering properties of apple tissues to the microstructural and mechanical characteristics, as affected by postharvest storage; (6) Investigate the feasibility of applying a two-layer diffusion model for determining the optical properties of fruit skin and flesh through solid model samples of known properties. 6   CHAPTER 2 LITERATURE REVIEW 2.1 Principles of Light Interaction with Biological Materials Light interaction with turbid biological materials is characterized by the fundamental optical properties of absorption, scattering and refraction, which are defined by the absorption coefficient, reduced scattering coefficient, and refractive index. These optical properties are present in the radiation transport equation governing light propagation through biological tissues. The definitions of these optical parameters are presented in the following subsections. 2.1.1 Absorption Coefficient Absorption coefficient (µa) quantifies the conversion of light energy into other forms of energy such as heat, electricity, or chemical energy. Absorption or the decrease in the amount of electromagnetic radiation is proportional to the incident light intensity and the distance over which the absorption takes place in an absorbing-only medium [Figure 2.1(a)]. According to the Beer-Lambert law, it can be expressed in the following equation: I  I 0 exp(   a d ) (2.1) where I is the intensity of the transmitted light, I0 is the incident light intensity, d is the thickness -1 -1 of the material in mm or cm, corresponding to the units of µa in mm or cm . Light absorption is mainly related to the chemical composition of the material. A relationship between the absorption and the chemical composition may be established, which could be used to evaluate quality, ripeness, and damage of agricultural and food products. Chemical bonds of the biological materials absorb light energy at specific wavelengths. The major absorbers for many horticultural products in the visible wavelength range of 400-770 nm 7   are pigments including chlorophylls, carotenoids, anthocyanin, lycopene and other color compounds, while water, fats, carbohydrates and proteins usually have absorption bands in the near-infrared wavelength region of 770-2,500 nm. Qin and Lu (2008) showed that as tomato fruit ripened from the green (unripe) to red (ripe) stage, their chlorophylls content decreased significantly. As a result, the absorption coefficient at 675 nm also decreased dramatically. An absorption peak due to water occurs to the µa spectrum of fruit at the wavelength of 975 nm (Nicolai et al., 2008; Valero et al., 2004). (a) Incident intensity I0 Length d Absorption µa Length d (b) Transmitted intensity I Incident intensity I0 (c) Scattering µs Transmitted intensity I (d) Incident light θ1  n1  Incident photon Specular reflection n2  θ2  Refraction Figure 2.1 Light interaction with matter: (a) absorption, (b) scattering, (c) scattering event with anisotropy factor g, and (d) refraction (For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this dissertation). 8   2.1.2 Scattering Coefficient Scattering (µs) is a physical process that takes place when light interacts with scattering media, and the travelling path of the photons is no longer direct as shown in Figure 2.1(b). µs quantifies the probability of photon scattering unit path length, and is the inverse of the average distance that light travels among scattering events. Similar to absorption, scattering can be described by the non-scattered transmitted intensity I, incident intensity I0, and the thickness of the medium with the form identical to equation (2.1). I  I 0 exp(   s d ) (2.2) where µs represents a probability per unit length of a photon being scattered, which has the same -1 units as µa in mm -1 or cm . Light scattering in the tissue depends on many variables including the size of scattering particles, the wavelength of the light, and the variation of the refractive indices of the various tissue components. In agricultural products, scattering is closely related to the cellular structures and characteristics, and therefore it could provide useful information about their condition or quality. 2.1.3 Anisotropy Factor Anisotropy factor is defined as a measure of the amount of photons retained in the forward direction after a single scattering event. When a photon is scattered by a particle, its trajectory is deflected by scattering angle θ as shown in Figure 2.1 (c). It has been proved that the HenyeyGreenstein function can be used to describe the probability density function for scattering, and it is given by (Tuchin, 2000) 9   p( )  1 g2 2(1  g 2  2 g cos  )3/2 (2.3) where g is called anisotropy factor, and it is dimensionless, which is the expected value of cosine of the scattering angle θ  g  cos    p( ) cos   2 sin  d 0 (2.4) The scattering phase function p(θ) is usually determined from goniophotometric measurements in relatively thin tissue samples (Tuchin et al., 1994). The value of g varies in the range from 0 to 1 with g = 0 for isotropic scattering and g = 1 for the total forward scattering. For most biological tissues, g  0.70  0.99 in the visible and near-infrared region (Vo-Dinh, 2003). 2.1.4 Reduced Scattering Coefficient The reduced scattering coefficient (µs') is a property incorporating the scattering coefficient (µs) and anisotropy factor ( g ), which is expressed as follows (Graaff et al., 1993):  s '   s (1  g ) (2.5) It is used to describe the diffusion of photons in a random walk with the step size of 1/ µs' (mm or cm), where each step involves isotropic scattering. In many biological materials, scattering is dominant during a light transport process, which is known as the diffusion regime. Because the photons encounter many scattering events in small steps before an absorption event takes place, the total scattering could be considered as isotropic. Hence the exact value of the anisotropy factor is no longer needed for the description of light propagation in the tissues. As a result, µa and µs' are the only optical parameters in the diffusion regime, which is common when visible 10   and near-infrared light propagates through biological tissues. Other relevant optical parameters can be derived from µa and µs', which include total attenuation coefficient (µt'), effective attenuation coefficient (µeff), transport mean free path ( mfp' ), and transport albedo ( a ' ) t '  a  s ' (2.6)  eff  3  a (  a   s ')  1/ 2 (2.7) mfp '  ( a   s ')1 (2.8) a '   s '/ (  a   s ') (2.9) Values of the reduced scattering coefficient for some biological materials have been reported in the literature. For instance, the µs' spectrum of fruit tissues in the visible and nearinfrared region is relatively flat compared with the µa spectrum (Cubeddu et al., 2001a; Qin and Lu, 2008). Empirical equations for determining µs' values of Intralipid, a scattering material commonly used for creating model samples, in the wavelength range of 400-1,100 nm have been derived by van Staveren et al. (1991). 2.1.5 Refractive Index Determination of the optical properties requires knowledge of the refractive index of the sample and the medium surrounding the sample. The refractive index is a measure of light refraction at the boundary between two materials [Figure 2.1(d)], and is governed by Snell’s law n1 sin 1  n2 sin  2 (2.10) where n1 and n2 are the refractive indices of the media, and θ1 and θ2 are the incident light angle and refraction angle for the transmitted light. 11   Since biological tissues are heterogeneous in composition, the refractive indices are different for different parts of a biological cell. It is common that only the average value of refractive index for the tissue is calculated, which is defined by the refractive indices of the scattering center and ground matter, and their volume fraction of the scatters (Tuchin, 2000). In studying light propagation in biological tissues, the overall refractive index of 1.35 is often used and it is considered constant in the visible and near-infrared region (Haskell et al., 1994). 2.2 Modeling of Light Transport in Tissues Light propagation in tissues can be described by the electromagnetic (EM) theory, i.e., Maxwell’s equation. However, it is not feasible to directly use the EM theory for describing light interaction with biological tissues, because of the complexity of cellular constituents and physical structures of biological tissues. When light is incident upon the surface of the tissue, it may be reflected, absorbed, scattered or transmitted, and the relative contribution of each phenomenon depends on the optical properties of the tissue. For biological materials, several simple approaches like the Kubelka-Munk model have been employed to describe these phenomena (Budiastra et al., 1998). However, they are only valid in several particular situations and cannot be generalized to most cases of practical interest. Mathematical models derived from the radiation transport equation include diffusion approximation, adding-doubling, and Monte Carlo simulation, which have been successfully used for modeling light transport in tissues and, hence, are reviewed in this section. 2.2.1 Radiation Transport Model The radiation transport (RT) equation is based on the assumption that only the flow of energy through the medium is considered, and migrating particles do not interact with each other. 12   In the case of tissue optics, the particles are photons (i.e., probability electromagnetic wave packets) interacting with cells structures. The RT equation can be derived by considering the radiant energy balance in a small volume of tissue. It is a balance equation related to the change of energy radiance in time and in energy flow, loss due to the absorption and scattering, and gain due to scattering sources and radiation sources. The time-dependent RT equation is given by (Ishimaru, 1978) ˆ 1 L(r, s, t ) ˆ ˆ ˆ ˆ ˆ ˆ ˆ    L(r, s, t ) s  ( a   s ) L(r, s, t )  s  L(r, s ', t ) p( s  s ')d  '  Q(r, s, t ) (2.11) 4 t c 2 ˆ where L(r , s, t ) is the radiance with the unit of W/(m sr), which represents the energy transfer ˆ per unit time per unit solid angle  in direction s through a unit area at position r and time t , ˆ ˆ ˆ sr is in the unit of solid angle, s is a unit vector pointing in the direction of interest, p(s  s ') is ˆ ˆ the phase function, which equals the probability distribution of the scattering angle from s to s ' , ˆ Q(r, s, t ) is the source term that represents power injected into a unit solid angle centered on the 8 ˆ direction s in a unit volume at r , and c = (3×10 )/n (m/s) is the velocity of light in the medium, where n is the refractive index of the medium. Although the RT equation has been successfully applied in various fields for the light propagation in absorbing/scattering media, the derivation of direct analytical solutions for equation (2.11) for many practical problems is difficult. Therefore, further approximations or numerical methods are needed for specific cases. A number of approximate and stochastic models of photon transport derived from the RT equation have been developed, of which diffusion model, adding-doubling, and Monte Carlo simulation are the most widely used. 13   2.2.2 Diffusion Model Diffusion approximation to the RT equation is the most widely used model to describe the light propagation in biological materials. It assumes that scattering is dominant (i.e., µs' » µa) and the source-detector distance is greater than the transport mean free path that defines the distance traveled before the direction of light propagation is randomized (Dogariu and Ellis, 2010). An ˆ additional assumption is that the source term Q(r, s, t ) is isotropic with an equal probability of scattering in all solid angles combined with a net flux. With these conditions, the RT equation is approximated to the diffusion equation given by (Haskell et al., 1994) 1 (r, t )  D2(r, t )  a(r, t )  S (r, t ) c t where (r, t )   4 (2.12) ˆ L(r, s, t )d  is the fluence rate, D  [3( a   s ')]1 is the diffusion ˆ coefficient, and S (r, t )   Q(r, s, t )d  represents an isotropic source. 4 It has been confirmed that the diffusion theory is accurate for describing photon migration in infinite, homogeneous turbid media. For studying biological materials, techniques based on the diffusion theory have been developed, which include time-resolved (TR), frequency-domain (FD), and spatial-resolved (SR) techniques. TR method uses a short pulse laser to illuminate the sample and then collects the time-delayed remitted light at a specific distance from the source, FD employs a sinusoidally modulated source to measure the phase and modulation, and SR technique adopts a continuous-wave point light beam incident on the sample surface with the remitted light collected at different source-detector distances. Different forms of the diffusion equation for TR, FD, and SR have been derived (Welch and van Gemert, 1995), and details for these measurement techniques are presented in Section 2.3. 14   Several analytical solutions derived from the diffusion equation for the homogenous media of particular geometries (e.g., semi-infinite, slab, cylindrical, and spherical geometries) are available (Farrell et al., 1992; Welch and van Gemert, 1995). Solutions for inhomogeneous and irregular samples are usually obtained using numerical methods (e.g., Monte Carlo, finite element, and finite difference). Moreover, solutions in analytical form for multilayer media with certain geometries have also been derived (Kienle et al., 1998b; Swartling et al., 2003b). These solutions provide an effective means for modeling light propagation in turbid media of specific geometries under given boundary conditions, and they could be used to extract the optical properties of biological materials by means of inverse algorithm. 2.2.3 Adding-doubling Model As a purely numerical method, the adding-doubling method provides another way to solve the RT equation with high accuracy and flexibility in modeling any medium with mismatched boundary conditions and anisotropic scattering. The inverse adding-doubling method has advantages for measuring the optical properties, because the only values needed are the reflectance and transmittance of the medium. The following assumptions are given for this modeling approach: no time-dependence of light distribution, homogenous optical properties for tissue layers, an infinite plane-parallel slab of samples, uniform illumination of collimated or diffuse light, internal reflection at boundaries governed by Fresnel’s law, and unpolarized light (Prahl et al., 1993). In the adding-doubling method, the reflectance function R(v, v ') and transmittance function T(v, v ') are defined as the radiance reflected or transmitted by the slab in direction v ' for light incident from the v direction (Figure 2.2). Based on the numerical integration of functions with 15   quadrature, the total reflectance and total transmittance for the diffuse irradiance can be expressed (Welch and van Gemert, 1995) Rd   1 1  0 0 R(v, v ')2vdv 2v ' dv ' 1 1 Td    T (v, v ')2vdv2v ' dv ' 0 0 (2.13) (2.14) Then, the transport albedo a ' and transport optical thickness  ' may be calculated by   1  4 R  T 2 d d 1      1 T d   a'    4  1  Rd  Td  2 1     9  1 T d     ln(Td ) ln(0.05)  ln( Rd )  '  15( Rd Td ) 2 if if Rd  0.1 1  Td (2.15) Rd  0.1 1  Td if Rd  0.1 (2.16) if Rd  0.1 where a' = µs'/(µa+µs'), τ' = d(µa+µs'), and d is the physical thickness of the slab. Thus, the two parameters, i.e., the absorption coefficient µa and the reduced scattering coefficient µs', can be determined for a combination of Rd and Td using equations (2.13) through (2.16). The procedures for implementing the adding-doubling method are as follows: 1) guess a set of optical properties, 2) calculate the reflectance and transmittance using equations (2.13) and (2.14), 3) compare the calculated values with the measured reflectance and transmittance, and 4) repeat until the calculated values and measured values are matched. 16   v R(v, v ') T (v, v ') Figure 2.2 Light incident at an angle v that is reflected and transmitted by a slab at an angle v ' in the adding-doubling method (Welch and van Gemert, 1995). The total reflectance and the total transmittance are usually measured using integrating spheres. The results obtained from the adding-doubling method are accurate for a wide range of optical properties but at the cost of increased computational time. Therefore, it is mostly used as a standard method for measuring the optical properties of biological materials (Lualdi et al., 2001; Saeys et al., 2008). 2.2.4 Monte Carlo Simulation Another powerful tool for modeling light transport in tissues is the Monte Carlo (MC) simulation method. MC simulation is a stochastic numerical method in which the expected value of a certain random variable is equivalent to the value of a physical quantity to be determined (Wang et al., 1995a). MC simulation offers a flexible and accurate approach for quantifying the optical features of light transport that are difficult to measure directly. It can be simply implemented, and has the ability to handle any complex geometries and multilayer media. However, the method is statistical in nature and requires computing the propagation of a large quantity of photons. Therefore, significant computational time is needed to achieve required precision and resolution. 17   The motivation for MC simulations is to predict radiant energy transport in tissues, while ignoring such features as phase and polarization. In MC simulation of light propagation, a photon is launched into the turbid medium, and it then propagates in the medium described by the optical properties (i.e., µa, µs, g, and n) of the investigated material. The photon transported in the medium is traced step by step, which is expressed as probability distributions to describe the step size of photon movement between the sites of photon-tissue interaction, and the angles of deflection in a photon’s trajectory when a scattering event occurs (Wang et al., 1995a). The simulation of light propagation is then achieved by tracking millions of photons, which can score multiple physical quantities simultaneously. Since MC simulation for the light distribution approaches an exact solution of the RT equation as the number of photons increases, it is also considered a standard method for describing light propagation in tissues. Forward and inverse MC models have been used to analyze the diffuse reflectance and extract the optical properties of biological materials. Public MC simulation programs known as ‘MCML’ and ‘CONV’ developed by Wang et al. (1995a) have been widely employed in the biomedical field for the study of light transport in tissues. Several methods to accelerate the MC simulation of photon migration in turbid media were proposed for slab geometries and inhomogeneous media. Alerstam et al. (2008) found that general-purpose computing on graphics processing units (GPGPU) can dramatically increase the speed of MC simulations. Zolek et al. (2008) corrected the calculation of polar deflection angle to improve the MC simulations of light transport in tissue. Other numerical approaches, such as finite difference, MC diffusion hybrid model, finite element, and asymptotic approximation (Alexandrakis et al., 2000; Cui and Ostrander, 1992; Deulin and L'Huillier, 2006; Seo et al., 2007), have also been used to model light propagation in 18   turbid media. Although those numerical methods are generally considered more accurate than the diffusion approximation model, they suffer the main drawback of computational complexity with long computational time, especially for solving the inverse problems to obtain the optical properties of biological materials by iterative fitting methods. 2.3 Techniques for Measuring Optical Properties of Biological Materials Various techniques have been developed for measuring the absorption and scattering properties of biological materials. Because of the obvious advantages of noninvasive or nondestructive measurements, determination of the optical properties of biological materials based on diffuse reflectance measurements has significant potential in food and agriculture. The section reviews optical techniques along with their solutions to the diffusion model, and these techniques may be divided into three main categories: time-resolved (TR) (Patterson et al., 1989), frequency-domain (FD) (Patterson et al., 1991), and spatially-resolved (SR) (Groenhuis et al., 1983b). The TR, FD, and SR techniques described in this section are for both single-layer and two-layer biological materials with a semi-infinite geometry. 2.3.1 Time-resolved Technique Time-resolved technique is based on measurement of the attenuation, broadening and delay of a short light pulse, caused by the absorption and scattering events during photon propagation in highly scattering media. Photon movement in a turbid medium is complicated due to the strong effect of light scattering at different wavelengths. Photons reflected or transmitted across the medium may have transported along many different paths within the medium. In a scattering dominant medium, a photon traveling path can be quite tortuous and may escape from the medium within several centimeters from the injection point. The photon moves in random 19   directions with respect to its trajectory at injection and continues a random walk within the tissue until it escapes at the tissue surface or is absorbed. Because the pathlength that a photon has taken within the tissue is proportional to the time, time-resolved technique seeks to use the pathlength information implied by the time of escape to specify the optical properties of tissue traveled by the photons [Figure 2.3(a)]. (b) Light source Intensity Detected light Light source Detected light Intensity (a) Time Time Figure 2.3 Measurement principle of (a) time-resolved and (b) frequency-domain techniques. Photon propagation measurements in the time domain depend on the ability to extract photon information encoded in the temporal distribution of the re-emitted light, following the injection of a short monochromatic pulse in a diffusive medium. Thus, temporal resolution and high sensitivity become two critical factors in designing a time-resolved system (Cubeddu et al., 2000). Temporal resolution is mainly affected by the width of the light pulse and by the response of the detection apparatus, and the power of the injected light pulse should be fixed at appropriate values to avoid possible damage or injury to the sample when considering the sensitivity. The safety regulations for the biological tissues give the maximum permissible value 2 of 2 mW/mm for laser pulses in the wavelength range of 600-1,000 nm. For most time-resolved systems, photons are recorded using the time-correlated single-photon counting (TCSPC) 20   technique within a short time period at a specified distance from the beam incident point (O'Connor and Philip, 1984). The diffusion model for the time-resolved signals has been developed, which describes light propagation in its time distribution and allows the quantification of optical coefficients µa and µs'. Consider a homogenous, semi-infinite turbid medium, the time-resolved reflectance function is given by (Patterson et al., 1989) R1layer (r, t )  (4 Dc) 3/2 1 5/2 (s ') t  r 2  ( ')2  s exp   exp(act ) 4Dct     (2.17) where D  [3( a   s ')]1 is the diffusion constant, and c is the velocity of light in the turbid medium. The term µac represents the absorption of the light, and the other terms are due to the early spreading of the light in the tissue, such as t -5/2 dominates the early dynamic of the reflectance signals. Kienle et al. (1998a), using the Fourier transform approach, derived an analytical solution from the diffusion equation using proper boundary conditions for a two-layer medium that has the first layer of thickness d and the semi-infinite second layer. It is assumed that the thickness of -1 the first layer is larger than one transport mean free path [d > (µa1+µs1') ]. The time-domain reflectance R2-layer(r, t) for the two-layer medium can be obtained by calculation of the real and imaginary parts of the reflectance in the frequency domain R2-layer(r, f), as detailed in the following section, at many frequencies and by fast Fourier transform of these data. As a novel nondestructive method, applications of time-resolved technique have been extensively reported in biomedical research for measuring the optical properties of human tissues 21   (Chernomordik et al., 2002; Kacprzak et al., 2007; Svensson et al., 2005). Its applications for fruits were also reported in several studies. Cubeddu et al. (2001b) developed a time-resolved spectroscopic measurement system based on the TCSPC technique to measure the optical properties of peaches, kiwifruit and apples. Nicolai et al. (2008) used the same system to predict soluble solids content and firmness of pear, although no satisfactory calibration models could be established between the optical properties from TR measurements and quality attributes. 2.3.2 Frequency-domain Technique Frequency-domain is another well-recognized technique for quantitatively measuring the optical properties of biological tissues, which could provide information equivalent to that obtained in the time domain. In the frequency domain, the propagation and measurement of light are accomplished through sinusoidally modulated sources. The photon flux at the detector will also be sinusoidal in time but the oscillation will be delayed in time relative to the source and reduced in amplitude relative to the average flux [Figure 2.3(b)]. In this case, the measured quantities are the phase angle between the detected and source signals and the amplitude of the oscillation relative to the DC level (Welch and van Gemert, 1995). The signal measured at the frequency-domain and time-domain can be transformed reciprocally by the Fourier transform (Arridge et al., 1992). However, in practice, there are specific advantages for frequency-domain technique compared to the time-resolved measurement, such as cheaper and simpler instrumentation, easier correction for instrument response, real-time measurement of phase and modulation, and capability for characterizing an entire spectrum from the measurement at a single frequency. 22   The optical parameters (µa & µs') could be obtained from the diffusion model based on the frequency-domain measurement of the modulation ( M ) and the phase angle. Patterson et al. (1991) derived the analytical solutions for the modulation and phase by the Fourier transform of equation (2.17) for a homogenous and semi-infinite medium, which are given as follows  r    1  i  1layer (r , f )   r  tan 1  (1   0 2  2 i )1/ 2 exp(   i ) M1 layer ( r , f )  1   (2.18) (2.19) where  0  {3(  a   s ')[ r 2  (  s ') 2 ][(  a c ) 2  (2 f ) 2 ]1/ 2 c 1}1/ 2 1  2 f  a c  r   0 sin  tan 1  2  1  2 f  a c  i   0 cos  tan 1  2             {3 a (  a   s ')[ r 2  (  s ') 2 ]}1/ 2 (2.20) (2.21) (2.22) (2.23) For the fixed values of distance r and modulation frequency f, optical parameters µa and µs' can be obtained graphically or using more accurate numerical methods to solve equations (2.18) and (2.19). The choice of the optimum frequency is very important, which depends on the optical properties of the measured materials. A solution to the frequency-domain for a two-layer medium was derived with the phase angle θ2-layer(r, f) and modulation M2-layer(r, f) calculated based on frequency-domain reflectance R2-layer(r, f), which is given by (Kienle et al., 1998a) 23   2layer (r, f )  tan1 Im  R2layer (r, f )   (2.24) Re  R2layer (r, f )   1/2 2 2    Im  R2 layer (r , f )   Re  R2 layer (r , f )      M 2 layer (r , f )    Re  R2 layer (r , f   2 layer / 2 )        (2.25) where R2 layer (r , f )  C11  r , z  0, f   C2 D1 i (r, z, f )   1  r , z , f  z z 0 exp( j 2 ft )  0 i ( z, f , s)sJ0 (sr )ds 2 (2.26) (2.27) C1 and C2 in equation (2.26) are two constants (see details in Chapter 3). Equation (2.27) is the fluence rate  i ( r , z , f ) in layer i , J0 is the zeroth-order Bessel function, and i ( z , f , s ) is the fluence rate of each layer obtained from solving diffusion equations in the frequency domain by Fourier transformation. In applications, the frequency-domain technique can be generally classified into ‘multidistance’ and ‘multi-frequency’ categories. The former acquires reflected light signals at several different tissue locations, and the latter employs a single source-detector position with multiple source modulation frequencies (Madsen et al., 1994; Tromberg et al., 1997). The light source such as lamps and lasers used in the frequency-domain is usually modulated by electro-optic and acousto-optic modulators, and the modulated beam will therefore contain higher harmonics, which is related to the selection of a detection method. Detectors to be used in frequency-domain technique require a bandwidth encompassing the desired modulation frequency. In summary, avalanche photodiodes offer the best frequency response followed by microchannel plate photomultipliers (MCP-PMTs) and conventional PMTs, which have been widely used in 24   biomedical applications (Coquoz et al., 2001; Spichtig et al., 2009; Welch and van Gemert, 1995). The early frequency-domain instrument was developed for measuring the decays of fluorescence intensity from the flexible molecules in proteins (Lakowicz et al., 1989). It has been also successfully applied to in vitro spectroscopy studies of turbid media for medical diagnosis (Banos et al., 2007; Tromberg et al., 2000; Tromberg et al., 1996). 2.3.3 Spatially-resolved Technique Spatially-resolved technique was first proposed by Reynolds et al. (1976) for understanding light propagation in the turbid media. Later, Langerholc (1982) and Marquet et al. (1995) suggested that spatially-resolved measurement could be used to determine the optical properties of biological tissues. In this method, a small continuous-wave light beam perpendicularly illuminates the sample’s surface, and the remitted light is measured at different distances from the light source (Figure 2.4). The optical coefficients µa and µs' of biological materials can then be extracted from the measured spatially-resolved reflectance profiles using an appropriate analytical solution of the diffusion equation coupled with an inverse algorithm. 25   Spatial profile R(r) r Incident light Absorption & scattering Turbid material (µa, µs') Figure 2.4 Measurement principle for spatially-resolved technique. For the case of steady-state spatially-resolved reflectance for a homogeneous semi-infinite turbid medium, Farrell et al. (1992) derived an analytical solution from the diffusion equation using the extrapolated boundary conditions, at which the fluence is forced to zero by introducing a negative ‘image source’. The diffuse reflectance from the medium is calculated as the current across the boundary, and it is originated from a single isotropic point source located at a depth of one transport mean free path in the medium. The final expression of the reflectance R1-layer at the surface of the semi-infinite turbid medium is R1 layer ( r )  a' 4  1 1 exp(   eff r1 ) 1 4A 1 exp(   eff r2 )  (  (  eff  ) )(  eff  )   r1 t ' 3 t ' r2 r12 r2 2  t '    (2.28) where r is the source-detector distance, a '   s '/ (  a   s ') is the transport albedo, eff  [3a (a  s ')]1/2 is the effective attenuation coefficient, and t '   a   s ' is the total  attenuation coefficient, r1  z02  r 2  1/2 1/2 and r2   z0  2 zb 2  r 2      26   are the distances from the observation point at the interface to the isotropic source and the image source, z0  (  a   s ') -1 , zb  2 AD , and A  0.2190 for n  1.35 is the internal reflection coefficient related to the relative index of the tissue-air interface n, which can be calculated from an empirical equation developed by Groenhuis et al. (1983a). Later, Kienle et al. (1997) proposed an improved analytical solution by expressing the reflectance as the integral of the radiance over the backward hemisphere based on the study of Haskell et al. (1994). In this case, the radiance can be expressed as the sum of isotropic fluence rate and the flux. The details of Kienle’s single-layer model are presented in Chapter 3. For a two-layer medium in the spatially-resolved condition, Kienle et al. (1998b) solved the diffusion equation under the same conditions as that for time-resolved and frequency-domain methods. The Kienle’s two-layer diffusion model is detailed in Chapter 8. In practice, spatially-resolved measurement employs a point light source or narrow collimated beam of the constant intensity and multiple detectors at different source-detector distances. Optical fiber arrays and non-contact reflectance imagery are two commonly used sensing configurations in spatially-resolved measurement systems (Doornbos et al., 1999; Fabbri et al., 2003; Jones and Yamada, 1998; Pilz et al., 2008). The former requires multiple spectrometers or a single imaging spectrometer to measure diffuse reflectance at different distances from the light incident point. Optical properties at multiple wavelengths or over a specific spectral region can be obtained using this method. Yet the measurements need good contact between the detecting probes and the sample, which may not be suitable for agricultural and food products. The second method usually uses a CCD (charge-coupled device) camera to acquire diffuse reflectance from the scattering medium generated by a point light beam. The measurement can be achieved without contacting the investigated medium, which is particularly 27   advantageous for food and agricultural products because of the safety and sanitation requirements. However, most research on non-contact reflectance imagery mode can only provide optical property information at single or several wavelengths. Compared with time-resolved and frequency-domain methods, the spatially-resolved technique has advantages of being relatively easier to use and less expensive in instrumentation for food and agricultural applications. Moreover, it is faster in measurement in the range of millisecond for a broad spectral region versus in seconds at single wavelengths with the timeresolved method (Cubeddu et al., 2001a; Qin and Lu, 2006). However, all these techniques still face considerable challenges in achieving accurate measurement of the optical properties due to experimental difficulties and complex mathematical modeling. Reported results on the measured optical properties of biological tissues are still not satisfactory, as shown in Table 2.1. Pifferi et al. (2007) used time-resolved method to measure the optical properties of Intralipid, and obtained the maximum error of 27% at 730 nm compared with the result obtained by the standard method. Spichtig et al. (2009) reported the results of measuring silicone-based model samples with errors of 10% for µa and 31% for µs' at 690 nm and 830 nm using a frequency-domain method. Pilz et al.’s (2008) results demonstrated 10% accuracy of µa and 5% of µs' for measuring tissue model samples using spatially-resolved technique at one single wavelength of 633 nm. Also, Pifferi et al. (2005) applied a general protocol for assessing several optical methods for determining optical properties, and they reported significant differences from different instruments in measuring the same model samples with the maximum discrepancy of 32% for µa and 41% for µs'. While improved accuracies of determining µa and µs' were reported in several recent studies (Martelli and Zaccanti, 2007; Xu and Patterson, 2006), large differences are present for the same 28   model samples from these studies. Moreover, most reported studies (Table 2.1) only tested a few model samples with a relatively small range of µa and µs' values for a narrow spectral region. Table 2.1 Summary of the performance of time-resolved, frequency-domain and spatiallyresolved measurements (TS: time-resolved, FD: frequency-domain, SR: spatially-resolved, POM: polyoxymethylene). Method TS Reference AnderssonEngels et al. (1993) Wavelength 450-660 nm Sample Porcine muscle Results Measured optical properties of tissue over the visible spectral range Swartling et al. (2003a) 660 and 785 nm Solid tissue model samples -1 µa=0.05-0.3 cm Absolute differences of less than 0.05 cm-1 for µa, and 10% difference µs'=9-20 cm for µs' Breast tissue Average µa of Variations below 40% for µa and 20% for µs' -1 Svensson et al. (2005) 660, 786, 916 and 974 nm -1 0.04 cm -1 and µs' of 8 cm Pifferi et al. (2007) 550-600 nm, and 610-700 nm Maximum error of 27% at 730 nm with the median error of 7% over the whole range of µa Intralipid µa =0.01-0.07 -1 cm -1 µs'=6-8 cm FD Tromberg et al. (1997) 674, 811, 849, and 957 nm Detected physiological changes in breast tissue Human breast tissue µa =0.02-0.16 -1 cm -1 µs'=7-11 cm Coquoz et al. (2001) 674, 811, and 849 nm Liquid model samples -1 µa =0.04-0.4 cm Errors of 10-15% for µa, and 5-10% for µs' -1 µs'=0.2-0.9 cm Xu and Patterson 750 nm Five IntralipidIndian ink model 29   Errors in µa of 1-27%, Table 2.1 (cont’d) (2006) Spichtig et al. (2009) samples µs'/µa=110 -1 µs'=20 cm 690 and 830 nm and µs' of 3-42% Silicone-based model samples µa = 0.10-0.15 Errors in µa of 10% and µs' of 31% -1 cm , µs'= 4-11 -1 cm SR Bays et al. (1996) 514 , and 630 nm POM samples Human esophagus -1 µeff = 0-6 cm , Errors of 21% and 32% for µeff and µs' -1 µs'=0.2-14 cm Doornbos et al. (1999) 600-900 nm Liquid and solid model samples in vivo tissue -1 µa =0-1.5 cm 40% error of µa for the model sample -1 µs'=5-20 cm Qin et al. (2006) 530-900 nm Intralipid-dye model samples -1 µa =0-0.8 cm Errors in µa of 16% and µs' of 11% -1 µs'=2.2-23.2 cm Martelli and Zaccanti (2007) 750 nm Intralipid-Indian ink model samples µa and µs' with standard deviation <2% , and no validation Pilz et al. (2008) 633 nm Tissue model samples -1 µa =0.03-1.0 cm Errors in µa of 10% and -1 µs'=5-20 cm 30   µs' of 5% 2.4 Approaches to Inverse Problems for Determining Optical Parameters Direct solutions to the diffusion equation, or numerical methods like Monte Carlo simulation for quantitative description of light propagation in the turbid medium, are referred to as forward problems. The corresponding inverse problems, also called inverse radiation transport problems, deal with estimating the optical properties of turbid media using the reflectance data obtained from the experimental measurement. While forward problems can be solved accurately because there often exists a unique solution based on complete description of a physical system, inverse problems are more complicated and may be difficult to solve because it is challenging to find at least one model system that is consistent with the observed data with unique model parameters. Therefore, in determining optical properties of biological materials, a suitable inverse algorithm should be developed, and information involving measurement data uncertainty, model efficiency, parameter characteristics, and curve fitting errors should be analyzed. 2.4.1 Sensitivity Analysis Sensitivity analysis determines how a result (or the dependent variable) changes when parameters that affect the result are changed. It is performed by calculating sensitivity coefficients, which are the first derivatives of the dependent variable with respect to the parameters (Beck and Arnold, 1977). Sensitivity coefficients describe the response of the change of the dependent variable to perturbations in the values of the parameters, and show the relationship between the dependent variable and the parameters. Sensitivity coefficients can help us determine if unique solutions for estimating all optical parameters exist, and/or identify those parameters that can be estimated when it is not possible to uniquely estimate all parameters from the measurements (Beck and Arnold, 1977; Taktak et al., 1993). They can be calculated by 31   X  R  (2.29) where β represents optical parameters µa or µs'. Sensitivity coefficients of parameters can be computed by multiple methods, such as finite different, complex step, differentiation of analytical solutions, and sampling methods. In many applications such as design/analysis studies, sensitivity analysis provides important information for characterization of the changes of a result with parameter variations, optimization and parameter estimation with gradient-based techniques, optimal experiment design, and uncertainty analysis (Blackwell, 2009). Thus, it is important and useful to conduct sensitivity analysis in the study of inverse light transport problems. 2.4.2 Inverse Algorithm Algorithms for solving inverse problems or parameter estimations have been developed with the advances of mathematics, statistics and computer technologies. Regression or curve fitting is the most used inverse method, which can be divided into linear and nonlinear. In the linear regression, the most commonly used methods are ordinary least squares (OLS), maximum likelihood (ML), and maximum a posteriori (MAP) (Beck and Arnold, 1977). Estimating optical properties from the diffusion model is considered a nonlinear regression problem, which can be solved using nonlinear least squares method by minimizing the general sum of squares (SS) function when the dependent variable is nonlinear in terms of the parameters. The extreme found for the sum of squares functions when the model is linear in the parameters can be proved to be the correct one and the unique minimum point exists. However, it is difficult to assure for the nonlinear cases since there may be more than one extremum. Therefore, it is recommended that contours or profiles of sum of squares be plotted in the region in which the solution is expected. 32   Another possibility is to start the iteration procedure with different sets of initial parameter values. All methods for nonlinear least squares optimization are iterative. From a starting point x0 , the method produces a series of vectors x1 , x2 , x3 ,... , which converges to x* , a local minimum for the given function. Most methods have measures to enforce the descending condition F ( xk 1 )  F ( xk ) (2.30) It prevents convergence to a maximizer and also makes it less possible that it converges towards a saddle point by using descent methods such as the steepest descent, conjugate gradient, line search, and trust region and damped methods (Tarantola, 2005). The Newton’s, Gauss-Newton and Levenberg-Marquardt methods are the most used to achieve optimization (Levenberg, 1944; Marquardt, 1963; Wilkinson, 1961). The Newton’s method provides a general framework for the least squares problems. However, the convergence is very slow or even fails if the hypothesis of quadratic convergence is not satisfied. The Gauss-Newton method is a very efficient method based on the first derivatives of the components of the vector function. It normally only has linear convergence, and it can give quadratic convergence in special cases. Considerable research has shown that it gives quite good performance in many applications (Aigner and Juttler, 2009; Baltagi, 1996; Dastidar, 2006). The Levenberg-Marquardt method is a modification of the Gauss-Newton method, which is a good choice for small to medium sized nonlinear least squares problems. 2.4.3 Data Transformation and Weighting Methods Before an inverse algorithm is implemented, data pre-check or pre-process is usually developed to analyze nonlinear data. Data transformation is a powerful method for analyzing 33   nonlinear regression models by linear approximation or changing data properties to obtain better estimations. In statistics, data transformation refers to the linear or nonlinear transformation of a mathematical function that each data point yi is replaced with the transformed value yTi  f ( yi ) . The linear transformed data (such as dividing, multiplying all Y values by a constant or subtracting a constant from all Y values) can be analyzed with linear regression, but the nature of the best-fit curve will not be changed. Nonlinear transformation (such as converting Y values to their logarithms, square roots, integrals, derivative, or reciprocals) will change the relative position of data points from the curve and produce a different curve to minimize the sum-ofsquares. It is usually applied when the raw data does not meet the statistical assumptions or the data range is at several orders of magnitude (Motulsky and Christopoulos, 2004). For the spatially-resolved diffuse reflectance profile, the reflectance data (dependent variable: Y) decreases dramatically along the source-detector distance (independent variable: X), and the statistical assumptions that the Gaussian distribution of errors of the Y-data and constant variance of errors along the X-axis are violated when nonlinear regression is applied. Therefore, nonlinear data transformation to the Y-data could be useful to make the scatter more Gaussian. The logarithm function converts the multiplicative pattern (proportional-variance) to an additive pattern (constant-variance), and the integral function transforms instantaneous data to accumulated data, which alters the relationship between the independent variable and the dependent variable, and could make the experimental data meet the assumptions. The implementations of the logarithm and integral transformation methods for the experimental data and Kienle’s model for single-layer turbid media are given below as an example. In the logarithm transformation, the natural logarithm of diffusion equation (3.4) (see Chapter 3), called logarithm-transformed diffusion model (LTDM), is given as 34   R1layer ,log (r )  log  R1layer (r )    (2.31) The integral of diffusion equation, defined as the integral-transformed diffusion model (ITDM), is calculated by (Gobin et al., 1999) r R1layer ,int (r )   R1layer (  ) d   C1  f (r , z0 )  f (r , z ')  C2  g (r , z0 )  g (r , z ') 0 (2.32) where z '  z 0  2 zb f (r, z)  g (r , z )   (2.33)  1 exp(eff z)  exp eff (r 2  z2 )1/2    4 Deff  (2.34)  1 exp(eff z )  z (r 2  z 2 )1/2 exp  eff (r 2  z 2 )1/2    4 (2.35) If the scatter of Y data is Gaussian and the variance of the scatter is the same at all values of X, then the correct parameters can be found by least squares estimates without any data transformation. However, in some cases, the variance of the scatter often increases as Y increases. With this type of data, the least squares method is inappropriate because it tends to give undue weight to points with large Y values on the sum-of-squares value and ignores points with small Y values. To overcome this problem, a proper weighting method can also be applied in the nonlinear regression. One common method is the relative weighting which is defined as minimizing the sum-of-squares of the relative distances of the data from the curve 2 (  (Yobs - Y pred ) / Yobs  , where Yobs is the experimental data, and Y pred is the predicted   35   reflectance from the diffusion model). Other weighting methods, such as Poisson weighting (weighting by 1 / Yobs ) and weighting by observed variability, are also used in nonlinear regression (Motulsky and Christopoulos, 2004). 2.5 Nondestructive Techniques for Measurement of Maturity and Quality Quality of fruits and vegetables is of primary concern to growers, packers or processors, and retailers throughout the production, postharvest handling and processing, and marketing steps. It is a subjective judgment related to the learned criteria for each consumer. With the advent of instrumentation and computers, a large number of sensing techniques and devices to measure quality or quality-related attributes have been developed. Traditionally, quality detection of fruits and vegetables is achieved using destructive techniques based on the chemical and mechanical properties, like refractometry for soluble solids content (SSC) and mechanical force/deformation techniques (e.g., Magness-Taylor test) for firmness measurement. Destructive techniques are confined to testing a small batch of samples. Nondestructive sensing techniques, on the other hand, have the potential for rapid and nondestructive quality evaluation, sorting and grading of individual product items. The increasing demands for all-year-round supply of high quality fresh produce by retailers are driving the development of new sensing technologies for quality evaluation of horticultural products. 2.5.1 Maturity and Quality of Fruits It is difficult to establish one single maturity index due to the complexity of physiological and/or biochemical processes taking place in fruit during growth and ripening as well as further complication by such factors as variety or cultivar, geographic location, season, and climate condition. Hence a combination of indices is commonly used to assess the maturity of fruit, 36   which include ground color changes, firmness, soluble solids content, titratable acidity, ethylene production, and starch content. In addition, the maturity stage at which fruits should be harvested depends, to a great extent, on the final use of the products. For instance, if fruit are to be sold for fresh eating immediately after harvest, they may be allowed to ripen fully on the tree. However, if fruit are to be kept in storage for a long time period, they should be picked prior to full maturity or ripeness in order to prolong the postharvest life and minimize physiological disorders or quality decay during storage. Fruit quality is determined by the maturity of the fruit at harvest. It can be characterized by appearance (color, size and/or shape), flavor, texture and aroma. A definition of quality used in postharvest quality evaluation was given by Kramer and Twigg (1970): “the composite of those characteristics that differentiate individual units of a product, and have significance in determining the degree of acceptability of that unit by the buyer.” The specific quality requirements for horticultural products also depend on the end use and cultivar. Minimum quality standards have been established for fruits in many countries, and they can be divided into external and internal (Wills et al., 2004). The external standards, characterized by appearance (color, size, shape, and/or surface defect), are one of the key factors for consumers in purchase of fresh fruit. Many fruits undergo color changes as part of the ripening process. In some cases, color is a strong indicator of eating quality and shelf life. Internal quality is dependent on the properties, physical structures, and chemical compositions (e.g., firmness and SSC), which are directly related to the texture and flavor of fruits. And it is critical in the repeat purchase of a particular product or cultivar. Ripening causes changes in the composition and cellular structures of fruit and thus the eating quality. The sweetness taste development is the result of increased gluconeogenesis, 37   hydrolysis of polysaccharides, and decreased acidity. Softness of fruit is due to the loss of firmness in texture caused by cell wall depolymerization and an increase in solubility of the middle lamella. Thus, cellular characteristics such as cell wall strength, cell turgor, number of cells and their shape and size, and intercellular spaces have a significant effect on the texture of the fruit, which is, in turn, directly related to the mechanical properties of fruit tissue. Mechanical properties of fruit tissue are influenced by ripening stage and water status, and also determine the susceptibility of fruit to mechanical damage during harvest, transport, and postharvest storage and handling. Thus, better understanding of the structural and mechanical properties of fruit tissue is needed in the development of new or improved techniques for fruit quality evaluation and monitoring, and better processing operations for the food industry. 2.5.2 Mechanical Techniques Most nondestructive mechanical techniques measure elastic properties (e.g., modulus of elasticity or Young’s modulus) at very small deformations, which are related to the texture (e.g., firmness) of fruits and vegetables. Modulus of elasticity measures the capacity of the material to take elastic deformation, and it is the stress/strain ratio, which can be determined from the slope of the compressive force/deformation curve for cylindrical tissue specimens or whole fruit at small levels of deformation prior to rupture (Figure 2.5) (Abbott, 1999). 38   Figure 2.5 Force/deformation (F/D) curves for cylindrical tissue specimens from firm and soft apples under constant strain-rate (0.42 mm/s) compression (Abbott, 1999). Both sonic and ultrasonic vibration techniques have been applied for evaluating fruit firmness. Ultrasonic waves can be reflected, transmitted, refracted or diffracted as they interact with the fruit, which are directly related to the mechanical properties and geometry (i.e., size and shape) of the fruit. Recent advances in ultrasonic technology for monitoring the quality of fruits and vegetables have been reviewed by Mizrach (2008). The technique has been used for measuring firmness and mealiness in some fruits, and different ultrasonic sources are needed for different fruits in order to achieve good penetration (Bechar et al., 2005; Galili et al., 1993; Kim et al., 2004; Mizrach et al., 2003). Sonic vibration method, on the other hand, measures vibrational responses of fruit to sonic signals (i.e., impulse or sinusoidal) in the audible frequency range for up to a few thousand Hz. Sonic firmness is calculated based on the mass of the fruit and the first or second resonance frequency (Langenakens et al., 1997), depending on 39   the measurement configuration. The technique has been used for fruit firmness or stiffness measurement (Chen and Debaerdemaeker, 1993; De Belie et al., 2000; Diezma-Iglesias et al., 2004; Wang et al., 2006). Low-mass impact method provides another means for determining fruit firmness, and the technique currently is commercially available for online sorting and grading of fruit (De Ketelaere et al., 2006; Diezma-Lglesias et al., 2006). Although numerous nondestructive mechanical methods have been reported, they are often slow in measurement and difficult to implement for online applications. Moreover, these mechanical methods are only limited to firmness evaluation. Therefore, new sensors that are more reliable and faster and able to measure multiple quality attributes are needed to improve industry profitability and meet the increasing demands from the consumer for better and more consistent fresh products. 2.5.3 Optical Techniques The electromagnetic spectrum from the shortest to longest wavelengths includes Gamma rays, X-rays, ultraviolet, visible, infrared and radio waves (Figure 2.6). Optical techniques based on the electromagnetic radiation for a specific spectral range have been developed for nondestructive quality measurement of fruits and vegetables because of the rich information obtained from light interaction with materials. For fruits and vegetables, the most useful spectral regions are visible (400-770 nm) and near-infrared (770-2500 nm), while other spectral regions such as ultraviolet and X-rays have also been found useful. Visible/near-infrared spectroscopic and imaging techniques are most popularly used in postharvest quality inspection. 40   Figure 2.6 Diagram of the whole electromagnetic spectrum (Mohsenin, 1984). Color is the basis for sorting many products into commercial grades, which is directly related to pigments (e.g., chlorophylls, carotenoids, and anthocyanins) in the products, and can be measured in the visible range. Near-infrared (NIR) spectroscopy is well suitable for the determination of chemical compositions. The NIR region involves the responses of the molecular bonds O-H, C-H, C-O and N-H. These bonds are subject to vibrational energy changes when irradiated by NIR frequencies, and two vibration patterns exist in these bonds including stretch vibration and bent vibration. The energy absorption of organic molecules in the NIR region occurs when molecules vibrate or are translated into an absorption spectrum (Cen and He, 2007; Williams and Norris, 2001). Chemical constituents, such as ethanol, water, sugars (glucose and fructose), organic acids (citric, lactic, tartaric and malic acids), phenolic compounds and oxidation in fruits and vegetables, either in combination or alone, have strong influences on the 41   quality evaluation of fruits and vegetables. In practice, NIR measurement can be classified into three different modes: reflectance, transmittance and interactance (Figure 2.7). Transmittance measurement is often applied for liquid samples, such as fruit juice, using glass or quartz chamber with different sizes, while diffuse reflectance or interactance measurement is popular for solid samples such as fruits (Abbott et al., 1997; Lu et al., 2000; McGlone et al., 2003). Over the past two decades, considerable research has been reported on NIR technique for nondestructive determination of firmness, soluble solids content, acidity, maturity, dry matter, moisture and other characteristics of fruits and vegetables, including apple, pear, cherry, melon, tomato, cucumber, potato, etc. (Carlomagno et al., 2004; Lammertyn et al., 1998; Pedro and Ferreira, 2005; Sirisomboon et al., 2007). (a) (b) Figure 2.7 Three different modes for NIR spectroscopic measurements: a) reflectance, b) transmittance and c) interactance.     42   Figure 2.7 (cont’d) (c) Imaging techniques are widely used for classification or sorting of agricultural products based on the recognition of shape, size, surface defect, and color (Brosnan and Sun, 2004; Chao et al., 1999; Chinchuluun et al., 2009; Singh and Delwiche, 1994). Conventional imaging (e.g., monochrome and color camera) is based on the analysis of spatial information, acquired by a digital imaging device, from the object, which is effective when quality attributes of a product are related to its extrinsic characteristics. Recently, multispectral and hyperspectral imaging techniques have been used for quality evaluation of horticultural and food products. Multispectral imaging usually generates a set of images at fewer than 10 discrete wavelengths, which can be obtained either by positioning a bandpass filter in front of a monochrome camera (Park et al., 2004), or by capturing a series of spectral image using an acousto-optic tunable filter or a liquid crystal tunable filter (Peng and Lu, 2006; Tran, 2000). Hyperspectral imaging combines conventional spectroscopy and imaging techniques to acquire both spectral and spatial information from an object simultaneously, enabling us to analyze product properties or characteristics more reliably and completely than existing imaging and spectroscopic techniques. The technique may be implemented by acquiring a sequence of narrow-band spectral images or 43   capturing scan line images with a full spectral range to create 3-D spatial-spectral data or hypercubes (Bernhardt, 1995; Park et al., 2002). Many studies on multispectral imaging (Blasco et al., 2009; Lu and Peng, 2007; Park and Chen, 1994) and hyperspectral imaging (ElMasry et al., 2007; Lawrence et al., 2003; Lu, 2003) have been reported on quality and safety inspection of food and agricultural products. Besides spectroscopic and imaging techniques for the visible and NIR region, other optical techniques have also been developed for quality inspection of agricultural and food products, such as X-ray (Brecht et al., 1991; Haff et al., 2006; Schatzki et al., 1997), fluorescence and delayed light emission (Abbott and Massie, 1985; Gunasekaran, 1990; Kolb et al., 2006; Noh and Lu, 2007; Seiden et al., 1996), magnetic resonance and magnetic resonance imaging (Barreiro et al., 2000; Kim et al., 1999; Koizumi et al., 2009). Techniques based on electrochemistry, such as electronic nose and electronic tongue (Buratti et al., 2006; Infante et al., 2008; Rudnitskaya et al., 2006; Zoltan et al., 2008), are used for measuring the concentration of volatiles of fruits or vegetables to evaluate the ripeness. Optical techniques are advantageous for nondestructive maturity and quality measurement of horticultural products since they are generally nondestructive or noninvasive, fast, and more easy to implement for offline and online applications. However, spectroscopic techniques such as NIR mostly rely on statistical methods to establish empirical relationships between experimental data and chemical compositions, while most imaging techniques are mainly restricted to quantifying the physical or surface characteristics of the samples. Due to inadequate understanding of, and technical challenges in measuring, light transfer in biological materials, many optical techniques still are of empirical nature in measurement, and thus cannot provide accurate, complete information about the optical characteristics or properties of food products. 44   Consequently, they have not been able to accurately measure fruit quality or establish reliable, robust relationships or models for predicting the quality of fruit. 45   CHAPTER 3 OPTIMIZATION OF INVERSE ALGORITHM FOR ESTIMATING OPTICAL PROPERTIES OF BIOLOGICAL MATERIALS 3.1 Introduction When light scattering is dominant, the radiation transport equation can be simplified to a diffusion approximation equation (Ishimaru, 1978). Direct solutions to the diffusion equation coupled with appropriate boundary conditions, referred to as forward problems, provide a quantitative description of light transport in the medium. The corresponding inverse problems, also called inverse radiation transport problems, deal with estimating the optical properties of turbid media. Different approaches and methods have been developed for determining or measuring the optical properties of turbid media. They may be divided into three groups: direct, empirical and inverse (Cheong et al., 1990). Direct methods, such as the Beer-Lambert law for measuring collimated light transmission, are based on some fundamental concepts and rules without using any models to obtain optical parameters from measurements. The methods are conceptually simple, but they can be used only for samples of simple geometry (e.g., thin slabs). Empirical methods, on the other hand, obtain optical properties from the measured transmittance and reflectance using mathematical models that are based on simplifications and restrictive boundary conditions [e.g., Kubelka-Munk theory (Edstrom, 2004)]. As a result, empirical methods lack accuracy and generality, and they are only suitable for samples of specific geometry. Unlike direct and empirical methods, inverse methods determine optical properties by solving the radiation transport equation coupled with appropriate boundary conditions and an inverse algorithm. Inverse methods are of great interest to researchers because they are applicable to a wide range of intact turbid materials with minimal or no sample preparation. Inverse methods 46   may be implemented using diffusion theory, inverse Monte Carlo (IMC) simulation, and inverse adding-doubling (IAD). Being purely numerical, IMC and IAD achieve relatively high accuracy by describing the light propagation more rigorously using the radiation transport equation (Palmer and Ramanujam, 2007; Prahl et al., 1993). But they suffer the main drawback of computational complexity with long computational time. In this work, a diffusion model was chosen because it is more efficient and accurate for describing the radiation transport in turbid media, needs less computational time, and is particularly useful for nondestructive measurement. For the case of steady-state spatially-resolved reflectance for a homogeneous semi-infinite turbid medium, Kienle’s model was used in this study (Kienle et al., 1996). Because of the complexity of the analytical solution and potential experimental difficulties/errors in measuring diffuse reflectance profiles from the medium, there exist considerable difficulties in accurate estimation of the optical parameters from spatially-resolved diffuse reflectance. In this work, the optical parameter estimation is formulated as a nonlinear least squares optimization problem, which is based on several important assumptions (i.e., constant variance errors, uncorrelated errors, and the Gaussian distribution of errors). The results will not be valid if these assumptions are violated and the data are not presented in appropriate format. Proper data transformation and weighting methods should be considered when some of the assumptions are violated. Moreover, to improve the accuracy of the parameter estimation, the inverse algorithm and the experimental design need to be optimized and the information involving measurement uncertainty, model efficiency, curve fitting errors, and parameter characteristics should be acquired and analyzed. Therefore, the overall objective of this part of the research was to optimize the inverse diffusion theory algorithm and instrumental measurement for accurate estimation of the 47   absorption and reduced scattering coefficients from spatially-resolved diffuse reflectance data. The specific objectives of this research were to: a) examine different data transformation and weighting methods for nonlinear least squares estimates; b) perform sensitivity analysis to determine an appropriate data transformation and/or weighting method to improve the inverse algorithm for estimating the optical properties of turbid media; and c) assess the robustness of the diffusion model and the inverse algorithm using statistical analysis. 3.2 Forward Problem 3.2.1 One-layer Diffusion Model for the Steady-state Case Consider a semi-infinite turbid medium which is impinged upon perpendicularly by an infinitely small light beam. Under steady-state condition, the diffusion equation (2.12) can be rewritten as D2(r)  a(r)  S (r) (3.1) where r  ( x, y, z) , D  1 / [3(  a   s ')] is the diffusion constant, and S (r) is the isotropic scattering source. The solution of the fluence rate at the surface (z = 0) using the extrapolated boundary in which the fluence rate is forced to be zero is given by (Farrell et al., 1992)  ( r , z  0)  1  exp(   eff r1 ) exp(   eff r2 )     4 D  r1 r2  (3.2) where eff  [3 a (  a   s ')]1/ 2 is the effective attenuation coefficient, r1  ( z02  r 2 )1/2 , r2  [( z0  2 zb )2  r 2 ]1/2 , r  ( x2  y 2 )1/2 is the source-detector distance, z0  ( a   s ')1 , and zb  2 D 1  R e ff 1  R e ff for the extrapolated boundary condition. Reff is the effective reflection 48   coefficient, which is equal to 0.4498 for the refractive index n = 1.35 (Haskell et al., 1994). The diffuse reflectance calculated as the flux across the boundary is expressed by R flux ( r )  1 4  1 exp(  eff r1 ) 1 exp( eff r2 )   ( z0  2 zb )( eff  )  z0 ( eff  )  r1 r2 r12 r22     (3.3) The final solution for the steady-state spatially-resolved diffuse reflectance for a single-layer material derived by Kienle and Patterson (1997) is given below R1layer (r )  C1(r , z  0)  C2 R flux (r ) where C1  (3.4) 1 3 2 2 [1  R fres ( )]cos d  and C2  4 2 [1  R fres ( )]cos  d  are constants 4 determined by the relative refractive index mismatch at the tissue-air interface, in which R fres ( ) is the Fresnel reflection coefficient for a photon with an incident angle  relative to the normal to the boundary, and  is the solid angle. Details on these parameters can be found in the literature (Haskell et al., 1994). For a refractive index n = 1.35, a typical value for many biological materials (Mourant et al., 1997), C1 and C2 are equal to 0.1277 and 0.3269, respectively. Hence, the diffuse reflectance R(r) at the surface of the turbid medium is a function of the source-detector distance (r) as well as two unknown optical parameters of the medium. 3.2.2 Monte Carlo Simulation For accurate estimation of the optical parameters of turbid media, the diffusion model and inverse algorithm were validated by MC simulations. The MC simulation assumes that an infinitely narrow photon beam is normally incident on the surface of the medium, which is considered infinitely wide compared to the spatial extent of photon distribution. For a 49   homogenous medium, the light transport is determined by its refractive index (n), absorption coefficient (µa), and reduced scattering coefficient (µs'). A publicly available MC simulation program for light propagation in biological materials was used (Wang et al., 1995b). The low-noise reflectance generated by MC simulation provided an ideal measurement without introducing experimental uncertainties. Several groups of absorption and reduced scattering coefficients that are similar to the optical properties of fruits were selected. 3.3 Inverse Problem 3.3.1 Nonlinear Least Squares Inverse Algorithm To extract the optical properties from the spatially-resolved diffuse reflectance for a turbid medium, a suitable inverse algorithm should be selected. Equation (3.4) shows that the two optical parameters may be estimated based on the measured reflectance when no a priori information related to these parameters is given. Nonlinear least squares method was used to find the minimum of the sum-of-squares of the difference between the true reflectance and predicted reflectance values with estimated parameters. A large-scale algorithm such as a subspace trust-region method based on the interiorreflective Newton approach was selected to achieve the algorithm optimization (Thomas and Li, 1994). It is defined by minimizing a quadratic function subject only to an ellipsoidal constraint, and is globally and locally quadratically convergent. The interior-reflective Newton approach can generate iterates in the strictly feasible region by using a new affine scaling transformation, and the speed of convergence is accelerated by following a reflective line search technique. A preconditioned conjugate gradients (PCG) method was used to produce an approximate solution 50   to a large linear system in each iteration (Thomas and Li, 1994). The upper bandwidth of preconditioner for PCG was 2, which was determined by the preliminary implementation of the algorithm to achieve a better quality step towards the solution and also by consideration of the number of PCG iterations to reduce the computational time. The inverse algorithm for parameter estimations was implemented using the Toolbox function ‘lsqcurvefit’ in Matlab 7.5 (The MathWorks, Inc., Natick, MA). 3.3.2 Data Transformation and Weighting Methods In this study, two data transformation methods including logarithm and integral transformation for the experimental data and the diffusion model were used before the curve fitting process. In addition, the diffusion model with relative weighting (RWDM) was also used for the nonlinear least squares curve fitting to extract the optical properties, and the results were compared with those obtained from the data transformation methods with absolute weighting  [  Robs  R pred 3.3.3  2 ] . The details of LTDM, ITDM, and RWDM are presented in Section 2.4.3. Sensitivity Analysis To better understand how the changes in the reflectance profiles of the diffusion model are attributed to the changes in the optical parameters, sensitivity analysis was performed before the inverse algorithm was implemented. The sensitivity coefficients for the original diffusion model (ODM), LTDM, ITDM, and RWDM were calculated by equations (3.5)-(3.8), which were derived from equation (2.29). X ori _    R  51   (3.5)  log( R)  (3.6) X int_    Rint  (3.7) X RW _    R R (3.8) X log_    The sensitivity coefficients defined in equations (3.5)-(3.8) have the same unit as the dependent variable and can be readily compared for different parameters. In these equations, the sensitivity coefficients are expressed in absolute values in order to facilitate the comparison of different parameters. In this study, sensitivity coefficients were first calculated for the absorption and reduced scattering coefficients using the finite difference method. They were examined against the Monte Carlo simulation data to determine if it is possible to estimate individual optical parameters before an inverse algorithm was developed. 3.3.4 Statistical Analysis and Model Assessment Statistical analysis was conducted to interpret the results of parameter estimation. It systematically considers the characteristics of the model and algorithm. Several statistical parameters used for assessing the mathematical model and analyzing the curve-fitting results are listed in this section. The equations given below are based on the ODM. Corresponding equations for the transformed models can be derived easily from these given equations. To compare the estimated optical parameters by fitting the MC simulation data to the diffusion model with the true values, the relative errors are calculated by ˆ   t  100% t 52   (3.9) ˆ ˆ ˆ where  is the estimated value for optical parameter  a or s ' , and βt is the true value of µa or µs'. The residual for reflectance is given by resi  Robs,i  R pred ,i (3.10) The sum-of-squares for the n data points are minimized by the nonlinear least squares regression to obtain the estimated values of µa and µs'. n SS   resi 2 (3.11) i 1 The variance-covariance matrix gives the information about the standard error   of the parameters. The vector of   is the square root of the corresponding diagonal of the symmetric parameter    diag ([ J T J ]1 s2 ) (3.12) where    [ a ,  s ' ] , in which   a and   s ' are the standard errors of µa and µs', J   R pred is the estimated Jacobian, and s 2  1 SS is the mean squared residual with p n p parameters. To compute the confidence intervals (CIs) for the two parameters, the ‘NLPARCI’ (beta, residuals, Jacobian) function in Matlab was used. 3.4 Simulation Experiments To validate the diffusion model and the inverse algorithm by using MC simulations, the medium was considered to be turbid and semi-infinite. Absorption and reduced scattering 53   coefficients are two input parameters for the MC simulation model, and their values determine the simulation results. Thirty-six different combinations of µa and µs' and their transport mean -1 free path [1mfp' = (µa+µs') ] were selected (Table 3.1), which span a large range of values: -1 -1 0.04  µa  8 cm , 4  µs'  40 cm and 5  µs'/ µa  100 . These values were chosen based on the published data for the optical properties of fruit and other food products (Budiastra et al., 1998; Cubeddu et al., 2001b; Qin and Lu, 2008). For groups 5, 6 and 12, MC simulations gave zero reflectance at large source-detector distances due to large values for µa and µs'. Therefore, these three groups are not included in the following discussion. The refractive index of the medium was assumed to be 1.35, similar to that of fruit (Mourant et al., 1997), and the refractive 6 index for the medium above was 1.0, equal to that for the air. A total of 3×10 photons and 0.1 mm spatial resolution of both radial distance and depth were used to produce the reflectance for the spatial distance of 0.1- 10 mm. The MC generated diffuse reflectance profiles were then fitted by the inverse algorithm for the diffusion model to deduce the optical properties of the media. MC simulations introduced small stochastic variability in the calculated reflectance for any given set of optical properties; in this study the average variability was 0.8-3.1%, generated from five simulations for each set of optical properties. Table 3.1 Thirty-six combinations of the absorption (µa) and reduced scattering coefficients (µs') and their corresponding transport mean free path ( mfp' ) used in Monte Carlo simulations (unit: -1 cm for µa & µs', and mm for mfp'). Group no. 1 2 µa µs'/ µa = 5 0.80 1.40 µs' 4 7 mfp' Group no. 2.08 1.19 19 20 54   µa µs'/ µa = 50 0.08 0.14 µs' mfp' 4 7 2.45 1.40 Table 3.1 (cont’d) 3 4 *5 *6 2.00 4.00 6.00 8.00 7 8 9 10 11 *12 µs'/ µa = 10 0.40 0.70 1.00 2.00 3.00 4.00 10 20 30 40 4 7 10 20 30 40 0.83 0.42 0.28 0.21 21 22 23 24 2.27 1.30 0.91 0.45 0.30 0.23 25 26 27 28 29 30 µs'/ µa = 20 13 0.20 4 2.38 31 14 0.35 7 1.36 32 15 0.50 10 0.95 33 16 1.00 20 0.48 34 17 1.50 30 0.32 35 18 2.00 40 0.24 36 * These groups were excluded in the data analysis. 0.20 0.40 0.60 0.80 10 20 30 40 0.98 0.49 0.33 0.25 µs'/ µa = 70 0.06 0.10 0.14 0.29 0.43 0.57 4 7 10 20 30 40 2.46 1.41 0.99 0.49 0.33 0.25 µs'/ µa = 100 0.04 0.07 0.10 0.20 0.30 0.40 4 7 10 20 30 40 2.48 1.41 0.99 0.50 0.33 0.25 3.5 Results and Discussion 3.5.1 Comparison between the Diffusion Model and Monte Carlo Simulations Figure 3.1 shows the spatially-resolved diffuse reflectance calculated from the diffusion model and MC simulations for the two sets of optical properties with µs'/ µa = 5 and µs'/ µa = 100. Overall, the reflectance profiles calculated from the diffusion model well match those from MC simulations, except for small distances close to the light source, where relatively large differences are observed. However, the differences of reflectance are affected by the ratio of µs'/ µa. The average differences for groups 1-4 with µs'/ µa = 5, groups 7-11 with µs'/ µa = 10, and groups 13-36 with 20  µs'/ µa  100 are 18.1%, 10.8%, and 4.3%, respectively, for the source55   detector distances of greater than 1.5 mm. With relatively large values for µs'/ µa (i.e., ≥10) the differences between the diffusion model and MC simulations become smaller, which can be explained by the scattering dominated condition for the diffusion model. Hence the optical properties with 10  µs'/ µa  100 (a total of 29 groups in Table 3.1) were used in the parameter estimation. Larger deviations for reflectance are observed in Figure 3.1 when the detector is close to the light source (less than one transport mean free path), because the diffusion approximation is not valid in this situation. These results indicate that overall the diffusion model accurately quantifies light propagation in the turbid media. (a) 0 Reflectance (mm-2) 10 -2 10 -4 10 0 2 4 6 Distance (mm) (b) 8 10 0 2 4 6 Distance (mm) 8 10 0 Reflectance (mm-2) 10 -5 10 Figure 3.1 Comparison of the spatially-resolved diffuse reflectance obtained from the diffusion model (symbols) and Monte Carlo simulation (solid lines): (a) µs'/ µa = 5, and (b) µs'/ µa = 100. 56   3.5.2 Sensitivity Coefficients The sensitivity coefficients for the absorption and reduced scattering coefficients for the ODM, LTDM, ITDM, and RWDM were calculated as functions of the source-detector distance (ranging between 1 mfp' and 10 mm) for the corresponding transformed diffuse reflectance profiles [equations (3.5)-(3.8)]. Figure 3.2 shows the sensitivity coefficients as well as the diffuse -1 -1 - reflectance for two sets of optical properties with µa = 0.06 cm , µs' = 4 cm , and µa = 0.57 cm 1 -1 , µs' = 40 cm . In all of the four models, the magnitudes of the reduced scattering sensitivity coefficients are generally closer to those of the corresponding values of R. Moreover, the shapes of the two sensitivity coefficients are quite different. These observations show that values for the sensitivity coefficient of µs' are ‘large’ (i.e. on the order of R) and uncorrelated to those for the sensitivity coefficient of µa (different shapes), which are desirable conditions for estimating µs' in all four models with the exception of the sensitivity coefficients of µs' at the distance larger than 3 mm in Figure 3.2(a2). Furthermore, in general, the sensitivity coefficients for µs' are closer to the values of reflectance R at small source-detector distances than those at far sourcedetector distances, while the situation is different for the sensitivity coefficient of µa. This is because under the scattering dominant condition (µs'  µa), diffuse reflectance close to and far from the incident point of the light source does not have equal sensitivity to the optical properties of the medium. Overall, signals close to the source depend strongly on the reduced scattering property, whereas those far from the source exhibit large dependence on the absorption effect. This information is useful for determination of the minimum and maximum source-detector 57   distances for the curve fitting, and also provides a guide for the selection of appropriate data transformation and weighting methods. (a2) -2 4 x 10 Sensitivity coefficient (mm ) -2 Sensitivity coefficient (mm ) (a1) -3 3 2 1 0 4 6 8 Distance (mm) 10 0.3 0.2 0.1 0 2 Sensitivity coefficient (mm ) -2 -2 Sensitivity coefficient (mm ) 6 4 2 6 8 Distance (mm) 4 6 8 Distance (mm) 10 (b2) 8 4 ) 0.4 (b1) 0 ( 10 20 15 10 5 0 2 4 6 8 Distance (mm) -1 10 - Figure 3.2 Scaled sensitivity coefficients of the optical parameters [µa = 0.06 cm & µs' = 4 cm 1 -1 -1 for the left pane of plots (a1, b1, c1, d1), and µa = 0.57 cm & µs' = 40 cm for the right pane of plots (a2, b2, c2, d2)] as functions of the source-detector distance for ODM (a1, a2), LTDM (b1, b2), ITDM (c1, c2), and RWDM (d1, d2). Solid curves are for R, dash curves for µa and dot curves for µs' (a.u. = arbitrary unit). 58   Figure 3.2 (cont’d) (c2) 0.08 0.06 0.04 0.02 0 4 6 8 Distance (mm) 10 Sensitivity coefficient (a.u.) Sensitivity coefficient (a.u.) (c1) 0.15 0.1 0.05 0 2 1.5 1 0.5 0 4 6 8 Distance (mm) 10 (d2) Sensitivity coefficient (a.u.) Sensitivity coefficient (a.u.) (d1) 4 6 8 Distance (mm) 10 6 4 2 0 2 4 6 8 Distance (mm) 10 Values of the sensitivity coefficient for µa are relatively smaller compared with those for µs' due to much smaller values of µa than those of µs', especially for the ODM Figure 3.2a. However, with the use of the data transformation or relative weighting methods, values for the sensitivity coefficient of absorption have increased, as shown in Figure 3.2b, 3.2c, and 3.2d. This means that better estimations of µa can be achieved with the data transformation or relative weighting methods. Since in most cases, values of the reduced scattering sensitivity coefficient are still 59   larger in magnitude than those of the absorption sensitivity coefficient, µs', on average, would be estimated more accurately than µa for the diffusion model. 3.5.3 Estimation of Optical Parameters based on Monte Carlo Simulation Data Figure 3.3 shows the relative errors of estimating the absorption and reduced scattering coefficients of the turbid media by fitting ODM, LTDM, ITDM, and RWDM to the MC simulation data with the true values of optical properties shown in Table 3.1. The same fitting range with 1.5 mm minimum source-detector distance and 10 mm maximum source-detector distance for each group was used, which covered 86 detector positions with the resolution of 0.1 mm. The average relative errors for the 29 sets of optical properties extracted from ODM, LTDM, ITDM, and RWDM are 16.8%, 10.4%, 10.7%, and 11.4% for µa, respectively, and 8.1%, 6.6%, 7.0%, and 7.1% for µs', respectively. Better estimations of µs' were obtained, which are in agreement with the sensitivity analysis. The patterns of the relative errors in estimating µa and µs' obtained from these models are similar for the 29 simulation groups. Larger errors of µa for ODM appear in the groups with the largest and smallest mfp' for the same ratios of µs'/ µa, with the exception of µs'/ µa = 5. This is probably because the preselected 1.5 mm minimum sourcedetector distance is too small for those groups with large mfp' compared with the optimal position of 1-2 mfp' , and because the diffusion theory is unable to account for the nondiffuse component of the reflectance that is encountered for these small source-detector separations. Relatively large errors of estimating µs' with ODM were also obtained for those groups with the small value of mfp' [Figure 3.3(b)], due to the removal of the reflectance data close to the light 60   source. Therefore, it is important to select an optimal source-detector distance range in the experimental design and curve fitting process. It is also observed in Figure 3.3 that the error curves of µa and µs' are strikingly symmetric, which indicates that they are highly correlated. Hence when one of the parameters is overestimated, the other is likely to be underestimated in the nonlinear least squares. (a) a Relative error of  (%) 60 40 20 0 -20 -40 10 15 20 Group no. ( ) 25 30 35 25 30 35 Relative error of   (%) s (b) 20 10 0 -10 -20 -30 10 15 20 Group no. Figure 3.3 Relative errors of estimating 29 groups of µa (a) and µs' (b) shown in Table 3.1 by the original model, and the data transformation and relative weighting methods: ODM (·), LTDM (◦), ITDM (*) and RWDM (∆). The errors of estimating µa and µs' dramatically decrease when the data transformation and weighting methods are applied to most of the data groups. The logarithm transformation of the 61   raw data gives the best results with the smallest average errors of estimating the two optical parameters. However, the data transformation and weighting methods do not yield smaller errors -1 -1 -1 -1 for group 11 (µa = 3.00 cm , µs' = 30 cm ) and group 18 (µa = 2.00 cm , µs' = 40 cm ), when compared with those obtained using the original model. For these two cases, large values for the absorption and reduced scattering coefficients do not allow photons to travel farther away from the incident beam, resulting in extremely small reflectance at large source-detector distances. With small reflectance, the variability from MC simulations rapidly increases at large distances, which is further magnified when the data transformation and weighting methods are implemented. This suggests that it is also important to consider the signal-to-noise ratio (SNR) of the measured data used for the parameter estimation. To investigate the effect of the SNR of the experimental data on the parameters estimation, we examined the optical parameters extracted from the LTDM-generated diffuse reflectance data added with different levels of noise. Figure 3.4 shows absolute values of the relative errors of -1 estimating µa = 0.40 cm -1 and µs' = 20 cm under different noise levels. The error of estimating the optical parameters increases linearly with the noise level. If the maximum acceptable error is 5% in the parameters estimation, the noise level should be controlled to less than 5%. The errors of estimating µa are consistently higher than those of estimating µs' at different noise levels, which is consistent with the sensitivity analysis that µs' can be estimated more accurately. As the noise level increases, uncertainties in the MC simulation results also increase, which is clearly manifested by an unusually high estimation error at the noise level of 14% (Figure 3.4). Similar results and patterns were also obtained for the other groups of optical properties shown in Table 3.1. 62   Absolute values of relative error (%) 25 20 15 10 5 0 0 5 10 15 Noise level (%) Figure 3.4 Absolute values of the relative errors of estimating the optical properties (µa = 0.40 -1 -1 cm , µs' = 20 cm ) for different noise levels (each noise level includes 10 replications). Solid circles represent the absorption coefficient, while open circles are for the reduced scattering coefficient. 3.5.4 Statistical Analysis Nonlinear least squares estimates are based on a set of statistical assumptions. The above results were further tested by examining the validity of these assumptions. In this study, independent variable (source-detector distance r) is known precisely and nonstochastically, and all the “error” comes from the dependent variable (reflectance R). For different data transformation/weighting methods, the runs test (change of sign plus 1) of groups 25 and 30 was conducted based on the plots of the residual versus source-detector distance shown in Figure 3.5. The average numbers of runs of these two groups for the ODM, LTDM, ITDM, and RWDM are 4.5, 42, 3, and 41, respectively. The expected number of runs for the independent disturbances along the source-detector distance is about half of the measurements. There are 86 observations 63   for each group, and the numbers of runs from ODM and ITDM are considerably smaller than 86 / 2  43 , while the numbers of runs from LTDM and RWDM are close to the expected numbers of runs. Similar numbers of runs were obtained for other cases. It indicates that the residuals in ODM and ITDM are highly correlated, which violates the assumption of uncorrelated errors. Hence it is not appropriate to select the original model or integral transformation for the parameters estimation. The frequencies of residuals for groups 25 and 30 from ODM, LTDM, ITDM, and RWDM are shown in Figure 3.6. The distributions of the residuals approximately follow a normal distribution except for the one shown in Figure 3.6(c2). Further examination of the distributions of residuals was conducted for other groups of optical parameters; they follow similar patterns as shown in Figure 3.6: normal distributions are always observed for LTDM and RWDM, while this normal distribution assumption is violated by using ODM and ITDM in several cases. Hence, it can be concluded that the logarithm transformation and relative weighting methods are more suitable and reliable for accurate estimation of the absorption and reduced scattering coefficients of turbid media. x 10 (a) (b) 1 0 -1 -2 -3 0.2 -2 Residual (mm ) -2 Residual (mm ) 2 -4 0 2 4 6 8 Distance (mm) 0.1 0 -0.1 -0.2 10 0 2 4 6 8 Distance (mm) 10 Figure 3.5 Residual plots for the reflectance data versus source-detector distances from (a) ODM, -1 -1 (b) LTDM, (c) ITDM, and (d) RWDM with µa = 0.06 cm & µs' = 4 cm (dot curve) and µa = -1 -1 0.57 cm & µs' = 40 cm (asterisk curve). 64   Figure 3.5 (cont’d) x 10 (d) 0.2 -2 Residual (mm ) -2 Residual (mm ) 2 (c) -4 1 0 -1 -2 0.1 0 -0.1 -0.2 0   2 4 6 8 Distance (mm) 10 Frequency Frequency 20 10 10 30 20 10 0 -1 -2 0 2 -2 -4 Residual (mm ) x 10 (b1) -0.5 0 0.5 1 -2 -4 Residual (mm ) x 10 (b2) 15 15 Frequency Frequency 4 6 8 Distance (mm) 40 30 10 5 0 -0.1 2 (a2) (a1) 0 -4 0 -0.05 0 0.05 -2 Residual (mm ) 10 5 0 -0.2 0.1 -0.1 0 0.1 -2 Residual (mm ) 0.2 Figure 3.6 Residual histograms for the reflectance data from ODM (a), LTDM (b), ITDM (c), -1 -1 and RWDM (d) with 1) µa = 0.06 cm & µs' = 4 cm for the left pane of plots (a1, b1, c1, d1), -1 -1 and 2) µa = 0.57 cm & µs' = 40 cm for the right pane of plots (a2, b2, c2, d2). 65   Figure 3.6 (cont’d) 10 30 Frequency Frequency (c2) (c1) 5 0 -2 20 10 0 -15 -1 0 1 2 -4 Residual (a.u.) x 10 (d1) (d2) 15 Frequency Frequency 10 5 0 -0.1 -10 -5 0 5 Residual (a.u.) x 10-5 -0.05 0 0.05 10 5 0 -0.2 0.1 -2 Residual (mm ) -0.1 0 0.1 -2 Residual (mm ) 0.2 As examples, the estimated optical parameters for groups 25 and 30 (Table 3.1) extracted using the nonlinear least squares for the LTDM method are shown in Table 3.2. Figure 3.7 shows the corresponding 3D plots of sum-of-squares, which provide information about the characteristics of the standard error, confidence intervals and relative ease of convergence of the nonlinear regression routine (Mishra et al., 2008). It is observed that the sum-of-squares surface in Figure 3.7(a) is much shallower along both µa-axis and µs'-axis than in Figure 3.7(b), causing the difficult convergence of group 25 to find the lowest point on that surface, while it is easy to visualize the optimization of µa and µs' in Figure 3.7(b) for group 30. For the two groups, the 66   surfaces of both plots are shallower along the µa-axis than along the µs'-axis, resulting in the standard errors and confidence intervals for µa to be larger than those for µs', which are consistent with the results shown in Table 3.2. Better convergence can be obtained if the curve shows steeper changes along both µa-axis and µs'-axis. Table 3.2 Statistical results for estimating the optical parameters using the LTDM method. Group parameter no. 25 30 True value -1 (cm ) 0.06 µa Estimated -1 value (cm ) 0.07 Standard -1 Error (cm ) 0.015 Relative Error (%) 95% asymptotic confidence interval 16.7 [0.069, 0.072] [3.77, 3.81] 4.0 µs' 0.23 - 5.0 0.59 0.129 3.5 [0.560, 0.616] 40.0 µa 3.8 0.57 µs' 38.8 6.60 - 3.0 [37.34, 40.02] Sum of Squares (mm -4) (a) -7 x 10 4 7.4 3 7.2 2 0.382 -3 x 10 7  (mm -1) a 0.38 -1 s  (mm ) 0.378 0.376 6.8 Figure 3.7 3-D plot of sum of squares for (a) group 25 and (b) group 30 using the LTDM method.   67   Figure 3.7 (cont’d) Sum of Squares (mm-4) (b) -6 x 10 2 1 0 3.7 0.055 0.06 3.8 3.9 -1 s  (mm ) -1 a (mm ) 4 4.1 0.065   3.6 Conclusions In this research, an inverse algorithm combined with data transformation and weighting methods was applied in solving the inverse radiation transport problem for estimating the optical properties of biological materials. The results for the estimation of µa and µs' from the Monte Carlo generated reflectance data show that the minimum and maximum source-detector distances and signal-to-noise ratio should be considered in order to obtain sufficient and accurate information to uniquely determine µa and µs'. Sensitivity analysis demonstrates that the reduced scattering coefficient can be estimated more accurately than the absorption coefficient, which is also validated by the simulation results. The logarithm and integral transformation of the original data and the relative weighting method greatly improve the estimations of the two optical parameters with the relative errors of 10.4%, 10.7%, and 11.4% for µa, and 6.6%, 7.0% and 7.1% 68   for µs', which are much smaller than those obtained from the original diffusion model. Further statistical analysis shows that the logarithm transformation and relative weighting methods can improve the inverse algorithm to obtain more reliable estimations of the two optical parameters. The present work on parameter estimation, combined with statistical analysis, provides an appropriate means and guide in quantifying estimation errors for the optical properties. It suggests that attention should be paid on studying the characteristics and complexity of a mathematical model when it is used in an inverse algorithm for estimating their parameters. The methods proposed in this work should help us better interpret light propagation in biological materials and more accurately determine the relevant optical properties. 69   CHAPTER 4 OPTIMIZATION OF THE HYPERSPECTRAL IMAGING-BASED SPATIALLY-RESOLVED SYSTEM 4.1 Introduction Measurement of the optical properties of biological materials can be achieved using either direct or indirect methods. Direct methods (e.g., total transmittance and reflectance) are easy to carry out, but they are destructive or invasive. In contrast, indirect methods can be performed on intact samples nondestructively but need sophisticated instrumentation and complex mathematical models. Therefore, recent research has been focused on indirect methods, including time-resolved (TR), frequency-domain (FD), and spatially-resolved (SR). Based on radiation transfer theory, these indirect methods allow noninvasive or nondestructive measurement of the optical properties from intact biological samples. While TR and FD techniques have been extensively researched in the biomedical area, they may not be most suitable for food and agriculture because of expensive instrumentation, slow speed in measurement, and the requirement of good contact between the sample and detector. In comparison, spatially-resolved technique is less expensive in instrumentation, and it is easier to use and faster in measurement. The technique is, thus, potentially more viable for food and agriculture applications. A hyperspectral imaging-based spatially-resolved technique was recently developed in our lab to measure the optical properties of biological materials for a broad spectral range (i.e., 5001,000 nm) simultaneously (Cen and Lu, 2009; Qin and Lu, 2007). While the technique is promising for optical characterization of food and agricultural products, several critical issues need to be resolved. First, a reference method is needed to provide true values of the optical properties of samples, against which the hyperspectral imaging system can be evaluated and 70   validated. This would require model samples with known optical properties for a specific range of values. Several reference methods including transmittance, integrating sphere and empirical equation have been reported (Prahl et al., 1993; van Staveren et al., 1991); however, they have not been accepted as ‘standard’. Therefore, it is necessary to test and cross-validate these methods in order to provide accurate and reliable references for the optical system. Second, the hyperspectral imaging system needs to be such designed that it fully meets the requirements of the diffusion theory that governs light transfer in turbid biological materials. The hyperspectral imaging-based spatially-resolved system uses a continuous-wave point light beam to illuminate the sample. The shape and size of the beam can directly affect measurement accuracies. It is therefore important to examine and optimize the light beam. Third, an appropriate sourcedetector distance is critical for accurate estimate of the optical properties. Additionally, optimization of the inverse algorithm for the diffusion theory model is equally important in extracting the optical parameters from the reflectance profiles, which has been described in Chapter 3. This research was therefore aimed at addressing key technical issues in the development of the hyperspectral imaging-based spatially-resolved technique, so that it can provide accurate and reliable measurement of the optical properties of food and agricultural products. The specific objectives were to: a) evaluate and cross-validate three reference methods (i.e., transmittance, integrating sphere, and empirical equation) for determining the true values of optical properties of model samples, so as to establish a standard procedure of evaluating the hyperspectral imaging system; b) optimize the design of the light beam and source-detector distance in the hyperspectral imaging system, using Monte Carlo (MC) simulation and experiments for model samples, to improve the measurement of the absorption and reduced scattering coefficients; and c) 71   assess the performance of the hyperspectral imaging system in terms of accuracy, precision/reproducibility, and sensitivity. 4.2 Materials and Methods 4.2.1 Hyperspectral Imaging-based Spatially-resolved Technique In the spatially-resolved method, a small continuous-wave light beam perpendicularly illuminates the sample’s surface, and the remitted light is measured at different distances from the light source to obtain a spatially-resolved reflectance profile (see Figure 2.4). A laboratory hyperspectral imaging system in line scan mode used for acquiring spatiallyresolved diffuse reflectance profiles from a sample is schematically shown in Figure 4.1(a). The system mainly consisted of a high performance 512×512 pixel camera with a back-thinned illuminated CCD detector (Model C4880-21-24A, Hamamatsu Photonics Systems, Bridgewater, NJ, USA) covering the wavelengths of 200-1,100 nm, an imaging spectrograph (ImSpector V10, Spectral Imaging Ltd., Oulu, Finland) covering the wavelength region of 400-1,000 nm, a zoom lens with the focal length of 11-110 (10×) mm (Zoom 7000, Navitar Inc., Rochester, NY, USA), a quartz tungsten halogen light source (Oriel Instruments, Stratford, CT, USA) with a feedback controller, an optical fiber coupled with a focusing lens for delivering a small beam to the sample, and a computer installed with a frame grabber board for controlling the camera and acquiring images. The hyperspectral imaging technique, combining imaging and spectroscopy methods, acquires both spectral and spatial information simultaneously, and it is, therefore, ideally suitable for measuring spatially-resolved diffuse reflectance profiles for a broad wavelength region. To improve the repeatability of measurements, 10 line scans were taken from each sample for every 1 mm horizontal displacement over a range of 9 mm as the sample was moving at the 72   predetermined velocity during the image acquisition [Figure 4.1(b)]. The sample movement and positioning were controlled by the two motorized stages, and the distance between the sample and the lens was 243 mm. Spectral and spatial calibrations and nonuniform instrument response corrections for the hyperspectral imaging system were performed by following the procedure described in the literature (Lu and Chen, 1998; Qin and Lu, 2007). (a) Figure 4.1 Hyperspectral imaging-based spatially-resolved method: (a) schematic showing the major components of the system; (b) top view of multi-line scanning mode for acquiring spatially-resolved reflectance profiles.           73   Figure 4.1 (cont’d) (b) Scan line Point light source Sample moving direction Figure 4.2(a) shows a typical hyperspectral reflectance image for a liquid model sample made up of Intralipid fat emulsion as scatterers and blue dye as absorbers (see more details in Section 4.2.2). Each horizontal line taken from the image represents one spatially-resolved reflectance profile at a specific wavelength, as shown in Figure 4.2(b). Hence the reflectance image which had a calculated spectral resolution of 4.55 nm by the spectral calibration, in effect, consisted of more than 100 spatially-resolved reflectance profiles for the wavelengths of 5001,000 nm. Since the spatially-resolved reflectance profiles were symmetric to the light incident point, the two sides were averaged in the extraction of optical properties. Each spatially-resolved reflectance profile was then fitted by the diffusion model with the inverse algorithm, from which the spectra of absorption and reduced scattering coefficients for 500-1,000 nm were obtained. 74   (a) (b) 12000 Intensity(CCD count) Wavelength (nm) 400 600 800 -5 0 5 10 Distance (mm) 700 nm 8000 6000 4000 570 nm 2000 0 -10 1000 -10 10000 -5 0 5 Distance(mm) 10 Figure 4.2 Hyperspectral reflectance image of a liquid model sample: (a) 2-dimensional (2-D) display with intensities being indicated by pseudo colors, and (b) original spatially-resolved reflectance profiles at 570 nm and 700 nm. 4.2.2 Model Samples Preparation To optimize and evaluate the hyperspectral imaging system, 12 liquid model samples were created with two different dyes (Direct Blue 71 and Naphthol Green B, Sigma-Aldrich Inc., St. Louis, MO, USA) and their mixture as absorbers and fat emulsion (Intralipid-20%, SigmaAldrich Inc., St. Louis, MO, USA) as scatters for different concentration levels. The three dye absorbers were used in the model samples to allow fine adjustment of the absorption coefficient and comprehensive investigation of the accuracy for estimating µa over the entire spectral region of 500-1,000 nm because each dye had high absorption for a specific wavelength range. For each absorber, aqueous absorbing stock solution with the concentration of 1 mg/mL was first prepared, and it was then added to 250 mL distilled water followed by adding different volumes of Intralipid (10-50 mL) into the absorbing solutions to form the model samples with specified values of µa and µs'. The ranges of µa and µs' values for these model samples were 0.0-0.934 -1 -1 cm and 7.0-39.9 cm , respectively, at the wavelength range of 500-1,000 nm. 75   The model samples were first measured using the hyperspectral imaging system [Figure 4.1(a)]. Each sample was scanned in a 600 mL cylindrical beaker with 80 mm diameter and the solution height of approximately 50 mm, which can be considered semi-infinite. Prior to imaging acquisition, the liquid model samples were gently stirred by a magnetic stirrer until the absorbers and scatters were distributed homogenously. The beaker containing the liquid model sample was placed on a flat platform fixed on the motorized stage, and the images were then acquired after the solution was stable. After the imaging, confocal laser scanning microscopy (CLSM) images of these model samples were acquired for gaining an insight of the distribution and interaction of the dye and Intralipid particles in the solution. This is because when the dye and Intralipid were mixed, the two components could interact with each other, which might in turn affect the stability of its mixture and thus of the optical properties (Ren et al., 2007). Finally, the true optical properties of these model samples were measured using the reference methods discussed below. 4.2.3 Reference Methods The true values of µa and µs' of the model samples were needed to validate the hyperspectral imaging system. However, in the real situation, it is impossible to obtain true µa and µs' of the model samples. Therefore, nominal true values (or reference values) were derived as reliable estimates of the true values by using three commonly used methods, and evaluation and crossvalidation of these methods were performed. A collimated transmittance method was first employed for the standard absorption measurement using a miniature fiber optic spectrometer (USB4000-VIS-NIR, Ocean Optics, Dunedin, FL, USA) with a 10-mm pathlength quartz cuvette sitting in a cuvette holder (Figure 76   4.3). The true values of µa at the wavelengths of 500-1,000 nm were calculated by the BeerLambert law  a   ln Is  Id Ir  Id (4.1) where Is is the sample intensity, Id is the dark intensity, and Ir is the reference intensity (empty cuvette). The accuracy of the transmittance measurement was experimentally quantified by testing In-Spec® certified UV and visible standards (Standard 0.10, 0.5, 0.8, and Background solution, GFS Chemicals, Inc., Columbus, OH, USA) with known µa at selected wavelengths provided by the manufacturer. Spectrometer Light source Cuvette holder with sample Computer Figure 4.3 Collimated transmittance measurement for absorption using a miniature fiber optic spectrometer.  For each absorber, eight concentration levels of pure absorbing solutions with the absorption -1 coefficient varied from 0.0 to 1.800 cm for the wavelength range of 500-1,000 nm were made for transmittance measurements, which covered the range of µa for the model samples. The linear relationship between the absorption coefficient and the concentration of each absorber was 77   established at each wavelength, and the absorption coefficient was then derived based on the concentration of the absorber in the model samples. A 152.4 mm diameter integrating sphere (RT-060-SF, Labsphere, Inc., North Sutton, NH, USA) with a 12.7 mm diameter detector port and a 25.4 mm sample port was used for collecting the total transmittance and reflectance to estimate true values of the absorption and reduced scattering coefficients of the model samples (Figure 4.4). The sphere was coated with ® Spectraflect diffuse reflectance materials with 98% reflectance at 500-1,000 nm. A quartz tungsten halogen light source was used and the light received from the sample port of the integrating sphere by a 400 µm diameter fiber was delivered to the miniature fiber optic spectrometer, the same one used in the collimated transmittance measurement (Figure 4.3). Each model sample held in a 2 mm pathlength cuvette was placed in the sample holder. For measuring total diffuse reflectance [Figure 4.4(a)], the sample was placed at the sample port and the light beam was arranged at the beam entrance port. For total transmittance measurements [Figure 4.4(b)], the sample was placed at the sample port with other ports covered by the standard masks, and all light transmitted through the sample was collected by the detector. The same beam size was used for both reflectance and transmittance measurements to ensure measurement accuracy. (b) Detector (a) Detector Rsample Sample Specular exclusion Tsample Light beam Light beam Sample Figure 4.4 Single integrating sphere configurations for: (a) total diffuse reflectance measurement, and (b) total transmittance measurement. 78   After completion of the total transmittance and reflectance measurements, an inverse adding-doubling algorithm developed by Prahl et al. (1993) was used to calculate µa and µs' values of the model samples based on the total reflectance and transmittance. Prahl et al. (1993) reported that the accuracy of the integrating sphere with the adding-doubling algorithm for determining the optical properties (µa & µs') was 2-3% for most reflectance and transmittance values. However, a 1% variation in the reflectance and transmittance values would result in errors of 0.4% for µs' and 17% for µa. This suggests that the integrating sphere method could be problematic in estimating µa. For estimating the µs' of Intralipid scatterers, an empirical equation derived by van Staveren et al. (van Staveren et al., 1991) is widely used, which has been experimentally proved to have accuracies better than 6% in the wavelength range of 400-1,100 nm. Therefore, µs' values obtained from the integrating sphere were also compared with those from the empirical equation given below s '  (928 1.4  160 2.4 )C % (4.2) -1 where µs' is in cm ,  is the wavelength in µm, and C% is the percent Intralipid concentration. For easy comparison of the results obtained from the hyperspectral imaging system and the three ∆ reference methods, the absorption and reduced scattering coefficients are denoted as µa & µs' * * o ∆ for the hyperspectral imaging system, µa & µs' for the integrating sphere, µa for the o transmittance measurement, and µs' for the empirical equation. 79   4.2.4 Optimization of the Light Beam and Source-detector Distance To achieve accurate measurement of the optical properties, optimization of the hyperspectral imaging system was performed for two key factors, i.e., light beam and source-detector distance. In the experimental setup, a finite size beam was used to illuminate the sample, which deviates from the solution of the diffusion model that was derived for an infinitely small beam. Therefore, the effect of the incident beam on the determination of the optical properties was investigated based on MC simulation, which offers a flexible and accurate approach for quantifying the optical features of light transport that are difficult to measure experimentally. In the simulation experiment, a semi-infinite medium was described by its refractive index (n), absorption coefficient (µa) and reduced scattering coefficient (µs'). The spatially-resolved reflectance generated by MC simulation for an infinitely small photon beam was first convolved for the finite size beam using the programs ‘MCML’ and ‘CONV’ developed by Wang et al. (1997, 1995b). Two specific beams, i.e., Gaussian and circularly flat, were investigated. Meanwhile, the error of estimating µa and µs' caused by the finite size Gaussian beam with the diameter of 0.1 ≤ d ≤ 2 mm was quantified. Six sets of µa and µs' with the ranges of 0.060 ≤ µa ≤ -1 2.000 cm -1 and 4.0 ≤ µs' ≤ 40.0 cm and with the ratios of µs'/µa = 20 and µs'/µa = 70 were investigated. In conjunction with the MC simulations, actual beam profiles for the hyperspectral imaging system were measured and characterized. This was accomplished by moving and imaging the beam for every 0.12 mm distance (determined by the width of scan line) in the direction perpendicular to the scan line. Twelve 512×512 pixel images were acquired with the resolution of 0.1 mm along the scan line direction, and the 3-D beam profiles were then reconstructed from 80   these images at each wavelength, from which the beam profile type (i.e., Gaussian and flat) and beam circularity were determined. The circularity of the beam was calculated based on the roundness (Rd), which is given by Rd  4 A / P 2 (4.3) where A and P are the beam area and perimeter, respectively, and Rd = 1 for circle. Accurate information on the source-detector distance, including minimum source-detector distance (rmin), maximum source-detector distance (rmax) and spatial resolution, is required for determining the range of the spatially-resolved reflectance profile. MC simulation experiments were first conducted to investigate the effects of the minimum and maximum source-detector -1 distances on estimating 29 sets of µa and µs' with 0.040 ≤ µa ≤ 3.000 cm and 4.0 ≤ µs' ≤ 40.0 -1 cm . In experimental measurements, signal-to-noise ratio (SNR) should be considered when determining the optimal maximum source-detector distance since SNR varies at different radiation locations and decreases with the decreasing signal. Ten images from a model sample were acquired at the same location. The mean profile divided by the standard deviation of the CCD count at each wavelength was considered as an estimate of SNR. Sufficient data points for the spatially-resolved reflectance profile are needed to obtain good estimates of µa and µs'. Hence the effect of spatial resolution on calculating the optical coefficients was also examined. The reflectance profile with the resolution of 0.01 mm for µa = -1 1.000 cm -1 and µs' = 20.0 cm was firstly obtained from the diffusion model. Seven new reflectance profiles with the spatial resolution ranging from 0.07 mm to 0.25 mm were then derived from the original profile by integration of the light intensity over the area covered 81   between the two adjacent data points divided by the corresponding spatial resolution. These seven profiles were fitted by the diffusion model to obtain the estimated µa and µs', and they were then compared with the true values of µa and µs' obtained from the reflectance profile with the resolution of 0.01 mm. 4.2.5 Assessment for Accuracy, Precision/Reproducibility, and Sensitivity Assessment of the performance of the system is an important step to determine its validity, reliability, and robustness. The hyperspectral imaging system was evaluated for accuracy, precision/reproducibility, and sensitivity. The accuracy of a measurement is defined as the closeness of agreement between the measured values and the nominal true values of parameter  (αmeas & αtrue) such as µa and µs' obtained from the reference methods under optimal experimental conditions. Accuracy can be quantified by the relative error of the measurement, given in equation (4.4)   meas   true  100%  true (4.4) Accuracy is important for the absolute measurement, i.e., for simulating light transport and deriving the physiological information about the biological tissue. Precision/reproducibility is a term frequently used to characterize the random error and describe the system consistency of measurement over a long time period. It quantifies how the system is self-consistent over different times, and it is particularly important for a long time experiment like fruit quality monitoring during postharvest. Repetition of the measurement on the same model samples under the same experiment conditions over a period of four days was 82   performed. The precision of the system was evaluated by calculating the coefficient of variation (CV) with respect to the average values calculated over the entire experiment. Sensitivity is described by the minimum detectable value under an acceptable noise level, and it determines the detection limit. Clearly, the noise level is affected by the signal intensities and detection efficiencies, and thus the sensitivity of the system depends on the amount of signal acquired at each measurement point. Sensitivity analysis is especially important in measuring biological materials with small values of absorption coefficient. The analysis was performed on a blue dye model sample which had minimal levels of the absorption coefficient for the wavelengths of 500-1,000 nm. 4.3 Results and Discussion 4.3.1 Accuracy of the Reference Methods Figure 4.5 shows the measured absorption spectra of the three standard solutions by transmittance as well as the true absorption values provided by the manufacturer. The average error of µa by the transmittance measurement was 3% at 500-850 nm, and less than 5% at 850-1 900 nm with the minimum detectable value of 0.050 cm . The performance of the method for measuring µa above 900 nm was not evaluated because the true absorption values of these standard solutions were not available. However, it was found that the water absorption at 970 nm -1 measured by the transmittance method was 0.414 cm , which is within 1.4% of the reported data (Hale and Querry, 1973). Hence, it was concluded that the transmittance method gave accurate measurement of the absorption coefficient with errors less than 5%, and it was therefore appropriate for determining µa of the model samples. 83   2.5 µa (cm-1) 2 1.5 1 0.5 0 500 600 700 800 Wavelength (nm) 900 Figure 4.5 Spectra of the absorption coefficient for the three standard solutions measured by the transmittance method [asterisks () for standard 0.1, triangles () for standard 0.5, and squares () for standard 0.8] and their true vlaues (solid lines). Integrating sphere measurements of µs' for the 12 model samples were compared with the calculated values from the empirical equation. It was found that the empirical equation, on average, overestimated µs' measured with the integrating sphere method by 9% for the wavelengths of 500-900 nm when the Intralipid concentration ranged between 0.8% and 3.2%. However, at low Intralipid concentrations (i.e., <1.2%), the measured values of µs' from the integrating sphere matched well the calculated values from the empirical equation with the differences being less than 4.5% (Figure 4.6). As the concentration of the Intralipid increased, the measured value deviated more from the calculated value. This may be because the relationship between the reduced scattering coefficient and the particle concentration at high density Intralipid solutions was no longer linear, as reported by Giusto et al. (2003) and Zaccanti et al. (2003). The empirical equation was derived based on the linear relationship between µs' and the particle concentration. In addition, variation in batches of Intralipid might have also 84   existed, which, in turn, could have introduced errors in calculating µs' values. In view of these results, it was concluded that the integrating sphere method would give more accurate measurement of µs' than the empirical equation for a wider range of µs', and it was therefore used Average difference of µs' (%) to determine the true µs' of the model samples. 20 16 12 8 4 0 0.5 1 1.5 2 2.5 3 3.5 Concentration of Intralipid (%) Figure 4.6 Average differences of the reduced scattering coefficient for the model samples at different concentrations of Intralipid for 500-900 nm obtained from the integrating sphere measurement and the empirical equation. 4.3.2 Light Beam Characteristics Since a Gaussian beam, like most optical beams, does not have sharp physical edges, the beam size is usually determined between the two points that contain a selected percentage of the 2 ‘useful’ energy. In this study, the size of Gaussian beam was defined by the 1/e diameter, where the beam’s power is at 13.5% of the maximum height. Figure 4.7 shows the spatially-resolved reflectance profiles obtained for an infinitely small beam and flat and Gaussian beams of different sizes through Monte Carlo simulations. The Gaussian beam and the flat beam nearly 85   gave the same results for the same beam diameter of d = 1 mm. Finite circular beams with different sizes generated reflectance similar to that produced by the infinitely small beam when the source-detector distance was larger than the beam size (r > d). The average errors of estimating six sets of µa and µs' from the reflectance generated by the Gaussian beam of different sizes are presented in Figure 4.8(a). The error produced by the finite beams relative to the infinitely small beam was less than 1% for the beam size of less than 0.5 mm; it increased linearly with the beam size larger than 0.5 mm. Generally, a 1-mm light beam introduces around 5% error of estimating µa and µs' compared to the infinitely small beam. These results suggest that the beam in the system should be less than 1 mm in size in order to control the error to within 5%. The error introduced by a finite beam was also influenced by the value of the optical properties, as shown in Figure 4.8(b) and (c). The error curves are relatively flat for small values of µa and µs', which indicates that beam size had less effect on the smaller µa and µs' than on large µa and µs'. 86   0 10 Infinitely small beam Flat beam with d=1mm Gaussian beam with d=1mm Gaussian beam with d=2mm Reflectance (mm-2) -1 10 -2 10 -3 10 -4 10 0 2 4 6 Distance (mm) 8 10 Figure 4.7 Comparison of spatially-resolved reflectance produced by an infinitely small beam, a flat beam with diameter of d = 1 mm, and a Gaussian beam with d = 1 and 2 mm. (a) Relative error (%) 15 a s  10 5 0 -5 0 0.5 1 1.5 Beam size (mm) 2 Figure 4.8 Error analysis for different beam sizes: (a) average errors of estimating six sets of -1 -1 optical properties with 0.060 ≤ µa ≤ 2.000 cm and 4.0 ≤ µs' ≤40.0 cm , and the ratios of µs'/µa = 20 and µs'/µa = 70; and (b) and (c) relative errors for three sets of optical properties with -1 -1 increased µa and µs' from A to C with µa = 0.006, 0.029, 0.057 cm and µs' = 0.4, 2.0, 4.0 cm . 87   Figure 4.8 (cont’d) (c) (b) 30 A B C 30 Relative error of  s  (% ) Relative error of  a (%) 40 20 10 0 0 0.5 1 1.5 Beam size (mm) 2 A B C 20 10 0 0 0.5 1 1.5 Beam size (mm) 2 Simulation results for the light beam showed that it is necessary to understand beam characteristics (i.e., beam profile, shape, and size) in the experimental measurements. Figure 4.9 gives an example of the measured 3-D beam profiles and 2-D contours at the wavelengths of 650 nm and 950 nm. The beam at the visible and short near-infrared region had a good Gaussian distribution and its shape was circular with the roundness Rd = 0.986 (  1 ). Based on the commonly accepted method for defining the size of Gaussian beam, the beam size in the current setup was 0.83 mm, which would have contributed to less than 4% error in estimating µa and µs' based on the simulation experiments. Although a smaller beam is preferred, other factors such as light intensity or throughput and measurement repeatability also need to be considered in the optical design. 88   (b) Intensity (CCD count) Intensity (CCD count) (a) 650nm 10000 5000 0 1 0.5 0 D2 (mm) -0.5 -1 0 -1 0.5 1 950nm 3000 2000 1000 0 1 0.5 0 D2 (mm) -0.5 -0.5 D1 (mm) (c) -1 -1 0.5 1 -0.5 D1 (mm) (d) 1 1 650nm 950nm 0.5 0.5 D2 (mm) D2 (mm) 0 0 -0.5 0 -0.5 -1 -1 -0.5 0 D1 (mm) 0.5 -1 -1 1 -0.5 0 D1 (mm) 0.5 1 Figure 4.9 three-dimensional profiles (a and b) and 2-dimensional (c and d) contours of the incident light beam at wavelengths of 650 nm and 950 nm, where D1 is the direction along the scan line and D2 is perpendicular to the scan line. 4.3.3 Source-detector Distance Two examples of the relative errors between the true (values selected for Monte Carlo -1 simulations) and estimated values of µa = 0.290 cm -1 µs' = 30.0 cm -1 -1 & µs' = 20.0 cm , and µa = 0.430 cm & at different ranges of the source-detector distance are shown in Figure 4.10. For each set, 23 fittings were performed with 0.2 ≤ rmin ≤ 4 mm & rmax = 10 mm, and 4 ≤ rmax ≤ 10 mm & rmin = 1.5 mm. The estimated parameters from the reflectance data containing the smallest minimum source-detector distance (rmin = 0.2 mm) resulted in the largest errors for 89   both µa and µs' [Figure 4.10(a1) and (a2)]. Values of the absorption coefficient were systematically underestimated, while those of the reduced scattering coefficient were overestimated. These results may indicate the failure of the diffusion model which does not account for the nondiffusive component of the reflectance at these small source-detector distances. As the minimum source-detector distance was increased, the errors decreased, reaching minimum for rmin  0.5 mm (1 mfp', i.e., transport mean free path) to 2 mm (4 mfp') [1 -1 mfp' = (µa+µs') ], depending on the values of µa and µs'. For the minimum source-detector distance from 2 to 4 mm, the errors were relative stable. Also, the symmetrical curves of the relative errors of µa and µs' indicated that these errors were highly correlated. Patterns similar to Figure 4.10 were obtained for a wide range of the optical parameters. 8 Relative error (%) Relative error (%) 10 (a1) 12 4 0 -4 0 1 2 3 4 -8 -12 Minimum source-detector distance (mm)   (a2) 6 2 -2 0 1 2 3 4 -6 -10 Minimum source-detector distance (mm)   Figure 4.10 Relative errors of estimating µa (squares) and µs' (asterisks) from the spatiallyresolved reflectance data generated by Monte Carlo simulations when using different minimum -1 -1 and maximum source-detector distances: (a1) and (b1) µa = 0.290 cm & µs' = 20.0 cm , and -1 -1 (a2) and (b2) µa = 0.430 cm & µs' = 30.0 cm .   90   Figure 4.10 (cont’d) 20 (b1) Relative error (%) Relative error (%) 8 4 0 4 5 6 7 8 9 10 -4 -8 Maximum source-detector distance (mm) 10 0 4 5 6 7 8 9 10 -10 -20   (b2) Maximum source-detector distance (mm)   The optimal minimum source-detector distance of rmin  1-4 mfp' found in this study is different from the study of Nichols et al. (1997) (0.75-1 mfp') but is consistent with the finding of Farrell et al. (1992) that the minimum source-detector distance should be larger than 1 mfp'. However, in practice, there is no a priori information about the exact values of µa and µs' over a specified wavelength range, and thus it is difficult to determine the optimal minimum sourcedetector distance based on these loose criteria. But considering the average errors of 29 sets of µa and µs' that were investigated in the simulation experiments, the smallest average estimated errors of µa and µs' were obtained when rmin = 1.5 mm. Experimentally, the minimum sourcedetector distance should also be larger than the incident beam size (rmin > 0.83 mm). Hence, the 1.5 mm minimum source-detector distance selected in the system was considered optimal. Similar error patterns were found for the maximum source-detector distance, as shown in Figure 4.10(b1) and (b2). With rmax changing from 4 mm to 10 mm, the error varies between 0.03%-19% for µa, and 0.09%-13% for µs'. The optimal maximum source-detector distance 91   varied with the values of µa and µs', approximately between 10-20 mfp', which is also in agreement with Nichols et al. (1997) and Pham et al. (2000). Hence, the maximum sourcedetector distance should be optimized to obtain accurate estimates of µa and µs'. In practice, the optimal maximum source-detector distance may be largely determined by the SNR. As shown in Figure 4.11, to control the minimum SNR of 20 (corresponding to 5% of the variation coefficient), the signal level should be greater than 150 CCD counts at 650 nm with the maximum source-detector distance around 7.7 mm. A similar threshold of the signal level was obtained at other wavelengths with 150±12 CCD counts, while the maximum source-detector distance varied at different wavelengths for this signal level. Therefore, the threshold of 150 CCD counts (SNR = 20) was used in the curve fitting to determine the maximum source-detector distance at each wavelength. (a) Reflectance (CCD count) 5 10 4 10 3 10 2 10 1 10 2 3 4 5 6 7 Distance (mm) 8 9 10 Figure 4.11 Signal-to-noise ratio measurement of the hyperspectral imaging system: (a) average spatially-resolved reflectance profile of 10 measurements of a model sample at 650 nm, and (b) signal-to-noise ratio of the measurements within 10 mm source-detector distance. 92   Figure 4.11 (cont’d) (b) Signal-to-noise ratio 500 400 300 200 100 0 2 3 4 5 6 7 Distance (mm) 8 9 10 Further, the change of spatial resolution from 0.07 mm to 0.25 mm introduced an insignificant error of 0.53% for µa and 0.36% for µs' relative to the optical properties obtained for the resolution of 0.01 mm, as shown in Figure 4.12. In the system setup, the highest spatial resolution that the imaging sensor can achieve was 0.1 mm without binning. To reduce measurement time and also enhance the signal-to-noise ratio, the 0.2 mm spatial resolution with 2×2 binning was selected in our measurements since high resolution is not critical in this case. Relative error (%) 0.6 0.4 µa µs' 0.2 0.0 0.05 -0.2 0.10 0.15 0.20 0.25 -0.4 -0.6 Spatial resolution (mm) -1 -1 Figure 4.12 Errors of estimating µa = 1.00 cm and µs' = 20.0 cm introduced by different spatial resolutions relative to the optical properties obtained for the resolution of 0.01 mm. 93   4.3.4 Accuracy, Precision/Reproducibility, and Sensitivity Figure 4.13 shows the measured absorption and reduced scattering spectra for three model samples over 500-1,000 nm, determined by the hyperspectral imaging system and the reference methods. As can be seen from the plots on the left pane of Figure 4.13 (i.e., a1, b1 and c1), the shape of the absorption spectra for the three dyes obtained by transmittance, integrating sphere and hyperspectral imaging system was similar. They all showed a blue dye absorption peak at 585 nm [Figure 4.13(a1)], green dye absorption peak at 714 nm [Figure 4.13(b1)], mixed dye absorption peak around 600 nm [Figure 4.13(c1)], and water absorption at 970 nm in all three * plots. However, the µa values measured by the integrating sphere (µa ) were far off compared ∆ with those from transmittance (µa°) and hyperspectral imaging (µa ). Similar results were reported by Saeys et al. (2008) when they compared the integrating sphere measurement with the time-resolved method. It could be due to the fact that light losses at the sample edge are not accounted for in the integrating sphere measurement, and also, µa is very sensitive to the measured reflectance and transmittance according to Prahl et al. (1993) and Chen et al. (2006). Hence, the accuracy of hyperspectral imaging measurements for µa was evaluated against the transmittance method. The average errors of µa were calculated at 530-700 nm for the blue dye, 600-850 nm for the green dye, and 530-800 nm for the mixed dye, because the transmittance method was unable to measure µa at very small levels for some wavelengths. Overall, the average error of estimating µa for all model samples was 23% at 530-850 nm, and 16%, 26% and 26% for samples with blue, green, and mixed dyes, respectively. For the three model samples with small amount of blue dye, the absorption caused by the blue dye was almost zero at 750 nm, 94   and thus the absorption at 750 nm for these samples was mainly attributed to the water. The -1 reported results show that water absorption at 750 nm is between 0.0247 to 0.0286 cm (Prahl, -1 1998), and our result of µa = 0.0275 ± 0.0006 cm falls in this range and is only 1% different -1 compared with the recent report from Martelli et al. (2007) with µa = 0.0278 cm . (a1) (a2) µa (cm-1) µs' (cm-1) 1.4 Transmittance 1.2 Hyperspectral imaging 1.0 Integrating sphere 0.8 0.6 0.4 0.2 0.0 -0.2 500 600 700 800 900 1000 Wavelength (nm) 20 18 16 14 12 10 8 6 Empirical equation Hyperspectral Imaging Integrating sphere 500 1.4 1.2 1 0.8 0.6 0.4 0.2 0 (b2) 20 18 16 14 12 10 8 6 Transmittance Hyperspectral imaging Integrating sphere Empirical equation Hyperspectral imaging Integrating sphere µs' (cm-1) µa (cm-1) (b1) 600 700 800 900 1000 Wavelength (nm) 500 500 600 700 800 900 1000 Wavelength (nm) 600 700 800 900 1000 Wavelength (nm) Figure 4.13 Spectra of absorption and reduced scattering coefficients of three model samples with (a) blue dye, (b) green dye, and (c) mixed dye as absorbers measured by the hyperspectral imaging and reference methods. 95   Figure 4.13 (cont’d) 1.4 1.2 1 0.8 0.6 0.4 0.2 0 (c2) 20 18 16 14 12 10 8 6 Transmittance Hyperspectral imaging Integrating sphere 500 Empirical equation Hyperspectral imaging Integrating sphere µs' (cm-1) µa (cm-1) (c1) 500 600 700 800 900 1000 Wavelength (nm) 600 700 800 900 1000 Wavelength (nm) For most samples, µa from the hyperspectral imaging measurement, on average, was underestimated. A possible explanation for the discrepancies is that the dye might have interfered with the Intralipid particles. The confocal laser scanning microscopy (CLSM) images of the Intralipid solutions with blue and green dyes in Figure 4.14 showed that some dye particles, which were much smaller than Intralipid particles, went into the lipid vesicles in the Intralipid, resulting in a low dye concentration in the solution. It was also observed during the CLSM experiment that the Intralipid particles were relatively stable in the blue dye sample, while they were actively moving in the green dye sample. Since unstable scattering particles in the model samples could have an effect on the measurement, this may explain why better estimations of µa were obtained from the blue dye samples compared with those from the green and mixed dye model samples. It was also found that the hyperspectral imaging-estimated absorption above 900 nm was consistently lower for all the samples, particularly at the water absorption band of 970 nm. This might have been caused by the beam position change due to the different diffraction of 96   the focusing lens in the near-infrared region, which needs to be confirmed in further research. Therefore, the absorption above 900 nm was not used in calculating the accuracy. (a) (b)     Figure 4.14 Confocal laser scanning microscopy images of Intralipid solutions with (a) blue dye, and (b) green dye. The spectra of µs' presented in Figure 4.13(a2), (b2) and (c2) show that the measured values ∆ * (µs' ) matched well with those obtained from the integrating sphere measurement (µs' ) at 500900 nm, and lower than those calculated from the empirical equation (µs'°). The average error of estimating µs' was 7% at 500-900 nm compared with the integrating sphere measurement. The divergence of µs' above 900 nm appeared between the integrating sphere method and the wavelength-dependent exponential decay function of µs'. It is difficult to know which method is more accurate for determining µs' at the wavelength above 900 nm because no published paper reported the µs' values of Intralipid in these wavelengths. However, this wavelength-dependent 97   exponential decay function of µs' could be violated with the variation of the refractive index at long wavelengths according to the study of Mourant et al. (1997). When comparing our measurement accuracy with other techniques or other spatiallyresolved methods, one should notice the differences of experimental conditions, model samples and wavelength ranges. Bays et al. (1996) reported 32% accuracy of µs' at 633 nm for polyoxymethlene model samples by using a spatially-resolved reflectometry. Pifferi et al. (2005) applied a general protocol for assessing several optical methods for determining optical properties, and they reported large differences from different instruments in measuring the same model samples with the maximum discrepancies of 32% for µa at 970 nm and 41% for µs' at 820 nm. Spichtig et al. (2009) used the frequency-domain technique to measure the optical properties of Intralipid model samples with accuracies being less than 10% for µa and 31% for µs' at eight wavelengths. While covering a wider wavelength range of optical properties, the hyperspectral imaging-based spatially-resolved method has achieved superior results for measuring µa and µs', compared with these reported studies that were only conducted at single or several wavelengths. The system precision or reproducibility was calculated by repeated measurements of a blue dye model sample at four different days under the same experimental conditions with a 30 min warm-up time and the same level of source power and acquisition time. The coefficient of variation (or CV) in the absorption peak at 585 nm was 2.4% for µa and 3.8% for µs' with the maximum discrepancies of 4.9% and 6.3%, respectively. However, the model sample could have undergone some small changes over the four days, thus inducing some variations in the optical 98   properties of the tested sample. Hence the system has demonstrated good reproducibility and consistency over a prolonged time period. Ten measurements were performed on a model sample made up with the blue dye as the absorber because its minimum absorption coefficient was almost zero for the wavelengths of 700-900 nm. The result of the sensitivity test on µa is presented in Figure 4.15, with the CV of µa for the 10 measurements being plotted as a function of µa for different wavelengths. The -1 minimum detectable value of µa was around 0.0082 cm for a noise level of 10%. The sensitivity of measuring µs', as determined by the CV values, was always less than 4% because -1 µs' was much larger than µa for the investigated range of 7.0 ≤ µs' ≤ 40.0 cm . 3 Coefficient of variation of a (%) 10 2 10 1 10 0 10 0 0.02 0.04 0.06 0.08 0.1 0.12 -1 a (cm ) Figure 4.15 Coefficient of variation versus reordered ascending absorption coefficients of a model sample at different wavelengths. 99   4.4 Conclusions This research examined critical factors in the development of a hyperspectral imaging-based spatially-resolved technique (i.e., methods for providing reference measures of the optical properties, light beam and source-detector distance, etc.) for determining the absorption and reduced scattering coefficients of biological materials over the wavelengths of 500-1,000 nm. Monte Carlo simulations, coupled with experiments for model samples, demonstrated that to achieve the best performance for the hyperspectral imaging system, the light beam should be of circular shape and Gaussian type with the diameter of less than 1 mm, the optimal minimum source-detector distance should be about 1.5 mm, and the optimal maximum source-detector distance should be equivalent to 10-20 mfp' or determined by the minimum signal-to-noise ratio of 20 (or 150 CCD counts for the system). Under these optimal conditions, the hyperspectral imaging system achieved average accuracies of 23% for the absorption coefficient at 530-850 nm and 7% for the reduced scattering coefficient at 500-900 nm, when it was evaluated using the model samples against the transmittance and integrating sphere methods. These results are better than, or similar to, those obtained with other spatially-resolved sensing configurations, and frequency-domain and time-resolved techniques for single or several wavelengths. The system also had good reproducibility and sensitivity with the minimum detectable µa value of 0.0082 -1 cm . The research provided a systematic guide for optimizing and evaluating the optical system and offered solutions to improve accuracy in measuring the absorption and reduced scattering coefficients, which will be valuable for nondestructive quality evaluation of food and agricultural products. 100   CHAPTER 5 DEVELOPMENT OF AN OPTICAL PROPERTY MEASUREMENT PROTOTYPE 5.1 Introduction A general-purpose bench-top optical property measurement prototype, which is named ‘Optical Property Analyzer’ (or OPA, as shown in Figure 5.1), was designed and built for optical property measurement and hyperspectral image acquisition. The OPA was developed based on the concept of hyperspectral imaging-based spatially-resolved technology and the optimization studies on the inverse algorithm and optical designs presented in Chapters 3 and 4. Several important factors were considered in designing and assembling the OPA. First, it should be a stand-alone portable instrument incorporating the optimal optical designs, which is easy to use and relatively low in cost. Second, it should be able to measure the optical properties of turbid foods and biological materials for the visible and short-wave near-infrared region of 500-1,000 nm since rich information related to the quality attributes of food and agricultural products can be obtained from this wavelength window. Third, the prototype should be such designed that it can accommodate different sizes of samples whose effect on measurement of the optical absorption and scattering properties would be minimized or eliminated. Fourth, the OPA should also have full functions for acquiring and preprocessing hyperspectral images, which would expand the scope of use for the instrument and thus reduce overall costs for the user. Finally, it is important to develop integrated software for the OPA to control the hardware including camera, light source and linear stage, and also be able to perform real-time data processing and analysis for optical properties computation. Based on these considerations, the first version of the OPA was designed and built. 101   Figure 5.1 Optical Property Analyzer (or OPA) for measuring optical properties and acquiring hyperspectral images. 5.2 System Setup and Software Development 5.2.1 Optical Property Analyzer Setup The hardware of the OPA includes three main parts: imaging, illumination, and sample positioning as shown in Figure 5.2. The software for hardware control, imaging acquisition, and data processing and analysis was developed using Microsoft Visual C#. Details of the hardware and software in the OPA prototype are discussed in this section. 102   Figure 5.2 Schematic of the Optical Property Analyzer (or OPA) for measuring the optical properties of biological materials. The imaging unit consists of a high performance 14-bit electron-multiplying CCD (EMCCD) camera (Luca EM R604, ANDOR TM Technology, South Windsor, Connecticut, USA), an enhanced imaging spectrograph attached to the camera (ImSpector V10E, Spectral Imaging Ltd., Oulu, Finland) and a prime lens. Two different prime lenses are used alternatively to meet the different field of view requirements for two types of measurement; one is for optical properties measurement (Xenoplan 1.9/35, Schneider Optics, Hauppauge, NY, USA), and the other for hyperspectral image acquisition (Cinegon 1.4/12, Schneider Optics, Hauppauge, NY, USA). The EMCCD camera covers the wavelength range from 400 nm to 1,000 nm with the peak quantum efficiency of 65% at 600 nm. It can operate with the EM gain off for conventional CCD operation under strong illumination conditions, and with the EM gain on when the illumination intensity is low. It utilizes a monochrome megapixel frame transfer EMCCD sensor with 1,004×1,002 pixels of 8×8 µm, and is thermoelectrically cooled down to -20 °C. The sensor has 103   - an active pixel well depth of 30,000 electrons, a gain register pixel well depth of 80,000 e , a - - system readout noise of 18 e and less than 1 e for typical use and electron multiplication, respectively. The imaging spectrograph attached to the camera covers the wavelength range of 400-1,000 nm. The ImSpector Enhanced imaging spectrograph V10E is based on the principle of transmission diffraction grating in a prism-grating-prism configuration which provides high diffraction efficiency and good spectral linearity. The imaging spectrograph has a greater spatial image width, high resolution (small spot size), low smile and keystone distortions at subpixel level (smile <1.5 µm, keystone <1 µm), and high throughput with telecentric input. The imaging spectrograph has the slit width of 30 µm and slit length of 14.2 mm with the optical resolution of 2.8 nm. Based on the grating equation, different wavelengths can be diffracted to the same place on the detector surface. Therefore, an order-blocking filter (OBF) is installed between the imaging spectrograph and the camera close to the detector sensor to prevent the second-order below 500 nm from overlapping with the long wavelength range of 800-1,000 nm. Two illumination systems including a point light source and a line light source are used in the OPA for optical properties measurement and hyperspectral image acquisition, respectively. Only the point light source is shown in Figure 5.2 since this study was focused on optical properties measurement. A tungsten halogen light bulb with the output of 20 watts and the variation of 0.5% (HL-2000-HP, Ocean Optics, Dunedin, FL, USA) is connected to a DC regulated controller chips (PT6204N, 12V, Texas Instruments Inc., Dallas, Texas, USA), to provide point light, which covers a broad wavelength range of 360-2,000 nm with a cooling fan installed to keep the light source from overheating, thus maintaining its stability. A specially designed optical fiber coupled with a focusing lens is used to deliver a circular beam with the diameter of 1 mm at the focal point. The incident beam is arranged 1.5 mm away from the 104   scanning line with the angle of 15° with respect to the vertical axis, and it is parallel to the scanning line, which allows to maintain the constant source-scanning line distance, thus eliminating or minimizing errors in acquiring spatially-resolved reflectance that would be caused by the sample position’s variation. The sample positioning unit is composed of a motorized linear horizontal stage (Twintrac, TSZ8020, US23T22104-8LS, US Automation, Laguna Hills, CA, USA) with the maximum linear speed of 203 mm/sec, the maximum dynamic loading of 91 kg and a positioning accuracy of 0.0006 mm/mm, a manually adjustable vertical stage and a sample holder for positioning samples to the predetermined position. The stage motor (US23T22104-8LS, US Automation, Laguna Hills, CA, USA) is connected to the bi-polar stepper motor driver (SideStep, Probotix, Peoria, IL, USA). The light source controller and the motor controller are connected to a microcontroller board based on the ATmega168 (Arduino Duemilanove, Arduino, open-source electronics prototyping platform), which is, in turn, connected to a Windows-based computer with a USB cable. A specially developed software program in Arduino programming code is uploaded to the microcontroller to accept serial commands from the host application software installed in the computer to control stage movement and the light sources. The EMCCD camera is connected to the computer through another USB cable. The application software controls the camera operation by using Software Developer Kit (SDK) provided by the camera's manufacture. 5.2.2 Optical Property Analyzer Software Development The OPA software contains three components, including main program, optical property computation program, and microcontroller program, which were developed in three different environments, as shown in Table 5.1. The main program was developed using Microsoft Visual 105   Studio and written in C#.NET language. The graphical user interface uses Multiple Document Interface (MDI) style. It uses a third party external library to plot chart (PlotLab, a freeware downloaded from http://www.mitov.com/html/plotlab.html) and also an SDK library (Atmcd32cs.dll, library for C#) provided by Andor Corporation, the camera manufacturer. Table 5.1 Components of the Optical Property Analyzer control software. Component Main program Development Environment Microsoft Visual Studio (C# language) Function Multiple-Document Interface (MDI) container to call dialogs and functions; interact with other components to perform tasks (image acquisition, motor control, illumination control, and optical property computations) Optical property computation program MATLAB Perform optical property computation Microcontroller program Arduino (C language) Control linear stage; control light sources The OPA control software is designed for the system control (light source, camera and stage), image acquisition, and optical property computation. The main window of the optical property analyzer software is shown in Figure 5.3, which consists of File, Edit, View, Acquisition, Tool, Windows, and Help menu, and each menu contains submenu related to a particular task. File, Edit, and View menus are used to open, save, view and simply edit existing images (images with .bil extension). 106   Figure 5.3 Main window of the Optical Property Analyzer software. Figure 5.4 shows the window for setting parameters related to the camera, stage, illumination and image saving functions for the optical properties measurement. This display window is obtained by Control menu in Acquisition menu. Before images are acquired from samples, the camera dark current images are collected first by clicking ‘Background’, ‘Back Sub’ in Acquisition will be enabled, and the background image is then subtracted from the sample images. Auto Increment mode can be used in the experiment, which saves images in multiband format (ENVI format) with the given File Name and automatically adds sample number. 107   Figure 5.4 Display window of the Optical Property Analyzer for the image acquisition setup. A user-friendly interface for calculating the spectra of absorption and reduced scattering coefficients was also developed in the software, which includes two diffusion models and corresponding inverse algorithms, preselected parameters, and data pre-processing methods (Figure 5.5). After all parameters in the left column of ‘Optical Properties Computation Dialog’ are set up, image files, output directory and output file name are defined in the right column of 108   the dialog, and optical properties computations can be triggered with a processing message window by clicking ‘Compute OP’. In addition, the software has Calibration and Diameter menus as shown in Figure 5.6 for storing the spectral and spatial calibration information, and sample size information, which can be recalled during the optical properties computation. Figure 5.5 Display window of the Optical Property Analyzer for optical properties computation.     109   (a) (b) Figure 5.6 Display windows of the Optical Property Analyzer for (a) calibration information and (b) sample size information. 5.3 System Calibration and Evaluation To accurately measure spatially-resolved diffuse reflectance profiles from a sample, the Optical Property Analyzer was calibrated spectrally and spatially by following the procedure described in Lu and Chen (1998). Spectral calibrations were performed using a xenon lamp (model 6033, Newport, Irvine, CA, USA), a mercury-argon lamp (model 6035, Newport, Irvine, CA, USA) and a laser at 905 nm (model PM15/4586, Power Technology Inc., Little Rock, AR, USA). Spatial calibrations were achieved using a telecentricity target with a series of parallel lines (model 58404, Edmund Optics Inc., Barrington, NJ, USA). Finally, the OPA acquired twodimensional (2-D) spatially-resolved reflectance images with a spectral resolution of 3.1 nm/pixel in the vertical direction and a spatial resolution of 0.2 mm/pixel in the horizontal direction. 110   The performance of the OPA was also evaluated for accuracy, stability, precision/ reproducibility, and sensitivity by using the homogenous liquid model samples described in Section 4.2.2 and following the same procedures presented in Section 4.2.5 of Chapter 4. The average estimated error for all model samples was 24% for µa and 7% for µs'. Figure 5.7 shows the system reproducibility (or precision), and the coefficient of variation (or CV) in the absorption peak at 555 nm was less than 10% for µa, and less than 4% for µs'. It has to be mentioned that the absolute values of µa were very small, thus causing the relatively large error of estimating µa compared with that of µs'. The main error sources for estimating µa and µs' could have come from the light beam, source-detector distance, and inverse algorithm, according to the optimization studies of the inverse algorithm and optical designs described in Chapters 3 and 4. Additionally, the result of the sensitivity measurement on µa is presented in Figure 5.8. -1 The minimum detectable value of µa was 0.0117 cm . The sensitivity of µs', as determined by the CV values, was always less than 3% because µs' was much larger than µa for the investigated -1 range of 7.0 ≤ µs' ≤ 39.9 cm . These results showed that the Optical Property Analyzer has achieved acceptable accuracies for measuring the absorption and reduced scattering coefficients, which are either comparable to or superior to other reported studies using time-resolved, frequency-domain, or other types of spatially-resolved instruments based on the discussion in Section 4.3.4 of Chapter 4. 111   a Relative variations of  (%) (a) 15 10 5 0 -5 -10 1 2 3 day 4 5 4 5 Relative variations of   (%) s (b) 4 2 0 -2 -4 -6 1 2 3 day Figure 5.7 Reproducibility of the Optical Property Analyzer for measuring absorption and reduced scattering coefficients of a model sample at 555 nm at each measurement day with respect to the average value calculated over 5 days. 112   3 Coefficient of variation of a (% ) 10 2 10 1 10 0 10 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 -1 a (cm ) Figure 5.8 Coefficient of variation versus reordered ascending absorption coefficients of a model sample at different wavelengths. 113   CHAPTER 6 MATURITY/QUALITY ASSESSMENT FOR APPLE AND PEACH FRUIT USING OPTICAL PROPERTIES 6.1 Introduction Fruit quality represents a combination of the attributes and properties that give them value in terms of consumer’s satisfaction. Maturity determines the final quality of harvested fruit and how fruit should be handled and marketed after harvest. The maturity/quality of fresh fruit is determined by multiple indices including firmness, sugar (or soluble solids content), starch level, and titratable acidity, as well as external appearance characteristics such as skin and flesh color (Crisosto, 1994). Firmness directly influences the texture, shelf life and consumer acceptance. Soluble solids content (SSC) and acidity determine the flavor of the fruit, while changes in the skin/flesh color during fruit ripening are mainly related to the breakdown of chlorophylls and the increase of other pigments. Methods commonly used for maturity/quality measurement include Magness-Taylor (MT) firmness tester, Brix refractometry for SSC, and titratable methods for acidity. These methods are destructive, time consuming and prone to operational error, and only can test a small percentage of products. Therefore, nondestructive sensing techniques are needed for analyzing and monitoring fruit maturity/quality. Nondestructive assessment of maturity/quality of fresh products like apple is challenging because their properties and characteristics can vary greatly within each fruit and between cultivars and are also affected by environmental factors, postharvest storage condition, etc. A number of nondestructive techniques based on measurement of the optical, mechanical, and electrical properties have been developed to evaluate fruit maturity/quality. In this study, an advanced optical technique for measuring absorption and scattering properties was used for assessing the quality of apple and peach fruit. Absorption is mainly related to the chemical 114   compositions of the sample, while scattering is influenced by the structural properties. Therefore, separate quantification of absorption and scattering properties could provide a new means for assessing fruit internal quality. The objective of this research was therefore to determine the absorption and reduced scattering coefficient spectra of peach and apple fruit using hyperspectral imaging-based spatially-resolved (HISR) technique, and to assess the maturity/quality of the fruit with the measured optical properties. 6.2 Materials and Methods 6.2.1 Fruit Samples To determine fruit maturity/quality using the optical properties, peaches and apples harvested from the experimental orchard of Michigan State University’s Clarksville Horticultural Experiment Station in Clarksville, Michigan were used. Samples were selected such that they had large variations in the maturity/quality parameters to be measured. Five hundred ‘Redstar’ peaches were hand harvested at three different dates within one week in 2010. Tests for these fruit were conducted in the same harvest day. ‘Golden Delicious’ (GD) and ‘Delicious’ (RD) apples were harvested once a week during six consecutive weeks for the harvest season of 2010. Eighty apples for each cultivar were tested one day after each harvest, and the remaining apples were kept in refrigerated air storage at 0°C. Tests for the stored apples were begun one week after the last harvest for up to six weeks. A total of 1,039 GD and 1,040 RD apples were used in the experiment. Optical properties and maturity/quality indices were measured from the fruit samples using the OPA and standard destructive methods respectively. 115   All laboratory measurements for peaches and apples were performed after the samples had reached room temperature (~ 24°C). 6.2.2 Optical Properties Measurement The OPA was used to measure the optical properties of peaches and apples (Figure 5.1). Binning operations of 4×4 were carried out during the image acquisitions to reduce the exposure time of 500 ms for each scan, which resulted in hyperspectral scattering images of 251×250 pixels from the 1004×1002 EMCCD camera. Each sample was first moved to the pre-determined height via the vertical stage, and it then began to move horizontally by the horizontal stage in synchronization with image acquisitions by the imaging system. To improve the repeatability of measurements, 10 scans were taken from each peach sample for every 1 mm horizontal displacement increment for a total of 9 mm distance during the image acquisition at the predetermined speed. The motorized stage was operated in the ‘go-stop’ mode, so that the fruit sample remained stationary during the image acquisition. From the peach study, large variation among the scans for individual samples was observed. Therefore, for the apple experiment, 19 scans were acquired from each apple sample for every 0.5 mm horizontal displacement over a 9 mm distance in order to further reduce the variation. Figure 6.1 shows a typical hyperspectral reflectance image for a peach fruit sample. A raw reflectance spectrum obtained by taking a vertical line at the center spatial position is presented in Figure 6.1(b), which had an absorption valley around 675 nm due to the presence of chlorophyll-a in the fruit. Other absorption valleys were primarily caused by the spectral response of the optical system, which would be eliminated through the normalization and nonlinear curve-fitting process. Each horizontal line taken from the image represents one 116   spatially-resolved reflectance profile at a specific wavelength [Figure 6.1(c)], and hence the entire image, in effect, consisted of ninety-eight spatially-resolved reflectance profiles for the wavelengths of 515-1,000 nm with the interval of 5 nm. Similar images were obtained for the apple fruit; there were 101 reflectance profiles for GD apples extracted from the scattering images for 500-1,000 nm and 98 for RD apples for 515-1,000 nm. (a) (b) (c) (d) Figure 6.1 Hyperspectral reflectance image and optical property spectra of a peach sample: (a) 2D display of the original reflectance image, (b) a raw spectrum extracted for a specific scattering distance, (c) a spatially-resolved reflectance profile extracted for a selected wavelength, (d) preprocessed or averaged spatially-resolved reflectance profile at the selected wavelength, (e) the spectrum of absorption coefficient (µa), and (f) the spectrum of reduced scattering coefficient (µs'). 117   Figure 6.1 (cont’d) (e) (f) The acquired hyperspectral reflectance images were then analyzed to obtain µa and µs' spectra. Since the spatially-resolved reflectance profiles were symmetric to the light incident point, the two sides were averaged in the extraction of optical properties [Figure 6.1(d)]. Smoothing and normalization using the peak value for each profile were also applied to the averaged profiles to reduce the noise and avoid the need for absolute reflectance measurement. Each pre-processed spatially-resolved reflectance profile for every wavelength was then fitted by the Kienle diffusion model [equation (3.4)] using the least squares inverse algorithm presented in Chapter 3, from which the spectra of absorption and reduced scattering coefficients were obtained. The final µa and µs' spectra were determined by averaging the spectra of optical parameters from 10 or 19 images for each peach (or apple) sample, which are shown in Figure 6.1(e) and (f). 6.2.3 Acoustic Measurement Since it was the first time that the OPA was used to assess peach maturity/quality, acoustic firmness measurements were also performed after the optical measurement and before the 118   destructive measurements for each peach sample to compare with the measurements from the OPA. A commercial desktop acoustic firmness sensor (AFS, AWETA, Nootdorp, Netherlands) was used in this study. The acoustic firmness sensor detects the vibration pattern of the acoustic wave travelling across the fruit generated by gently tapping the fruit by a plastic probe, which is 2 2/3 related to the overall firmness of the fruit. The acoustic firmness index (FI = f m 6 /10 , where f is the first resonance frequency and m is the mass of the fruit) was obtained from each measurement. For each sample, three replicated measurements were performed at the same imaging area. Only average values were recorded for further analysis. 6.2.4 Destructive Maturity/Quality Measurement Peach maturity/quality measurements included firmness, SSC and skin and fresh color, whereas only firmness and SSC were measured for the apple fruit, by using standard destructive methods. The firmness of peaches and apples was measured from the same imaging area after removal of the fruit skin, using a Texture Analyzer (Model TA.XT2i, Stable Micro Systems, Inc., Surrey, UK) equipped with an 11-mm Magness-Taylor (MT) probe at a loading speed of 2 mm/s. Force/displacement curves were recorded for a penetration depth of 9 mm, and the following parameters were then determined from the curves: maximum force (FM in N), the area under the 2 force/displacement curve corresponding to maximum force (AM in mm ), and slope (force/deformation, SM in N/mm) measured between the point of initial contact and the point corresponding to maximum force. Soluble solids content (SSC in °Brix) was determined from the juice of the fruit samples using a hand-held digital refractometer (Model PR-101, ATAGO CO, Tokyo, Japan). A chroma meter (CR-400, Konica Minolta Sensing, Inc., Japan) was used 119   * * * * for color readings of the peach skin and flesh with the CIELAB color space (L , a , and b ). L * is a measure of lightness on a scale from 0 to 100, and can be easily explained. However, a and * b are merely coordinates and difficult to interpret separately (McGuire, 1992). Therefore, they * * * were converted to hue angle [h°, arc tan (b /a )] and chroma (C , a*2  b*2 ). Hue angle is a more appropriate measure of color. The color goes from red to bluish/green with the hue angle from 0° to 180°. Chroma is the degree of departure from gray to white towards the pure hue color or pure chromatic color (McGuire, 1992). 6.2.5 Data Analysis Statistical methods were applied for assessment of the peach and apple maturity/quality using the optical data. For the peach study, simple linear regression analyses were first performed on the relationship among the fruit maturity/quality parameters (firmness, SSC, and skin and flesh color). Prediction models for peach firmness, SSC and color parameters based on spectra of the absorption and reduced scattering coefficients and their three combinations (i.e., µa, 1/2 µs', µa&µs', µa×µs', µeff = [3µa(µa+µs')] , where µa&µs' refers to the combination of the two coefficient spectra in sequence, µa×µs' is the multiplication of the two coefficients, and µeff is the effective attenuation coefficient which defines the ability of light to penetrate the tissue) were developed using partial least squares (PLS). Simple regressions were performed between the MT firmness and acoustic/impact data. The 500 peach samples were first sorted in ascending order * * for firmness (FM, AM, and SM), SSC, and skin and flesh color (L , h°, and C ), respectively, and they were then divided into two groups with 75% for the calibration set and 25% for the 120   validation set. A four-fold cross-validation strategy by repeating the above subsampling procedure for four times was used to ensure all the samples in the dataset were eventually used for both calibration and validation. With this sampling technique and large data size for modeling, the effect caused by the sample selection on the calibration and validation results could be reduced or eliminated. The same modeling method and sampling technique were also used for the prediction of apple firmness and SSC. Partial least squares (PLS) regression is a bilinear modeling method where the original independent information is projected onto a small number of latent variables (LVs) to simplify the relationship between independent and dependent variables for prediction with the smallest number of LVs. A leave-one-out cross-validation method was applied on the calibration set to determine the optimal number of LVs based on the predicted residual error sum of squares (PRESS) value, which could avoid the overfitting problem. The calibration models were then evaluated for the validation set of samples using the statistical parameters, including standard error of prediction (SEP) calculated as the standard deviation of the prediction residuals, and correlation coefficient (r) between the predicted and the measured values of fruit maturity/quality parameters. 6.3 Results and Discussion 6.3.1 Statistical Data of Fruit Maturity/Quality The statistical data of peach firmness, SSC, and skin and flesh color for all tested samples are summarized in Table 6.1. The MT firmness distribution of the tested peaches, as measured by each of the three extracted parameters, was relatively uniform. While the SSC of the samples had a narrow range from 6.7 to 13.0 °Brix, with the mean of 10.0 °Brix and the standard deviation of 121   1.1 °Brix. Skin color variations in terms of lightness, hue angle, and chroma for the samples were greater than those for the fruit flesh. Table 6.1 Statistics of fruit maturity/quality parameters for 500 ‘Redstar’ peaches measured by standard methods. Maturity/quality parameter Magness- Maximum Taylor force, FM (N) firmness Area, AM (mm2) Slope, SM (N/mm) Mean 42.4 Standard deviation 26.2 Maximum Minimum 109.9 3.7 146.8 70.8 355.4 11.4 9.2 4.2 20.7 0.8 Sugar content Soluble solids content, SSC (°Brix) 10.0 1.1 13.0 6.7 Skin color Lightness, L* Hue angle, h° Chroma, C* 24.4 7.8 45.1 7.9 25.1 2.9 41.4 18.0 19.7 7.8 45.5 7.5 Lightness, L* Hue angle, h° Chroma, C* 49.4 5.5 57.4 16.2 24.4 2.2 32.5 17.0 55.2 6.2 66.3 16.6 Flesh color Table 6.2 shows the statistical data of fruit firmness and SSC for the two apple cultivars for three test groups (freshly harvested, after storage, and combined). The firmness of GD and RD apples decreased with storage time, whose mean values changed from 78.6 N and 79.7 N for the freshly harvested group to 48.4 N and 57.0 N for the after-storage group, respectively. The mean values of SSC, however, only changed slightly with storage time. 122   Table 6.2 Statistics of the firmness (i.e., maximum Magness-Taylor force or FM) and soluble solids content (SSC) of ‘Golden Delicious’ (GD) and ‘Delicious’ (RD) apples for the freshly harvested, after-storage, and combined groups. GD Mean Standard deviation Maximum Minimum Total samples Freshly harvested Firmness SSC (N) (°Brix) 78.6 12.9 8.3 1.1 105.1 16.9 50.2 10.1 479 RD Mean Standard deviation Maximum Minimum Total samples 79.7 13.4 114.2 32.2 480 6.3.2 11.4 1.3 15.2 8.6 After-storage Firmness SSC (N) (°Brix) 48.4 13.6 8.8 1.3 76.4 17.9 29.4 9.7 560 Combined Firmness (N) 62.3 17.3 105.1 29.4 1039 57.0 13.9 92.8 16.8 560 67.5 17.7 114.2 16.8 1040 13.3 1.0 17.0 9.7 SSC (°Brix) 13.3 1.3 17.9 9.7 12.4 1.5 17.0 8.6 Peach Maturity/Quality Evaluation 6.3.2.1 Absorption and Scattering Spectra of Peach Fruit Figure 6.2 shows the absorption and reduced scattering coefficients spectra of peach fruit with different skin colors (h° = 23.7, 21.9, and 20.4) and firmness levels (FM = 72, 47, and 21 N) for the spectral region of 515-1,000 nm. The µa spectra of the three peaches with different skin colors had two prominent peaks around 675 nm and 970 nm, which were caused by the absorption of chlorophyll-a and water in the fruit tissue, respectively [Figure 6.2(a)]. Values of the absorption coefficient for 720-850 nm were relatively small and consistent, compared with the rest of the spectral region. Since the color of peaches at different maturity stages varies greatly due to the changes of major pigments (i.e., chlorophyll and anthocyanin), the absorption spectrum for the dark red peach presented another peak around 525 nm due to the large amount of anthocyanin. For the peaches with light red and red colors, smaller absorption peaks from 123   chlorophyll-b were also observed around 620 nm. In addition, it was found that the absorption peak around 675 nm decreased with the increase of absorption at 525 nm, as the value of hue angle decreased with the fruit color that changed from light red to dark red (Figure 6.2). This is because during the fruit ripening, the chlorophyll content of peaches is expected to decrease, while anthocyanin is produced and then becomes the dominant pigment. Large changes in the pattern of the absorption spectra demonstrate that the absorption properties can be useful for evaluating the maturity/quality of peach fruit. (a) 0.5 Peach (light red, ho =23.7) 0.4 -1 0.3  a (cm ) Peach (red, ho =21.9) 0.2 Peach (dark red, ho =20.4) 0.1 0 550 600 650 900 950 1000 Peach (hard, 72N) Peach (medium, 47N) Peach (soft, 21N) 15 10 ) -1 850 (b) 20 s  (cm 700 750 800 Wavelength (nm) 5 550 600 650 700 750 800 Wavelength (nm) 850 900 950 1000 Figure 6.2 Spectra of (a) absorption coefficient for three peaches with different skin colors (light red, red, dark red), and of (b) reduced scattering coefficient for three peaches at different firmness levels (hard, medium, soft). 124   Compared to the absorption spectra, the reduced scattering coefficient spectra were relatively flat, with few features, for the spectral region of 515-1,000 nm [Figure 6.2(b)]. For most of the tested samples, their µs' values decreased steadily with the increasing wavelength. This pattern of changes is consistent with Mie scattering theory and other reported studies that scattering is wavelength-dependent (Keener et al., 2007; Michels et al., 2008). Considerable variations in the µs' spectra exist for different firmness levels. Hard peaches had consistently higher µs' values than the softer ones throughout the entire wavelength range. The softening of fruit is primarily due to the loss of firmness in texture caused by cell wall depolymerization and an increase in solubility of the middle lamella, leading to cell separations. Eventually, the scattering particles’ density in the fruit tissue would decrease during the softening process. Moreover, the particle size can also change due to the solubilization of pectic substances. These physical properties are related to the ability of the scattering particles in the fruit tissue to scatter -b light, according to the empirical equation (µs' = aλ ) for describing the relationship between µs' and wavelength λ, in which a is proportional to the density of the scattering particles, and b depends on the particle size (Mourant et al., 1997). For a few very soft peaches (with MT firmness below 16 N), an opposite trend was observed for the scattering spectra that the µs' value increased with the increasing wavelength. Dramatic drops of the µs' value for soft fruit could have invalidated the scattering-dominant assumption and, thus, the diffusion model, which could, in turn, have resulted in erroneous µs' and µa spectra. This issue warrants further investigation. The optical properties of fruits measured using different techniques have been reported; however the reported values of µa and µs' are quite different even for fruit of the same variety or 125   cultivar. When comparing the results from different studies, one should notice different experimental setups, fruit conditions (i.e., maturity level, postharvest storage, etc.), and other factors such as measurement techniques, each of which could affect the optical properties measurement. To roughly compare the results of this study with other studies, the distributions of µa and µs' values for all tested samples at four wavelengths (525 nm for anthocyanin, 620 nm for chlorophyll-b, 675 nm for chlorophyll-a, and 970 nm for water) are plotted in Figure 6.3. The µa -1 values for most of the samples at 525 nm, 620 nm, 675 nm, and 970 nm were 0.030 cm with -1 -1 -1 -1 the range of 0.002-0.550 cm , 0.090 cm with the range of 0.005-0.464 cm , 0.090 cm with -1 -1 the range of 0.0003-0.375 cm , and 0.450 cm -1 with the range of 0.108-1.144 cm , respectively. Large variations of µa values existed due to different maturity levels of the tested fruit. However, the µa values for most of the tested peach fruit at the four wavelengths obtained from this study are comparable with those of stone fruit from other reported studies. Cubeddu et -1 al. (2001a) showed µa values of one peach fruit ranging approximately from 0.020 to 0.450 cm -1 at 650-1,000 nm, with 0.085 cm -1 at 675 nm and 0.45 cm at 970 nm. Tijskens et al. (2007) reported µa values of several groups of nectarines at the commercial harvest stage for 670 nm -1 with the range from 0.05 to 0.35 cm . The variations of µs' values for all tested samples were -1 relatively consistent at the four wavelengths ranging approximately from 3 cm -1 -1 to 23 cm , and the µs' values for most of the samples were around 11 cm , which are different from the µs' of -1 ~20-24 cm at 700 nm for a peach fruit reported by Cubeddu et al. (2001a). 126   (a) Sample number 120 100 80 120 525 nm 620 nm 675 nm 970 nm Sample number 140 60 40 (b) 100 80 60 525 nm 620 nm 675 nm 970 nm 40 20 20 0 0 0 0.25 0.5 0.75 1 µa at four wavelengths (cm-1) 0 5 10 15 20 25 -1) µs' at four wavelengths (cm Figure 6.3 Distributions of (a) absorption (µa) and (b) reduced scattering coefficient (µs') for 500 ‘Redstar’ peaches at four wavelengths (525 nm for anthocyanin, 620 nm for chlorophyll-b, 675 nm for chlorophyll-a, and 970 nm for water). 6.3.2.2 Prediction of Peach Maturity/Quality The PLS prediction results of ‘Redstar’ peach firmness, SSC, and skin and flesh color using µa, µs' and their combinations (µa & µs', µa × µs' and µeff) are shown in Table 6.3. Individual µa and µs' spectra showed correlations with various maturity/quality parameters with values of the correlation coefficient (r) varying from 0.420 to 0.855 for µa and from 0.204 to 0.840 for µs'. In most cases, prediction results obtained from µa spectra were better than those from µs' spectra, which is similar to that for apples reported by Qin et al. (2009). However, the PLS prediction models developed for µs' are simpler with fewer factors and faster computation, compared to those for µa. The physiological process of fruit ripening is accompanied with simultaneous changes in the chemical compositions, which are related to µa, and the cellulosic structures, which influence µs'. Hence the combinations of µa and µs' (µa&µs', µa×µs', and µeff) could yield 127   better predictions of peach maturity/quality attributes. As shown in Table 6.3, using the combined data of µa and µs' generally improved prediction results with r = 0.724 (SEP = 18.13 N) for firmness parameter FM, r = 0.458 (SEP = 0.96 °Brix) for SSC, r = 0.893 (SEP = 3.54) for the * * skin color parameter L , and r = 0.722 (SEP = 3.32) for the flesh color L . It was observed that most of the improved results were obtained using the effective attenuation coefficient (µeff), which is reciprocal of light penetration depth, and describes how easily a medium can be penetrated by the light. Hence, it seems reasonable to have obtained better predictions for the peach maturity/quality parameters using the µeff values (except for AM and SSC). The PLS results for prediction of peach firmness, SSC, skin and flesh color demonstrated that the optical properties are useful for assessing multiple peach maturity/quality attributes simultaneously. 128   Table 6.3 Partial least squares (PLS) prediction results of ‘Redstar’ peach maturity/quality parameters using absorption coefficient (µa), reduced scattering coefficient (µs') and their three combinations (µa & µs', µa × µs', and µeff).* Maturity/ quality Fact. Firmness 10 FM r µa µs' µa & µs' SEP Fact. r SEP µa × µs' Fact. r SEP Fact. r µeff SEP Fact. r SEP 0.713 18.47 6 0.494 10 12 0.720 18.38 11 0.722 18.43 10 0.724 18.13 AM 10 0.547 59.90 5 0.383 9 11 0.582 58.51 8 0.570 59.07 9 0.553 59.92 SM 14 0.665 3.23 6 0.559 14 11 0.711 3.05 12 0.685 3.13 14 0.712 3.01 Soluble solids content SSC 15 0.420 0.99 5 0.315 13 12 0.458 0.96 14 0.451 0.96 13 0.419 0.99 Skin color * 13 L h° 16 * 14 C 0.855 0.717 0.847 4.07 2.82 4.19 6 6 6 0.837 0.795 0.840 13 13 13 15 16 0.884 0.765 3.65 2.58 13 0.883 3.66 14 16 13 0.881 0.738 0.882 3.70 2.70 3.69 13 13 13 0.893 0.778 0.886 3.54 2.50 3.63 Flesh color * 16 L h° 8 * 25 C 0.660 0.457 0.580 3.62 1.99 5.33 6 4 6 0.630 0.204 0.575 21 6 20 15 10 0.674 0.433 3.63 2.03 10 0.640 4.87 18 7 19 0.711 0.447 0.640 3.40 1.99 4.86 21 6 20 0.722 0.462 0.645 3.32 1.98 4.84 * See Table 6.1 for the explanation of individual maturity/quality parameters and SEP = standard error of prediction. 129    Further analysis showed that firmness predictions (r = 0.724 and SEP = 18.13 N) by µeff based on the PLS model were better than those obtained using the acoustic method with r = 0.639 (Figure 6.4). These results suggest that there is good potential of using HISR technique for peach firmness assessment. Further research for other peach cultivars is needed to determine the accuracy and reliability of the new technique for firmness measurement. The prediction results for the three flesh color parameters were poorer than those for the skin * * color (r = 0.722, 0.462, and 0.645 for L , h°, and C , with the corresponding SEP = 3.32, 1.98, and 4.84, respectively). The presence of fruit skin and smaller variations in flesh lightness, hue angle and chroma compared with the fruit skin for the test samples (Table 6.1) could be the possible reasons for lower correlations for flesh color prediction. For SSC, the prediction results were not as good as those for firmness and color. Previous studies for apple quality assessment using optical properties also showed poorer SSC predictions than firmness predictions (Qin et al., 2009). 100 r=0.724 SEP=18.13N 80 60 40 20 0 0 20 40 60 80 100 120 Actual MT firmness, FM (N) (b) Acoustic firmness, FI (Hz2 g2/3) Optical prediction, FM (N) (a) 30 y = 0.0941x + 9.8821 r = 0.639 25 20 15 10 5 0 0 20 40 60 80 100 120 Actual MT firmness, FM (N) Figure 6.4 (a) Prediction of Magness-Taylor (MT) firmness (FM) of ‘Redstar’ peaches using effective attenuation coefficient (µeff) with partial least squares (PLS) and (b) correlation between acoustic data and MT firmness (FM). 130    6.3.3 Apple Internal Quality Evaluation 6.3.3.1 Absorption and Scattering Spectra of Apple Fruit Typical absorption and reduced scattering coefficient spectra for the two apple cultivars for the spectral region of 500-1,000 nm for GD and 515-1,000 nm for RD are shown in Figure 6.5. The absorption spectra of apples showed a pattern similar to that of peach fruit in Figure 6.2. The µa spectra had two prominent peaks around 675 nm and 970 nm, which were caused by the absorption of chlorophyll-a and water in the fruit tissue, respectively [Figure 6.5(a)]. Since chlorophyll in the apple tissue decreased during the fruit ripening, the values of µa at 675 nm for -1 all samples ranged from 0.028 cm -1 -1 -1 to 0.841 cm for GD and from 0.041 cm to 0.762 cm for RD. Similar to those of peaches, the µs' spectra of apples were also relatively flat with fewer features [Figure 6.5(b)]. For most of the tested samples, their µs' values decreased steadily with -1 the increasing wavelength. The µs' values ranged between 7-20 cm -1 for GD and 4-26 cm for RD, which are at least one order in magnitude greater than the µa values. 20 (a) (b) GD RD -1 -1 a (cm 0.4 0.2 0 500 600 700 800 900 Wavelength (nm) GD RD 15 s  (cm ) 0.6 ) 0.8 10 1000 5 500 600 700 800 900 Wavelength (nm) 1000 Figure 6.5 Spectra of (a) absorption and (b) reduced scattering coefficients for four ‘Golden Delicious’ (GD) apples and four ‘Delicious’ (RD) apples. 131    6.3.3.2 Prediction of Apple Firmness and Soluble Solids Content Table 6.4 presents prediction results for the firmness of GD and RD apples using µa and µs', and their combinations (µa&µs', µa×µs', and µeff) for the freshly harvested, after-storage and combined groups. Likewise, the SSC prediction results using the optical properties for each group are summarized in Table 6.5. Both µa and µs' for each cultivar were correlated with the apple firmness and SSC. Firmness prediction results obtained from the µa spectra (r = 0.6870.885 for GD, and r = 0.744-0.844 for RD) were better than those from the µs' spectra (r = 0.6300.793 for GD, and r = 0.702-0.768 for RD) for each of the three sample groups. This was also true for the SSC evaluation. Fewer factors were used in the PLS prediction models with µs', because its spectra had fewer features. Using the three combinations of µa and µs' generally improved the prediction of apple firmness and SSC. With the best combination of µa and µs', the correlations for the firmness of RD and GD apples for the three test groups were 0.692-0.892 and 0.788-0.863, respectively, and they were 0.741-0.791 and 0.536-0.842 for SSC, respectively. Table 6.4 Firmness prediction results for ‘Golden Delicious’ and ‘Delicious’ apples for the freshly harvested, after-storage and combined groups.* Optical parameter Freshly harvested µa Golden Delicious Factors r 16 0.687 SEP 6.10 Delicious Factors r 20 0.785 SEP 8.43 0.630 6.49 7 0.764 8.68 µa&µs' 12 0.651 6.37 18 0.809 7.91 µa×µs' 17 0.692 6.06 19 0.813 7.83 µeff Afterstorage 6 16 0.687 6.10 22 0.822 7.67 µa 17 0.726 6.11 25 0.744 9.36 µs' 7 0.689 6.44 7 0.702 9.90 µs' 132    Table 6.4 (cont’d) µa&µs' 14 0.712 6.24 17 0.788 8.57 µa×µs' 17 0.734 6.02 20 0.765 9.00 µeff 18 0.730 6.07 29 0.762 9.10 µa 34 0.885 8.14 34 0.844 9.56 µs' 9 0.793 10.60 9 0.768 11.35 µa&µs' 38 0.881 8.29 19 0.852 9.31 µa×µs' 39 0.892 7.89 35 0.857 9.12 33 0.884 µeff * SEP = standard error of prediction. 8.16 38 0.863 8.94 Combined Table 6.5 Prediction results for the soluble solids content of ‘Golden Delicious’ and ‘Delicious’ apples for the freshly harvested, after-storage and combined groups.* Optical parameter Freshly harvest µa Golden Delicious Factors r 16 0.787 SEP 0.70 Delicious Factors r 17 0.823 SEP 0.76 6 0.489 0.99 7 0.784 0.84 12 0.781 0.70 17 0.842 0.73 17 0.791 0.69 20 0.816 0.78 µeff 16 0.777 0.72 22 0.821 0.77 µa 18 0.713 0.95 14 0.533 0.88 7 0.561 1.12 6 0.460 0.92 14 0.741 0.92 13 0.518 0.89 18 0.732 0.93 11 0.502 0.90 µeff 19 0.726 0.94 18 0.536 0.87 µa 20 0.760 0.84 23 0.812 0.88 8 0.544 1.09 8 0.750 0.99 16 0.768 0.83 19 0.825 0.85 23 0.778 0.82 21 0.804 0.89 0.773 0.83 27 0.805 0.89 µs' µa&µs' µa×µs' Afterstorage µs' µa&µs' µa×µs' Combined µs' µa&µs' µa×µs' 22 µeff *SEP = standard error of prediction. For the freshly harvested and after-storage groups, the correlations were relatively low because of smaller firmness variations for each group (Table 6.2). However, when the data from 133    these two groups were pooled, improved correlations for the firmness prediction of both cultivars were obtained, although the SEP values increased slightly due to larger firmness variations in the combined group. For SSC, the best predictions were achieved for the freshly harvested group with r = 0.791 (SEP = 0.69 °Brix) for GD and r = 0.842 (SEP = 0.73 °Brix) for RD. Comparable results were obtained for the combined group because the SSC in the apple fruit did not change significantly during storage. Figure 6.6 shows the firmness and SSC predictions for GD and RD apples obtained with the best combination of µa and µs' for the combined group. Better predictions of firmness for GD and RD apples with the correlation of 0.892 and 0.863, respectively, were obtained than for SSC predictions with r = 0.778 and 0.825. These results are better than, or comparable with, other reported studies using hyperspectral scattering technique 100 R = 0.892 SEP = 7.89 N 80 (a) Predicted SSC (°Brix) Predicted firmness (N) (Mendoza et al., 2010; Qin et al., 2009). GD 60 40 16 R = 0.778 SEP = 0.82 °Brix 14 (b) GD 12 10 8 20 20 40 60 80 Actual firmness (N) 100 8 10 12 14 Actual SSC (°Brix) 16 Figure 6.6 Prediction of fruit firmness (a, c) and soluble solids content or SSC (b, d) using the best combinations of µa and µs' for the combined group of ‘Golden Delicious’ (GD) and ‘Delicious’ (RD) apples.   134    100 80 R = 0.863 SEP = 8.94 N (c) Predicted SSC (°Brix) Predicted firmness (N) Figure 6.6 (cont’d) RD 60 40 20 20 40 60 80 Actual firmness (N) 100 16 R = 0.825 SEP = 0.85 °Brix 14 (d) RD 12 10 8 8 10 12 14 Actual SSC (°Brix) 16 6.4 Conclusions The spectra of absorption and reduced scattering coefficients for ‘Redstar’ peaches for the spectral region of 515-1,000 nm, ‘Golden Delicious’ apples for 500-1,000 nm, and ‘Delicious’ apples for 515-1,000 nm were measured using the Optical Property Analyzer. Significant changes in the absorption and reduced scattering spectra were observed for peaches with different levels of maturity. The absorption spectra were shaped by major pigments (i.e., chlorophyll and anthocyanin) and water in the fruit, while the reduced scattering spectra decreased in value during the fruit softening. Similar absorption and scattering features were also found for the apple samples. Both µa and µs' were correlated with the peach firmness, SSC, and skin and flesh color and with the apple firmness and SSC; however, better correlations were obtained using the absorption spectra than using the scattering spectra. The combinations of µa and µs' (i.e., µa&µs', µa×µs', and µeff) based on PLS models gave better predictions of firmness (FM), SSC, and skin * and flesh color parameter L for peach fruit, with values of the correlation coefficient (r) being 135    0.724 (SEP = 18.13 N), 0.458 (SEP = 0.96 °Brix), 0.893 (SEP = 3.54), and 0.722 (SEP = 3.32), respectively. For the apple study, the best firmness prediction was obtained with r = 0.892 (SEP = 7.89 N) for GD and r = 0.863 (SEP = 8.94 N) for RD. For SSC evaluation, better results were achieved for the freshly harvested group with r = 0.791 (SEP = 0.69 °Brix) for GD and r = 0.842 (SEP = 0.73 °Brix) for RD. These two studies demonstrated that the hyperspectral imaging-based spatially-resolved technique is useful for measuring the optical properties of fruit, and it provides a new means of assessing fruit maturity/quality.   136    CHAPTER 7 RELATIONSHIP BETWEEN THE OPTICAL AND MECHANICAL/STRUCTURAL PROPERTIES OF APPLE TISSUE 7.1 Introduction The texture of fresh fruit depends on variety, maturity, climate condition during growth, and postharvest storage condition, and it is largely determined by the structural and mechanical properties of fruit tissue. The softening and its resultant change in the texture of apples during the ripening mainly results from cell wall depolymerization and the dissolving of the middle lamella (Wood et al., 2009). Accompanied with the tissue softening and the weakened adhesion between cells are changes in the cell wall strength, cell turgor and the number, size, and shape of intercellular spaces, which, in turn, affect the fruit firmness or rigidity. It is thus critical to understand how the structural and/or mechanical properties of cells and their changes affect the textural properties and ultimately the quality of fruit in storage. Considerable studies have been reported on the relationship between fruit quality and physiological changes. Both destructive and nondestructive techniques that are based on quasistatic, impact, vibration, and acoustic principles have been used for measuring the mechanical and/or textural properties (e.g., firmness and strength or elasticity) of fruit. Electron microscopy was also utilized to evaluate the structural changes of fruit during ripening and storage, and to examine or quantify the physiological activities in the fruit. While it is well known that the optical absorption and scattering properties are related to the composition and cellular structure of fruit tissue, no research has been attempted to elucidate and quantify their interrelationships. Better understanding of these relationships can help us develop better techniques for fruit quality inspection, monitoring and control. 137    The research reported in this chapter was therefore aimed at gaining a fundamental, quantitative understanding of the relationships between the optical and mechanical/structural properties of apple fruit. The specific objectives of the research were to: 1) measure the absorption and reduced scattering coefficient spectra and the mechanical and structural properties of apple tissues for ‘Golden Delicious’ and ‘Granny Smith’ cultivars over time during high temperature/high humidity storage, and 2) relate the optical measurements to the elastic modulus (Young’s modulus) of flesh tissues and cell structures, including cell size/shape and morphological characteristics. 7.2 Materials and Methods 7.2.1 Apple Samples Two apple cultivars, ‘Golden Delicious’ (GD) and ‘Granny Smith’ (GS), were chosen in this study. The apples were harvested at the commercial maturity stage in 2010 from the experimental orchards of Michigan State University’s Horticultural Teaching and Research Center in Holt, MI and Clarksville Horticultural Experiment Station in Clarksville, MI. GD apples were harvested three weeks earlier, and they were kept in a controlled atmosphere (CA) environment (2% O2 and 3% CO2 at 0 °C), until after GS apples were received. The test apples were visually inspected to ensure they were absent of damage or blemishes, and they were of relatively uniform in size with the equatorial diameter of 73-78 mm for GD and 77-85 mm for GS. The apples were randomly grouped into five groups of 10 fruit each for each cultivar. They were kept at ~ 22 °C and ~ 95% relative humidity for 1, 8, 15, 22, and 30 days, respectively, to accelerate the changes in fruit quality and condition. 138    For each storage time, one group of apples was taken out from the storage. Measurements were first performed of the optical properties of all apple fruit. Thereafter, compression tests (see the following subsection for details) were conducted for each apple sample. Finally, microscopic image analyses of fruit tissue specimens were performed. Because of the extensive time needed, and the difficulty, in preparing tissue specimens for micro-structural analysis, only three specimens for each group of 10 apple samples were prepared for confocal laser scanning microscopy (CLSM) and two specimens for scanning electron microscopy (SEM). To monitor the apple firmness and weight change with storage time, another subset of 10 intact fruit samples for each cultivar were kept in the same storage condition for up to 30 days. Acoustic firmness and optical properties measurements were taken for these samples at each of the five test dates. 7.2.2 Optical Properties Measurement The Optical Property Analyzer (OPA) (see Chapter 5) was used to measure the spectral absorption and reduced scattering coefficients of GD and GS apples for the spectral region of 500-1,000 nm, following the procedure described in Section 6.2.2 of Chapter 6. Mean spectra of absorption and reduced scattering coefficients were calculated for the 10 fruit of each group for each storage time. 7.2.3 Acoustic Firmness Measurement and Compression Test for Measuring Mechanical Properties The commercial desktop firmness sensing system (AFS, AWETA, Nootdorp, Netherlands), described in Chapter 6, was used to monitor the apple firmness and weight change during storage. The fruit was gently tapped by a plastic rod in the instrument. Four parameters, including the 139    2 2/3 mass of the fruit (m), first resonance frequency (f), the acoustic firmness index (FI=f m 6 /10 ), and impact firmness (IF), were obtained from each measurement. For each sample, three replicate measurements were performed at the same equatorial location where the optical property measurements had been performed earlier. Only the average values were recorded for further analysis. Compression tests were conducted following the optical properties and acoustic measurements. The test apples were first cut into two halves parallel to the longitudinal axis as shown in Figure 7.1(a). Cylindrical specimens of 13.9 mm in diameter and 13 mm in length were then taken from the fruit flesh of one apple half 5 mm beneath the skin, using a cork borer and a double-blade knife [Figure 7.1(b) and (c)]. The specimen size was determined based on the literature (Abbott and Lu, 1996) and preliminary study of the effect of specimen size on the mechanical properties measurement. Each flesh specimen was compressed between two parallel plates using a Texture Analyzer (Model TA.XT2i, Stable Micro Systems, Inc., Surrey, UK) at a loading speed of 0.42 mm/s [Figure 7.1(d)]. Force/displacement curves were recorded for up to 4 mm of displacement; and failure stress (σ), failure strain (ε), and Young’s modulus (or apparent modulus of elasticity, E) were then calculated from these curves. (a) (b) 13.9 mm 13 mm 5 mm Figure 7.1 Experimental procedures for compression test of apple tissue specimens: (a) cylindrical specimen taken from the fruit flesh of the apple; (b) cork borer and double-blade knife used for cutting tissue specimens; (c) final tissue specimen with 13.9 × 13 mm (D × L); and (d) compression test setup with the Texture Analyzer. 140    Figure 7.1 (cont’d) (d) (c) 7.2.4 Confocal Laser Scanning Microscopy and Scanning Electron Microscopy for Microstructural Analysis Confocal laser scanning microscopy (CLSM, Olympus FluoView 1000 LSM, Olympus America Inc., Center Valley, PA, USA) and scanning electron microscopy (SEM, JEOL 6400V, Japan Electron Optics Laboratories, Tokyo, Japan) were conducted on the apple tissue specimens to quantify their microstructural changes for different storage times ranging from 1 to 30 days. CLSM enables us to view a single layer of cells from the tissue and reduces noise from adjacent layers of cells, and it was used for quantitative analysis of apple tissues for cell information (e.g. cell size and shape parameters). SEM allows surface examination to obtain cellular characteristics at higher resolutions. For CLSM, tissue specimens from the outer cortex of apples about 10 mm away from the skin were cut by a cork borer with a diameter of 13.9 mm, using the same sampling procedure as shown in Figure 7.1(a) for the compression test except for the sampling location or depth. A razor blade was then used to obtain cylindrical specimens of 1 mm thickness. Thereafter, the specimens were soaked with 0.1% Cango Red (Sigma-Aldrich Co. LLC, St. Louis, MO) for 10 min and then were washed by water several times prior to CLSM imaging. 141    For SEM, tissue specimens with the same dimensions as those for CLSM were cut from the fruit. The specimen preparation procedure was according to that described by Kim et al. (1990). The specimens were first fixed for 2 h in 2% glutaraldehyde/0.1M phosphate buffer at 0°C. The fixed specimens were then washed with 0.1M phosphate buffer for at least 30 min at 0°C, and post-fixed in 1% phosphate buffered osmium tetroxide (OsO4) solution for 1 h at 0°C. After another wash with phosphate buffer, they were dehydrated with a series of ethanol concentrations for various time periods as follows: 40% (5 min), 50% (5 min), 60% (5 min), 70% (15 min), 80% (10 min), 90% (15 min), and 100% (three times for 10 min each). The fixed and dehydrated specimens were quench-cooled in nitrogen slush under vacuum within the preparation chamber of the SEM (< 210°C). Each specimen was fractured using a movable blade, sublimated for 15 min at -90°C, and coated by gold/palladium. The specimens were then viewed by the SEM with 10 kV accelerating voltage at < 130°C. In order to obtain more comprehensive information on the cellular structure change during storage, SEM images were acquired at two different resolutions (magnifications of 50× and 1500×). 7.2.5 Data Analysis To quantify the optical property changes of apple tissue with storage time and their relationship with the mechanical/structural properties, image analysis and statistical analysis were carried out using Matlab programs. Analysis of variance with the storage time was first applied to the mean values of optical properties, acoustic firmness data, and compression test data for each storage time. Young’s modulus calculated as the gradient of the initial linear slope between 5.8% and 9.8% strain from the stress/strain curves obtained by the compression test was 142    used for characterizing the mechanical properties of fruit tissues and for understanding the structural changes in apples during storage. Cell morphological parameters including area, equivalent diameter, major and minor length, 2 perimeter, elongation (major length/minor length) and roundness (4 × π × area/perimeter ) were extracted from the CLSM images of apple tissues. Six images acquired from three specimens (two images for two locations from each specimen) for each cultivar were quantified using image segmentations to obtain average cell parameters for each storage time. Analysis of variance (ANOVA) and least significant difference (LSD) test were conducted to evaluate the cell changes among the storage days. The procedures of image processing are shown in Figure 7.2. First, the original images were preprocessed by converting the RGB images into binary images. Second, morphological operations on the binary images were implemented using the function ‘bwmorph’ in Matlab Imaging Processing Toolbox by bridging unconnected pixels, removing spur pixels, and performing dilation using the structuring element ‘ones (3)’. Third, the watershed transformation was performed on the images after the morphological operations, and morphological structuring elements were created. Finally, the segmented cells were labeled, and their morphological parameters were calculated using Matlab function ‘regionprops’. Original Preprocessed Clean-bridge Figure 7.2 Procedures of image processing for extracting cell parameters. 143    Figure 7.2 (cont’d) Clean-spur Clean-dilate Watershed Dilate-disk Overlay Label The SEM images of apple tissues acquired at the low resolution were analyzed to further quantify their microstructure using histogram probability and fractal analysis methods. The histogram of an image is a graphic description of the distribution of the gray levels, which are the gray level values versus the number of pixels at that value. It provides information related to the nature and characteristics of the image. Statistical-based histogram features including energy and entropy were extracted from the SEM images, where the histogram is used as a model of the probability distribution of the gray levels. The energy describes the distribution of the gray levels of an image and is defined as (Umbaugh, 2005). energy  n 1 2   p( g ) (7.1) g 0 where p( g )  N (g) is the first-order histogram probability, in which M is the pixel number in M the image, N(g) is the pixel number at gray level g, and n is the total number of gray levels. The entropy tells us how many bits we need to code the image data, and is given by 144    entropy  L 1  g 0 p( g ) log 2  p( g ) (7.2) The value of entropy increases when the pixel values in the image are distributed among more number of different gray levels. A complex image has higher entropy than a simple image, and the energy tends to vary inversely. Fractal analysis was used for description and quantification of surface topographies, providing texture information of the SEM image. Three parameters were extracted from the fractal analysis, which are fractal dimension at low scale (D) related to the surface roughness, the minimum elemental cell (dmin) related to the cell size, and the periodical region (dper) at high scale related to the distance between cells. The procedure to compute the fractal dimension (D) and parameters dper and dmin are based on the variogram of the intensity pixels in grayscale images and the method was proposed and successfully used in scanning electron microscopy images by Bonetto and Ladaga (1998), Bianchi and Bonetto (2001) and Quevedo et al. (2002). The variogram method has been detailed by Bianchi and Bonetto (2001) as follows. It calculates the fractal dimension from the V variance of the grey level distribution of the image corresponding to the studied surface sample V  z  s 2H (7.3) where the angle brackets <> denote the expectation value, Δzi is the difference of the gray level between two different positions in the image for a step i of length s measured in pixels. H value can be obtained by means of the graph slope in logarithmic scale from V versus s, which is related to the fractal dimension D by D  3 H 145    (7.4) The variogram presents a fractal behavior at a certain scale and simultaneously it produces minimums at other scales. Such minimums result from the presence of a periodical structure in the gray levels. This periodic region is described as a spectrum of space frequencies represented by “k wave vectors”. kxi and kyi can be obtained by projecting the “ki wave vector” on the pair of orthogonal coordinate axes, and therefore 2 2  k 2  k x    k y  (7.5) The periodic region parameter is then defined as d per  2 2 2  kx    k y  (7.6) The parameter dper is a measurement of the average diameter between the nearest elementary cells of enough statistic weight so as to produce periods. The parameter dmin is the intersection between the fractal and periodic zones that explains the highest value of the smallest cell size of enough statistic weight. These three parameters (D, dmin, and dper) extracted from fractal analysis have provided good characterization of image texture (Bonetto and Ladaga, 1998). Finally, general linear regression analysis was performed on correlating acoustic firmness, Young’s modulus and cell structural parameters to the optical properties of apple fruit. 7.3 Results and Discussion 7.3.1 Changes in the Optical Properties of Apples during Storage Figure 7.3 shows the mean absorption (µa) and reduced scattering coefficients (µs') spectra over 500-1,000 nm for 10 apples each of GD and GS for each test date. While the general 146    patterns of µa and µs' spectra for GD and GS apples were similar, their absolute values were different. (a1) µa (cm-1) 0.6 0.4 15 Day 1 Day 8 Day 15 Day 22 Day 30 0.2 13 12 11 10 9 0 0.8 (b1) 0.6 0.4 500 600 700 800 900 1000 Wavelength (nm) Day 1 Day 8 Day 15 Day 22 Day 30 µs' (cm-1) 500 µa (cm-1) (a2) 14 µs' (cm-1) 0.8 0.2 0 500 600 700 800 900 1000 Wavelength (nm) 15 14 13 12 11 10 9 600 700 800 900 1000 Wavelength (nm) (b2) 500 600 700 800 900 1000 Wavelength (nm) Figure 7.3 Mean spectra of 10 apple samples each for the five times of storage: (a1) absorption and (a2) reduced scattering coefficients of ‘Golden Delicious’ (GD) apples; and (b1) absorption and (b2) reduced scattering coefficients of ‘Granny Smith’ (GS) apples. There were two prominent peaks around 675 nm and 970 nm for the µa spectra of both GD and GS apples, which were caused by the absorption of chlorophyll-a and water in the fruit tissue, respectively. There were also small chlorophyll-b absorption peaks around 620 nm and water absorption peaks around 750 nm. Other absorption peaks were also observed in the µa spectra of GD and GS apples around 525 nm due to the presence of anthocyanin in the apple 147    skin. Moreover, the absorption peak around 675 nm for GD apples decreased consistently with -1 storage time; the average value of µa changed from 0.324 cm -1 to 0.103 cm over 30 days of storage, and dramatic decreases of µa values took place after eight days of storage. On the other hand, absorption around 525 nm increased slightly over time, with the mean values of µa -1 changing from 0.179 cm -1 to 0.207 cm during the 30 days of storage. This is because during storage, the chlorophyll content of GD apples was expected to decrease as the apple skin color changed from green to yellow and more anthocyanin was produced, which then became the dominant pigment. Similar changes also happened to the absorption spectra for GS apples. -1 However, the decrease of µa values around 675 nm from 0.248 cm -1 to 0.235 cm for GS apples was much smaller than that for GD apples, whereas the increase of µa values around 525 nm -1 -1 from 0.150 cm to 0.306 cm was much greater than that for GD apples. Overall, the change in absorption for the spectral region of 750-1,000 nm was much smaller during the entire storage period, compared to that for the visible spectral region of 500-750 nm. Compared to the absorption spectra, the reduced scattering coefficient spectra exhibited a different pattern for the spectral region of 500-1,000 nm [Figure 7.3(a2) and (b2)]. The µs' spectra were relatively flat, with few features. For most of the tested fruit samples of both cultivars, their µs' values decreased steadily with the increasing wavelength throughout the spectral region of 500-1,000 nm. Like the absorption coefficient, the reduced scattering coefficient also exhibited a consistent pattern of change during storage and its values for the fresh apples (i.e., day 1) were considerably higher than those for the stored apples for all other 148    test dates, except for GS for the region of 900 nm or higher due to the large variation of GS apples with small change of their µs' values. The maximum value for the mean reduced scattering coefficient for the entire wavelength range decreased with storage time consistently -1 from 13.4 cm -1 -1 to 11.4 cm for GD and from 12.4 cm -1 to 10.0 cm for GS over 30 days of storage. 7.3.2 Changes in Acoustic Firmness and Mechanical Properties during Storage Acoustic and impact firmness (e.g., FI and IF) for both cultivars decreased steadily with storage time [Figure 7.4(a)]. FI and IF decreased by 42.3% and 21.2% for GD, and 16.9% and 6.7% for GS, respectively, over the storage period of 30 days. GS apples had a lower rate of loss in firmness during the storage. The weight change of apples due to water loss can affect the mechanical and microstructural properties of the fruit. Figure 7.4(b) shows that there were 3.1% and 1.9% decreases in the average weight of GD and GS apples, respectively, but the decreases were not statistically significant. 149    45 105 90 35 75 25 60 15 45 5 350 Weight (g) GD_FI GS_FI GD_IF GS_IF (a) Impact firmness, IF Acoustic firmness, FI (Hz2 kg2/3) 55 GD GS (b) 300 250 30 200 150 100 50 1 8 15 22 30 15 22 30 Storage days Storage days Figure 7.4 Changes in a) acoustic firmness (FI) and impact firmness (IF), and (b) fruit weight, for the same 10 fruit each of ‘Golden Delicious’ (GD) and ‘Granny Smith’ (GS) measured at five storage times (the vertical bars denote two standard deviations). 1 8 The stress/strain curves for a GD apple and a GS apple measured for each of the five storage times are shown in Figure 7.5. For both cultivars, the stress/strain curve at day 1 had a sharp peak at the yield point with a sudden drop in the failure stress because the flesh tissue was still crispy or crunchy. As the storage time increased, the stress/strain curve became smoother, and it did not show a prominent yield point at days 8-30 for GD and at days 15-30 for GS. Day 1 Day 8 Day 15 Day 22 Day 30 (a) 0.5 0.4 0.6 Stress, σ (MPa) Stress, σ (MPa) 0.6 0.3 0.2 0.1 (b) 0.5 0.4 0.3 0.2 0.1 0 0 5 10 15 20 Strain, ε (%) 25 0 30 0 5 10 15 20 25 Strain, ε (%) 30 Figure 7.5 Stress/strain curves obtained from the compression test for the flesh specimens of (a) five ‘Golden Delicious’ apples and (b) five ‘Granny Smith’ apples, one for each storage time. 150    Figure 7.6 shows the mean values and statistical differences of Young’s modulus for GD and GS apples at five storage times. GS apples had significantly higher values of Young’s modulus compared with GD apples, which was in agreement with the relative values of acoustic firmness measurements for the two cultivars (Figure 7.4). The mean value of Young’s modulus decreased by 45.6% for GD and 34.2% for GS within 30 days of storage. This indicated that the rigidity of the cellular structures had decreased significantly over 30 days of storage time, and thus the fruit became soft at the end of storage. A dramatic decrease (significant at the 0.05 level) in Young’s modulus for GD apples occurred between day 1 and day 8, while a significant decrease in Young’s modulus for GS apples was only observed for day 22 and 30. 6 a a GD GS Young's modulus (MPa) a 5 4 b A b 3 BC B C 2 D 1 Day1 Day8 Day15 Stroage days Day22 Day30   Figure 7.6 Mean values of Young’s modulus for 10 ‘Golden Delicious’ (GD) (columns with capital letters) and ‘Granny Smith’ (columns with small-cap letters) apples for each of the five storage times (columns with different letters for the same cultivar are significantly different with p < 0.05). 151    7.3.3 Microstructural Changes in Apple Tissue during Storage ‘Golden Delicious’ apples are crispy, whereas GS fruit are firm, crunchy and juicy. These textural differences are related to the different cellular structures of the two cultivars. The CLSM images showed that the cells of GS tissues [Figure 7.7(b1)] were more rounded compared with those of GD tissues [Figure 7.7(a1)]. Closer contact between the cells of GS tissues was observed, while there was more space between the cells of GD tissues. Weaknesses at the cell conjunction may indicate less firmness or mealiness of the texture. During the storage, the cellto-cell adhesion was reduced and cell separation along the middle lamella increased. More broken cells were observed for both cultivars at day 30, as shown in Figure 7.7 (a2) and (b2). (a1) (a2) 100 µm  Figure 7.7 Confocal laser scanning microscopic images of a cross-section of apple flesh tissue 10 mm beneath the skin for ‘Golden Delicious’ apple at day 1 (a1) and day 30 (a2), and for ‘Granny Smith’ apple at day 1 (b1) and day 30 (b2) (‘→’ denotes the connection between cells, and ‘×’ denotes the broken cells).     152    Figure 7.7 (cont’d) (b2) (b1) 100 µm  Table 7.1 summarizes the cell morphological parameters for GD and GS apple tissue specimens, extracted from the CLSM images. The average cell area of GD apples (19.0-15.7 3 2 3 2 ×10 µm for day 1-30) was smaller than the average cell area of GS (28.3-22.2 ×10 µm for day 1-30). The same trend was observed for the equivalent diameter, major and minor length, and perimeter. The cell size of both cultivars decreased with storage time from day 1 to day 30. While significant differences in each of the cell size parameters for the two cultivars were found between day 1 and day 30, differences between two adjacent test days were not always significant. For the shape-related parameters (i.e., elongation and roundness), no consistent trend with the storage time was found for the two cultivars. During the storage, the breakdown of pectin in the weak intercellular lamella and the water redistribution in the cells and cellular space would have caused changes in the cell size and shape. But the cell shape change was difficult to quantify based on the two-dimensional CLSM images, since the cell deformation could have happened in all directions. 153    Table 7.1 Mean and standard error of the cell size/shape parameters for ‘Golden Delicious’ and ‘Granny Smith’ apple tissues for different storage days.* Storage day Cell size parameters Equiv. diameter Major length Minor length Perimeter (µm) (µm) (µm) (µm) Area (µm2) Golden Delicious a Day 1 18990±698 150.2±2.9 Day 8 17677±658 143.7±2.9 Day 15 16939±550 Day 22 16478±613 ac bc bc 141.9±2.5 139.2±2.7 a ac bc bc b b Day 30 15670±527 135.5±2.6 Granny Smith a a Day 1 28332±1120 184.3±3.8 b b Day 8 24292±1098 167.8±4.1 Day 15 21648±930 Day 22 23645±991 b b 158.2±3.6 166.0±3.8 b b 181.7±3.5 175.6±3.6 173.4±3.2 174.4±3.6 161.6±3.0 223.0±4.7 200.9±5.0 192.7±4.6 197.6±4.5 a a a a b a b b b 129.8±2.8 123.3±2.6 121.0±2.3 116.5±2.5 117.8±2.4 157.6±3.6 144.9±3.8 134.7±3.2 144.8±3.7 a ac bc bc bc a b c bc Cell shape parameters RSE of Elongation Roundness area (%) 545.2±10.4 524.3±10.5 515.2±9.2 1.60 ac 1.426±0.032 bce 1.527±0.044 de 705.0±15.5 603.7±14.7 570.0±13.0 595.1±13.5 a b b b 1.438±0.037 1.486±0.035 bc 510.9±10.1 485.1±9.1 a 1.420±0.028 1.92 1.408±0.033 1.396±0.023 1.470±0.036 1.429±0.054 ab ab ab a b a a a a 0.766±0.0061 0.754±0.0060 0.763±0.0053 0.757±0.0060 0.773±0.0054 0.765±0.0064 0.774±0.0059 0.777±0.0062 0.778±0.0087 ab ab a b ab a a a a b b b bc b b a Day 30 22181±910 161.9±3.6 188.0±4.2 143.4±3.4 576.0±12.7 1.346±0.020 0.789±0.0040 Values in the same column with different letters are significantly different at the 95% confidence interval (p<0.05); Standard error = standard deviation / cell number ; RSE (relative standard error) = standard error/mean  100%. 154  To further quantify the changes of tissue cells during the storage, the cell connections, denoted by the red lines in Figure 7.8(a1), were manually counted and measured for each scanning electron microscopic (SEM) image acquired at the low resolution for day 1 and day 30. The average connection length and standard error (standard deviation divided by the number of the cell connection) had decreased approximately from 87.80 ± 1.02 µm to 62.39 ± 0.52 µm for GD apples, and from 99.03 ± 0.71 µm to 78.29 ± 1.11 µm for GS apples for 30 days of storage. The high-resolution SEM images of GD and GS apples for 1 and 30 days of storage have presented significant differences in the cellular structures, as shown in Figure 7.8(b1)-(b2) and Figure 7.9(b1)-(b2). For both cultivars, cell separation in the middle lamella was observed after 30 days of storage, which was associated with pectin solubilization and the subsequent reduction in the connection between cells. This was also consistent with the changes observed from the CLSM images. At day 1, the cells of GD had a smoother external surface with some exudates [Figure 7.8(b1)], which indicated that the cell structure degradation might have begun. For day 30, more exudates were observed on the surface of cell walls [Figure 7.8(b2)]. Similar changes were also found in the cells of GS specimens. The cellular materials had appeared to be pulled away in the cells of GS even at day 1 [Figure 7.9(b1)]. Although attention was paid to the initial selection of fruit with the same level of maturity for both cultivars, it was still difficult to achieve this by visual inspection of the external appearance of each fruit. Therefore, image textural analysis was performed on the low-resolution SEM images to obtain quantitative information for the cellular structures. 155  a1 a2 500 µm 500 µm (b1) (b2) 20 µm 20 µm Figure 7.8 Scanning electron microscopy images of one cross-section of apple flesh tissue 10 mm beneath the skin for ‘Golden Delicious’ at day 1 (a1) and day 30 (a2) acquired at low resolution for an overview of the cell structures, and at day 1 (b1) and day 30 (b2) acquired at high resolution for a more detailed view of the cell-to-cell connection [the red lines in (a1) denotes five cell connections counted manually as examples].         156  (a1) (a2) 500 µm 500 µm (b2) (b1) 20 µm 20 µm Figure 7.9 Scanning electron microscopy images of one cross-section of apple flesh tissue 10 mm beneath the skin for ‘Granny Smith’ at day 1 (a1) and day 30 (a2) acquired at low resolution for an overview of the cell structures, and at day 1 (b1) and day 30 (b2) acquired at high resolution for a more detailed view of the cell-to-cell connection [the red lines in (a1) denote five cell connections counted manually as examples]. Table 7.2 shows the morphological features obtained from the SEM images of apple tissues. The energy and entropy changes for the two cultivars were not consistent over the storage time, and so was the change in fractal dimension D. There were differences in the microstructural characteristics, cell topology and physiological activity during storage as well as the level of maturity at harvest for GD and GS apples. These differences have thus resulted in different 157  morphological features, as described by their kinetic changes. A general tendency of decrease in dmin and increase in dper during storage was observed. A greater value of dmin or dper corresponds to a larger cell size or a longer distance between cells. Since only two specimens were used for each cultivar for each storage time in SEM imaging, no statistical test could be done for the five storage times. Compared with the cell size parameters obtained from the CLSM images, the change pattern for the morphological features extracted from the SEM images for different storage times was not apparent. It has to mention that only the intact cells in the CLSM images were segmented for extracting the cell size parameters; whereas for the analysis of SEM images, the entire images including both intact cells and void structures (i.e., pores or broken cells) were considered, which would have posed more difficulty for the differentiation and characterization of morphological features from the SEM images. Further, large variations in the SEM images were observed because of the small imaging area. It is, therefore, suggested that more specimens for the same apple fruit should be considered in future experiment in order to better assess the biological changes of apple fruit during storage. Table 7.2 Morphological features extracted from the scanning electron microscopic (SEM) images acquired at the low resolution for ‘Golden Delicious’ and ‘Granny Smith’ apple tissues.* Storage day First order histogram probability Fractal analysis parameters Energy Entropy D dmin (pixel) dper (pixel) 4.942 5.188 5.300 4.974 5.276 2.589 2.579 2.571 2.546 2.561 40.68 35.44 29.77 31.38 33.51 159.4 162.2 164.8 161.4 173.7 5.084 5.239 5.216 5.283 2.577 2.622 2.581 2.608 40.00 42.09 38.73 42.60 169.5 182.3 173.7 168.4 Golden Delicious Day 1 20352 Day 8 17048 Day 15 18098 Day 22 17415 Day 30 16800 Granny Smith Day 1 15607 Day 8 15936 Day 15 13975 Day 22 16998 158  Table 7.2 (cont’d) Day 30 18832 5.310 2.578 39.50 189.1 *D is the fractal dimension at low scale, dmin is the parameter of the minimum elemental cell, and dper is the parameter of the periodical region at high scale. 7.3.4 Correlations between Optical Properties and Mechanical/Structural Properties Table 7.3 shows correlation coefficients for the four optical parameters (µa, 675 nm, µs', 675 nm, aµs' and bµs') with the acoustic and impact firmness (FI and IF) and Young’s modulus (E) for GD and GS apples, where µa, 675 nm and µs', 675 nm represent the absorption and reduced scattering coefficients at the wavelength of 675 nm, which is related to the chlorophyll absorption, and aµs' and bµs' are the parameters related to the µs' spectra (i.e., µs' = aµs'λ -bµs' , where λ is the wavelength in nm). Optical parameter µa, 675 nm was highly correlated with the acoustic and impact firmness (r=0.870 and 0.918, respectively) for GD; however, the correlations for GS were much lower (r=0.334 and 0.421 respectively). Scattering parameter µs', 675 nm generally had better correlation with acoustic and impact firmness for both cultivars (r=0.9030.993). Both µa, 675 nm and µs', 675 nm had much lower correlations with Young’s modulus for the two cultivars, although better correlations were obtained for GD than for GS (Table 7.3), which could have been due to the relatively small change of chlorophyll in GS apples during the storage. It should also be noted that the measurement of Young’s modulus was subject to error due to variations in the dimension and shape, and sampling position of tissue specimens, whereas acoustic measurements were taken for the entire fruit and, thus, were more consistent. 159  Table 7.3 Correlations of selected optical parameters with acoustic and impact firmness and Young’s modulus for ‘Golden Delicious’ and ‘Granny Smith’ apples.* Optical parameters µa, 675 nm µs', 675 nm aµs' bµs' Golden Delicious FI IF E 0.870 0.918 0.585 Granny Smith FI IF 0.334 0.421 E 0.292 0.903 0.932 0.766 0.993 0.992 0.694 0.948 0.939 0.947 0.902 0.941 0.584 0.924 0.938 0.804 0.974 0.974 0.620 *µa, 675 nm and µs', 675 nm are absorption and reduced scattering coefficients at 675 nm, and aµs' - bµs' , where λ is the wavelength and bµs' are the parameters related to the µs' spectra (µs' = aµs'λ in nm); FI and IF are the acoustic firmness index and impact firmness index, respectively, and E is Young’s modulus. Since aµs' and bµs' are the two parameters in the wavelength-dependent exponential function for quantifying the entire µs' spectra, they should be more useful to characterize the main features of µs' spectra, compared with using single wavelengths. As expected, higher correlations between these two optical parameters and acoustic/impact firmness (r=0.902-0.974) were obtained. Likewise, better correlations of aµs' and bµs' with the Young’s modulus were also obtained, compared with µa, 675 nm and µs', 675 nm for GD, and compared with µa, 675 nm for GS. Correlations between the optical properties including transport mean free path (mfp) and tissue structural parameters (i.e., cell size parameters from CLSM images and fractal analysis parameters from SEM images) are shown in Table 7.4 and Table 7.5. The transport mean free path [or mfp = 1/(µa+µs')] represents the average path distance by moving photons between successive interactions with the fruit tissue. There were positive correlations between the optical parameters (i.e., µa, 675 nm, µs', 675 nm, aµs' and bµs', mfp) and cell size parameters (i.e., area, 160  equivalent diameter, major length, minor length and perimeter) for both cultivars. Values of the five optical parameters decreased with storage time, as the cell size decreased due to the disintegration and/or separation of cells. Better correlations (r = 0.843-0.955) were obtained between µa, 675 nm and cell size parameters for GD apples, which could have been attributed to the larger change in the absorption coefficient at 675 nm from day 1 to day 30, compared with GS apples. For GS apples, aµs' had better correlation (r = 0.876-0.906) with cell size parameters among the five optical parameters. Table 7.4 Correlations of optical parameters with the cell size parameters extracted from the confocal laser scanning microscopic images of ‘Golden Delicious’ and ‘Granny Smith’ apples. Area Optical parameter Golden Delicious 0.934 µa, 675 nm Equivalent diameter Major length Minor length Perimeter 0.941 0.912 0.843 0.955 µs', 675 nm 0.777 0.657 0.845 0.686 0.842 aµs' 0.767 0.774 0.670 0.754 0.730 bµs' 0.660 0.657 0.627 0.600 0.646 mfp'675nm 0.830 0.872 0.890 0.728 0.891 Granny Smith µa, 675 nm 0.630 0.607 0.593 0.613 0.547 µs', 675 nm 0.581 0.572 0.716 0.372 0.681 aµs' 0.903 0.906 0.891 0.876 0.902 bµs' 0.850 0.845 0.796 0.861 0.796 mfp'675nm 0.704 0.690 0.826 0.499 0.779 Table 7.5 Correlation of optical parameters with the fractal analysis parameters extracted from the scanning electron microscopic images of ‘Golden Delicious’ and ‘Granny Smith’ apples. Optical parameter µa, 675 nm Golden Delicious Granny Smith dmin 0.512 dmin 0.740 dper -0.820 161  dper -0.706 Table 7.5 (cont’d) µs', 675 nm 0.276 -0.755 0.517 -0.556 aµs' 0.821 -0.572 0.244 -0.689 bµs' 0.770 -0.582 0.549 -0.612   The four optical parameters (not including mfp) were positively correlated with dmin,but were negatively correlated with dper (Table 7.5). dmin decreased with storage time due to the change of the cell size, while dper increased as a result of the rupture and clustering of cells, causing separation between the cell structures. These microstructural changes were primarily responsible for changes in the optical properties. However, correlations between the optical parameters and fractal analysis parameters from the SEM images were relatively low due to the high variability in the SEM images and fewer specimens used for each test date. Hence in further research, more specimens with greater physiological variations should be considered in order to more accurately quantify the relationship between the optical and microstrucultural parameters. 7.4 Conclusions This research, for the first time, provided quantitative information on how the optical absorption and scattering properties were related to the elastic properties and the cellular structure of apples stored in a high temperature (~22 C)/high humidity (95%) environment. The µa value at 675 nm for both ‘Golden Delicious’ and ‘Granny Smith’ cultivars decreased consistently over 30 days of storage due to the decrease in chlorophyll-a, while the µa value increased around 525 nm, likely resulting from the increase in anthocyanin during the fruit storage. The reduced scattering coefficient of apple fruit also exhibited a consistent pattern of 162  change during storage; its value for the fresh apples was considerably higher than that for the stored apples for all test dates. The changes or decreases in the optical properties were accompanied with decreases in the acoustic/impact firmness, Young’s modulus, and cell size parameters with storage time. These changes were directly related to changes in the cellular structures/properties, primarily caused by the degradation of pectin and decrease in the molecular weight distribution of hemicelluloses (Harker et al., 1997). The optical parameters (i.e., µa, 675 nm, µs', 675 nm, aµs', bµs', and mfp) had various degrees of correlation  with acoustic/impact firmness, Young’s modulus, and cell morphological parameters for both cultivars. These findings suggest the potential of using optical properties to study or monitor the mechanical properties and microstructural changes of apples during storage and/or in the softening or ripening process. 163  CHAPTER 8 DETERMINATION OF THE OPTICAL PROPERTIES OF TWO-LAYER TURBID MATEIRALS 8.1 Introduction Previous studies have shown that hyperspectral imaging-based spatially-resolved technique is useful for determining the optical properties of fruits and food products that are approximately homogeneous. It is, however, desirable to consider fruit to be composed of two homogenous layers of skin and flesh in order to better quantify light propagation in the fruit and more accurately assess fruit quality attributes. Light transfer in two-layer media is more complex than in homogeneous media. Hence, determination of the optical properties of two-layer media entails great challenges both mathematically and experimentally. A two-layer diffusion model has four (or five) property parameters (i.e., µa and µs' of each layer, plus the unknown thickness of the first layer). Parameter estimation accuracy, computational time, and model efficiency must be considered in the development of an inverse algorithm for layered media. Several studies have provided analytical solutions to the diffusion equations for layered media (Alexandrakis et al., 2001; Hollmann and Wang, 2007; Kienle et al., 1998b; Schmitt et al., 1990). Numerical methods (e.g., Monte Carlo, asymptotic approximation, and finite element) have also been used to extract optical properties from the diffusion models (Gonzalez-Rodriguez and Kim, 2008; Schweiger and Arridge, 1997; Seo et al., 2007). Kienle et al. (1998b) derived an analytical form of the twolayer diffusion model, which enables fast forward computation and inverse algorithm implementation to estimate the optical properties. The increased number of free parameters in the two-layer model can dramatically increase computational time, exacerbate the process of 164  estimating the optical parameters, and/or cause convergence problems for the inverse algorithm. It is thus important to understand the intrinsic properties of the diffusion model and evaluate the feasibility of estimating all parameters simultaneously prior to implementing an inverse algorithm. Therefore, the overall objective of this research was to explore the feasibility of determining the optical properties of the surface layer (skin) and subsurface layer (flesh) for those fruits and food products that can be modeled as two-layer homogeneous media. The specific objectives were to: 1) perform sensitivity coefficient analysis and develop an inverse algorithm for a twolayer diffusion model to estimate the optical properties of each homogeneous layer; 2) create solid model samples and measure their optical properties using an integrating sphere method; and 3) validate the diffusion model and the inverse algorithm using Monte Carlo simulations and experimental data acquired from the model samples using hyperspectral imaging-based spatiallyresolved technique. 8.2 Materials and Methods 8.2.1 Two-layer Diffusion Model Consider a two-layer turbid medium which is impinged upon perpendicularly by an infinitely small light beam. The thickness of the first layer (d) is assumed to be larger than one transport mean free path of the first layer [ z0  1/ (a1  s1 ') , in which µa1 and µs1' are absorption and reduced scattering coefficients of the first layer]. Under steady-state condition, the fluence rate, defined as the irradiance that is incident from all solid angles onto a small sphere in a two-layer turbid medium, is given by (Kienle et al., 1998b) D121(r )  a11(r )    x, y, z  z0  165  0 zd (8.1) D2 2 2 (r )  a 2 2 (r )  0 dz (8.2) where  i is the fluence rate of layer i , and Di  1 / [3(  ai   si ')] is the diffusion constant. Equations (8.1) and (8.2) can be converted to ordinary differential equations using a twodimensional Fourier transform (1998b). The following boundary conditions are applied for obtaining solutions for the fluence rates [ i ( z , s ) ] of each layer on the frequency domain (Haskell et al., 1994; Kienle et al., 1998b): the same refractive index for the first layer and second layer, the zero fluence rate at the extrapolated boundary, a finite photon fluence rate as z   , and the continuity of fluence rate at the boundary between the first and second layer. After the differential equations in the frequency domain are solved, the two-dimensional inverse Fourier transform is implemented to obtain the following fluence solution on the Cartesian coordinate system:  i (r , z )  1  2    1   ( z , s ) exp  i  s1 x  s2 y  ds1ds2  i ( z , s ) sJ 0 ( sr ) ds   2   i 2 0    (8.3) where r  ( x2  y 2 )1/2 , and J 0 is the zeroth-order Bessel function. This inverse transform is achieved by numerically evaluating the integral using an adaptive Gauss-Kronrod quadrature (Shampine, 2008). Spatially-resolved diffuse reflectance R(r) is obtained from the integration of radiance over the solid angle accepted by the fiber. Then, the final solution of the two-layer diffusion model is given as (Kienle et al., 1998b) R2layer (r )  C11  r , z  0   C2 D1  1  r , z  z z 0 (8.4) The details related to the parameters C1 and C2 can be found in Chapter 3. Based on equation (8.4), the diffuse reflectance R(r) at the surface of the two-layer turbid medium is a function of the source-detector distance (r) as well as five unknown optical parameters of the two layers and 166  the thickness of the first layer (µa1, µs1', µa2, µs2', and d, where µa2 and µs2' are absorption and reduced scattering coefficients of the second layer). Researchers have used the two-layer model [equation (8.4)] to model photon propagation in turbid media under steady-state conditions (Hollmann and Wang, 2007; Kienle et al., 1998b; Tseng et al., 2008). 8.2.2 Model Samples Preparation and Measurement To investigate the feasibility of the two-layer diffusion model and test the performance of the hyperspectral imaging system for measuring the optical properties of two-layer materials, the two-layer model samples were created, which were made up of two-component roomtemperature-vulcanizing silicone (part A + part B, ELASTOSIL®RT 604, Wacker Chemie AG, Munich, Germany) as the host, aluminum oxide particles (Al2O3, Duke Scientific, Fremont, California, USA) as scattering materials, and blue dye (Direct Blue 71, Sigma-Aldrich, St. Louis, MO, USA) as absorbers. The stock solution of blue dye was prepared in ethanol with the concentration of 0.6 mg/L. The solution of the absorber was mixed with 360 g part A of silicone with a stirrer and heated gently for several hours to remove the ethanol. The scattering materials were dispersed in 40 g part B of silicone and put in the sonicator for several minutes to break up large particles. Then, these components were mixed together and the bubbles produced in the process of mixing were removed in the vacuum chamber. Two silicone mixtures were made by adding different absorbing stock solutions (34 mL and 27 mL) and scattering materials (8 g and 5 g). Two cylindrical model samples with a diameter (D) of 81 mm and a height (h) of 57 mm and two disk samples of 81 mm diameter and a height of 2.64 mm and 3.29 mm respectively were created by pouring the mixtures into the plastic cylindrical and disk modes. One disk sample corresponded to one cylindrical sample with the same optical properties. Finally, the two-layer 167  model samples were obtained by attaching a disk sample with different optical properties onto a cylindrical sample. These two model samples would satisfy the assumption for the two-layer diffusion model with a limited thickness of the first layer and semi-infinite second layer. An integrating sphere system, which has been used as a standard method for measuring the optical properties of solid turbid samples in many reported studies (Prahl et al., 1993; Saeys et al., 2008), was used to obtain reference values of the absorption and reduced scattering coefficients for each layer of the two model samples. Total diffuse reflectance and transmittance measurements were performed on the two model disk samples of 3.29 mm and 2.64 mm thick, respectively, using the same configurations shown in Figure 4.4 in Chapter 4. An inverse addingdoubling algorithm was used to calculate the scattering and absorption coefficients of the model samples based on the total reflectance and transmittance. A detailed description of the method can be found in Prahl et al. (1993). 8.3 Validation of the Two-layer Diffusion Model and Inverse Algorithm Like the single-layer diffusion model described in Chapter 3, the sensitivity analysis was also performed, and sensitivity coefficients were first calculated for each of the four optical parameters (i.e., the absorption and reduced scattering coefficients for each layer) by the following equations, which were derived from equation (2.29): R ai  ai R ai R si '   si ' R  si ' where µai and µsi' are the absorption and reduced scattering coefficients of layer i. 168  (8.5) (8.6) For accurate determination of the optical properties of fruit and food samples, the diffusion model and the inverse algorithm were then validated by Monte Carlo simulations. A publicly available MC simulation program for multi-layered media was used (Wang et al., 1995b). The low-noise reflectance generated by MC simulation provided an ideal measurement without 6 experimental uncertainties. A total of 3×10 photons were used to produce the reflectance for a spatial distance of 0.1-10 mm at 0.1 mm spatial resolution for both radial distance and depth. Four different combinations of µa1, µs1', µa2, and µs2' were selected, which span a large range of -1 values: 0.05 ≤ µa ≤ 2.60 cm -1 and 12 ≤ µs' ≤ 36 cm . These values were selected based on the published data for fruit skin and flesh (Budiastra et al., 1998; Saeys et al., 2008). The diffusion model was first compared with MC simulation for generating spatially-resolved reflectance profiles for a two-layer medium with the selected values for the four optical parameters. Then, the MC generated reflectance profiles were fitted by the inverse algorithm of the two-layer diffusion model, from which the optical properties for each layer were determined. MC simulations of light propagation in a two-layer turbid medium introduced small stochastic variability in the calculated reflectance for any given set of optical properties. The resulting variability was less than 1.7% in 20 simulations for each set of optical properties. The diffusion model and inverse algorithm were further validated with experimental data obtained from the two-layer model samples using the laboratory hyperspectral imaging system (Figure 4.1), following the same procedures for liquid model samples described in Section 4.2.1. The time required for each scan of the two-layer mode sample was 100-200 ms. Two cases were considered: 1) only two unknown optical parameters for one of the two layers with known optical properties for the other layer, and 2) four unknown optical properties for both layers. 169  Finally, the optical properties of the model samples estimated by the inverse algorithm were compared with those measured using the integrating sphere method. 8.4 Results and Discussion 8.4.1 Monte Carlo Simulations of Spatially-resolved Diffuse Reflectance Figure 8.1 compared diffuse reflectance from the two-layer diffusion model and MC simulations for the four sets of optical properties (µa1/µa2 = 0.50 and µs1'/µs2' = 0.86, µa1/µa2 = 6.50 and µs1'/µs2' = 1.80, µa1/µa2 = 0.80 and µs1'/µs2' = 1.58, µa1/µa2 = 2 and µs1'/µs2' = 0.75). The reflectance profiles calculated from the diffusion equation matched those from the MC simulations. The differences for the given values of optical properties used in Figure 8.1(a), (c) and (d) were less than 6% for source-detector distances greater than 1.5 mm. For case 2 with µa1/µa2 = 6.50 and µs1'/µs2' = 1.80 [Figure 8.1(b)], the average difference was 12%, greater than those for the other three cases. This was likely due to low levels of diffuse reflectance resulting from high absorption of both layers. Larger deviations of reflectance were observed when the detector was close to the light source (less than one transport mean free path), because the diffusion approximation is not valid in this situation. These results indicate that overall the diffusion model accurately quantifies light propagation in two-layer turbid media. 170  (a) 10 10 -2 -4 -6 0 2 0 4 6 8 Distance (mm) (c) 10 10 10 10 10 (b) 0 -2 -4 -6 0 2 4 6 8 Distance (mm) (d) 10 0 2 4 6 8 Distance (mm) 10 0 -2 -2 Reflectance (mm ) 10 10 Reflectance (mm ) 10 -2 Reflectance (mm ) 0 -2 Reflectance (mm ) 10 10 10 10 -2 -4 -6 0 2 4 6 8 Distance (mm) 10 10 10 -2 -4 Figure 8.1 Comparison of spatially-resolved diffuse reflectance obtained from the two-layer diffusion model (symbols) and Monte Carlo simulations (solid lines): (a) µa1/µa2 = 0.50 and µs1'/µs2' = 0.86 (0.05, 12, 0.10, 14, 1), (b) µa1/µa2 = 6.50 and µs1'/µs2' = 1.80(2.60, 36, 0.40, 20, 0.5), (c) µa1/µa2 = 0.80 and µs1'/µs2' = 1.58 (0.32, 19, 0.40, 12, 1), and (d) µa1/µa2 = 2 and µs1'/µs2' = 0.75 (1.00, 15, 0.50, 20, 1). Values in the parentheses are µa1, µs1', µa2, µs2', d with -1 the units of cm 8.4.2 for optical properties and mm for the thickness of the first layer. Sensitivity Coefficient Analysis The scaled sensitivity coefficients for the absorption and reduced scattering coefficients of the two layers were calculated as functions of the source-detector distance for the corresponding diffuse reflectance profiles (Figure 8.2). For the four combinations of optical properties discussed in this study, the magnitudes of the reduced scattering sensitivity coefficients are 171  almost equal to those of R, while the sensitivity coefficients of absorption are relatively smaller, especially in Figure 8.2(a) because of the small values of the selected absorption coefficients in this case. For the other three combinations, the sensitivity coefficients of absorption are still on the same order of R at the large source-detector distances. Also it is noted that the shapes of the four sensitivity coefficients are quite different. These observations show that the four sensitivity coefficients of µa1, µs1', µa2, and µs2' are ‘large’ (i.e., on the order of R) and uncorrelated (different shapes), which are desirable conditions for estimating these parameters. Since the reduced scattering sensitivity coefficients are larger in magnitude than the absorption sensitivity coefficients in Figure 8.2, on average µsi' should be estimated more accurately than µai for such 10 10 10 (a) 0 -2 10 Scaled Sensitivity(mm ) -2 Scaled Sensitivity(mm ) model. -2 -4 -6 0 2 4 6 8 Distance (mm) 10   10 10 10 10 (b) 0 -2 -4 -6 0 2 4 6 8 Distance (mm) 10 Figure 8.2 Scaled sensitivity coefficients as functions of source-detector distances for four combinations of optical properties (a) µa1/µa2 = 0.50 and µs1'/µs2' = 0.86 (0.05, 12, 0.10, 14, 1), (b) µa1/µa2 = 6.50 and µs1'/µs2' = 1.80 (2.60, 36, 0.40, 20, 0.5), (c) µa1/µa2 = 0.80 and µs1'/µs2' = 1.58 (0.32, 19, 0.40, 12, 1), and (d) µa1/µa2 = 2 and µs1'/µs2' = 0.75 (1.00, 15, 0.50, 20, 1), and -1 the values in the bracket are µa1, µs1', µa2, µs2', d with the unit cm for optical properties and mm for the thickness of the first layer (‘-’, ‘·’, ‘+’, ‘ₒ’, and ‘∆’denote reflectance, and scaled sensitivity coefficients of µa1, µs1', µa2, and µs2', respectively).   172  10 10 10 (c) 0 -2 10 Scaled Sensitivity(mm ) -2 Scaled Sensitivity(mm ) Figure 8.2 (cont’d) -2 -4 -6 0 2 4 6 8 Distance (mm) 10 10 10 10 10 (d) 0 -2 -4 -6 0 2   4 6 8 Distance (mm) 10   The effect of the thickness of the first layer on the sensitivity coefficients of optical properties is presented in Figure 8.3 for four different thicknesses of the first layer with the same optical properties. Values for the four sensitivity coefficients vary with the thickness (d) that changes from 0.85 mm to 6 mm. For the thin thickness of the first layer, it is easy to estimate optical properties of the second layer. While for the larger d, it is more difficult to obtain accurate estimations of optical properties of the second layer due to small sensitivity coefficients for µa2 and µs2'. This may be explained by the fact that only a small number of photons 10 10 10 (a) 0 -2 10 Scaled Sensitivity(mm ) -2 Scaled Sensitivity(mm ) propagating through the second layer are remitted when the first layer is relatively thick. -2 -4 -6 0 2 4 6 8 Distance (mm) 10 10 10 10 10 (b) 0 -2 -4 -6 0 2 4 6 8 Distance (mm) 10 Figure 8.3 Scaled sensitivity coefficients as functions of source-detector distances for µa1 = 0.05 -1 -1 -1 -1 cm , µs1' = 12 cm , µa2 = 0.10 cm , µs2' = 14 cm , and the thickness of the first layer: (a) d = 173  0.85 mm, (b) d = 2 mm, (c) d = 4 mm, and (d) d = 6 mm (‘-’, ‘·’, ‘+’, ‘ₒ’, and ‘∆’denote reflectance, and scaled sensitivity coefficients of µa1, µs1', µa2, and µs2', respectively). 8.4.3 10 10 10 (c) 0 -2 10 Scaled Sensitivity(mm ) -2 Scaled Sensitivity(mm ) Figure 8.3 (cont’d) -2 -4 -6 0 2 4 6 8 Distance (mm) 10 10 10 10 10 (d) 0 -2 -4 -6 0 2 4 6 8 Distance (mm) 10 Extraction of Optical Properties from Monte Carlo Simulation Data When the optical properties of the second layer and the thickness of the first layer are treated as known, the inverse algorithm gives almost zero errors in estimating absorption and reduced scattering coefficients for the first layer from the noiseless reflectance data produced by the diffusion model. These results indicate that it is possible to accurately estimate two optical parameters in the diffusion model using nonlinear least squares method. Figure 8.4 shows the estimated absorption and reduced scattering coefficients of the first layer from the two-layer diffusion model fitted to the MC data that has a noise level of 2%. For -1 case 1 [Figure 8.4(a)], the actual value of µa1 varies between 0.20 to 1.20 cm -1 -1 with µs1' = 12 -1 cm , µa2 = 0.50 cm , µs2' = 11 cm , and d = 2 mm. The average errors of estimating µa1 and µs1' are 4.2% and 4.1% compared with the true values. While for case 2 [Figure 8.4(b)], µs1' -1 changes in the range of 15-45 cm with the fixed values of µa1, µa2, µs2' and L. The average 174  errors are 12.7% for µa1 and 5.8% for µs1'. The optical parameters of the second layer were also estimated with the optical parameters and the thickness of the first layer being treated as known. For the same combination of the optical parameters as in case 1, the average errors are 10.0% for µa2, and 9.3% for µs2', higher than those for case 1. The statistic uncertainty of MC simulations reduces the accuracy of the inverse algorithm for parameter estimations compared with the fitting results using the noiseless diffusion model-generated data. (a) 0.14 Estimated  a1 (mm -1) 0.12 0.1 0.08 0.06 0.04 0.02 0 0.02 0.04 0.06 0.08 True a1 (mm 0.1 0.12 -1 ) Figure 8.4 Estimated absorption and reduced scattering coefficients of the first layer from fitting -1 the diffusion model to the Monte Carlo simulation data: (a) µa1 varies 0.20 - 1.20 cm with µs1' -1 -1 -1 -1 = 12 cm , µa2 = 0.50 cm , µs2' = 11 cm , and d = 2 mm, and (b) µs1' varies 15-45 cm -1 -1 -1 µa1= 0.10 cm , µa2 = 0.05 cm , µs2' = 11 cm , and d = 2 mm. 175  with Figure 8.4 (cont’d) (b) 4.5 Estimated s1 (mm -1) 4 3.5 3 2.5 2 1.5 1.5 2 2.5 3 True s1 (mm 3.5 4 4.5 -1 ) The four optical parameters of µa1, µs1', µa2, µs2'  were also estimated simultaneously with the known d, and the same four sets of optical properties (µa1/µa2 = 0.50 and µs1'/µs2' = 0.86, µa1/µa2 = 6.50 and µs1'/µs2' = 1.80, µa1/µa2 = 0.80 and µs1'/µs2' = 1.58, µa1/µa2 = 2 and µs1'/µs2' = 0.75) shown in Figure 8.1 were used. The average errors of estimating µa1, µs1', µa2, and µs2' from the reflectance generated by the diffusion model are 1.4%, 0.01%, 0.4% and 0.1%, respectively, which shows that the inverse algorithm can accurately estimate four parameters based on the noiseless data. Mean and maximum errors in the estimated absorption and reduced scattering coefficients of both layers from the MC generated reflectance are presented in Figure 8.5. The mean errors for estimating the first layer’s optical properties are 14.5% and 8.1% for the absorption and reduced scattering coefficients, respectively. For the second layer, the mean 176  errors are 29.3% and 9.1% for the absorption and reduced scattering coefficients, respectively. The relatively smaller errors are obtained for the optical parameters of the first layer compared with those in the second layer, and the errors of estimating the reduced scattering coefficients are smaller than those for the absorption coefficients. The maximum error for µa2 reaches 40%, two to four times greater than that for the other three parameters. The relative scales of the estimation errors for the four optical parameters are in agreement with the findings from the sensitivity analysis. Comparison of the results for estimation of two and four parameters demonstrates that the errors increased with the number of parameters to be estimated in the two-layer model. Also, to obtain a good estimation, the signal-to-noise ratio must be considered in the measurement of the reflectance. 45 mean max 40 Relative error (%) 35 30 25 20 15 10 5 0 a1 a2 s1 s2 Optical properties   Figure 8.5 Relative errors of estimating the optical properties of two layers from the Monte Carlo diffuse reflectance for four combination sets of optical properties: µa1/µa2 = 0.50 and µs1'/µs2' = 177  0.86 (0.05, 12, 0.10, 14, 1), µa1/µa2 = 6.50 and µs1'/µs2' = 1.80 (2.60, 36, 0.40, 20, 0.5), µa1/µa2 = 0.80 and µs1'/µs2' = 1.58 (0.32, 19, 0.40, 12, 1), and µa1/µa2 = 2 and µs1'/µs2' = 0.75 (1.00, 15, -1 0.50, 20, 1). Values in the parentheses are µa1, µs1', µa2, µs2', d with the units of cm properties and mm for the thickness of the first layer. 8.4.4 for optical Optical Properties of Model Samples Measured with Integrating Sphere Values of the absorption and reduced scattering coefficients for the two model disk samples measured using the integrating sphere at wavelengths of 500-1,000 nm are shown in Figure 8.6. Light absorption by the blue dye absorber is observed at 535 nm for both samples. The absorption peak around 910 nm was caused by the silicone. The reduced scattering spectra are relatively flat compared with the absorption spectra of the model samples. The standard method has good repeatability of measuring the model samples, whose variation coefficients for the absorption and reduced scattering coefficients from four replicated measurements are 2.4% and 2.3%, respectively. The optical properties of the two-layer model samples at 535 nm and 700 nm are listed in Table 8.1; they were used as the true optical values when compared with those obtained from the spatially-resolved technique with the inverse algorithm. (a) 0.6 -1 a (cm ) Disk sample 1 Disk sample 2 0.4 0.2 0 500 600 700 800 Wavelength (nm) 900 1000 Figure 8.6 Absorption and reduced scattering spectra of the homogenous model disk samples measured by the integrating sphere. 178  Figure 8.6 (cont’d) (b) 14 -1 s (cm ) 12 10 Disk sample 1 Disk sample 2 8 6 500 600 700 800 Wavelength (nm) 900 1000   Table 8.1 Optical properties of two two-layer model samples at the wavelengths of 535 nm and 700 nm determined by the integrating sphere and the adding-doubling method (Prahl et al., 1993). Model sample Wavelength (nm) Sample 1 µa1 -1) µs1' -1) µa2 -1) µs2' -1) 535 700 (cm 0.51 0.18 (cm 12.97 11.72 (cm 0.45 0.13 (cm 8.08 7.29 535 700 0.45 0.13 8.08 7.29 0.51 0.18 12.97 11.72 d (mm) 2.64 Table 8.1 (cont’d) Sample 2 3.29   8.4.5 Optical Properties of Model Samples Determined from Hyperspectral Imaging Measurements Figure 8.7 shows the measured reflectance profiles of the two two-layer model samples at wavelengths of 535 nm and 700 nm, and also the theoretical profiles calculated from the twolayer diffusion model using the optical parameters given in Table 8.1. These two particular wavelengths were selected in this study because they represented the lower and upper range of values for the absorption and scattering coefficients for the model samples over the spectral range of 500-1,000 nm (Figure 8.6). Hence the estimation errors for the two wavelengths should 179  be representative for the other wavelengths. Because the diffuse reflectance profiles from the experiment and the diffusion model were not on the same scale, the profiles were normalized at a distance of 1.5 mm. As shown in Figure 8.7, the experimental measurements covering the source-detector distance from 1.5 mm to 10 mm are in good agreement with the theoretical predictions. However, differences between the measurements and predictions become larger at distances close to the source due to the violation of diffusion theory assumptions. (a) (b) 4 4 10 Reflectance (a.u.) Reflectance (a.u.) 10 2 10 0 10 -2 10 4 6 8 Distance (mm) 10 10 (c) 4 2 4 6 8 Distance (mm) 10 (d) 4 10 Reflectance (a.u.) 10 Reflectance (a.u.) 0 10 -2 2 2 10 0 10 2 10 0 10 -2 -2 10 2 10 10 2 4 6 8 Distance (mm) 10 2 4 6 8 Distance (mm) (d) 10 Figure 8.7 Comparison of diffuse reflectance from measurements (symbols) and the two-layer diffusion model (solid lines) for sample 1 at 535 nm (a) and 700 nm (b), and sample 2 at 535 nm (c) and at 700 nm (d) (See Table 8.1 for the two samples). Figure 8.8 shows estimated values for the absorption and reduced scattering coefficients of the first layer when the optical properties of the second layer and the thickness of the first layer 180  are considered as known parameters. The errors of estimating µa1 and µs1' of the model samples are 11.3-23.0% and 3.8-18.4%, respectively. Overall, relatively small values of the optical coefficients produce larger errors in the estimation of the optical properties (values of µa1 and µs1' at 700 nm are less than those at 500 nm). While the results are worse than those obtained from the MC simulations, they are better than those reported by Kienel et al. (1998b) for twolayer turbid media with the errors in µa of 6-24% and in µs' of 0-88%. 25 535 nm 700 nm Relative error (%) 20 15 10 5 0 S1-a1 S1-s1 S2-a1 Optical properties S2-s1   Figure 8.8 Relative errors of the estimated optical properties for two model samples (denoted as S1 and S2) at wavelengths of 535 nm and 700 nm (see Table 8.1 for the optical property data for the two model samples). When the optical properties of the second layer are considered as only unknown and those of the first or upper layer are treated as known, the inverse algorithm does not give reasonable estimates. This may have been attributed to the effect of the relatively thick first layer (2.64 mm 181  and 3.29 mm). These results suggest that the two-layer model coupled with the hyperspectral imaging technique can be used to obtain the optical properties of fruit skin with the thickness similar to those studied in this research. However, further study is needed to determine if the method is also suitable for thin-skin fruit like apple. Unsatisfactory results for estimating the optical properties of the lower or second layer for the model samples suggest that further improvements in the method are needed in order to measure the optical properties of both fruit skin and flesh simultaneously. 8.5 Conclusions This research has demonstrated a methodology for determining the optical properties of twolayer turbid materials, using a steady-state diffusion model for spatially-resolved reflectance acquired by a hyperspectral imaging system. The diffusion model accurately described light propagation in two-layer turbid media. Accurate estimations for two and four parameters were achieved based on the noiseless data from the diffusion model. The inverse algorithm also gave good estimates of two unknown optical parameters of the first layer or the second layer for MC generated reflectance data, with errors of 8.5% for µa1 and 5.0% for µs1', 10.0% for µa2 and 9.3% for µs2'. The errors, however, were greater in estimating four optical parameters. Experimental validation with the model samples showed larger errors (11.3-23.0% for µa1 and 3.8-18.4% for µs1') in estimating the optical properties of the first layer. Satisfactory results can be achieved when two unknown parameters for the first or upper layer are to be determined. While the results for determining the two optical parameters of the second or lower layer of the model samples are still unsatisfactory, Monte Carlo simulations for 182  ideal samples suggest that it is possible to accurately estimate the optical properties of both layers simultaneously. Hence further study should be focused on improving the reflectance measurement by the hyperspectral imaging system and the inverse algorithm. To accurately evaluate the diffusion model for determining the optical properties of two-layer samples, it is important that standard model samples with known optical properties be made available for calibrating or validating the inverse algorithm. Efforts should also be made to reduce the computational complexity and time for estimating the four optical parameters plus the thickness of the first layer simultaneously. 183  CONCLUSION AND FUTURE RESEARCH This dissertation research has made significant progress in the development of new methods and techniques for optical characterization and nondestructive quality evaluation of horticultural products. A hyperspectral imaging-based spatially-resolved technique, along with a bench-top optical property measuring prototype, was developed, through optimization of the inverse algorithm and optical designs, to accurately and reliably measure the optical absorption and scattering properties of horticultural and food products. The instrument was used for measuring the optical properties of apple and peach and for assessing their maturity/quality attributes. Moreover, research was conducted in quantification of the relationships of the optical properties with the mechanical properties and microstructural characteristics of apples that were subject to accelerated softening in a high temperature (~22 ºC)/high humidity (95%) environment. Further study was also carried out to evaluate the feasibility of determining the optical properties of twolayer media representing fruit skin and flesh. The following major conclusions were drawn from the research: (1) An inverse algorithm for estimating the absorption and reduced scattering coefficients from spatially-resolved reflectance profiles based on a single-layer diffusion model was developed and optimized by using data transformation and weighting methods. The logarithm and integral transformation of the original data and the relative weighting method greatly improved the estimations of the two optical parameters with the relative errors of 10.4%, 10.7%, and 11.4% for µa, and 6.6%, 7.0% and 7.1% for µs' for liquid model samples made up of dye and intralipid, which were much smaller than those obtained from the original diffusion model. Further statistical analysis showed that the logarithm transformation and 184  relative weighting methods improved the inverse algorithm for more reliable estimation of the two optical parameters. (2) Main optical design factors, including light beam and source-detector distance, in the development of hyperspectral imaging-based spatially-resolved technique were examined and optimized through Monte Carlo (MC) simulation and experimental studies. To achieve the best performance of the system, the light beam should be of circular shape and Gaussian type with the diameter of less than 1 mm, the optimal minimum source-detector distance should be about 1.5 mm, and the optimal maximum source-detector distance should be equivalent to 10-20 mfp' or determined by the minimum signal-to-noise ratio of 20 (or 150 CCD counts for the system used in this research). (3) A new multi-purpose Optical Property Analyzer (OPA) prototype was designed and constructed, based on the optimization studies, for optical characterization and hyperspectral imaging of biological materials for the spectral region of 500-1,000 nm. Software incorporating two popular diffusion models (Farrell et al., 1992; Kienle and Patterson, 1997) with the corresponding inverse algorithms was developed for the system control and automatic imaging acquisition and optical properties computation. The OPA was extensively tested and evaluated for measurement accuracy, precision/reproducibility, and sensitivity using liquid model samples with known optical properties. The test results showed that the OPA has achieved acceptable accuracies for the absorption and reduced scattering coefficients, which are either comparable or superior to other reported studies using more sophisticated time-resolved, frequency-domain, or other types of spatially-resolved instruments. 185  (4) Spectra of the absorption and reduced scattering coefficients for one cultivar of peach and two cultivars of apple were determined for the spectral region of 500-1,000 nm using the Optical Property Analyzer. The absorption spectra of fruit were shaped by major pigments (i.e., chlorophyll and anthocyamin) and water in the fruit tissue, while the reduced scattering spectra decreased steadily with the increasing wavelength for most of the fruit samples. The best prediction results for peach fruit were obtained with the combinations of absorption and scattering spectra, with values of the correlation coefficient (r) for Magness-Talyor firmness, * * soluble solids content (SSC), the skin color parameter (L ) and the flesh color parameter (L ) being 0.724 (SEP = 18.13 N), 0.458 (SEP = 0.96 °Brix), 0.893 (SEP = 3.54), and 0.722 (SEP = 3.32), respectively. For the apple study, the best firmness prediction was obtained with r = 0.892 (SEP = 7.89 N) for ‘Golden Delicious’ (GD) and r = 0.863 (SEP = 8.94 N) for ‘Delicious’ (RD). The best result for SSC was achieved for the freshly harvested group with r = 0.791 (SEP = 0.69 °Brix) for GD and r = 0.842 (SEP = 0.73 °Brix) for RD. It demonstrated that measurement of optical properties provides a new means for quantifying the quality of fruit. (5) The absorption and reduced scattering coefficients, acoustic/impact firmness, Young’s modulus, and cell size for GD and ‘Granny Smith’ apples all decreased with time when they were subject to accelerated softening at high temperature(~22 C)/high humidity (95%). The decreases were caused by changes in the cellular structure and properties, primarily resulting from the degradation of pectin and decrease on the molecular weight distribution of hemicelluloses (Harker et al., 1997). Various correlations for the optical parameters with acoustic/impact firmness, Young’s modulus, and cell size and morphological parameters were found. These findings suggest the potential of using optical properties to study or 186  monitor changes in the mechanical properties and physiological condition of apple fruit during storage or in the softening process. (6) The inverse algorithm of the two-layer diffusion model gave good estimates of two unknown optical parameters of the first layer or the second layer for MC generated reflectance data, with errors of 8.5% for µa1 and 5.0% for µs1' for the first or top layer, and 10.0% for µa2 and 9.3% for µs2' for the second or substrate layer. The errors, however, were considerably greater in estimating the four optical parameters. Experimental validation with the solid model samples showed larger errors (11.3-23.0% for µa1 and 3.8-18.4% for µs1') in estimating the optical properties of the first layer. The results and findings from this research provide a fundamental understanding of optical characterization of biological materials, and a systematic guide in designing spatially-resolved optical systems for nondestructive quality measurement of horticultural and food products. The new optical property measuring instrument developed in this research can have applications for a wide range of food and biological materials. The following are recommendations for future work: (1) It is much more difficult to tackle inverse light transport problems than forward problems. Therefore, attention should be paid to studying the characteristics and complexity of the diffusion model when it is used in an inverse algorithm for estimating optical parameters. This research was based on the diffusion model with several assumptions (i.e., scatteringdominant, isotropic source, and source-detector distance greater than one transport mean free path). Further work is needed to extend the diffusion model to different biological materials and expand its scope of application. 187  (2) Although the optical system and algorithm of the OPA has been carefully evaluated and optimized using the model samples, other factors such as roughness and possible geometric irregularities on the surface of fruit should be studied for their effect on the optical properties measurement and on maturity/quality assessment of fruit. (3) This research has shown that the optical properties were correlated with the cell size and cellular structures of fruit tissue. To better quantify the interaction of light with plant tissue, it would be useful to consider a multi-scale approach (Kienle et al., 2007) for describing the light propagation in different volumes of biological tissue at macroscopic and microscopic scales, respectively. (4) Based on the preliminary study of the two-layer diffusion model for determining optical properties of fruit skin and flesh, further efforts should be made to develop light propagation models and optical techniques for measuring the optical properties of heterogeneous or multilayer biological materials, as well as to reduce the computational complexity and time for estimating the absorption and reduced scattering coefficients. (5) The OPA is a general-purpose optical instrument designed for measuring a wide range of biological materials. To expand its application to other areas and other agricultural and food products, the OPA should be further evaluated for its stability, accuracy, and measurement limit. 188  BIBLIOGRAPHY 189  BIBLIOGRAPHY Abbott, J. A. 1999. Quality measurement of fruits and vegetables. Postharvest Biology and Technology 15(3): 207-225. Abbott, J. A., and R. Lu. 1996. Anisotropic mechanical properties of apples. Transactions of the ASABE 39(4): 1451-1459. Abbott, J. A., R. Lu, B. L. Upchurch, and R. L. Stroshine. 1997. Technologies for nondestructive quality evaluation of fruits and vegetables. In Horticultural Reviews, 20: 1-120. J. Janick, ed. New York, N.Y.: John Wiley and Sons. Abbott, J. A., and D. R. Massie. 1985. Delayed light-emission for early detection chilling in cucumber and bell pepper fruit. Journal of the American Society for Horticultural Science 110(1): 42-47. Aigner, M., and B. Juttler. 2009. Distance regression by Gauss-Newton-type methods and iteratively re-weighted least-squares. Computing 86(2-3): 73-87. Alerstam, E., T. Svensson, and S. Andersson-Engels. 2008. Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration. Journal of Biomedical Optics 13(6): 060504. Alexandrakis, G., D. R. Busch, G. W. Faris, and M. S. Patterson. 2001. Determination of the optical properties of two-layer turbid media by use of a frequency-domain hybrid Monte Carlo diffusion model. Applied Optics 40(22): 3810-3821. Alexandrakis, G., T. J. Farrell, and M. S. Patterson. 2000. Monte Carlo diffusion hybrid model for photon migration in a two-layer turbid medium in the frequency domain. Applied Optics 39(13): 2235-2244. Anderson, E. R., D. J. Cuccia, and A. J. Durkin. 2007. Detection of bruises on Golden Delicious apples using spatial-frequency-domain imaging. In Proc. SPIE 6430: 64301O. T. VoDinh, W. S. Grundfest, D. A. Benaron, G. E. Cohn, and R. Raghavachari, eds. San Jose, California. Andersson-Engels, S., R. Berg, A. Persson, and S. Svanberg. 1993. Multispectral tissue characterization with time-resolved detection of diffusely scattered white-light. Optics Letters 18(20): 1697-1699. Arridge, S. R., M. Cope, and D. T. Delpy. 1992. The theoretical basis for the determination of optical pathlengths in tissue-temporal and frequency-analysis. Physics in Medicine and Biology 37(7): 1531-1560. Baltagi, B. H. 1996. Testing for random individual and time effects using a Gauss-Newton regression. Economics Letters 50(2): 189-192. 190    Banos, A., M. Wolf, C. Grawe, M. Stahel, D. Haensse, D. Fink, and R. Hornung. 2007. Frequency domain near-infrared Spectroscopy of the uterine cervix during cervical ripening. Lasers in Surgery and Medicine 39(8): 641-646. Barreiro, P., C. Ortiz, M. Ruiz-Altisent, J. Ruiz-Cabello, M. E. Fernandez-Valle, I. Recasens, and M. Asensio. 2000. Mealiness assessment in apples and peaches using MRI techniques. Magnetic Resonance Imaging 18(9): 1175-1181. Bays, R., G. Wagnieres, D. Robert, D. Braichotte, J. F. Savary, P. Monnier, and H. van den Bergh. 1996. Clinical determination of tissue optical properties by endoscopic spatially resolved reflectometry. Applied Optics 35(10): 1756-1766. Bechar, A., A. Mizrach, P. Barreiro, and S. Landahl. 2005. Determination of mealiness in apples using ultra-sonic measurements. Biosystems Engineering 91(3): 329-334. Beck, J., and K. Arnold. 1977. Parameter Estimation in Engineering and Science. New York: Wiley. Bernhardt, P. A. 1995. Direct reconstruction methods for hyperspectral imaging with rotational spectrotomography. Journal of the Optical Society of America a-Optics Image Science and Vision 12(9): 1884-1901. Bianchi, F. D., and R. D. Bonetto. 2001. FERImage: An interactive program for fractal dimension, d(per) and d(min) calculation. Scanning 23(3): 193-197. Blackwell, B. 2009. Sensitivity analysis, uncertainty propagation and validation for computational models. In Inverse Problems Symposium. Michigan State University, East Lansing, MI, USA. Blasco, J., N. Aleixos, J. Gomez-Sanchis, and E. Molto. 2009. Recognition and classification of external skin damage in citrus fruits using multispectral data and morphological features. Biosystems Engineering 103(2): 137-145. Blasco, J., N. Aleixos, and E. Molto. 2003. Machine vision system for automatic quality grading of fruit. Biosystems Engineering 85(4): 415-423. Bonetto, R. D., and J. L. Ladaga. 1998. The variogram method for characterization of scanning electron microscopy images. Scanning 20(6): 457-463. Brecht, J. K., R. L. Shewfelt, J. C. Garner, and E. W. Tollner. 1991. Using X-ray-computed tomography to nondestructively determine maturity of green tomatoes. HortScience 26(1): 45-47. Brosnan, T., and D. W. Sun. 2004. Improving quality inspection of food products by computer vision - a review. Journal of Food Engineering 61(1): 3-16. 191    Budiastra, I. W., Y. Ikeda, and T. Nishizu. 1998. Optical methods for quality evaluation of fruits (part 1)- optical properties of selected fruits using the Kubelka-Munk theory and their relationships with fruit maturity and sugar content. Journal of JSAM 60(20): 117-128. Buratti, S., A. Rizzolo, S. Benedetti, and D. Torreggiani. 2006. Electronic nose to detect strawberry aroma changes during osmotic dehydration. Journal of Food Science 71(4): E184-E189. Bykov, A. V., M. Y. Kirillin, A. V. Priezzhev, and R. Myllyla. 2006. Simulations of a spatially resolved reflectometry signal from a highly scattering three-layer medium applied to the problem of glucose sensing in human skin. Quantum Electronics 36(12): 1125-1130. Carlomagno, G., L. Capozzo, G. Attolico, and A. Distante. 2004. Non-destructive grading of peaches by near-infrared spectrometry. Infrared Physics & Technology 46(1-2): 23-29. Cen, H., and R. Lu. 2009. Quantification of the optical properties of two-layer turbid materials using a hyperspectral imaging-based spatially-resolved technique. Applied Optics 48(29): 5612-5623. Cen, H. Y., and Y. He. 2007. Theory and application of near infrared reflectance spectroscopy in determination of food quality. Trends in Food Science & Technology 18(2): 72-83. Chao, K., Y. R. Chen, H. Early, and B. Park. 1999. Color image classification systems for poultry viscera inspection. Applied Engineering in Agriculture 15(4): 363-369. Chen, C., J. Q. Lu, H. F. Ding, K. M. Jacobs, Y. Du, and X. H. Hu. 2006. A primary method for determination of optical parameters of turbid samples and application to intralipid between 550 and 1630nm. Optics Express 14(16): 7420-7435. Chen, H., and J. Debaerdemaeker. 1993. Effect of apple shape on acoustic measurements of firmness. Journal of Agricultural Engineering Research 56(3): 253-266. Cheong, W. F., S. A. Prahl, and A. J. Welch. 1990. A review of the optical properties of biological tissues. IEEE Journal of Quantum Electronics 26(12): 2166-2185. Chernomordik, V., D. W. Hattery, D. Grosenick, H. Wabnitz, H. Rinneberg, K. T. Moesta, P. M. Schlag, and A. Gandjbakhche. 2002. Quantification of optical properties of a breast tumor using random walk theory. Journal of Biomedical Optics 7(1): 80-87. Chinchuluun, R., W. S. Lee, and R. Ehsani. 2009. Machine vision system for determining citrus count and size on a canopy shake and catch harvester Applied Engineering in Agriculture 25(4): 451-458. Coquoz, O., L. O. Svaasand, and B. J. Tromberg. 2001. Optical property measurements of turbid media in a small-volume cuvette with frequency-domain photon migration. Applied Optics 40(34): 6281-6291. 192    Crisosto, C. H. 1994. Stone fruit maturity indices: a descriptive review. Postharvest News and Information 5(6): 65N-68N. Cubeddu, R., C. D'Andrea, A. Pifferi, P. Taroni, A. Torricelli, G. Valentini, C. Dover, D. Johnson, M. Ruiz-Altisent, and C. Valero. 2001a. Nondestructive quantification of chemical and physical properties of fruits by time-resolved reflectance spectroscopy in the wavelength range 650-1000 nm. Applied Optics 40(4): 538-543. Cubeddu, R., C. D'Andrea, A. Pifferi, P. Taroni, A. Torricelli, G. Valentini, M. Ruiz-Altisent, C. Valero, C. Ortiz, C. Dover, and D. Johnson. 2001b. Time-resolved reflectance spectroscopy applied to the nondestructive monitoring of the internal optical properties in apples. Applied Spectroscopy 55(10): 1368-1374. Cubeddu, R., A. Pifferi, and A. Torricelli. 2000. Measuring fresh fruit and vegetable quality: advanced optical methods. In Fruit and Vegetable Processing: Improving Quality, 150168. W. Jongen, ed. Boca Raton, FL, USA: CRC Press LLC. Cui, W. J., and L. E. Ostrander. 1992. The relationship of surface reflectance measurements to optical-properties of layered biological media. IEEE Transactions on Biomedical Engineering 39(2): 194-201. Dahm, D. J., and K. D. Dahm. 2001. The physics of near-infrared scattering. In Near-Infrared Technology in the agricultural and Food Industries, 1-17. P. Williams and K. Norris, eds. St. Paul, MN. Dastidar, S. G. 2006. Gompertz: A Scilab program for estimating Gompertz curve using GaussNewton method of least squares. Journal of Statistical Software 15(12). De Belie, N., S. Schotte, J. Lammertyn, B. Nicolai, and J. De Baerdemaeker. 2000. Firmness changes of pear fruit before and after harvest with the acoustic impulse response technique. Journal of Agricultural Engineering Research 77(2): 183-191. De Ketelaere, B., M. S. Howarth, L. Crezee, J. Lammertyn, K. Viaene, I. Bulens, and J. De Baerdemaeker. 2006. Postharvest firmness changes as measured by acoustic and lowmass impact devices: a comparison of techniques. Postharvest Biology and Technology 41(3): 275-284. Deulin, X., and J. P. L'Huillier. 2006. Finite element approach to photon propagation modeling in semi-infinite homogeneous and multilayered tissue structures. European Physical Journal Applied Physics 33(2): 133-146. Diezma-Iglesias, B., M. Ruiz-Altisent, and P. Barreiro. 2004. Detection of internal quality in seedless watermelon by acoustic impulse response. Biosystems Engineering 88(2): 221230. Diezma-Lglesias, B., C. Valero, F. J. Garcia-Ramos, and M. Ruiz-Altisent. 2006. Monitoring of firmness evolution of peaches during storage by combining acoustic and impact methods. Journal of Food Engineering 77(4): 926-935. 193    Dogariu, A., and J. Ellis. 2010. Volume Scattering in Random Media. In Handbook of Optics, Volume II - Physical Optics. M. Bass, and V. N. Mahajan, eds. Columbus, OH: McGrawHill. Doornbos, R. M. P., R. Lang, M. C. Aalders, F. W. Cross, and H. Sterenborg. 1999. The determination of in vivo human tissue optical properties and absolute chromophore concentrations using spatially resolved steady-state diffuse reflectance spectroscopy. Physics in Medicine and Biology 44(4): 967-981. Durduran, T., A. G. Yodh, B. Chance, and D. A. Boas. 1997. Does the photon-diffusion coefficient depend on absorption? Journal of the Optical Society of America a-Optics Image Science and Vision 14(12): 3358-3365. Edstrom, P. 2004. Comparison of the DORT2002 radiative transfer solution method and the Kubelka-Munk model. Nordic Pulp & Paper Research Journal 19(3): 397-403. ElMasry, G., N. Wang, A. ElSayed, and M. Ngadi. 2007. Hyperspectral imaging for nondestructive determination of some quality attributes for strawberry. Journal of Food Engineering 81(1): 98-107. Fabbri, F., M. A. Franceschini, and S. Fantini. 2003. Characterization of spatial and temporal variations in the optical properties of tissuelike media with diffuse reflectance imaging. Applied Optics 42(16): 3063-3072. Farrell, T. J., M. S. Patterson, and B. Wilson. 1992. A diffusion-theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical-properties in vivo. Medical Physics 19(4): 879-888. Galili, N., A. Mizrach, and G. Rosenhouse. 1993. Ultrasonic testing of whole fruit for nondestructive quality evaluation. In ASABE Paper No. 936026. St. Joseph, Mich.: ASABE. Giusto, A., R. Saija, M. A. Iati, P. Denti, F. Borghese, and O. I. Sindoni. 2003. Optical properties of high-density dispersions of particles: application to intralipid solutions. Applied Optics 42(21): 4375-4380. Gobin, L., L. Blanchot, and H. Saint-Jalmes. 1999. Integrating the digitized backscattered image to measure absorption and reduced-scattering coefficients in vivo. Applied Optics 38(19): 4217-4227. Gonzalez-Rodriguez, P., and A. D. Kim. 2008. Light propagation in two-layer tissues with an irregular interface. Journal of the Optical Society of America a-Optics Image Science and Vision 25(1): 64-73. Graaff, R., J. G. Aarnoudse, F. F. M. Demul, and H. W. Jentink. 1993. Similarity relations for anisotropic scattering in absorbing media. Optical Engineering 32(2): 244-252. 194    Groenhuis, R. A. J., H. A. Ferwerda, and J. J. Tenbosch. 1983a. Scattering and absorption of turbid materials determined from reflection measurements.1. theory. Applied Optics 22(16): 2456-2462. Groenhuis, R. A. J., J. J. Tenbosch, and H. A. Ferwerda. 1983b. Scattering and absorption of turbid materials determined from reflection measurements. 2. Measuring method and calibration. Applied Optics 22(16): 2463-2467. Gunasekaran, S. 1990. Delayed light-emission as a means of quality evaluation of fruits and vegetables. Critical Reviews in Food Science and Nutrition 29(1): 19-34. Haff, R. P., D. C. Slaughter, Y. Sarig, and A. Kader. 2006. X-ray assessment of translucency in pineapple. Journal of Food Processing and Preservation 30(5): 527-533. Hale, G. M., and M. R. Querry. 1973. Optical Constants of Water in 200 nm to 200 µm Wavelegnth region. Applied Optics 12(3): 555-563. Haskell, R. C., L. O. Svaasand, T. T. Tsay, T. C. Feng, and M. S. McAdams. 1994. Boundaryconditions for the diffusion equation in radiative-transfer. Journal of the Optical Society of America a-Optics Image Science and Vision 11(10): 2727-2741. Hollmann, J. L., and L. V. Wang. 2007. Multiple-source optical diffusion approximation for a multilayer scattering medium. Applied Optics 46(23): 6004-6009. Infante, R., M. Farcuh, and C. Meneses. 2008. Monitoring the sensorial quality and aroma through an electronic nose in peaches during cold storage. Journal of the Science of Food and Agriculture 88(12): 2073-2078. Ishimaru, A. 1978. Wave Propagation and Scattering in Random Media, Vol. 1: Single Scattering and Transport Theory. Academic Press, New York, NY. Jones, M. R., and Y. Yamada. 1998. Determination of the asymmetry parameter and scattering coefficient of turbid media from spatially resolved reflectance measurements. Optical Review 5(2): 72-76. Kacprzak, M., A. Liebert, P. Sawosz, N. Zolek, and R. Maniewski. 2007. Time-resolved optical imager for assessment of cerebral oxygenation. Journal of Biomedical Optics 12(3): 034019. Keener, J. D., K. J. Chalut, J. W. Pyhtila, and A. Wax. 2007. Application of Mie theory to determine the structure of spheroidal scatterers in biological materials. Optics Letters 32(10): 1326-1328. Kienle, A., T. Glanzmann, G. Wagnieres, and H. van den Bergh. 1998a. Investigation of twolayered turbid media with time-resolved reflectance. Applied Optics 37(28): 6852-6862. Kienle, A., L. Lilge, M. S. Patterson, R. Hibst, R. Steiner, and B. C. Wilson. 1996. Spatially resolved absolute diffuse reflectance measurements for noninvasive determination of the 195    optical scattering and absorption coefficients of biological tissue. Applied Optics 35(13): 2304-2314. Kienle, A., R. Michels, J. Schafer, and O. Fugger. 2007. Multiscale Description of Light Propagation in Biological Tissue. In Advances in Medical Engineering, 355-360. Springer Berlin Heidelberg. Kienle, A., and M. S. Patterson. 1997. Improved solutions of the steady-state and the timeresolved diffusion equations for reflectance from a semi-infinite turbid medium. Journal of the Optical Society of America a-Optics Image Science and Vision 14(1): 246-254. Kienle, A., M. S. Patterson, N. Dognitz, R. Bays, G. Wagnieres, and H. van den Bergh. 1998b. Noninvasive determination of the optical properties of two-layered turbid media. Applied Optics 37(4): 779-791. Kim, K. B., H. M. Jung, M. S. Kim, and G. S. Kim. 2004. Evaluation of fruit firmness by ultrasonic measurement. In Key Engineering Materials, 1049-1054. S. S. Lee, D. J. Yoon, J. H. Lee, and S. Lee, eds. Switzerland: Trans Tech Publications. Kim, N. K., and Y. C. Hung. 1990. Microscopic measurement of apple bruise. Food Structure 9(2): 97-104. Kim, S. M., P. Chen, M. J. McCarthy, and B. Zion. 1999. Fruit internal quality evaluation using on-line nuclear magnetic resonance sensors. Journal of Agricultural Engineering Research 74(3): 293-301. Koizumi, M., S. Naito, H. Kano, and T. Haishi. 2009. Examination of the Tissue Water in Cucumber Fruit by Small Dedicated Magnetic Resonance Imaging with a 1-T Permanent Magnet. Journal of the Japanese Society for Food Science and Technology-Nippon Shokuhin Kagaku Kogaku Kaishi 56(3): 146-154. Kolb, C. A., E. Wirth, W. M. Kaiser, A. Meister, M. Riederer, and E. E. Pfundel. 2006. Noninvasive evaluation of the degree of ripeness in grape berries (Vitis vinifera L. cv. Bacchus and Silvaner) by chlorophyll fluorescence. Journal of Agricultural and Food Chemistry 54(2): 299-305. Kramer, A., and B. A. Twigg. 1970. Quality Control for the Food Industry. AVI Publishing, Westport, Connecticut, USA. Lakowicz, J. R., G. Laczko, I. Gryczynski, H. Szmacinski, W. Wiczk, and M. L. Johnson. 1989. Frequency-domain fluorescence spectroscopy- principles, biochemical applications and future-developments. Berichte Der Bunsen-Gesellschaft-Physical Chemistry Chemical Physics 93(3): 316-327. Lammertyn, J., B. Nicolai, K. Ooms, V. De Smedt, and J. De Baerdemaeker. 1998. Nondestructive measurement of acidity, soluble solids, and firmness of Jonagold apples using NIR-spectroscopy. Transactions of the ASABE 41(4): 1089-1094. 196    Langenakens, J. J., X. Vandewalle, and J. DeBaerdemaeker. 1997. Influence of global shape and internal structure of tomatoes on the resonant frequency. Journal of Agricultural Engineering Research 66(1): 41-49. Langerholc, J. 1982. Beam broadening in dense scattering media. Applied Optics 21(9): 15931598. Lawrence, K. C., W. R. Windham, B. Park, and R. J. Buhr. 2003. A hyperspectral imaging system for identification of faecal and ingesta contamination on poultry carcasses. Journal of Near Infrared Spectroscopy 11(4): 269-281. Levenberg, K. 1944. A method for the solution of certain non-linear problems in least squares. Quarterly Journal of Applied Mathmatics II(2): 164-168. Lu, R. 2003. Detection of bruises on apples using near-infrared hyperspectral imaging. Transactions of the ASABE 46(2): 523-530. Lu, R., H. Cen, M. Huang, and D. P. Ariana. 2009. Optical properties of bruised apple tissue. Lu, R., and Y. R. Chen. 1998. Hyperspectral imaging system for safety inspection of food and agricultural products. In Proc. SPIE 3544, 121. Boston, MA. Lu, R., and Y. Peng. 2007. Development of a multispectral imaging prototype for real-time detection of apple fruit firmness. Optical Engineering 46(12): 123201. Lu, R., and Y. Peng. 2006. Hyperspectral scattering for assessing peach fruit firmness. Biosystems Engineering 93(2): 161-171. Lu, R. F., D. E. Guyer, and R. M. Beaudry. 2000. Determination of firmness and sugar content of apples using near-infrared diffuse reflectance. Journal of Texture Studies 31(6): 615-630. Lualdi, M., A. Colombo, B. Farina, S. Tomatis, and R. Marchesini. 2001. A phantom with tissuelike optical properties in the visible and near infrared for use in photomedicine. Lasers in Surgery and Medicine 28(3): 237-243. Madsen, S. J., P. Wyss, L. O. Svaasand, R. C. Haskell, Y. Tadir, and B. J. Tromberg. 1994. Determination of the optical-proeprties of the human uterus using frequency-domain photon migration and steady-state techniques. Physics in Medicine and Biology 39(8): 1191-1202. Marquardt, D. W. 1963. An algorithm for least-squares estimation of nolinear parameters. Journal of the Society for Industrial and Applied Mathematics 11(2): 431-441. Marquet, P., F. Bevilacqua, C. Depeursinge, and E. B. Dehaller. 1995. Determination of reduced scattering and absorption-coefficients by a single charge-coupled-device array measurement. 1. comparison between experiments and simulations. Optical Engineering 34(7): 2055-2063. 197    Martelli, F., and G. Zaccanti. 2007. Calibration of scattering and absorption properties of a liquid diffusive medium at NIR wavelengths. CW method. Optics Express 15(2): 486-500. McGlone, V. A., D. G. Fraser, R. B. Jordan, and R. Kunnemeyer. 2003. Internal quality assessment of mandarin fruit by vis/NIR spectroscopy. Journal of near Infrared Spectroscopy 11(5): 323-332. McGuire, R. G. 1992. Reporting of objective color measurements. HortScience 27(12): 12541255. Mendoza, F. A., R. Lu, D. P. Ariana, H. Cen, and B. B. Bailey. 2011. Integrated spectral and image analysis of hyperspectral scattering data for prediction of apple fruit firmness and soluble solids content. Postharvest Biology and Technology 62(2): 149-160. Michels, R., F. Foschum, and A. Kienle. 2008. Optical properties of fat emulsions. Optics Express 16(8): 5907-5925. Mishra, D. K., K. D. Dolan, and L. Yang. 2008. Confidence intervals for modeling anthocyanin retention in grape pomace during nonisothermal heating. Journal of Food Science 73(1): E9-E15. Mizrach, A. 2008. Ultrasonic technology for quality evaluation of fresh fruit and vegetables in pre- and postharvest processes. Postharvest Biology and Technology 48(3): 315-330. Mizrach, A., A. Bechar, Y. Grinshpon, A. Hofman, H. Egozi, and L. Rosenfeld. 2003. Ultrasonic classification of mealiness in apples. Transactions of the ASABE 46(2): 397-400. Mohsenin, N. N. 1984. Electromagnetic Radiation Properties of Foods and Agricultural Products. Gordon and Breach, Science Publishers, Inc., New York, NY. Motulsky, H., and A. Christopoulos. 2004. Fitting models to biological data using linear and nonlinear regression. Oxford University Press, New York. Mourant, J. R., T. Fuselier, J. Boyer, T. M. Johnson, and I. J. Bigio. 1997. Predictions and measurements of scattering and absorption over broad wavelength ranges in tissue phantoms. Applied Optics 36(4): 949-957. Nichols, M. G., E. L. Hull, and T. H. Foster. 1997. Design and testing of a white-light, steadystate diffuse reflectance spectrometer for determination of optical properties of highly scattering systems. Applied Optics 36(1): 93-104. Nicolai, B. M., E. Lotze, A. Peirs, N. Scheerlinck, and K. I. Theron. 2006. Non-destructive measurement of bitter pit in apple fruit using NIR hyperspectral imaging. Postharvest Biology and Technology 40(1): 1-6. Nicolai, B. M., B. E. Verlinden, M. Desmet, S. Saevels, W. Saeys, K. Theron, R. Cubeddu, A. Pifferi, and A. Torricelli. 2008. Time-resolved and continuous wave NIR reflectance 198    spectroscopy to predict soluble solids content and firmness of pear. Postharvest Biology and Technology 47(1): 68-74. Noh, H. K., and R. F. Lu. 2007. Hyperspectral laser-induced fluorescence imaging for assessing apple fruit quality. Postharvest Biology and Technology 43(2): 193-201. O'Connor, D. V., and D. Philip. 1984. Time-correlated Single Photon Counting. Academic Press, London. Palmer, G. M., and N. Ramanujam. 2007. Monte Carlo-based inverse model for calculating tissue optical properties. Part I: Theory and validation on synthetic phantoms. Applied Optics 46(27): 6847-6847. Park, B., and Y. R. Chen. 1994. Intensified multispectral imaging-system for poultry carcass inspection. Transactions of the ASABE 37(6): 1983-1988. Park, B., K. C. Lawrence, W. R. Windham, and R. J. Buhr. 2002. Hyperspectral imaging for detecting fecal and ingesta contaminants on poultry carcasses. Transactions of the ASABE 45(6): 2017-2026. Park, B., K. C. Lawrence, W. R. Windham, and D. P. Smith. 2004. Multispectral Imaging system for fecal and ingesta detection on poultry carcasses. Journal of Food Process Engineering 27(5): 311-327. Patterson, M. S., B. Chance, and B. C. Wilson. 1989. Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical-properties Applied Optics 28(12): 2331-2336. Patterson, M. S., J. D. Moulton, B. C. Wilson, K. W. Berndt, and J. R. Lakowicz. 1991. Frequency-domain reflectance for the determination of the scattering and absorption properties of tissue. Applied Optics 30(31): 4474-4476. Pedro, A. M. K., and M. M. C. Ferreira. 2005. Nondestructive determination of solids and carotenoids in tomato products by near-infrared spectroscopy and multivariate calibration. Analytical Chemistry 77(8): 2505-2511. Peng, Y., and R. Lu. 2006. An LCTF-based multispectral imaging system for estimation of apple fruit firmness: Part I. Acquisition and characterization of scattering images. Transactions of the ASABE 49(1): 259-267. Pham, T. H., F. Bevilacqua, T. Spott, J. S. Dam, B. J. Tromberg, and S. Andersson-Engels. 2000. Quantifying the absorption and reduced scattering coefficients of tissuelike turbid media over a broad spectral range with noncontact Fourier-transform hyperspectral imaging. Applied Optics 39(34): 6487-6497. Pifferi, A., A. Torricelli, A. Bassi, P. Taroni, R. Cubeddu, H. Wabnitz, D. Grosenick, M. Moller, R. Macdonald, J. Swartling, T. Svensson, S. Andersson-Engels, R. L. P. van Veen, H. Sterenborg, J. M. Tualle, H. L. Nghiem, S. Avrillier, M. Whelan, and H. Stamm. 2005. 199    Performance assessment of photon migration instruments: the MEDPHOT protocol. Applied Optics 44(11): 2104-2114. Pifferi, A., A. Torricelli, P. Taroni, D. Comelli, A. Bassi, and R. Cubeddu. 2007. Fully automated time domain spectrometer for the absorption and scattering characterization of diffusive media. Review of Scientific Instruments 78(5): 053103. Pilz, M., S. Honold, and A. Kienle. 2008. Determination of the optical properties of turbid media by measurements of the spatially resolved reflectance considering the point-spread function of the camera system. Journal of Biomedical Optics 13(5): 054047. Prahl, S. 1998. Optical absorption of water. http://omlc.ogi.edu/spectra/water/index.html. Prahl, S. A., M. J. C. Vangemert, and A. J. Welch. 1993. Determining the optical-properties of turbid media by using the adding-doubling method. Applied Optics 32(4): 559-568. Qin, J., and R. Lu. 2006. Hyperspectral diffuse reflectance imaging for rapid, noncontact measurement of the optical properties of turbid materials. Applied Optics 45(32): 83668373. Qin, J., and R. Lu. 2007. Measurement of the absorption and scattering properties of turbid liquid foods using hyperspectral imaging. Applied Spectroscopy 61(4): 388-396. Qin, J., and R. Lu. 2008. Measurement of the optical properties of fruits and vegetables using spatially resolved hyperspectral diffuse reflectance imaging technique. Postharvest Biology and Technology 49(3): 355-365. Qin, J., R. Lu, and Y. Peng. 2009. Prediction of apple internal quality using spectral absorption and scattering properties. Transactions of the ASABE 52(2): 499-507. Quevedo, R., L. G. Carlos, J. M. Aguilera, and L. Cadoche. 2002. Description of food surfaces and microstructural changes using fractal image texture analysis. Journal of Food Engineering 53(4): 361-371. Ren, T. Z., Z. Y. Yuan, and B. L. Su. 2007. Encapsulation of direct blue dye into mesoporous silica-based materials. Colloids and Surfaces a-Physicochemical and Engineering Aspects 300(1-2): 79-87. Reynolds, L., C. Johnson, and A. Ishimaru. 1976. Diffuse reflectance from a finite blood medium- applications to modeling of fiber optic catheters. Applied Optics 15(9): 20592067. Rudnitskaya, A., D. Kirsanov, A. Legin, K. Beullens, J. Lammertyn, B. M. Nicolai, and J. Irudayaraj. 2006. Analysis of apples varieties - comparison of electronic tongue with different analytical techniques. Sensors and Actuators B: Chemical 116(1-2): 23-28. 200    Saeys, W., M. A. Velazco-Roa, S. N. Thennadil, H. Ramon, and B. M. Nicolai. 2008. Optical properties of apple skin and flesh in the wavelength range from 350 to 2200 nm. Applied Optics 47(7): 908-919. Schatzki, T. F., R. P. Haff, R. Young, I. Can, L. C. Le, and N. Toyofuku. 1997. Defect detection in apples by means of x-ray imaging. Transactions of the ASABE 40(5): 1407-1415. Schmitt, J. M., G. X. Zhou, E. C. Walker, and R. T. Wall. 1990. Multilayer model of photon diffusion in skin. Journal of the Optical Society of America A 7(11): 2141-2153. Schweiger, M., and S. R. Arridge. 1997. The finite-element method for the propagation of light in scattering media: Frequency domain case. Medical Physics 24(6): 895-902. Seiden, P., R. Bro, L. Poll, and L. Munck. 1996. Exploring fluorescence spectra of apple juice and their connection to quality parameters by chemometrics. Journal of Agricultural and Food Chemistry 44(10): 3202-3205. Seo, I., J. S. You, C. K. Hayakawa, and V. Venugopalan. 2007. Perturbation and differential Monte Carlo methods for measurement of optical properties in a layered epithelial tissue model. Journal of Biomedical Optics 12(1): 014030. Shampine, L. F. 2008. Vectorized adaptive quadrature in MATLAB. Journal of Computational and Applied Mathematics 211(2): 131-140. Singh, N., and M. J. Delwiche. 1994. Machine vision methods for defect sorting stonefruit. Transactions of the ASABE 37(6): 1989-1997. Sirisomboon, P., M. Tanaka, S. Fujita, and T. Kojima. 2007. Evaluation of pectin constituents of Japanese pear by near infrared spectroscopy. Journal of Food Engineering 78(2): 701707. Spichtig, S., R. Hornung, D. W. Brown, D. Haensse, and M. Wolf. 2009. Multifrequency frequency-domain spectrometer for tissue analysis. Review of Scientific Instruments 80(2): 024301. Svensson, T., J. Swartling, P. Taroni, A. Torricelli, P. Lindblom, C. Ingvar, and S. AnderssonEngels. 2005. Characterization of normal breast tissue heterogeneity using time-resolved near-infrared spectroscopy. Physics in Medicine and Biology 50(11): 2559-2571. Swartling, J., J. S. Dam, and S. Andersson-Engels. 2003a. Comparison of spatially and temporally resolved diffuse-reflectance measurement systems for determination of biomedical optical properties. Applied Optics 42(22): 4612-4620. Swartling, J., A. Pifferi, A. M. K. Enejder, and S. Andersson-Engels. 2003b. Accelerated Monte Carlo models to simulate fluorescence spectra from layered tissues. Journal of the Optical Society of America a-Optics Image Science and Vision 20(4): 714-727. 201    Taktak, R., J. V. Beck, and E. P. Scott. 1993. Optimal experimental -design for estimating thermal-properties of composite -materials. International Journal of Heat and Mass Transfer 36(12): 2977-2986. Tarantola, A. 2005. Inverse Problem Theory and Methods for Model Parameter Estimation. Society for Industrial and Applied Mathematics, Philadelphia, PA. Thomas, F. C., and Y. Y. Li. 1994. On the convergence of interiro-reflective newton methods for nonlinear minimization subject to bounds. Mathematical Programming 67(2): 189-224. Tijskens, L. M. M., P. E. Zerbini, R. E. Schouten, M. Vanoli, S. Jacob, M. Grassi, R. Cubeddu, L. Spinelli, and A. Torricelli. 2007. Assessing harvest maturity in nectarines. Postharvest Biology and Technology 45(2): 204-213. Tran, C. D. 2000. Visualising chemical composition and reaction kinetics by the near infrared multispectral imaging technique. Journal of near Infrared Spectroscopy 8(2): 87-99. Tromberg, B. J., O. Coquoz, J. Fishkin, T. Pham, E. R. Anderson, J. Butler, M. Cahn, J. D. Gross, V. Venugopalan, and D. Pham. 1997. Non-invasive measurements of breast tissue optical properties using frequency-domain photon migration. Philosophical Transactions of the Royal Society of London Series B-Biological Sciences 352(1354): 661-668. Tromberg, B. J., N. Shah, R. Lanning, A. Cerussi, J. Espinoza, T. Pham, L. Svaasand, and J. Butler. 2000. Non-invasive in vivo characterization of breast tumors using photon migration spectroscopy. Neoplasia 2(1-2): 26-40. Tromberg, B. J., L. O. Svaasand, M. K. Fehr, S. J. Madsen, P. Wyss, B. Sansone, and Y. Tadir. 1996. A mathematical model for light dosimetry in photodynamic destruction of human endometrium. Physics in Medicine and Biology 41(2): 223-237. Tseng, S. H., A. Grant, and A. J. Durkin. 2008. In vivo determination of skin near-infrared optical properties using diffuse optical spectroscopy. Journal of Biomedical Optics 13(1): 014016. Tuchin, V. 2000. Tissue Optics: Light Scattering Methods and Instruments for Medical Diagnosis. SPIE PRESS, Bellingham, Washington, USA. Tuchin, V. V., S. R. Utz, and I. V. Yaroslavsky. 1994. Tissue optics, light-distribution, and spectroscopy. Optical Engineering 33(10): 3178-3188. Umbaugh, S. E. 2005. Computer Imaging: Digital Image Analysis and Processing. CRC Press, Boca Raton, FL. Valero, C., M. Ruiz-Altisent, R. Cubeddu, A. Pifferi, P. Taroni, A. Torricelli, G. Valentini, D. S. Johnson, and C. J. Dover. 2004. Detection of internal quality in kiwi with time-domain diffuse reflectance spectroscopy. Applied Engineering in Agriculture 20(2): 223-230. 202    van Staveren, H. J., C. J. M. Moes, J. Vanmarle, S. A. Prahl, and M. J. C. Vangemert. 1991. Light-scattering in intralipid-10% in the wavelength range of 400-1100 nm. Applied Optics 30(31): 4507-4514. Vanstaveren, H. J., C. J. M. Moes, J. Vanmarle, S. A. Prahl, and M. J. C. Vangemert. 1991. Light-scattering in intralipid-10% in the wavelength range of 400-1100 nm. Applied Optics 30(31): 4507-4514. Vo-Dinh, T. 2003. Biomedical Photonics Handbook. CRC Press LLC, Boca Raton, FL. Wang, J., A. H. Gomez, and A. G. Pereira. 2006. Acoustic impulse response for measuring the firmness of mandarin during storage. Journal of Food Quality 29(4): 392-404. Wang, L. H., S. L. Jacques, and L. Q. Zheng. 1997. CONV - convolution for responses to a finite diameter photon beam incident on multi-layered tissues. Computer Methods and Programs in Biomedicine 54(3): 141-150. Wang, L. H., S. L. Jacques, and L. Q. Zheng. 1995a. MCML - Monte-carlo modeling of light transport in multilayered tissues. Computer Methods and Programs in Biomedicine 47(2): 131-146. Wang, L. H., S. L. Jacques, and L. Q. Zheng. 1995b. MCML - Monte Carlo modeling of light transport in multilayered tissues. Computer Methods and Programs in Biomedicine 47(2): 131-146. Welch, A. J., and M. J. C. van Gemert. 1995. Optical-Thermal Response of Laser-Irradiated Tissue. Plenum Press, New York, NY. Welzel, J., C. Reinhardt, E. Lankenau, C. Winter, and H. H. Wolff. 2004. Changes in function and morphology of normal human skin: evaluation using optical coherence tomography. British Journal of Dermatology 150(2): 220-225. Wilkinson, G. 1961. Statistical estimations in enzyme kinetics. Biochemical Journal 80(2): 324332. Williams, P., and K. Norris. 2001. Near-infrared Technology in the Agricultural and Food Industries. American Association of Cereal Chemists, St. Paul, Minnesota. Wills, R., B. McGlasson, D. Graham, and D. Joyce. 2004. Postharvest: An Introduction to the Physiology & Handling of Fruit, Vegetables & Ornamentals. 4th ed. Hyde park Press, Adelaide, South Australia. Wilson, B. C., and M. S. Patterson. 2008. The physics, biophysics and technology of photodynamic therapy. Physics in Medicine and Biology 53(9): R61-R109. Wood, D. F., S. H. Imam, W. J. Orts, and G. M. Glenn. 2009. Fresh fruit-microstructure, texture and quality. In Proc. SPIE 7378, M. T. Postek, D. E. Newbury, S. F. Platek, and D. C. Joy, eds. Monterey, CA. 203    Xu, H. P., and M. S. Patterson. 2006. Determination of the optical properties of tissue-simulating phantoms from interstitial frequency domain measurements of relative fluence and phase difference. Optics Express 14(14): 6485-6501. Zaccanti, G., S. Del Bianco, and F. Martelli. 2003. Measurements of optical properties of highdensity media. Applied Optics 42(19): 4023-4030. Zolek, N. S., S. Wojtkiewicz, and A. Liebert. 2008. Correction of anisotropy coefficient in original Henyey Greenstein phase function for Monte Carlo simulations of light transport in tissue. Biocybernetics and Biomedical Engineering 28(4): 59-73. Zoltan, K., K. D. Balazs, and F. Andras. 2008. Qualitative analysis of fruit juice by electronic tongue. Elelmiszervizsgalati Kozlemenyek 54(3): 151-162. 204