Some Matlab Task With Biomedical Science. Hello,These are some of matlab tasks and related to bio medical science as well.I have attached the assignment and the text book as well.you can check before applying for the job.I need two soluiton of this assignment with four questions for two students.solutions should be litltle differnt so professor cant judget is done by same person.Our deadline is 2 days .
Let me know if any body is sure to do it. Please no time wasters.
B I C E N T C N N I A I
|| 1 8 O 7 r|®WILEY ÜJ2 0 0 7
B I C E N T E N N I A L
T H E WILEY BICENTENNIAL-KNOWLEDGE FOR GENERATIONS
//~)ach generation has its unique needs and aspirations. When Charles Wiley first opened his small printing shop in lower Manhattan in 1807, it was a generation of boundless potential searching for an identity. And we were there, helping to define a new American literary tradition. Over half a century later, in the midst of the Second Industrial Revolution, it was a generation focused on building the future. Once again, we were there, supplying the critical scientific, technical, and engineering knowledge that helped frame the world. Throughout the 20th Century, and into the new millennium, nations began to reach out beyond their own borders and a new international community was born. Wiley was there, expanding its operations around the world to enable a global exchange of ideas, opinions, and know-how.
For 200 years, Wiley has been an integral part of each generation’s journey, enabling the flow of information and understanding necessary to meet their needs and fulfill their aspirations. Today, bold new technologies are changing the way we live and learn. Wiley will be there, providing you the must-have knowledge you need to imagine new worlds, new possibilities, and new opportunities.
Generations come and go, but you can always count on Wiley to provide you the knowledge you need, when and where you need it!
iLLu^M. (2, WILLIAM J . P E S C E
PRESIDENT A N D CHIEF EXECUTIVE OFFICER PETER B O O T H WILEY
CHAIRMAN DF T H E BOARD
BIOMEDICAL OPTICS
PRINCIPLES AND IMAGING
Lihong V. Wang
Hsin-i Wu
B I C E N T E N N I A L ·
; | 1 8 0 7
©WILEY 2 OO 7
B I C E N T E N N I A L
WILEY-INTERSCIENCE A John Wiley & Sons, Inc., Publication
Copyright © 2007 by John Wiley & Sons, Inc. All rights reserved.
Published by John Wiley & Sons, Inc., Hoboken, New Jersey. Published simultaneously in Canada.
No part of this publication may be reproduced, stored in a retrieval system, or transmitted in any form or by any means, electronic, mechanical, photocopying, recording, scanning, or otherwise, except as permitted under Section 107 or 108 of the 1976 United States Copyright Act, without either the prior written permission of the Publisher, or authorization through payment of the appropriate per-copy fee to the Copyright Clearance Center, Inc., 222 Rosewood Drive, Danvers, MA 01923, (978) 750-8400, fax (978) 750-4470, or on the web at www.copyright.com. Requests to the Publisher for permission should be addressed to the Permissions Department, John Wiley & Sons, Inc., 111 River Street, Hoboken, NJ 07030, (201) 748-6011, fax (201) 748-6008, or online at http://www.wiley.com/go/permission.
Limit of Liability/Disclaimer of Warranty: While the publisher and author have used their best efforts in preparing this book, they make no representations or warranties with respect to the accuracy or completeness of the contents of this book and specifically disclaim any implied warranties of merchantability or fitness for a particular purpose. No warranty may be created or extended by sales representatives or written sales materials. The advice and strategies contained herein may not be suitable for your situation. You should consult with a professional where appropriate. Neither the publisher nor author shall be liable for any loss of profit or any other commercial damages, including but not limited to special, incidental, consequential, or other damages.
For general information on our other products and services or for technical support, please contact our Customer Care Department within the United States at (800) 762-2974, outside the United States at (317) 572-3993 or fax (317) 572-4002.
Wiley also publishes its books in a variety of electronic formats. Some content that appears in print may not be available in electronic formats. For more information about Wiley products, visit our web site at www.wiley.com.
Wiley Bicentennial Logo: Richard J. Pacifico
Library of Congress Cataloging-in-Publication Data:
Wang, Lihong V. Biomedical optics : principles and imaging / Lihong V. Wang, Hsin-i
Wu. p. ; cm.
Includes bibliographical references and index. ISBN: 978-0-471-74304-0 (cloth) 1. Imaging systems in medicine. 2. Lasers in medicine. 3. Optical
detectors. I. Wu, Hsin-i. II. Title. [DNLM: 1. Optics. 2. Diagnostic Imaging—methods. 3. Light.
4. Models, Theoretical. 5. Tomography, Optical—methods. WB 117 W246b 2007] R857.06W36. 2007 616.0754-dc22
2006030754
Printed in the United States of America.
109 8 7
To our families, mentors, students, and friends
CONTENTS
Preface xiii
1. Introduction 1
1.1. Motivation for Optical Imaging 1
1.2. General Behavior of Light in Biological Tissue 2
1.3. Basic Physics of Light-Matter Interaction 3
1.4. Absorption and its Biological Origins 5
1.5. Scattering and its Biological Origins 7
1.6. Polarization and its Biological Origins 9
1.7. Fluorescence and its Biological Origins 9
1.8. Image Characterization 10
Problems 14
Reading 15
Further Reading 15
2. Rayleigh Theory and Mie Theory for a Single Scatterer 17
2.1. Introduction 17
2.2. Summary of Rayleigh Theory 17
2.3. Numerical Example of Rayleigh Theory 19
2.4. Summary of Mie Theory 20
2.5. Numerical Example of Mie Theory 21
Appendix 2A. Derivation of Rayleigh Theory 23
Appendix 2B. Derivation of Mie Theory 26
Problems 34
Reading 35
Further Reading 35
vii
VIII CONTENTS
Monte Carlo Modeling of Photon Transport in Biological Tissue
3.1. Introduction
3.2. Monte Carlo Method
3.3. Definition of Problem
3.4. Propagation of Photons
3.5. Physical Quantities
3.6. Computational Examples
Appendix 3A. Summary of MCML
Appendix 3B. Probability Density Function
Problems
Reading
Further Reading
Convolution for Broadbeam Responses
4.1. Introduction
4.2. General Formulation of Convolution
4.3. Convolution over a Gaussian Beam
4.4. Convolution over a Top-Hat Beam
4.5. Numerical Solution to Convolution
4.6, Computational Examples
Appendix 4A. Summary of CONV
Problems
Reading
Further Reading
Radiative Transfer Equation and Diffusion Theory
5.1. Introduction
5.2. Definitions of Physical Quantities
5.3. Derivation of Radiative Transport Equation
5.4. Diffusion Theory
5.5. Boundary Conditions
5.6. Diffuse Reflectance
37
37
37
38
39
50
55
58
60
60
62
62
67
67
67
69
71
72
77
77
80
81
81
83
83
83
85
88
101
106
CONTENTS IX
5.7. Photon Propagation Regimes 114
Problems 116
Reading 117
Further Reading 118
Hybrid Model of Monte Carlo Method and Diffusion Theory 119
6.1. Introduction 119
6.2. Definition of Problem 119
6.3. Diffusion Theory 119
6.4. Hybrid Model 122
6.5. Numerical Computation 124
6.6. Computational Examples 125
Problems 132
Reading 133
Further Reading 133
Sensing of Optical Properties and Spectroscopy 135
7.1. Introduction 135
7.2. Collimated Transmission Method 135
7.3. Spectrophotometry 139
7.4. Oblique-Incidence Reflectometry 140
7.5. White-Light Spectroscopy 144
7.6. Time-Resolved Measurement 145
7.7. Fluorescence Spectroscopy 146
7.8. Fluorescence Modeling 147
Problems 148
Reading 149
Further Reading 149
Ballistic Imaging and Microscopy 153
8.1. Introduction 153
8.2. Characteristics of Ballistic Light 153
X CONTENTS
8.3. Time-Gated Imaging 154
8.4. Spatiofrequency-Filtered Imaging 156
8.5. Polarization-Difference Imaging 157
8.6. Coherence-Gated Holographic Imaging 158
8.7. Optical Heterodyne Imaging 160
8.8. Radon Transformation and Computed Tomography 163
8.9. Confocal Microscopy 164
8.10. Two-Photon Microscopy 169
Appendix 8A. Holography 171
Problems 175
Reading 177
Further Reading 177
9. Optical Coherence Tomography 181
9.1. Introduction 181
9.2. Michelson Interferometry 181
9.3. Coherence Length and Coherence Time 184
9.4. Time-Domain OCT 185
9.5. Fourier-Domain Rapid-Scanning Optical Delay Line 195
9.6. Fourier-Domain OCT 198
9.7. Doppler OCT 206
9.8. Group Velocity Dispersion 207
9.9. Monte Carlo Modeling of OCT 210
Problems 213
Reading 215
Further Reading 215
10. Mueller Optical Coherence Tomography 219
10.1. Introduction 219
10.2. Mueller Calculus versus Jones Calculus 219
10.3. Polarization State 219
10.4. Stokes Vector 222
CONTENTS XI
10.5. Mueller Matrix 224
10.6. Mueller Matrices for a Rotator, a Polarizer, and a Retarder 225
10.7. Measurement of Mueller Matrix 227
10.8. Jones Vector 229
10.9. Jones Matrix 230
10.10. Jones Matrices for a Rotator, a Polarizer, and a Retarder 230
10.11. Eigenvectors and Eigenvalues of Jones Matrix 231
10.12. Conversion from Jones Calculus to Mueller Calculus 235
10.13. Degree of Polarization in OCT 236
10.14. Serial Mueller OCT 237
10.15. Parallel Mueller OCT 237
Problems 243
Reading 244
Further Reading 245
11. Diffuse Optical Tomography 249
11.1. Introduction 249
11.2. Modes of Diffuse Optical Tomography 249
11.3. Time-Domain System 251
11.4. Direct-Current System 252
11.5. Frequency-Domain System 253
11.6. Frequency-Domain Theory: Basics 256
11.7. Frequency-Domain Theory: Linear Image Reconstruction 261
11.8. Frequency-Domain Theory: General Image Reconstruction 267
Appendix 11 A. ART and SIRT 275
Problems 276
Reading 279
Further Reading 279
12. Photoacoustic Tomography 283
12.1. Introduction 283
12.2. Motivation for Photoacoustic Tomography 283
XII CONTENTS
12.3. Initial Photoacoustic Pressure 284
12.4. General Photoacoustic Equation 287
12.5. General Forward Solution 288
12.6. Delta-Pulse Excitation of a Slab 293
12.7. Delta-Pulse Excitation of a Sphere 297
12.8. Finite-Duration Pulse Excitation of a Thin Slab 302
12.9. Finite-Duration Pulse Excitation of a Small Sphere 303
12.10. Dark-Field Confocal Photoacoustic Microscopy 303
12.11. Synthetic Aperture Image Reconstruction 307
12.12. General Image Reconstruction 309
Appendix 12A. Derivation of Acoustic Wave Equation 313
Appendix 12B. Green Function Approach 316
Problems 317
Reading 319
Further Reading 319
13. Ultrasound-Modulated Optical Tomography 323
13.1. Introduction 323
13.2. Mechanisms of Ultrasonic Modulation of Coherent Light 323
13.3. Time-Resolved Frequency-Swept UOT 326
13.4. Frequency-Swept UOT with Parallel-Speckle Detection 329
13.5. Ultrasonically Modulated Virtual Optical Source 331
13.6. Reconstruction-Based UOT 332
13.7. UOT with Fabry-Perot Interferometry 335
Problems 338
Reading 339
Further Reading 339
Appendix A. Definitions of Optical Properties 343
Appendix B. List of Acronyms 345
Index 347
PREFACE
Biomedical optics is a rapidly growing area of research. Although many universi- ties have begun to offer courses on the topic, a textbook containing examples and homework problems has not been available. The need to fill this void prompted us to write this book.
This book is based on our lecture notes for a one-semester (45 lecture hours) entry-level course, which we have taught since 1998. The contents are divided into two major parts: (1) fundamentals of photon transport in biological tissue and (2) optical imaging. In the first part (Chapters 1 -7) , we start with a brief introduc- tion to biomedical optics and then cover single-scatterer theories, Monte Carlo modeling of photon transport, convolution for broadbeam responses, radiative transfer equation and diffusion theory, hybrid Monte Carlo method and diffusion theory, and sensing of optical properties and spectroscopy. In the second part (Chapters 8-13), we cover ballistic imaging, optical coherence tomography, dif- fuse optical tomography, photoacoustic tomography, and ultrasound-modulated optical tomography.
When the book is used as the textbook in a course, the instructor may request a solution manual containing homework solutions from the publisher. To ben- efit from this text, students are expected to have a background in calculus
® and differential equations. Experience in MATLAB or C/C++ is also helpful. Source codes and other information can be found at ftp://ftp.wiley.com/public/ scLtech_med/biomedical_optics.
Although our multilayered Monte Carlo model is in the public domain, we have found that students are able to better grasp the concept of photon transport in biological tissue when they implement simple semiinfinite versions of the model. For this reason, we encourage the use of simulations whenever appropriate.
Because a great deal of material beyond our original lecture notes has been added, two semesters are recommended to cover the complete textbook. Alterna- tively, selected chapters can be covered in a one-semester course. In addition to serving as a textbook, this book can also be used as a reference for professionals and a supplement for trainees engaged in short courses in the field of biomedical optics.
We are grateful to Mary Ann Dickson for editing the text and to Elizabeth Smith for redrawing the figures. We appreciate Sancy Wu’s close reading of
xiii
XIV PREFACE
the manuscript. We are also thankful to the many students who contributed to the homework solutions. Finally, we wish to thank our students Li Li, Manojit Pramanik, and Sava Sakadzic for proofreading the book.
LIHONG V. WANG, PH.D. HSIN-I Wu, PH.D.
CHAPTER 1
Introduction
1.1. MOTIVATION FOR OPTICAL IMAGING
The most common medical imaging modalities include X-ray radiography, ultra- sound imaging (ultrasonography), X-ray computed tomography (CT), and mag- netic resonance imaging (MRI). The discovery of X rays in 1895, for which Roentgen received the first Nobel Prize in Physics in 1901, marked the advent of medical imaging. Ultrasonography, which is based on sonar, was introduced into medicine in the 1940s after World War II. The invention of CT in the 1970s, for which Cormack and Hounsfield received the Nobel Prize in Medicine in 1979, initiated digital cross-sectional imaging (tomography). The invention of MRI, also in the 1970s, for which Lauterbur and Mansfield received the Nobel Prize in Medicine in 2003, enabled functional imaging with high spatial resolution. Optical imaging, which is compared with the other modalities in Table 1.1, is currently emerging as a promising new addition to medical imaging.
Reasons for optical imaging of biological tissue include
1. Optical photons provide nonionizing and safe radiation for medical appli- cations.
2. Optical spectra—based on absorption, fluorescence, or Raman scatter- ing—provide biochemical information because they are related to molec- ular conformation.
3. Optical absorption, in particular, reveals angiogenesis and hyperme- tabolism, both of which are hallmarks of cancer; the former is related to the concentration of hemoglobin and the latter, to the oxygen satura- tion of hemoglobin. Therefore, optical absorption provides contrast for functional imaging.
4. Optical scattering spectra provide information about the size distribution of optical scatterers, such as cell nuclei.
5. Optical polarization provides information about structurally anisotropic tissue components, such as collagen and muscle fiber.
Biomedical Optics: Principles and Imaging, by Lihong V. Wang and Hsin-i Wu Copyright © 2007 John Wiley & Sons, Inc.
1
INTRODUCTION
TABLE 1.1. Comparison of Various Medical Imaging Modalities
Characteristics
Soft-tissue contrast Spatial resolution Maximum imaging depth Function Nonionizing radiation Data acquisition Cost
X-ray Imaging
Poor Excellent Excellent None No Fast Low
Ultrasonography
Good Good Good Good Yes Fast Low
MRI
Excellent Good Excellent Excellent Yes Slow High
Optical Imaging
Excellent Mixed” Good Excellent Yes Fast Low
“High in ballistic imaging (see Chapters 8-10) and photoacoustic tomography (see Chapter 12); low in diffuse optical tomography (see Chapter 11).
6. Optical frequency shifts due to the optical Doppler effect provide infor- mation about blood flow.
7. Optical properties of targeted contrast agents provide contrast for the molecular imaging of biomarkers.
8. Optical properties or bioluminescence of products from gene expression provide contrast for the molecular imaging of gene activities.
9. Optical spectroscopy permits simultaneous detection of multiple contrast agents.
10. Optical transparency in the eye provides a unique opportunity for high- resolution imaging of the retina.
1.2. GENERAL BEHAVIOR OF LIGHT IN BIOLOGICAL TISSUE
Most biological tissues are characterized by strong optical scattering and hence are referred to as either scattering media or turbid media. By contrast, optical absorption is weak in the 400-1350-nm spectral region. The mean free path between photon scattering events is on the order of 0.1 mm, whereas the mean absorption length (mean path length before photon absorption) can extend to 10-100 mm.
Photon propagation in biological tissue is illustrated in Figure 1.1. The light source is spatially a pencil beam (an infinitely narrow collimated beam) and temporally a Dirac delta pulse. The optical properties (see Appendix A) of the tissue include the following: refractive index n — 1.37, absorption coefficient [ia =z 1.4 cm- 1 , scattering coefficient [is = 350 cm”1, and scattering anisotropy g = 0.8. The mean free path equals 28 μιτι, corresponding to a propagation time of 0.13 ps. The transport mean free path equals 140 μηι, corresponding to a propagation time of 0.64 ps. Note how widely the photons spread versus time in relation to the two time constants mentioned above. This diffusion-like behavior of light in biological tissue presents a key challenge for optical imaging. Various techniques have been designed to meet this challenge.
BASIC PHYSICS OF LIGHT-MATTER INTERACTION 3
Air I Laser beam
Tissue 1 x
ΙΟΟμπι/div ~ 1 ί ί ί ί 1 1 1 1 i
Geometry
0.15 ps
0.55 ps
1.55 ps
0.05 ps
0.35 ps
1.05 ps ‘.,, _ JULU ,
/
2.05 ps
Figure 1.1. Snapshots of the simulated photon density distribution in a piece of biological tissue projected along the y axis, which points out of the paper.
1.3. BASIC PHYSICS OF LIGHT-MATTER INTERACTION
Absorption of a photon can elevate an electron of a molecule from the ground state to an excited state, which is termed excitation. Excitation can also be caused by other mechanisms, which are either mechanical (frictional) or chemical in nature. When an electron is raised to an excited state, there are several possi- ble outcomes. The excited electron may relax to the ground state and give off luminescence (another photon) or heat. If another photon is produced, the emis- sion process is referred to as fluorescence or phosphorescence, depending on the lifetime of the excited electron; otherwise, it is referred to as nonradiative relax- ation. Lifetime is defined as the average time that an excited molecule spends in the excited state before returning to the ground state. The ratio of the number of photons emitted to the number of photons absorbed is referred to as the quantum yield of fluorescence. If the excited molecule is near another molecule with a sim- ilar electronic configuration, the energy may be transferred by excitation energy transfer—the excited electron in one molecule drops to the ground state while the energy is transferred to the neighboring molecule, raising an electron in that molecule to an excited state with a longer lifetime. Another possible outcome is photochemistry, in which an excited electron is actually transferred to another
4 INTRODUCTION
Excited state
hvA w +
Absorption
Internal conversion
Nonradiative relaxation
I
Vibrational energy levels
Intersystem crossing Excited triplet state
Fluorescence
Virtual state N
Phosphorescence
hvA hvR
Raman
Ground state
Figure 1.2. Jablonski energy diagram showing excitation and various possible relaxation mechanisms. Each hv denotes the photon energy, where subscripts A, F, P, and R denote absorption, fluorescence, phosphorescence, and Raman scattering, respectively.
molecule. This type of electron transfer alters the chemical properties of both the electron donor and the electron acceptor, as in photosynthesis.
A Jablonski energy diagram describing electronic transitions between ground states and excited states is shown in Figure 1.2. Molecules can absorb photons that match the energy difference between two of their discrete energy levels, provided the transitions are allowed. These energy levels define the absorption and the emission bands.
Fluorescence involves three events with vastly different timescales. Excita- tion by a photon takes place in femtoseconds (1 fs = 10-15 s, about one optical period). Vibrational relaxation (also referred to as internal conversion) of an excited-state electron to the lowest vibrational energy level in the excited state lasts for picoseconds (1 ps = 10~12 s) and does not result in emission of a photon (nonradiative). Fluorescence emission lingers over nanoseconds (1 ns = 10~9 s). Consequently, fluorescence lifetime is on the order of a nanosecond.
Phosphorescence is similar to fluorescence, but the excited molecule further transitions to a metastable state by intersystem crossing, which alters the electron spin. Because relaxation from the metastable state to the ground state is spin- forbidden, emission occurs only when thermal energy raises the electron to a state where relaxation is allowed. Consequently, phosphorescence depends on temperature and has a long lifetime (milliseconds or longer).
Two types of photon scattering by a molecule exist: elastic and inelastic (or Raman) scattering. The former involves no energy exchange between the photon and the molecule, whereas the latter does. Although both Raman scattering and fluorescence alter the optical wavelength, they have different mechanisms. In
ABSORPTION AND ITS BIOLOGICAL ORIGINS 5
Raman scattering, the molecule is excited to a virtual state; in fluorescence, the molecule is excited to a real stationary state. In both cases, the excited molecule relaxes to an energy level of the ground state and emits a photon. The molecule may either gain energy from, or lose energy to, the photon. If the molecule gains energy, the transition is known as a Stokes transition. Otherwise, the transition is known as an anti-Stokes transition. The scattered photon shifts its frequency accordingly since the total energy is conserved. Raman scattering can reveal the specific chemical composition and molecular structure of biological tissue, whereas elastic scattering can reveal the size distribution of the scatterers.
1.4. ABSORPTION AND ITS BIOLOGICAL ORIGINS
The absorption coefficient μα is defined as the probability of photon absorp- tion in a medium per unit path length (strictly speaking, per unit infinitesimal path length). It has a representative value of 0.1 cm – 1 in biological tissue. The reciprocal of ia is referred to as the mean absorption length.
For a single absorber, the absorption cross section oa, which indicates the absorbing capability, is related to its geometric cross-sectional area og through the absorption efficiency Qa : σα = Qaog. In a medium containing many absorbers with number density Na, the absorption coefficient can be considered as the total cross-sectional area for absorption per unit volume:
ν>α = Νασα. (1.1)
Here, absorption by different absorbers is considered to be independent. According to the definition of the absorption coefficient, light attenuates as it
propagates in an absorbing-only medium according to the following equation:
— = – μ β Λ χ , (1.2)
where / denotes the light intensity and x denotes the distance along the light propagation direction. Equation (1.2) means that the percentage of light being absorbed in interval (JC, X -f- dx) is proportional to the product of ia and dx; the negative sign is due to the decrease of / as x increases. Integrating Eq. (1.2) leads to the well-known Beer law
I(x) = / 0 βχρ( -μ β χ ) , (1.3)
where IQ is the light intensity at x — 0. Beer’s law actually holds even for a tortuous path. The transmittance is defined by
rw = T ,̂ (1.4)
which represents the probability of survival after propagation over x.
6 INTRODUCTION
Wavelength (nm)
Figure 1.3. Molar extinction coefficients of oxy- and deoxyhemoglobin versus wave- length.
Optical absorption in biological tissue originates primarily from hemoglobin, melanin, and water. Hemoglobin has two forms: oxygenated and deoxygenated. Figure 1.3 shows the molar extinction coefficients—the extinction coefficient divided by ln(10) (see Section 7.3) per unit molar concentration—of oxy- and deoxyhemoglobin as a function of wavelength, where the extinction coefficient is defined as the probability of photon interaction with a medium per unit path length. Although extinction includes both absorption and scattering, absorption dominates scattering in hemoglobin. The molar extinction spectra of oxy- and deoxyhemoglobin are distinct but share a few intersections, termed isosbestic points. At these points, the absorption coefficient of an oxy- and deoxyhe- moglobin mixture depends only on the total concentration, regardless of the oxygen saturation.
The absorption coefficients of some primary absorbing biological tissue com- ponents are plotted as a function of wavelength in Figure 1.4. Melanin absorbs ultraviolet (UV) light strongly but longer-wavelength light less strongly. Even water can be highly absorbing in some spectral regions. At the 2.95-μπι water absorption peak, the penetration depth is less than 1 μτη since xa — 12,694 cm- 1 .
The absorption coefficients of biological tissue at two wavelengths can be used to estimate the concentrations of the two forms of hemoglobin based on the following equations:
μα(λ,) = 1η(10)εοχ(λ1)Γοχ + 1η(10)εα6(λ1)Γα6, (1.5)
μα(λ2) = ln(10)8ox(X2)Cox + ln(10)ede(X2)Cde. (1.6)
Here, λι and λ2 are the two wavelengths; εοχ and ede are the known molar extinction coefficients of oxy- and deoxyhemoglobin, respectively; Cox and Cde
SCATTERING AND ITS BIOLOGICAL ORIGINS 7
W 10 3 104
Wavelength (nm)
105
Figure 1.4. Absorption coefficients of primary biological absorbers.
are the molar concentrations of oxy- and deoxyhemoglobin, respectively, in the tissue. Once Cox and Cde are obtained, the oxygen saturation (SO2) and the total concentration (CHÖ) of hemoglobin can be computed as follows:
C0 so2 =
CfJb = Cox + Cde-
(1.7)
(1.8)
This principle provides the basis for pulse oximetry and functional imaging. Angio- genesis can increase Cnb, whereas tumor hypermetabolism can decrease SO2.
1.5. SCATTERING AND ITS BIOLOGICAL ORIGINS
Scattering of light by a spherical particle of any size can be modeled exactly by the Mie theory, which reduces to the simpler Rayleigh theory if the spherical particle is much smaller than the wavelength. In a scattering medium containing many scatterers that are distributed randomly in space, photons usually encounter multiple scattering events. If scatterers are sparsely distributed (where the mean distance between particles is much greater than both the scatterer size and the wavelength), the medium is considered to be loosely packed. In this case, scat- tering events are considered to be independent; hence, single-scattering theory applies to each scattering event. Otherwise, the medium is considered to be densely packed. In this case, scattering events are coupled; thus, single-scattering theory does not apply. In this book, we consider only loosely packed scattering media. Keep in mind that one must differentiate a single coupled-scattering event (which involves multiple particles) from successive independent-scattering events (each of which involves a single particle).
8 INTRODUCTION
The scattering coefficient [is is defined as the probability of photon scattering in a medium per unit path length. It has a representative value of 100 cm- 1 in bio- logical tissue. The reciprocal of μν is referred to as the scattering mean free path.
For a single scatterer, the scattering cross section σν, which indicates the scattering capability, is related to its geometric cross-sectional area og through the scattering efficiency Qs : as = Qsog. For a medium containing many scatterers with number density Ns, the scattering coefficient can be considered as the total cross-sectional area for scattering per unit volume:
is = Nsas. (1.9)
The probability of no scattering (or ballistic transmittance T) after a photon propagates over path length x can be computed by Beer’s law:
T(x) = ε χ ρ ( – μ , χ ) . (1.10)
Optical scattering originates from light interaction with biological structures, which range from cell membranes to whole cells (Figure 1.5). Photons are scat- tered most strongly by a structure whose size matches the optical wavelength and whose refractive index mismatches that of the surrounding medium. The indices of refraction of common tissue components are 1.35-1.36 for extracellular fluid, 1.36-1.375 for cytoplasm, 1.38-1.41 for nuclei, 1.38-1.41 for mitochondria and organelles, and 1.6-1.7 for melanin. Cell nuclei and mitochondria are primary scatterers. The volume-averaged refractive index of most biological tissue falls within 1.34-1.62, which is greater than the refractive index of water (1.33).
The extinction coefficient μ,, also referred to as the total interaction coefficient, is given by
μ, = μ« + μ*. ( 1 . 1 1 )
The reciprocal of μ, is the mean free path between interaction events.
Cells
Nuclei
Mitochondria
Lysosomes, vesicles
Striations in collagen fibrils Macromolecular aggregates
Membranes
ΙΟμπι
1 μηι
0.1 μπ
0.01 μπι
Figure 1.5. Biological structures of various sizes for photon scattering.
FLUORESCENCE AND ITS BIOLOGICAL ORIGINS 9
1.6. POLARIZATION AND ITS BIOLOGICAL ORIGINS
Linear birefringence (or simply birefringence), which is also known as double refraction, is the most important polarization property. A linearly birefringent material has dual principal indices of refraction associated with two linear polar- ization states of light (orientations of the electric field). The index of refraction for light polarization that is parallel with the optical axis of the material (e.g., the orientation of collagen fibers) is commonly denoted by ne, while the light is referred to as the extraordinary ray. By contrast, the index of refraction for light polarization that is perpendicular to the optical axis is commonly denoted by n0, while the light is referred to as the ordinary ray. If ne > n0, the birefringence is said to be positive. Conversely, if ne < n0, the birefringence is said to be negative.
Similarly, a circularly birefringent material has dual principal indices of refrac- tion associated with the left and the right circular polarization states of light; as a result, it can rotate a linear polarization. The amount of rotation depends on the properties and the concentration of the active material, the optical wavelength, and the path length. If the other parameters are known, the amount of rotation can reveal the concentration.
Collagen, muscle fibers, myelin, retina, and keratin have linear birefringence. Collagen I is intensely positively birefringent, whereas collagen III is weakly negatively birefringent. Amino acids and glucose have circular birefringence; amino acids are levorotatory (exhibit left rotation) to linearly polarized light, whereas glucose is dextrorotatory (exhibits right rotation).
1.7. FLUORESCENCE AND ITS BIOLOGICAL ORIGINS
Fluorescence has the following characteristics:
1. Fluorescence light is red-shifted (wavelength is increased or frequency is reduced) relative to the excitation light; this phenomenon is known as the Stokes shift. The primary origins include the initial vibrational relaxations and the subsequent inclined fluorescence transitions to higher vibrational energy levels of the ground state. Other origins include excited-state reac- tions, complex formations, and resonance energy transfers.
2. Emission wavelengths are not only longer than but also independent of the excitation wavelength. Although the initial excited state is related to the excitation wavelength, a vibrational relaxation to the same intermediate state terminates the memory of such a relationship.
3. Fluorescence light is incoherent even if the excitation light is coherent because the uncertain delays in the vibrational relaxations spread over more than one light period.
4. Fluorescence spectrum, when plotted against the frequency, is generally a mirror image of the absorption spectrum for the following reasons: (a) before excitation, almost all the molecules are at the lowest vibra- tional energy level of the ground state; (b) before emission, almost all the
10 INTRODUCTION
molecules are at the lowest vibrational energy level of the first excited state; (c) the least photon energy for excitation equals the greatest emis- sion photon energy; (d) the vibrational energy levels in the ground and first excited states have similar spacing structures; and (e) the probability of a ground-state electron excited to a particular vibrational energy level in the first excited state is similar to that of an excited electron returning to a corresponding vibrational energy level in the ground state.
The properties of some endogenous fluorophores are listed in Table 1.2 (where a denotes maximum absorption wavelength; ε denotes molar extinction coefficient; x denotes maximum excitation wavelength; m denotes maximum emission wavelength; Y denotes quantum yield of fluorescence). Fluorescence can provide information about the structure, dynamics, and interaction of a bioassembly. For example, mitochondrial fluorophore NADH (nicotinamide adenine dinucleotide, reduced form) is a key discriminator in cancer detection; it tends to be more abundant in cancer cells owing to their higher metabolic rate. NAD(P)H (nicoti- namide adenine dinucleotide phosphate, reduced form) has a lifetime of 0.4 ns when free but a longer lifetime of 1-3 ns when bound.
1.8- IMAGE CHARACTERIZATION
Several parameters are important in the characterization of medical images. In this section, the discussion is limited primarily to two-dimensional (2D) images, but the principles involved can be extended to one-dimensional (ID) or three- dimensional (3D) images.
When a high-contrast point target is imaged, the point appears as a blurred blob in the image because any practical imaging system is imperfect. The spatial distri- bution of this blob in the image is referred to as the point spread function (PSF). The PSF is sometimes called the impulse response (or Green’s function) because
TABLE 1.2. Properties of Endogenous Fluorophores at Physiologic pH
Fluorophore
Ceroid
Collagen, elastin FAD Lipofuscin
NAD+
NADH
Phenylalanine Tryptophan Tyrosine
Xö(nm)
—
— — —
260 260 340 260 280 275
e(cm_1M_1)
—
— — —
18 x 103
14.4 x 103
6.2 x 103
0.2 x 103
5.6 x 103
1.4 x 103
λ*(ηηα)
340-395
325 450
340-395
— 290 340 — 280 —
Xw(nm)
430-460 540-640
400 515
430-460 540-540
— 440 450 280 350 300
Y
—
— — —
— — —
0.04 0.2 0.1
IMAGE CHARACTERIZATION 11
a geometric point can be represented by a spatial Dirac delta function (an impulse function). When two point targets are too close to each other, the combined blob in the image can no longer be clearly resolved into two entities. The full width at half maximum (FWHM) of the PSF is often defined as the spatial resolution. Even though an ideal geometric point target cannot be constructed or detected in reality, a point target needs only to be much smaller than the spatial resolution.
Sometimes, a line spread function (LSF), which is the system response to a high-contrast geometric line, is measured instead of a PSF. For a linear system, an LSF can be related to a PSF on the (JC, y) plane by
LSF( – /
PSF(;t, y)dy. (1.12)
Likewise, an edge spread function (ESF), which is the system response to a high- contrast semiinfinite straight edge, can be measured as well. For a linear system, an ESF can be related to an LSF as follows (Figure 1.6):
ESF(JC) -Γ J-c
LSF(x)dx
LSF(x) = —ESF(x). dx
(1.13)
(1.14)
In a linear, stationary, and spatially translation-invariant system, image function i(r) equals the convolution of object function o(r) with point spread function PSF(r):
/ ( r ) = o ( r ) * * P S F ( r ) , (1.15)
1
0.8
0.6
0.4
0.2
n
–
—r
J – ** S
__ —
/ /
/ /
/ /
1
™ “” ” Lol’ ESF 1
H
, x» -“””
1 ^ ” * * *
-0.5 0.5
Figure 1.6. Illustration of an LSF and an ESF.
12 INTRODUCTION
where r — (x, y) and ** represents 2D spatial convolution. Equation (1.15) can be expressed in several forms:
i(r)= if o(r’)PSF(? -r’)dr’
o(x, / ) P S F ( J C -x’,y- y’)dx’dy’ (1.16)
o(r-r”)PSF(r”)dr”.
Taking the 2D Fourier transformation of Eq. (1.15) yields
/(ρ,ξ) = 0(ρ,ξ)//(ρ,ξ). (1.17)
Here, p and ξ represent the spatial frequencies; / represents the image spectrum; O represents the object spectrum; and H represents the PSF spectrum, which is the system transfer function (STF). The amplitude of the STF is referred to as the modulation transfer function (MTF):
ΜΤΡ(ρ,ξ) = |//(ρ,ξ)|. (1.18)
Similarly, for an LSF, the MTF is based on the ID Fourier transformation:
(1.19) /
+oo exp(-j2npx)[LSF(x)]dx
-oo
Most imaging systems act as lowpass filters, resulting in blurring of the fine structures.
The visibility of a structure in an image depends on, among other factors, the contrast C:
Δ / C = — . (1.20)
(/)
While (/) is the average background image intensity, Δ / is the intensity variation in the region of interest (Figure 1.7).
Contrast does not represent a fundamental limitation on visualization since it can be artificially enhanced by, for example, subtracting part of the background (thresholding) or raising the intensity to some power. Statistical noise does, however, represent a fundamental limitation. The signal-to-noise ratio (SNR) is defined as
SNR= — , (1.21)
IMAGE CHARACTERIZATION 13
/ ▲
L-JT” τ
Figure 1.7. Illustration of image contrast.
where σ/ denotes the standard deviation of the background intensity, that is, the noise representing the root-mean-squared (rms) value of the intensity fluctuations.
Ultimately, the ability to visualize a structure depends on the contrast-to-noise ratio (CNR), which is defined as
Δ / C N R = — , (1.22)
07
which can be rewritten as
CNR = C SNR. (1.23)
The field of view (FOV) in an image refers to the extent of the image field that can be seen all at once. A tradeoff often exists between FOV and spatial resolution. For example, “zooming in” with a camera compromises the FOV for resolution.
The maximum imaging depth in tomography is the depth limit at which the SNR or the CNR is acceptable. A tradeoff often exists between maximum imaging depth and depth resolution. The ratio of maximum imaging depth to depth reso- lution, referred to as the depth-to-re solution ratio (DRR), represents the number of effective pixels in the depth dimension. A DRR of 100 or greater is considered to indicate high resolution in terms of pixel count.
The frame rate is defined as the number of frames of an animation that are displayed per second, measured in frames per second (fps); it measures how rapidly an imaging system produces consecutive 2D images. At or above the video rate (30 fps), the human eye cannot resolve the transition of images; hence, the animation appears smooth.
In this book, the object to be imaged is typically a scattering medium, which can be a biological tissue phantom, a sample (specimen) of biological tissue, or an insitu or in vivo biological entity. Sometimes, “sample” refers broadly to the object to be imaged.
14 INTRODUCTION
Example 1.1. Derive Eq. (1.13).
On the basis of ID convolution followed by a change of variable, we derive
E S F U ) = / LSFU – x’) dx = / LSF(x”)d(-x”)
= j LSF(x’)dx’. (1.24) J — oo
PROBLEMS
1.1 Derive the following relationship between electromagnetic wavelength λ in the unit of μπ and photon energy hv in electron volts (eV): hv — 1.24, where h denotes the Planck constant and v denotes the electromagnetic frequency.
1.2 In a purely absorbing (nonscattering) medium with absorption coefficient μα, what percentage of light is left after a lightbeam propagates a length of L? Plot this percentage as a function of L in MATLAB.
1.3 In a purely absorbing (nonscattering) medium with absorption coefficient μα, derive the average length of survival of a photon.
1.4 In a purely scattering (nonabsorbing) medium with scattering coefficient μs, what percentage of light has not been scattered after the original light- beam propagates a length of L?
1.5 In a purely scattering (nonabsorbing) medium with scattering coefficient μν, derive the average length of survival of a photon.
1.6 In a scattering medium with absorption coefficient μα and scattering coef- ficient μ5, what percentage of light has survived scattering and absorption after the original lightbeam propagates a length of L? Of the percentage that has been absorbed and scattered, what is the percentage that has been absorbed?
1.7 In MATLAB, draw a 2D diagram to simulate a random walk by follow- ing the subsequent steps: (1) start the point at (0,0); (2) sample a random number x that is evenly distributed in interval (0,1]; (3) determine a step size by s = lOOlnQcj); (2) sample a random number X2 that is evenly dis- tributed in interval (0,1]; (4) determine an angle by α = 2π*2; (5) move the point by step size s along angle a; (6) repeat steps 2-5 20 times to obtain a trajectory; (7) repeat steps 1-6 3 times to trace multiple trajec- tories.
1.8 Derive the oxygen saturation SO2 and the total concentration of hemo- globin CHÖ based on Eqs. (1.5) and (1.6).
FURTHER READING 15
1.9 Download the data for the molar extinction coefficients of oxy- and deoxy- hemoglobin as a function of wavelength from the Web (URL: http://omlc. ogi.edu/spectra/) and plot the two curves in MATLAB.
1.10 Download the data for the molar extinction coefficients of oxy- and deoxy- hemoglobin as a function of wavelength from the Web (URL: http://omlc. ogi.edu/spectra/). Download the data for the absorption coefficient of pure water as a function of wavelength as well. Using physiologically repre- sentative values for both oxygen saturation SO2 and total concentration of hemoglobin CHÖ* compute the corresponding absorption coefficients. Plot the three absorption spectra on the same plot in MATLAB. Identify the low-absorption near-IR window that provides deep penetration.
READING
Drezek R, Dunn A, and Richards-Kortum R (1999): Light scattering from cells: Finite- difference time-domain simulations and goniometric measurements, Appl. Opt. 38: 3651-3661. [See Section 1.5, above (in this book).]
Jacques SL (2005): From http://omlc.ogi.edu/spectra/ and http://omlc.ogi.edu/classroom/. (See Sections 1.5 and 1.6, above.)
Richards-Kortum R and Sevick-Muraca E (1996): Quantitative optical spectroscopy for tissue diagnosis, Ann. Rev. Phys. Chem. 47: 555-606. (See Section 1.7, above.)
Wang LHV and Jacques SL (1994): Animated simulation of light transport in tissues. Laser-tissue interaction V, SPIE 2134: 247-254. (See Section 1.2, above.)
Wang LHV, Jacques SL, and Zheng LQ (1995): MCML—Monte Carlo modeling of pho- ton transport in multi-layered tissues, Comput. Meth. Prog. Biomed. 47(2): 131-146. (See Section 1.2, above.)
Wang LHV (2003): Ultrasound-mediated biophotonic imaging: A review of acousto- optical tomography and photo-acoustic tomography, Disease Markers 19(2-3): 123-138. (See Section 1.1, above.)
FURTHER READING
Hecht E (2002): Optics, Addison-Wesley, Reading, MA. Lakowicz JR (1999): Principles of Fluorescence Spectroscopy, Kluwer Academic/Plenum,
New York. Macovski A (1983): Medical Imaging Systems, Prentice-Hall, Englewood Cliffs, NJ. Mourant JR, Freyer JP, Hielscher AH, Eick AA, Shen D, and Johnson TM (1998): Mech-
anisms of light scattering from biological cells relevant to noninvasive optical-tissue diagnostics, Appl Opt. 37(16): 3586-3593.
Shung KK, Smith MB, and Tsui BMW (1992): Principles of Medical Imaging, Academic Press, San Diego.
Tuchin VV (2000): Tissue Optics: Light Scattering Methods and Instruments for Medical Diagnosis, SPIE Optical Engineering Press, Bellingham, WA.
Welch AJ and van Gemert MJC (1995): Optical-Thermal Response of Laser-Irradiated Tissue, Plenum Press, New York.
CHAPTER 2
Rayleigh Theory and Mie Theory for a Single Scatterer
2.1. INTRODUCTION
Both the Rayleigh and the Mie theories, which are based on the Maxwell equations, model the scattering of a plane monochromatic optical wave by a single particle. Even if the particle size is much greater than the optical wave- length, the wave is diffracted by the particle with an effective cross section that is usually different from the geometric cross section. The Rayleigh theory is applicable only to particles that are much smaller than the optical wavelength, whereas the Mie theory is valid for homogeneous isotropic spheres of any size. The Mie theory reduces to the Rayleigh theory when the particle is much smaller than the wavelength.
2.2. SUMMARY OF RAYLEIGH THEORY
The Rayleigh theory (Appendix 2A) models the scattering of light by particles that are much smaller than the optical wavelength. Figure 2.1 shows the spherical polar coordinates used for light scattering. The incident light propagates along the z axis; the scatterer is located at the origin; field point P is located at (r, θ, φ). The distribution of the scattered light intensity for unpolarized incident light is given by
( l+cos29)fc4 |q | 2 , /<>, Θ) = —2 /0. (2.1)
Here, a denotes the polarizability of the particle; /o denotes the incident light intensity; k denotes the propagation constant (also referred to as the magnitude of the wavevector or the angular wavenumber) in the background medium. We have
k=2-^, (2.2)
Biomedical Optics: Principles and Imaging, by Lihong V. Wang and Hsin-i Wu Copyright © 2007 John Wiley & Sons, Inc.
17
18 RAYLEIGH THEORY AND MIE THEORY FOR A SINGLE SCATTERER
* y
Figure 2.1. Spherical polar coordinates used for light scattering, where Θ denotes the polar angle and φ denotes the azimuthal angle.
where nt, denotes the refractive index of the background medium and λ denotes the wavelength in vacuum. Substituting Eq. (2.2) into Eq. (2.1), we obtain /(r, Θ) a l/λ4. This strong wavelength dependence explains the blue sky in broad day- light because blue light is scattered much more strongly than red light.
The scattering cross section is given by
Snk4a2 (2.3)
The polarizability of a sphere with radius a is given by
4i + 2 (2.4)
where nre is the relative refractive index of the particle:
“rel = — · nb
(2.5)
Here, ns is the refractive index of the sphere and nt, is the refractive index of the background. Substituting Eq. (2.4) into Eq. (2.3), we obtain
8πα x 2 V 4 a. =
<ι – 1 Λ«ι + 2
(2.6)
where size parameter x is defined as
x = ka. (2.7)
NUMERICAL EXAMPLE OF RAYLEIGH THEORY 19
Substituting Eqs. (2.7) and (2.2) into Eq. (2.6), we obtain os oc a6/4. From Eq. (2.6), the scattering efficiency is given by
0* = ΊΓ Sx4 ‘ * ‘ l2 »rel ~ 1
i L + 2 ‘rel (2.8)
The scattering efficiency depends only on x and ηκ. If nre is close to unity, Eq. (2.8) reduces to
32JC4
Qs = ^ – K e i – l | 2 · (2.9)
Note that rcrei can be complex, in which case the imaginary part is responsible for absorption.
2.3. NUMERICAL EXAMPLE OF RAYLEIGH THEORY
The Rayleigh theory can compute scattering cross section os and scattering effi- ciency Qs. As an example, the following parameters are given:
1. Diameter of sphere: 2a = 20 nm 2. Wavelength in vacuum: λ = 400 nm 3. Refractive index of sphere: ns = 1.57 4. Refractive index of background: n^ — 1.33 5. Specific weight of sphere: ρ̂ = 1.05 g/cm3
6. Specific weight of background: p^ = 1.00 g/cm3
7. Concentration of spheres in background by weight: Cwt = 1 x 10~5
We compute in SI units as follows:
1. Propagation constant in background medium: k = Inn^/X = 2.09 x 107 m-1.
2. Relative refractive index of sphere: nre — ns/ni, — 1.18. 3. Size parameter: x = ka = 0.209. 4. From Eq. (2.6), σν = 2.15 x 10~20 m2. 5. From Eq. (2.8), Qs = 6.83 x 10~5. 6. Compute the number density of scatterers Ns. For a sphere, the mass density is
Ps = ms/Vs, where ms denotes the mass and Vs denotes the volume—Vs = (4/3)πα3. For the background, the mass density is p^ = my, j V&, where m^ is the total mass of the background and V& is the total volume of the background.
2 0 RAYLEIGH THEORY AND MIE THEORY FOR A SINGLE SCATTERER
The concentration by weight is Cwt — (NsVb)ms/mi,. Therefore, we have Ns = Cwt(9h/ps)/Vs = 2.27 x 1018 m~3.
7. From Section 1.5, μ5 = Nsos = 0.0488 m _ l .
We can implement the Rayleigh theory with the following MATLAB script:
% Rayleigh scattering % Use SI units
diameter=input(‘Diameter of sphere (e.g., 20 nm):’)*1e-9; radius=diameter/2; lambda=input(‘Wavelength (e.g., 400 nm):’)*1e-9; n_sphere=input(‘Refractive index of sphere (e.g., 1.57):’); n_background=input(‘Refractive index of background (e.g., 1.33):’); w_sphere=input(‘Specific weight of sphere(e.g., 1.05 g/cc):’)*1e3; w_background=input(“Specific weight of background(e.g., 1 g/cc):’)*1e3; concentration=input(‘Concentration by weight (e.g., 1e-5):’);
k=2*pi*n_background/lambda x=k*radius n_rel=n_sphere/n_background
Qs = 8*x A 4 /3*abs( (n_rer2 – 1 ) / (n_re l A 2 + 2) )Λ2 sigma_s=Qs*pi*radiusΛ2
vol_sphere = 4*p i /3* rad ius”3 N_s=concentration*w_background/(vol_sphere*w_sphere) mu_s=N_s*sigma_s
% Output resu l t s {‘wavelength(nm) ‘ ,’Qs ( – ) V m u s ( /cm) ‘ ; lambda*1e9, Qs, mu_s/1e2}
2.4. SUMMARY OF MIE THEORY
The Mie theory (Appendix 2B) models the scattering of light by a spherical particle of any radius a. The sphere is made of homogeneous and isotropic material and is irradiated by a plane monochromatic wave. In practice, we can treat the incident wave as a plane wave if the wavefront is much wider than both the wavelength and the particle size.
Application of the Mie theory is straightforward. The scattering efficiency Qs and the scattering anisotropy g (defined by g = (cosG)) can be computed as follows:
2 °°
Qs = -2 Σ ( 2 / + 1}(ι^ι2 + ι^/i2)’ (2·10)
^ ^ Ϊ Σ [ Τ Τ Γ R e ^ > + * ‘ * ‘ + > ) + mh^01^} (2·η)
NUMERICAL EXAMPLE OF MIE THEORY 21
Size parameter x = ka. Coefficients ai and b[ are given by
Ψ/ΟΟΨ/00 – Ππ;ΐΨ/Ο0Ψ/(*) 0/ =
bi =
ΨϋΟζ/(*)-Λι*ιΨ/ϋΟζ(*) (2.12)
πι*ιΨί(3θζ/(*)-Ψ/ϋΟζί(*) ‘
where superscript prime denotes first-order differentiation and size parameter y is defined by
2nnsa y = nrdx = . (2.13)
The Riccati-Bessel functions are defined by
‘ π ζ ΐ / 2
< π ζ ΐ / 2
Ψ/(ζ) – zji(z) = ( y ) ‘ 7/+i/2(z) – St(z), (2.14)
X/(z) = -zyiiz) = – ( y ) F/+i/2(z) = Q(z), (2.15)
ζ/ω = Ψ/ω + ιχ/ω = zApfe) = ( y ) , / 2 #/(+ί/2ω. (2·16> Here, / and I + are the orders; ji( ) and y/( ) denote the spherical Bessel functions of the first and second kind, respectively; 7/( ) and K/( ) denote the (2) Bessel functions of the first and second kind, respectively; h) {) denotes the spherical Hankel function of the second kind; //, }{ ) denotes the Hankel function of the second kind; and 5/() and C/() are alternative symbols that are commonly used. Note that
A/(2)() = 7/( ) – i > ( ) , (2.Π)
/ / / 2 ) ( ) = J / O – i T / O . (2.18)
If ttrei is complex, the extinction instead of scattering efficiency that also contains a component representing absorption can be computed.
2.5. NUMERICAL EXAMPLE OF MIE THEORY
For a spherical particle of any size, the Mie theory can compute the scattering efficiency Qs, the scattering anisotropy g, and the scattering cross section os. For a scattering medium, we can further compute the scattering coefficient μ5 and the reduced scattering coefficient μ^. When the Mie theory is implemented in MATLAB or another high-level computer language, the following derivative
2 2 RAYLE1GH THEORY AND MIE THEORY FOR A SINGLE SCATTERER
identities are useful:
J,'(z) = –Ji(z) + Ji-dz). (2 z
y;(z) = –Yi(z) + Yi-(z)· (2
An example MATLAB script is shown below:
% Mie theory % Use SI units
diameter = input(‘Diameter of sphere (e.g., 579 nm):’)*1e-9; radius = diameter/2; lambda = input(‘Wavelength (e.g., 400 nm):’)*1e-9; n_s = input(‘Refractive index of sphere (e.g., 1.57):’); n_b = input(‘Refractive index of background (e.g., 1.33):’); w_s = input(‘Specific weight of sphere(e.g., 1.05 g/cc):’)*1e3; w_b = input(‘Specific weight of background(e.g., 1.0 g/cc):’)*1e3; concentration = input(‘Concentration by weight (e.g., 0.002):’);
k = 2*pi*n_b/lambda x = k*radius n_rel = n_s/n_b y = n_rel*x
% Calculate the summations err = le-8; Qs = 0; gQs = 0; for n = 1:100000
Snx = sqrt(pi*x/2)*besselj(n+0.5,x); Sny = sqrt(pi*y/2)*besselj(n+0.5,y); Cnx = -sqrt(pi*x/2)*bessely(n+0.5,x); Zetax = Snx+i*Cnx;
% Calculate the first-order derivatives Snx_prime = – (n/x)*Snx+sqrt(pi*x/2)*besselj(n-0.5,x); Sny_prime = – (n/y)*Sny+sqrt(pi*y/2)*besselj(n-0.5,y); Cnx_prime = – (n/x)*Cnx-sqrt(pi*x/2)*bessely(n-0.5,x); Zetaxprime = Snx_prime + i*Cnx_prime;
an_num = Sny_prime*Snx-n_rel*Sny*Snx_prime; an_den = Sny_prime*Zetax-n_rel*Sny*Zetax_prime; an = an_num/an_den;
bn_num = n_rel*Sny_prime*Snx-Sny*Snx_prime; bn_den = n_rel*Sny_prime*Zetax-Sny*Zetax_prime; bn = bn_num/bn_den;
Qs1 = (2*n+1)*(abs(an)^2+abs(bn)A2); Qs = Qs+Qs1;
APPENDIX 2A. DERIVATION OF RAYLEIGH THEORY 2 3
if Π > 1 gQsl = (n-1)*(n+1)/n*real(an_1*conj(an)+bn_1*conj(bn))…
+(2*n-1)/((n-1)*n)*real(an_1*conj(bn_1)); gQs = gQs+gQs1;
end
an_1 = an; bn_1 = bn;
if abs(Qs1)<(err*Qs) & abs(gQs1)<(err*gQs) break;
end end
Qs = (2/xA2)*Qs; gQs = (4/x”2)*gQs; g = gQs/Qs;
vol_s = 4*pi/3*radiusA3 N_s = concentration*w_b/(vol_s*w_s) sigma_s = 05*ρί*Γ3αίυ8Λ2; mu_s = N_s*sigma_s
mu_s_prime = mu_s*(l-g);
% Output resul ts {‘wavelength(nm)’,’Qs ( – ) V g ( – ) V m u s ( /cm)’ , ‘mus_prime(/cm)’ ; . . .
lambda*1e9,Qs,g,mu_s*1e-2,mu_s_prime*1e-2}
Below, we present a numerical example with the following inputs:
1. Diameter of sphere: 2a = 579 nm 2. Wavelength: λ = 400 nm 3. Refractive index of sphere: ns == 1.57 4. Refractive index of background: η, = 1.33 5. Specific weight of sphere: ps = 1.05 g/cm3
6. Specific weight of background: p^ = 1.0 g/cm3
7. Concentration of spheres in the solution by weight: Cwt = 0.002
The MATLAB script gives the following outputs: Qs = 2.03, g = 0.916, is = 100 cm”1, and μ̂ = 8.40 cm”1.
In Figure 2.2, Qs and g, which are calculated using a modified version of the MATLAB program presented above, are plotted against ka, where ns = 1.40 and nh = 1.33. Note that g is less than unity even for large x values.
APPENDIX 2A. DERIVATION OF RAYLEIGH THEORY
The Rayleigh theory is derived here. The polarizability a is defined as the pro- portionality constant between the induced oscillating dipole moment ρβχρ(ιωί)
2 4 RAYLEIGH THEORY AND MIE THEORY FOR A SINGLE SCATTERER
(b) x – ka
Figure 2.2. (a) Qs versus ka, where the dashed line represents the asymptote from the Rayleigh theory; (b) anisotropy g versus ka.
and the electric field of the incident linearly polarized wave Εο^χρ(ίωί), where ω denotes angular frequency and t denotes time:
p = VLEQ. (2.21)
On the basis of the dipole radiation theory for a <^ λ, the electric field of the scattered wave in the far field (r > λ) is
k2puxvy ikr t, = e (2.22)
APPENDIX 2A. DERIVATION OF RAYLEIGH THEORY 2 5
where γ is the angle between the directions of the scattered light propagation and the dipole oscillation. From Eq. (2.21) and (2.22), the scattered light intensity is
9 £4|/?|2sin2 γ £4 |a |2sin2y 7 &4|al2sin2y / = E2 = — ^ – Y- = —L» L|£0|2 = _Ji L/0f ( 2 2 3 )
where IQ is the incident light intensity. The incident light is assumed to propagate in the positive z direction. Thus, its
electric field lies on the xy plane. We express both the unit vector of polarization p and the unit vector that points to the field point P from the scatterer at the origin r in terms of the unit vectors of the Cartesian coordinates, (ex, ey, ez):
p = ex cos ψ + ey sin ψ, (2.24)
r = ex sin Θ cos φ -f ey sin Θ sin φ -f ez cos Θ. (2.25)
Thus, we obtain
cosy = p · r = sinθcos(φ — ψ), (2.26)
which leads to
sin2 Y = 1 – cos2 γ = 1 – sin2 Θ c o s 2 ^ – ψ). (2.27)
If the incident light is unpolarized, sin2 y needs to be averaged over angle ψ:
(sin2 y) = 1 – – sin2 Θ = – (1 + cos2 Θ). (2.28)
Hence
£ 4 | a | 2 ( l+cos 2 0) /(r, Θ) – 2 r 2 -/o, (2.29)
which is Eq. (2.1). The total scattering cross section σ̂ is defined as
os = – I 7(Γ,θ)Γ2ί/Ω. (2.30)
Evaluating Eq. (2.30) with differential solid angle element dQ = sin θί/θ^φ, we obtain
a , f l + c o s 2 0 8π^4|α|2 σ, = 2πΓ |α | ζ / sinOJO = — , (2.31) t4|a|2 Γ Jo
which is Eq. (2.3).
2 6 RAYLEIGH THEORY AND MIE THEORY FOR A SINGLE SCATTERER
APPENDIX 2B. DERIVATION OF MIE THEORY
The Mie theory is an exact solution of the Maxwell equations for a plane monochromatic electromagnetic wave scattered by a homogeneous sphere of radius a with an isotropic relative index of refraction nrt. An abbreviated deriva- tion of the Mie theory is presented here. The general idea is to (1) solve the Maxwell equations inside and outside the sphere with undetermined coefficients in the solution expansions and (2) determine these coefficients by applying bound- ary conditions on the spherical surface. The Mie theory for a cylindrical scatterer is beyond the scope of this book.
We assume that all conditions for the following Maxwell equations are met:
V · E = 0, (2.32)
V · B = 0, (2.33)
r) R
S/xE = , (2.34) dt
V x * = ( – ) – . (235)
Here, E and B are the electric and magnetic fields, respectively; n is the refractive index of the medium; c is the speed of light in vacuum.
2B.1. Vector and Scalar Wave Equations
The following vector wave equation for both E and B can be obtained from Eqs. (2.32)-(2.35):
/n2 32A
= Ü w ,236) V2A where A represents either E or B. Each vector component satisfies the following scalar wave equation:
ν2ψ = (ö) w i237)
Example 2.1. Derive Eq. (2.36) for E.
Operating Vx from the left on both sides of Eq. (2.34), we obtain V x (V x E) = -3 (V x B)/dt. From V x (V x E) = V(V · E) – V2E = – V 2 £ , we get
9 – dV x B V2E = . (2.38) dt
Substituting Eq. (2.35) into Eq. (2.38), we obtain Eq. (2.36) for E.
APPENDIX 2B. DERIVATION OF MIE THEORY 2 7
2B.2. Solution of Scalar Wave Equation
The standard procedure for solving Eq. (2.37) is by separation of the variables. Let Ψ(χ, 0 = X(x)T(t), where x represents the spatial coordinates. Substituting this expression into Eq. (2.37) and dividing both sides by XT, we obtain
V2X /n2 1 /d2T
The left-hand side of this equation is a function of only x, whereas the right-hand side is a function of only t thus, the two sides must equal a constant, which is termed a separation variable constant’, this constant is denoted by — β2. The time-dependent part can be expressed as
(d2T + ω2Τ = 0, (2.40)
where ω = ßc/n. Since ω is the angular frequency of the wave, we have β = k = nko, where k and ko are the propagation constants in the medium and in vacuum, respectively. The solution of Eq. (2.40) is
C O S ( 0 / ) . (2.41) sinooi J
The pair of braces represents a linear combination of the two functions inside. The spatially dependent part is a scalar Helmholtz equation
V2X + k2X = 0, (2.42)
which can be expressed in spherical coordinates as
– x — [r2— + – (sinO— I + = X + k2X = 0. r2dr drj r2 sind dd 3Θ/ r2 sin2 Θ d<>2 J
(2.43)
Letting X = /?(λ-)Θ(θ)Φ(φ), substituting it into Eq. (2.43), and then dividing the equation by /?(Γ)Θ(Θ)Φ(Φ), we obtain
1 a / 2dR 1 3 / . dS 1 32Φ f 2 Λ -Ϊ r2— ) + -T ( s inG— 1 + = T + * = 0. r2Rdr dr ) r 2 0 s in0aOV 3Θ / Γ2Φδίη2θ3φ2
(2.44)
Multiplying Eq. (2.44) by r2 sin2 Θ, we see that the third term on the left-hand side depends only on φ, whereas the other terms are independent of φ. Thus, we let
1 32Φ φ ” ^
2 constant = — m . (2.45)
2 8 RAYLEIGH THEORY AND MIE THEORY FOR A SINGLE SCATTERER
Thus, we have
92Φ 3φ2
which has the following solutions:
+ m2cD = 0, (2.46)
Φ cos(m(|)) sin(m$) I * (2.47)
From condition Φ(0) = Φ(2π), we have m = 0, ± 1 , ±2, The r0-dependent part from Eq. (2.44) is
l a / 7dR i a / . a © —^ r2— I + -^ sinG— ) – – r2Rdr dr r 2 0 s i n 0 a 0 V 9Θ / r 2sin29
+ r = 0. (2.48)
Multiplying this equation by r2, we see that the first and fourth terms depend only on r, whereas the other terms are independent of r. Letting
i a 0sin0äÖ
we have
s i n0— –
r2 dr dr
—7Γ- = constant = —/(/ + 1), sin20
+ .2 K / + D
R = 0
(2.49)
(2.50)
and
d-z2)- d2@ Λ ^ Θ Γ / ( /+ ! ) – m dz2 dz -z2
where z = cos θ. The solutions of Eq. (2.50) are
1θ = 0, (2.51)
R jl(kr) | _ | ji(nk0r) yiihr) I I yi(nkQr)
(2.52)
Here, / is an integer; ji(x) and yi(x) are the spherical Bessel functions of the first and second kinds, respectively, given by
JiM = J^JIM/2)W,
yiM = J 2^Yi+(/2)(x)-
(2.53)
(2.54)
APPENDIX 2B. DERIVATION OF MIE THEORY 2 9
The solutions of Eq. (2.51) are
Θ ~ ( ^ ( ; 0 8 θ > } , (2.55) l Ö/,m(cOS0) j
where P/m (cosΘ) and Qm (cosΘ) are the associated Legendre polynomials of the first and second kinds, respectively.
Recalling that Ψ ( ί , 0 = X(x)T(t) and X = / ? (Γ)Θ(Θ)Φ(Φ) , we obtain the solution of the scalar wave equation [Eq. (2.37)] in spherical polar coordinates (see Figure 2.1):
ψ _ | cos(ü>f) ] ί cos(/w<|>) I | ji(nk0r) J J P/,m(cosO) | { sin(ü>r) J [ sin(/w<|>) J { yi(nk0r) J { ß/,m(cos9) J
We use expO’cof) from the linear combination of the time-dependent factors for the phasor expressions of the waves. Of course, we can also use exp(—ι’ωί) consistently. Function Qi m (cosO) is dropped off because it has singularities at Θ € {0, π}. For a wave inside the sphere, function ylinker) is dropped off because ylinker) -» —oo when r -> 0. For an outgoing spherical wave outside
(2) the sphere, the spherical Hankel function of the second kind h {nk^r) is chosen because of its asymptotic behavior:
(2)
h) (nkor) — exp(—ink$r). (2.57) nkor
Therefore, we use
*~«HÄ)H*?w>},,-<a,,’)· (“8) If exp(—ΐωί) is used, the spherical Hankel function of the first kind h) nkor) should be used instead for the outgoing spherical wave.
2B.3. Theorem Relating Solutions of Scalar and Vector Wave Equations
From Eq. (2.37), the time-independent scalar wave equation (the scalar Helmholtz equation) is given by
V2X + rc2^X = 0, (2.59)
which has the following solutions:
cos(m(|)) 1 [ ji(nkor) } ( 4 ^ ( 2 · 6 0 ) sin(ra<t>) I I hi(nk0r)
3 0 RAYLEIGH THEORY AND MIE THEORY FOR A SINGLE SCATTERER
Similarly, from Eq. (2.36), the time-independent vector equation (the vector Helmholtz equation) is given by
V2A+n2klA = 0. (2.61)
The solution to the vector Helmholtz equation can be found from the following theorem. If X satisfies the scalar wave equation [Eq. (2.59)], vectors Μχ and Νχ, defined by
Mx = V x (rX) and nk0Nx = V x Mx, (2.62)
must satisfy the vector Helmholtz equation [Eq. (2.61)], where Μχ and Νχ are related by
nk0Mx = V x Nx. (2.63)
The full components of Μχ and Νχ are
Μθ
Μφ
nkoNr
nk0Ne
nk0N^
Equations (2.67) can also be expressed as
1 d / dX 1 d2X nk0Nr = — – s i n 0 — – — – j – – ^ . (2.70)
If u and v are two solutions to the scalar wave equation, and Mu, Nu, Mv, and Nv are the derived vector fields, the spatial components of the solutions of the vector wave equations are
E = Mv + iNu and B = -(-Mu + iNv). (2.71) c
o, 1 d(rX)
rsinG 3φ 1 d(rX)
~~r 9Θ ‘ i a * 2HrX)~ r
1 32(rX)
~r 3r39 ‘ 1
r sin Θ d2(rX)
+ n2k^r X,
(2.64)
(2.65)
(2.66)
(2.67)
(2.68)
(2.69)
APPENDIX 2B. DERIVATION OF MIE THEORY 3 1
The components of E and B can thus be written in terms of the scalar solutions u and v and their first and second derivatives.
Example 2.2. Show that vector Μχ defined by Eq. (2.62) satisfies the vector wave equation given by Eq. (2.61).
Multiplying r and then operating Vx from the left through Eq. (2.59), we obtain
V x (?V2X) + n2klMx = 0. (2.72)
From the vector identity
V[V · (VX x r)] = r(V x VX) – VX(V x r) = 0, (2.73)
we derive
V x V x (VX x r) = V[V · (VX x r)] – V2(VX x r)
= -V 2 (VX x r) . (2.74)
Also
V x [V x (VX x r)] = V x [(r · V)VX – ?V2X – (VX · V)r + VX(V · r)]
= V x [(r · V)VX] – V x [rV2X] – V x [(VX · V)r]
+ V x [ V X ( V . r ) ] , (2.75)
The first, third, and fourth terms on the right-hand side can be evaluated as follows:
r · V = r— => V x [(r · V)VX] = V x ( r —VX ) = 0, (2.76) Or dr )
(VX .V)? = V X 4 V x [(VX · V)?] = V x VX = 0, (2.77)
V · r = 3 => V x [VX(V · r)] = 0. (2.78)
Therefore, Eq. (2.75) becomes
V x [V x (VX x r)] = – V x (rV2X). (2.79)
From Eqs. (2.74) and (2.79), we obtain
V2(VX x r) = V x (?V2X). (2.80)
We also obtain
V2MX = V2[V x (rX)] = V2[VX x r 4- X(V x r)] = V2(VX x r). (2.81)
3 2 RAYLEIGH THEORY AND MIE THEORY FOR A SINGLE SCATTERER
Merging Eqs. (2.72), (2.80), and (2.81) yields
V2Mx+n 2klMx = 0, (2.82)
which shows that Μχ is a valid solution of the vector wave equation.
Example 2.3. Verify the relationships in Eqs. (2.64)-(2.66).
From the vector operation of curl in spherical coordinates
V X V :
we obtain
f
r sin0 9 9Ve — (snGVd.) 9Θ φ 9φ
φ Γ 9 dVr~ + – – ( r V e ) – – ^ r _dr 9Θ _
Θ + –
r
‘
M x = V x (rX) ■■ Θ
sin 3X
Θ9φ
” 1 3Vr _sin9 9φ
.ax – Φ — ,
– έ(^φ)] (2.83)
(2.84)
which can be rewritten as Eqs. (2.64)-(2.66).
2B.4 Solution for Coefficients from Boundary Conditions
The origin of the coordinates is set at the center of the spherical scatterer. The positive z axis is set along the propagation direction of the incident wave. The x axis is set in the plane of electric vibration of the linearly polarized incident wave. Solutions u and v are chosen in the following forms:
1. For the incident plane wave outside the spherical particle, we have
°° 2/ + 1 u =η^6χρ( /ωΟοο8φ^(- / ) / —-—-P/ , i (cose)7 / ( / : r ) , (2.85)
/= i /(/ + 1)
2/ + 1 v = η^χρ^ωΟύη^Υ^^ί)1 -—— P/j(cosO)y/(*r), (2.86)
/=!
where k represents the propagation constant of the background. 2. For the scattered wave outside the spherical particle, we have
°° 2/ + 1 u = rib exp(io)t) cosφ 22, ~~ai^~^ fyiicosOJAj (fcr), (2.87)
/=! / ( /+!) 00 2/4-1
v = n^exp(i^0 s i n φ ^ ] – & / ( – / ) ‘ PlA[cosü)hfkr), (2.88) / = i
APPENDIX 2B. DERIVATION OF MIE THEORY 3 3
where a and b are coefficients to be determined. 3. For the wave inside the spherical particle, we have
00 2/ + 1 u = η5εχρ(ίωί)οο$φΣαι(-ί)
1 ——–Pij (cos d)jt(nrekr), (2.89)
°° 2/ + 1 υ = π,βχρίιωθβίηφ Υ ] φ(-ί)1 ——— P/,i(cos0).//(ftrei*r), (2.90)
where Q and d/ are also coefficients to be determined.
To determine these coefficients, we substitute solutions E and B from Eqs. (2.85) -(2.90) into the following boundary conditions:
än x (E0 – Ei) = unx (B0 – Bi) = 0. (2.91)
Here, un is the unit vector perpendicular to the boundary surface; subscripts o and i represent the outside and inside, respectively. The obtained coefficients a and b are shown in Eq. (2.12).
2B.5 Scattering Efficiency and Anisotropy
Substituting the asymptotic expression of the spherical Hankel function of the second kind into Eqs. (2.87) and (2.88) leads to
exp(—ikr -f ίωί) v-> 2/ -f 1 u = -i-^—, c o s φ Γ ο , — — Ρ,,i(cosΘ), (2.92)
kr ^ /(/ + 1)
exp(—ikr + ιωί) ^ 2/ -f- 1 υ = – ί F . – s i ^ V ^ / — — f t , , ( c o s e ) . (2.93)
fcr ^ /(/ + 1) The resulting field components are
exO(—ikr + ιωί) ΕΘ = Β^ = ~i
FV y – cos φ52(θ), (2.94)
exp(—/fcr -f /ωί) – £ φ = £θ = – t y βίηφΑΚΘ). (2.95)
kr The amplitude functions are given by
2/4-1 5ι(β) = ΣΖ 77ϊ^-ΤΤ[α’π'(°08θ> + *ix/(cose)], (2.96)
S2Q) = Z^JTTT—[&/^(cose)+a,T/(cose)], (2.97)
3 4 RAYLEIGH THEORY AND MIE THEORY FOR A SINGLE SCATTERER
where
π,(cosΘ) = J . Λ , (2.98) sinO
x/(cose) = ^ P / j ( c o s e ) . (2.99) αυ
The scattering efficiency Qs, defined as the ratio of the scattering cross-sectional area as to the physical cross-sectional area πα2, can be expressed in terms of the amplitude functions:
Qs = ^ = – ^ f ( |5,(θ) |^θ82φ + |52(θ)|28ίη2φ)^Ω. (2.100)
Following the integration over φ, Eq. (2.100) becomes
Qs = f (|5,(θ)|2 + |52(θ)|2)8ίηθ^θ. (2.101) x2 Jo
Likewise, the scattering anisotropy g — (cos Θ) can be evaluated by
gQs = f (|51(0)|2 + |52(0)|2)cosesineJ6. (2.102)
The integrations over Θ in Eqs. (2.101) and (2.102) can be completed using the orthogonality relations of π/ and x/, which yields Eqs. (2.10) and (2.11).
PROBLEMS
2.1 Show that
2.2 Show that
2.3 Show that
Sji(kr) k . 0 = -zr—rVji-dkr) – (I + l)ji+(kr)]. or 21 + 1
= —ψ/(ζ) + ψ/-ι(ζ). dz z
— — = — ζ / ( ζ ) + ζ / – ι ( ζ ) . dz z
2.4 Derive the net radiation force exerted by light on a spherical particle. (Hint: The photon momentum equals the photon energy divided by the speed of light.)
2.5 Plot the angular distribution of the scattered light in the Rayleigh scattering regime and calculate the scattering anisotropy.
FURTHER READING 3 5
2.6 Implement the Mie theory using an alternative program to calculate the scattering efficiency, scattering anisotropy, scattering coefficient, and reduced scattering coefficient of spherical particles suspended in a back- ground medium. (a) Use the example in Section 2.5 to verify your program. (b) Duplicate Figure 2.2. (c) Summarize the asymptotic dependence of the scattering cross section
on the particle size as ka varies. 2.7 Derive coefficients a and b shown in Eqs. (2.12) by completing the
derivations in Appendix 2B. 2.8 (a) Extend the Mie theory to absorbing scatterers.
(b) Extend the Rayleigh theory to absorbing scatterers.
READING
Van de Hülst HC (1981): Light Scattering by Small Particles, Dover Publications, New York. (See Sections 2.2, 2.4, 2.5, and 2.5, above.)
FURTHER READING
Bohren CF and Huffman DR (1983): Absorption and Scattering of Light by Small Particles, Wiley, New York.
Ishimaru A (1997): Wave Propagation and Scattering in Random Media, IEEE Press, New York.
CHAPTER 3
Monte Carlo Modeling of Photon Transport in Biological Tissue
3.1. INTRODUCTION
Photon transport in biological tissue can be numerically simulated by the Monte Carlo method. The trajectory of a photon is modeled as a persistent random walk, with the direction of each step depending on that of the previous step. By contrast, the directions of all the steps in a simple random walk are independent. By tracking a sufficient number of photons, we can estimate physical quantities such as diffuse reflectance.
3.2. MONTE CARLO METHOD
Although widely used in various disciplines, the Monte Carlo method defies a succinct definition. Here, we adopt the description provided by Lux and Koblinger (1991):
In all applications of the Monte Carlo method, a stochastic model is constructed in which the expected value of a certain random variable (or of a combination of several variables) is equivalent to the value of a physical quantity to be determined. This expected value is then estimated by averaging multiple independent samples representing the random variable introduced above. For the construction of the series of independent samples, random numbers following the distribution of the variable to be estimated are used.
It is important to realize that the Monte Carlo method estimates ensemble- averaged quantities.
The Monte Carlo method offers a flexible yet rigorous approach for simulating photon transport in biological tissue. An ensemble of biological tissues is mod- eled for the averaged characteristics of photon transport; the ensemble consists of all instances of the tissues that are microscopically different but macroscop- ically identical. Rules are defined for photon propagation from the probability
Biomedical Optics: Principles and Imaging, by Lihong V. Wang and Hsin-i Wu Copyright © 2007 John Wiley & Sons, Inc.
37
3 8 MONTE CARLO MODELING OF PHOTON TRANSPORT IN BIOLOGICAL TISSUE
distributions of, for example, the angles of scattering and the step sizes. The statistical nature requires tracking a large number of photons, which is compu- tationally time-consuming. Multiple physical quantities can be simultaneously estimated, however.
In this chapter, photons are treated as waves at each scattering site but as clas- sical particles elsewhere. Coherence, polarization, and nonlinearity are neglected. Structural anisotropy—not to be confused with scattering angular anisotropy—in tissue components, such as muscle fibers or collagens, is neglected as well.
3.3. DEFINITION OF PROBLEM
The problem to be solved begins with an infinitely narrow photon beam, also referred to as a pencil beam, that is perpendicularly incident on a multilay- ered scattering medium (Figure 3.1); various physical quantities are computed as responses. The pencil beam can be represented by an impulse (Dirac delta) function of space, direction, and time; thus, the responses are termed impulse responses or Green’s functions. The layers are infinitely wide and parallel to each other. Each layer is described by the following parameters: thickness d, refractive index n, absorption coefficient μβ, scattering coefficient μ.ν, and scat- tering anisotropy g. The top and the bottom ambient media are each described by a refractive index. Although never infinitely wide in reality, a layer can be so treated if it is much wider than the photon distribution.
Three coordinate systems are defined. A global Cartesian coordinate system (JC, y, z) is used for tracking photons (Figure 3.1); the xy plane coincides with the surface of the scattering medium; the z axis is along the pencil beam.
A global cylindrical coordinate system (r, φ’, ζ), which shares the z axis with the Cartesian coordinate system, is used for recording photon absorption as a function of r and z. The photon absorption distribution has cylindrical symmetry
Photon beam
Layer 1
Layer 2
Layer Nj
Figure 3.1. Schematic of a scattering medium with Ni layers. The y axis points out of the paper.
PROPAGATION OF PHOTONS 3 9
because of the homogeneity of each layer and the axial alignment of the pencil beam. The r coordinate is also used for recording both the diffuse reflectance and the diffuse transmittance that are functions of r as well as a, where a is the polar angle of the propagation direction of a reemitted photon with respect to the normal vector of the exit surface of the scattering medium (—z axis for the top surface and +z axis for the bottom surface). One can further resolve reemitted photons with the azimuthal angle φ’.
A local moving spherical coordinate system whose z axis is dynamically aligned with the propagation direction of the photon is used for sampling the scattering angles. Once the polar angle Θ and the azimuthal angle φ are sampled, they are converted to direction cosines in the global Cartesian coordinate system.
The physical quantities to be computed include relative specific absorption, rel- ative fluence, relative diffuse reflectance, and relative diffuse transmittance, all of which are relative to the incident energy. The relative specific absorption A(r, z) represents the probability of photon absorption per unit volume by the medium. From A(r, z), the relative fluence F(r,z)—which is the probability of photon flow per unit area—can be computed. The unscattered absorption from the first photon interaction events, which always take place on the z axis, is recorded sepa- rately. The relative diffuse reflectance Rd(r, a) for the top surface and the relative diffuse transmittance 7^(r, a) for the bottom surface—collectively referred to as the relative diffuse reemittance—are defined as the probability of photon reemis- sion from the surfaces per unit area at r per unit solid angle around a, where a solid angle has the unit of steradians (sr). Like unscattered absorption, specularly reflected and unscattered transmitted photons are recorded separately. Physical quantities of lower dimensions can be computed from those of higher dimen- sions. For brevity, relative physical quantities are written as physical quantities in this chapter unless otherwise noted.
Simulated physical quantities are represented in grids on the coordinate sys- tems. For photon absorption, a 2D homogeneous grid system is set up in the r and z directions. The grid element sizes in the r and z directions are ΔΓ and Δζ, respectively, and the total numbers of grid elements are Nr and Nz, respec- tively. For reemitted photons, a ID grid system in the a direction is further set up with Na grid elements. Since a has a range of π/2, the grid element size is Δα = π/(2Να). For convenience, the grid elements that should appear after a physical quantity are sometimes represented in this chapter by coordinates.
For consistency, centimeters (cm) are used as the basic unit of length through- out the simulation. For example, the thickness of each layer and the grid element sizes in the r and z directions are measured in cm. The absorption and scattering coefficients are measured in reciprocal centimeters (cm-1).
3 A PROPAGATION OF PHOTONS
This section presents the rules that govern photon propagation. Figure 3.2 shows a basic flowchart for the photon tracking part of the Monte Carlo simulation of
4 0 MONTE CARLO MODELING OF PHOTON TRANSPORT IN BIOLOGICAL TISSUE
/ Launch photon Udimensionless s_ = 0y
Figure 3.2. Flowchart for tracking photons in a scattering medium by the Monte Carlo method, where s_ denotes the dimensionless step size (to be discussed) and Y and N represent yes and no, respectively.
light transport in multilayered scattering media. The Monte Carlo simulation was written in ANSI (American National Standards Institute) Standard C as a software package entitled MCML (Appendix 3A). This software can be executed on any computer platform that supports ANSI Standard C. The following subsections implement many of the boxes in the flowchart.
PROPAGATION OF PHOTONS 41
3.4.1. Sampling of a Random Variable
The Monte Carlo method relies on the sampling of random variables from their probability distributions. The probability density function (PDF) /?(χ) defines the distribution of χ over interval (a, b). The interval can also be closed or half- closed in some cases, which usually makes no practical difference. For readers unfamiliar with PDF, a brief review is given in Appendix 3B.
To sample χ, we choose a value repeatedly based on its PDF. First, a pseu- dorandom number ξ that is uniformly distributed between 0 and 1 is generated by computer. Then, χ is sampled by solving the following equation:
/ ‘ Ja
ρ ( χ ) ^ χ = ξ. (3.1)
Since the left-hand side represents the cumulative distribution function (CDF) ^ (x ) , we have
P(lO = *. (3.2)
This equation means that if Ρ(χ) is sampled uniformly by ξ between 0 and 1, the inverse transformation correctly samples χ as illustrated in Figure 3.3:
χ = />”ι(ξ). (3.3)
This sampling method, referred to as the inverse distribution method (IDM), is invoked repeatedly below.
Figure 3.3. Illustration of the inverse distribution method (IDM) for sampling a random variable χ based on a random number ξ uniformly distributed between 0 and 1.
4 2 MONTE CARLO MODELING OF PHOTON TRANSPORT IN BIOLOGICAL TISSUE
Example 3.1. Prove Eq. (3.1)
We show that when χ is sampled according to Eq. (3.1), χ follows its intended CDF Ρ(χ) .
If χ is sampled according to Eq. (3.1), the probability that a selected χ is less than xi can be expressed by
/>ιθΜ{χ<Χι} = />{ξ<ξι} · (3.4)
Here, PIDM{ } denotes the probability in the IDM of the event in the braces; P{} denotes the true probability of the event in the braces based on the CDF of the random variable; χι is related to ξι through Eq. (3.1):
*;■= Γ%(χ)<ίχ. (3-5) Ja
Because ξ is equidistributed between 0 and I, we have Ρ{ξ < ξι} = ξ. Thus, we obtain
^ I D M { X < X I } = [ ‘ PdOdm. (3.6) Ja
Since the right-hand side is the CDF Ρ(χ) , we have
A D M ( X < X I } – / > ( X I ) , (3.7)
which means that the sampled χ indeed follows its intended CDF.
3.4.2. Representation of a Photon Packet
A simple variance reduction technique, referred to as implicit photon capture, is used to improve the efficiency of the Monte Carlo simulation. This technique enables many photons to propagate as a packet of equivalent weight W along a particular trajectory.
Related parameters are structured in C to make the program easier to write, read, and modify. Thus, parameters for a photon packet are grouped into a single structure defined by
typedef struct { double x, y ,z ; /* Cartesian coordinates of photon packet. */ double ux, uy, uz ; / * d i rect ion cosines of photon propagation. */ double w; /* weight of photon packet. */ char dead; /* 0 i f photon is propagating, 1 i f terminated. */ double s_; /* dimensionless step size. */ long scatters; /* number of scatter ing events experienced. */ short layer; /* index of layer where photon packet resides. */
} PhotonStruct;
PROPAGATION OF PHOTONS 4 3
Structure members x, y, and z represent the coordinates of a photon packet, (JC, y, z), respectively. Structure members ux, uy, and uz represent the direction cosines of the propagation direction of the photon packet, {μχ, μ>;, μ^}, respectively. Structure member w represents the weight of the photon packet, W.
Structure member dead represents the status of the photon packet and is initialized to zero when the photon packet is launched. If the photon packet has reemitted from the scattering medium or has not survived a Russian roulette (to be discussed), this structure member is set to unity, which signals the program to stop tracking the current photon packet.
Structure member s_ represents the dimensionless step size, which is defined as the integration of extinction coefficient μ, over the trajectory of the photon packet. In a homogeneous medium, the dimensionless step size is simply the physical path length multiplied by the extinction coefficient.
Structure member s c a t t e r s stores the number of scattering events expe- rienced by the photon packet. If this structure member is zero when the weight of a photon packet is recorded, the weight contributes to the unscat- tered physical quantities.
Structure member l aye r is the index of the layer in which the photon packet resides; it is updated only when the photon packet leaves the current layer.
3.4.3. Launching of a Photon Packet
One photon packet is launched at a time orthogonally onto the scattering medium. For each new photon packet, the coordinates (x, y, z) are initialized to (0,0,0), the direction cosines {μ*, μγ,μζ] are initialized to (0,0,1), and the weight W is ini- tialized to unity. Several other structure members—including dead, s c a t t e r s , and layer—are also initialized.
If the upper ambient medium and the first layer have mismatched indices of refraction (no and n, respectively), specular (Fresnel) reflection occurs. Specular reflectance^ in normal incidence is given by
K S P = ( ^ V . 0.8)
If the first layer of refractive index n is a clear medium, it causes multiple spec- ular reflections and transmissions. If interference effect is neglected, an effective specular reflectance can be computed by
Asp = flspl + – j · (3.9)
Here, Rsp and Rsp2 are the specular reflectances on the two boundaries of the first layer:
*spi = (**^λ) , (3.10)
4 4 MONTE CARLO MODELING OF PHOTON TRANSPORT IN BIOLOGICAL TISSUE
( Π — ηι *sp2= – ~ – – , (3.11)
where n2 denotes the refractive index of the second layer. When the photon packet enters the scattering medium, the weight is decreased
by #sP:
W=l- Rsp. (3.12)
If Eq. (3.8) is applicable, structure member l ayer is set to the first layer. If Eq. (3.9) is applicable, it is set to the second layer, and member z is set accord- ingly.
Example 3.2. Compute the specular reflectance in normal incidence between air and water, glass, and soft tissue, respectively.
For water, n = 1.33. For an air-water interface, Rsp = 2.0%, where HQ = 1 for air.
For glass, n — 1.5. For an air-glass interface, Rsp — 4.0%. For soft tissue, n — M.37. For an interface between air and soft tissue,
Rsp = -2 .4%.
Example 3.3. Compute the effective specular reflectance in normal incidence for a glass slide that is sandwiched between air and water.
For air, no = 1. For glass, n — 1.5. For water, n2 = 1.33. Thus, R$p = 0.04, and /?sp2 = 0.0035. The effective specular reflectance Rsp = 0.0432. In this case, Rsp is quite close to the direct sum of Rsp and /?sp2<
3.4.4. Stop Size of a Photon Packet
The step size of a photon packet is sampled by the IDM [Eq. (3.3)], based on the PDF of the dimensional free-path length s (0 < s < oo)) of a photon. We first consider a homogeneous medium. According to the definition of the extinction coefficient μ,, we have
~dP{s > s’) 1 . . . . . μ — (3.13) ^ P{s>s’} dsf
where P{] denotes the probability of the event in the braces. The first frac- tion on the right-hand side represents the probability of interaction in interval (s s’ + ds’), and the second fraction represents normalization by the path length. Rearranging Eq. (3.13) yields
d[n(P{s > s’})] = – μ , ds’. (3.14)
PROPAGATION OF PHOTONS 4 5
Integrating this equation over s’ in interval (0, s), we obtain
P{s > 5ι} = βχρ(-μ , ί ι ) . (3.15)
This is a form of the well-known Beer law. From Eq. (3.15), the probability that an interaction occurs within s is given by
P{s < s] — 1 – exp(—its). (3.16)
The corresponding PDF is given by dP{s <s]
P(s) = = μ, cxp(-itsi). (3.17) ds
According to the IDM, the CDF of s in Eq. (3.16) can be equated to ξ to yield the sampling equation for the step size
ln(l ~ %) n 15n s = , (3.18)
where 1 — ξ can be replaced by ξ for simplicity because ξ is uniformly distributed between 0 and 1:
„ = – ί ^ . (3.19) μ*
We then consider a multilayered medium where the photon packet may experience a free path over multiple layers before an interaction occurs. In this case, the counterpart of Eq. (3.15) becomes
-Y^VtiSi . P{s > st) = exp I – 2 ^ VtiSi 1 . (3.20)
Here, the summation is over all the segments that the photon packet has traversed before an interaction occurs; μπ is the extinction coefficient for the ith segment, Si is the length of the ith segment, and st is the total step size:
st Σ*· <3·21)
Equating Eq. (3.20) to ξ, we obtain the sampling equation
£ > „ * , = – 1 η ( ξ ) , (3.22)
which is a generalized form of Eq. (3.19). The left-hand side of Eq. (3.22) is the total dimensionless step size. Note that photon paths in a clear medium do not change the dimensionless step size because the extinction coefficient is zero.
4 6 MONTE CARLO MODELING OF PHOTON TRANSPORT IN BIOLOGICAL TISSUE
Equation (3.22) is used to sample the step size in MCML, where the dimen- sionless step size s_ is initialized to — 1η(ξ). A photon packet may travel multiple substeps of size s, in a multilayered scattering medium before reaching an interac- tion site. Only when the photon packet has completed — 1η(ξ) in dimensionless step size does an interaction occur. In an interaction event, the entire photon packet must experience both absorption and scattering. Since the step size is modeled, the simulation is intrinsically time-resolved.
3.4.5. Movement of a Photon Packet
Once the dimensional substep size s; is determined, the photon packet is moved. The coordinates of the photon packet are updated by
x^-x + [ixSj, y<^y + VLySi, ζ « – ζ + μζ5/ , (3.23)
where the arrows indicate quantity substitutions. The variables on the left-hand side have the new values, and the variables on the right-hand side have the old values. In C/C 4- +, an equal sign (=) is used for this purpose.
3.4.6. Absorption of a Photon Packet
Once the photon packet reaches an interaction site, a fraction of the weight (AW) is absorbed:
AW = — W. (3.24)
If the photon packet has not been scattered, AW is recorded into an array for unscattered absorption. Otherwise, it is recorded into A(r, z) at the local grid element:
A(r,z) <~ A(r,z) + AW. (3.25)
The weight must be updated by
W+-W-AW. (3.26)
The photon packet with the new weight then undergoes scattering at the interac- tion site.
3.4.7. Scattering of a Photon Packet
For scattering, the polar angle Θ (0 < θ < π) and the azimuthal angle φ (0 < φ < 2π) are sampled by the IDM. The probability distribution of cos Θ is commonly given by the Henyey-Greenstein phase function that was originally proposed for galactic scattering:
P(cos Θ) – * – ^ . (3.27) 2(1 + gl — 2 gcosO)-5/2
PROPAGATION OF PHOTONS 4 7
The anisotropy g, defined as (cos0), has a value between —1 and 1. A value of zero indicates isotropic scattering, and a value close to unity indicates dom- inantly forward scattering. For most biological tissues, g is M).9. In addition to the Henyey-Greenstein phase function, the Mie theory can also provide a phase function. Note that the phase function is irrelevant to the phase of an electromagnetic wave.
Applying the IDM [Eq. (3.1)] to Eq. (3.27), we sample cosG as follows:
COS0 = (ΐΐ4ξ)} lf * * ° . (3.28) 2% – 1 if g = 0
The azimuthal angle φ, which is assumed to be uniformly distributed over interval [0, 2π), is sampled with another independent pseudorandom number ξ:
φ = 2πξ. (3.29)
Once the local polar and azimuthal angles are sampled, the new propagation direction of the photon packet can be represented in the global coordinate sys- tem as
/ s in0^j^zcos( | ) — μν sinc|)) μ* = , l· μ* cosO,
iz cos φ 4-μν = p = h μν cos Θ, (3.30) ι2
μζ = -J — μ2 8 ΐηθ^8φ + μ^θ8θ .
If the photon direction is sufficiently close to the z axis (e.g., |μζ | > 0.99999), the following formulas are used instead so that division by a small number is avoided:
μχ = sin Θ cos φ,
μ^ = sin Θ sin φ, (3.31)
μ̂ = s g n ^ ) c o s 0 ,
where sign function sgn^ z ) returns 1 when μ̂ is positive and —1 when μ̂ is negative.
Note that the direction cosines are in the global Cartesian coordinate system, whereas the polar and azimuthal angles are in the local spherical coordinate system. Since trigonometric operations are computationally intensive, alternative algebraic operations are encouraged in the program when possible.
4 8 MONTE CARLO MODELING OF PHOTON TRANSPORT IN BIOLOGICAL TISSUE
Example 3.4. Derive Eq. (3.28) for g φ 0 using the IDM [Eq. (3.1)].
Here, χ = cos6 and χ e [—1, 1]. Therefore, we have
/ x /»cose
ρ(χ) d = / p(COS0′) dcOSd’
= izi! / » L·^ 2g V A / I + ^ – ^ ^ C O S G 1 + g /
which leads to Eq. (3.28) for g φ 0.
3.4.8. Boundary Crossing of a Photon Packet
During a step of dimensionless size s_, the photon packet may hit a boundary of the current layer. Several steps are involved in boundary crossing.
Step 1. The distance df, between the current location (x, y, z) of the photon packet and the boundary of the current layer in the photon propagation direction is computed by
db =
zo-z if μζ < 0
oo if μ, = 0 , (3.32) Ζ — Ζ
μ< if μζ > 0
where zo and z are the z coordinates of the upper and lower boundaries of the current layer. If μζ approaches zero, the distance approaches infinity, which is represented by constant DBL-MAX in C. Step 2. The dimensionless step size s_ and the distance db are compared as follows:
db[Lt < * – , (3-33)
where μ, is the extinction coefficient of the current layer. If Eq. (3.33) holds, the photon packet is moved to the boundary; s_ is updated by s_ <— s_ — ^ μ , ; the simulation proceeds to step 3. Otherwise, the photon packet is moved by s-/it to reach the interaction site, s- is set to zero to signal the generation of the next dimensionless step size, and the photon packet experiences absorption and scattering. Since division by μ, is avoided, Eq. (3.33) is applicable to clear layers (μ, = 0) as well. Step 3. If the photon packet hits a boundary, the specular reflectance is computed. The angle of incidence a, of the photon packet is first calculated by
a , = c o s – V z l ) . (3.34)
PROPAGATION OF PHOTONS 4 9
The angle of transmission at is then computed by Snell’s law
rii sin a,· = nt sin a,, (3.35)
where Π[ and nt represent the refractive indices of the incident and transmit- ted media, respectively. If n, > nt and a, is greater than the critical angle sin_1(nr/n/), the local reflectance /?,((*,) equals unity. Otherwise, /?,·(α,·) is cal- culated by Fresnel’s formula:
*/(««)= 2 sin2 (a,· — a,) tan2 (a,; — at) sin2 (a, + ar) tan2 (a,· + a , )
(3.36)
which is an average of the reflectances for two orthogonal linear polarization states because light is assumed to be randomly polarized. Step 4. To determine whether the photon packet is reflected or transmitted, a pseudorandom number ξ is generated. If ξ < /?/(α,·), the photon packet is reflected; otherwise, it is transmitted. If reflected, the photon packet stays on the boundary and its direction cosines are updated by reversing the z component:
{μχ, μ^, μζ} <r- [[ix, μν, – μ , } . (3.37)
If transmitted into a neighboring layer, the photon packet continues to propagate with an updated direction and step size. The new direction cosines are
s i n o t ‘ n M v [i =μχ- , (3.38)
sin a,· S i n 0 t ‘ / – J i m
μ Y = Vy- , (3.39) y sin a,
μζ = s g n ^ ) cos a,. (3.40)
On the basis of Snell’s law, Eqs. (3.38) and (3.39) can be computed more effi- ciently by
μχ=μχ —, (3.41) nt
f Hi
μ ^ μ ^ —. (3.42)
If transmitted into an ambient medium, the photon packet contributes to reemit- tance. If the photon packet has not been scattered, its weight is recorded into the unscattered reemittance; otherwise, its weight is recorded into either the diffuse reflectance Rd(r, at) or the diffuse transmittance 7^(r, at):
Rd(r,at) *-Rd(r,at) + W if z = 0; (3.43)
7^(r, at) <- Td(r, at) -f W if z is at the bottom of the medium.
5 0 MONTE CARLO MODELING OF PHOTON TRANSPORT IN BIOLOGICAL TISSUE
When reemitted from the scattering medium, the photon packet completes its history (or Markov chain).
Reemission at an interface can be modeled alternatively. Rather than making the reflection of the photon packet an all-or-none event, the photon packet can be partially reflected and partially transmitted. A fraction 1 – /?/((*,) of the current weight of the photon packet is reemitted from the scattering medium and recorded to the local reemittance. Then, the weight is updated by W <— W/?/(a;). The photon packet with the new weight is reflected and further propagated.
3.4.9. Termination of a Photon Packet
A photon packet can be terminated from the scattering medium automatically by reemission as discussed above. If the weight of a photon packet has been sufficiently decreased by many interaction events, further propagation of the photon packet yields little useful information unless interest is focused on a late stage of photon propagation. However, photon packets must be properly terminated so that energy is conserved.
A technique called Russian roulette is used in MCML to terminate a photon packet when the weight falls below a threshold W^ (e.g., W^ — 0.0001). This technique gives the photon packet one chance in m (e.g., m — 10) of surviving with a weight of m W. In other words, if the photon packet does not survive the Russian roulette, it is terminated with the weight set to zero; otherwise, the photon packet increases in weight from W to mW. This technique is mathematically summarized as
ί mW if ξ < 1 , W <- m (3.44)
0 if %> -,
where ξ is a uniformly distributed pseudorandom number (0 < ξ < 1). This method terminates photons in an unbiased manner while conserving the total energy.
3.5. PHYSICAL QUANTITIES
In this section, the process of obtaining physical quantities is discussed in detail. The units for some of the physical quantities are given at the end of their respec- tive expressions.
Physical quantities are stored in arrays. Although photon packets propagate in infinite continuous space (limited only by computational precision), weights are recorded in finite discrete space (limited by grid element sizes). When a photon packet is recorded, its physical location may not fit into the grid system. In this case, the last grid element in the direction of the overflow collects the weight. Therefore, the last grid elements in the r and z directions do not reflect the actual values at the corresponding locations. Angle a, however, does not overflow.
PHYSICAL QUANTITIES 51
As a rule of thumb for most problems, each spatial grid element should mea- sure about 10% of the penetration depth or the transport mean free path. If the grid elements are too small, the relative errors—which are determined by the number of events occurring in each grid element—will be too large. If the grid elements are too large, the dependence of the physical quantities will not be represented with good resolution.
3.5.1. Reflectance and Transmittance
The diffuse reflectance and the diffuse transmittance are represented in MCML by two arrays, Rd-ralhJa] and Td_r0L[ir, ia], respectively, where ir and ia are the indices of r and a (0 < ir < Nr — 1 and 0 < ia < Na — 1), respectively. The unscattered reflectance and the unscattered transmittance are stored in Rd_r[— 1] and Td_r[— 1], respectively.
It can be shown that the optimal coordinates of the simulated physical quanti- ties for the grid elements in the radial and angular directions are as follows (see Problem 3.1):
r(ir) = H) 1 + Ar (cm), (3.45) 12(/r + i ) αθ’α) = ( ι’α + – j Δα + 1 – -Aaco t l – A a j cot L· + – j Aa (rad).
(3.46)
The deviation of the optimized point from the center is 25% for the first radial grid element but decreases as the index of the grid element increases. Because the optimized coordinates are computed only after all photon packets are simulated, this optimization does not increase simulation time but improves accuracy.
After multiple photon packets (N) have been tracked, raw data Rd^rAh, ia] and Td_ra[ir, ia] represent the total accumulated weights in each grid element. To compute the integrated weights in each direction of the 2D grid system, we sum the 2D arrays in the other dimension:
Na-
Rd-rUA = Σ Rd-rAir, l a] , (3.47) ία=0
Nr-
Rd^aUa] = Σ Rd-raUr, I 0 ] , (3 .48) i r=0
Na-
Td_rVr]= £ Td_ra[ir,ial (3.49) ί’α=0
Nr-
Td-oilia] = Σ Td-raUrJal (3-50) ,7=0
5 2 MONTE CARLO MODELING OF PHOTON TRANSPORT IN BIOLOGICAL TISSUE
To compute the total diffuse reflectance and the total diffuse transmittance, we sum the ID arrays:
Nr-
Rd = Σ Rd-r[ir]> ( 3 5 1 )
/r=0
N,–
Τα=Σ Td-rVr]. (3.52) /r=0
These arrays describe the total weights in each grid element on the basis of N initial photon packets of unit weight. Raw data /?</_ra[/r, *<*] and 7j_m[/ r , ia] are converted into probabilities of reemission per unit cross-sectional area per unit solid angle as follows:
Rd-ra[h,ia] <- – -prr: (cm sr ), (3.53)
Td_r*[ir,i*+- J d~rAlrl^AT ( cm-V” 1 ) . (3.54)
AacosaAilN The area Aa and the solid angle ΔΩ are given by
= 2π(,ν + 1)
ΔΩ = 4π sin I ia + ~ I Δα sin I – Δ α 1 (sr). (3.56)
Raw Rd-rUr] and Td_r[iA are converted to probabilities of reemission per unit area as follows:
* r f – r [ i V ] ^ % ^ (cm”2), (3.57) N Aa
7 i _ r [ i V ] ^ – % x ^ (cm-2). (3.58) N Aa
Raw Rd_r[— 1] and 7^_Γ[—1] are converted to total unscattered reflectance and total unscattered transmittance, respectively, by dividing them by N. Then, Rd-r[—l] is augmented by the specular reflectance or the effective specular reflectance. Raw Rd^aUa] and 7^_α[/α] are converted to the probabilities of reemission per unit solid angle as follows:
Ärf-J/α] ^ – £ χ 7 ~ (sr”‘), (3.59) ΝΑΩ
Td^Va)^7^^ ( s r 1 ) . (3.60) ΝΑΩ
Aa = 2π ( ir + – ) (Ar)2 (cm2), (3.55)
PHYSICAL QUANTITIES 5 3
Raw Rj and 7^ are converted to probabilities as follows:
Rd Rd < (dimensionless), (3.61)
N T
Td < (dimensionless). (3.62) N
3.5.2. Absorption and Fluence
At each interaction, the absorbed weight is stored in the specific absorption array Arz[ir, iz], where ir and iz are the indices of grid elements in the r and z directions, respectively (0 < ir < Nr — 1 and 0 < iz < Nz — 1). The unscattered absorption is stored in Arz[—l,iz]. Whereas the optimal coordinate for ir is shown in Eq. (3.45), the coordinate for iz is simply
A) = (/z + ) zdz) = iz + 2 ) Az’ ( 3 ‘ 6 3 )
Raw Arz[ir, iz] represents the total accumulated weight in each grid element. The total weight in each grid element in the z direction is computed by summing the 2D array in the r direction:
Nr-
Az[iz]= J^Arz[ir,izl (3.64) i r = 0
Next, the total weight absorbed in each layer, A/[//], can be computed by
A/[//] = ^ A J i J , (3.65)
where // is the index of a layer, and the summation includes any iz that leads to a z coordinate in layer //. Further, the total weight absorbed in the scattering medium A can be computed by
Nz- Α = Σ A*vJ’ ( 3 · 6 6 )
i z=0
Then, these raw quantities are converted into the final physical quantities as follows:
Arz[irJz]+- A’z[h’*z] (cm”3), (3.67)
NAaAz
Ard-hit]*- An[J*’il] (cm-3), (3.68)
NAz
AzVz]^^r ( c m _ l ) ‘ ( 3 6 9 > NAz
5 4 MONTE CARLO MODELING OF PHOTON TRANSPORT IN BIOLOGICAL TISSUE
Mh] AiUi] < (dimensionless), (3.70)
N A
A < (dimensionless). (3.71)
The ID array Λ/[//] represents the probability of photon absorption in each layer. The quantity A represents the probability of photon absorption by the entire scattering medium. The 2D array Arz[ir, iz] represents the probability of photon absorption per unit volume, which can be converted into fluence Frz[ir,iz] as follows:
FrzUrJz)^ a U (cm”2), (3.72)
where ia denotes the local absorption coefficient. This equation breaks down in a non-absorbing medium.
The ID array Az[iz] represents the probability of photon absorption per unit length in the z direction, which can be converted to a dimensionless quantity Fz[iz] as follows:
M’z] Fz[iz] = (dimensionless). (3.73) μ«
Fz[iz] represents the internal fluence as a function of z apart from a constant factor. This equation also breaks down in a non-absorbing medium.
Example 3.5. Show the equivalence of Eq. (3.73) to the convolution over an infinitely wide uniform beam (see Chapter 4).
According to Eqs. (3.64), (3.67) and (3.69), the final converted Az[iz] and Arz [ir, iz] have the following relationship:
Nr-
Az[iz]= Σ Α Γ Ζ [ Ι Γ , / 2 ] Δ Ο ( Ι Γ ) , (3.74) ir=0
where Aa(ir) is computed by Eq. (3.55). Employing Eqs. (3.72) and (3.73), we convert Eq. (3.74) to
Nr- F^z]= £ F r z [ / r , / J A a ( / r ) . (3.75)
/r=o
This equation is a discrete representation of the following integral: poo
Fz(z)= / Frz(r,z)2nrdr. (3.76) Jo
This integration is essentially the convolution over an infinitely wide uniform beam of unit fluence. Therefore, Fz[iz] represents the fluence distribution along the z axis in response to such a beam.
COMPUTATIONAL EXAMPLES 5 5
3.6. COMPUTATIONAL EXAMPLES
As computational examples, MCML simulation results are compared with results from other theories and other investigators’ Monte Carlo simulations. These com- parisons partially validate the MCML.
3.6.1. Total Diffuse Reflectance and Total Transmittance
The total diffuse reflectance Rd and the total transmittance Tt (sum of both the unscattered and the scattered transmittances) are computed for a slab of scatter- ing medium with the following optical properties: relative refractive index nre = 1, absorption coefficient μα = 10 cm- 1 , scattering coefficient μ5 = 90 cm- 1 , anisotropy g = 0.75, and thickness d = 0.02 cm. The relative refractive index is defined as the ratio of the refractive index of the scattering medium to that of the ambient medium. If nTe = 1, the boundaries are termed refractive-index- matched. After 10 Monte Carlo simulations of 50,000 photon packets each are completed, the averages and the standard errors (standard deviations of the averages) of the total diffuse reflectance and total transmittance are computed and compared with the data from van de Hulst’s (1980) table and Prahl et al.’s (1989) Monte Carlo simulations (Table 3.1). Because the unscattered trans- mittance is exp[—(μα + [Ls)d] = e~2 = 0.13534, the total diffuse transmittance equals 0.66096-0.13534 = 0.52562. All results are in good agreement. It is worth noting that standard errors are expected to decrease proportionally with the square root of the number of photon packets tracked owing to the central limit theorem.
For a semiinfinite scattering medium that has a refractive-index-mismatched boundary (nre φ 1), the total reflectance is computed similarly and compared in Table 3.2 with the data from Giovanelli (1955) and Prahl et al.’s (1989) Monte Carlo simulations. The scattering medium has the following optical properties: nTe = 1.5, ia = 10 cm”1, [is = 90 cm- 1 , g = 0 (isotropic scattering). Then, 10 Monte Carlo simulations of 5000 photon packets each are completed to compute the average and the standard error of the total diffuse reflectance.
TABLE 3.1. Comparison of Results from MCML with van de Hulst’s Table and the Monte Carlo Data of Prahl et al. (1989)α
Source
MCML van de Hust (1980) Prahl et al. (1989)
Rd Average
0.09734 0.09739 0.09711
Rd Error
0.00035
0.00033
Tt Average
0.66096 0.66096 0.66159
Tt Error
0.00020
0.00049 a”Rd average” and “/?</ error” list the averages and the standard errors of the total diffuse reflectance, respectively. Columns ‘Ύ, average” and ‘T, error” list the averages and the standard errors of the total transmittance, respectively.
5 6 MONTE CARLO MODELING OF PHOTON TRANSPORT IN BIOLOGICAL TISSUE
TABLE 3.2. Comparison of Total Reflectance0 from MCML with Data from Giovanelli (1955) and Prahl et al. (1989)
Source Rj Average Rj Error
MCML 0.25907 0.00170 Giovanelli 1995 0.2600 — Prahl et al. (1989) 0.26079 0.00079
“The total reflectance includes both the specular reflectance (0.04) and the diffuse reflectance
3.6.2. Angularly Resolved Diffuse Reflectance and Transmittance
The angularly resolved diffuse reflectance and transmittance are also computed for a slab of scattering medium with the following optical properties: nre — 1, xa = 10 cm- 1 , [is = 90 cm- 1, g = 0.75, and d = 0.02 cm. In this simulation, 500,000 photon packets are tracked, and the number of angular grid elements is 30. The results from MCML are compared in Figure 3.4 with the data from van de Hülst’s (1980) table. Because van de Hülst used a different definition of reflectance and transmittance and also used an incident flux of π, his data are multiplied by cos a and then divided by π before the comparison.
3.6.3. Depth-Resolved Fluence
As an example, the depth-resolved internal fluence Fz[iz] is simulated by MCML for two semiinfinite scattering media with refractive-index-matched and refractive-index-mismatched boundaries, respectively (Figure 3.5). The optical parameters are nrf> = 1.0 or 1.37, μα = 0.1 cm- 1 , μ5 = 100 cm- 1, and g = 0.9. One million photon packets are tracked in each simulation. The grid element size and the number of grid elements in the z direction are 0.005 cm and 200, respectively.
As can be seen, the fluence near the surface is greater than unity because scattered light augments the fluence. Furthermore, the internal fluence in the scattering medium with a refractive-index-mismatched boundary is greater than that in the medium with a refractive-index-matched boundary, because photons can be bounced back into the scattering medium by the mismatched boundary.
When z is greater than the transport mean free path lr diffusion theory predicts that the internal fluence distribution is
F(z) = KF0exp(~y (3.77)
Here, K is a scaling factor that depends on the relative index of refraction, F() is the incident irradiance (unity in MCML), and δ is the penetration depth. Consequently, the two curves in the diffusive regime should be separated by a constant factor only, which means that the tails of the two curves should be
COMPUTATIONAL EXAMPLES 57
0.025
0.02
a: 0.015
C 0.01
Q 0.005
– D D D O
^ n n ° · α ~ D ° u □ a D D Rfc
D MCML • van de Hülst
D
0π Ο.ΐπ (a)
0.2π 0.3π Exit angle a (rad)
0.4π 0.5π
0.8
0.7
‘fe 0.6
£ 0 . 5
0.4
1 0.3
0.2
0.1
m Γ D
r D
i ■ ■ ■■ T ■ l
ü MCML 1 1 • van de Hülst H
J I D J
I D J D
Γ D Ί D
l · -J Γ D | D D D D D * ^ D
Οπ Ο.ΐπ (b)
0.2π 0.3π Exit angle a (rad)
0.4π 0.5π
Figure 3.4. Angularly resolved (a) diffuse reflectance Rd(a) and (b) diffuse transmittance Td(OL).
parallel in a log-linear plot (Figure 3.5). The two curves shown here are parallel when z > l’t = [μα + μ*(1 – g)]~l ^ 0.1 cm.
We fit the parallel portions of the two curves with exponential functions. The decay constants for the curves are approximately 0.578 cm for the refractive- index-matched boundary and 0.575 cm for the refractive-index-mismatched boundary. Both values are close to the one predicted by the diffusion theory
δ = 1
ν 3 μ α [ μ 0 + μ 5 ( 1 – g)] = 0.574 cm, (3.78)
which is independent of the relative index of refraction.
5 8 MONTE CARLO MODELING OF PHOTON TRANSPORT IN BIOLOGICAL TISSUE
10
c 3
1
0 0.2 0.4 0.6 0.8 1
z (cm)
Figure 3.5. Comparison of the internal fluence as a function of depth z for two semiin- finite scattering media with a refractive-index-matched boundary and a refractive-index- mismatched boundary, respectively. Minus (—) represents a dimensionless unit.
APPENDIX 3A. SUMMARY OF MCML
The entire source code for MCML can be found on the Web at ftp://ftp.wiley.com/ public/sci_tech_med/biomedical_optics. The whole program for MCML is divided into several files. Header file mcml. h defines the data structures and some con- stants. File mcmlmain.c contains the main function and the status-reporting function. File mcmlio.c deals with reading and writing data. File mcmlgo.c contains most of the photon-tracking code. File mcmlnr. c contains several func- tions for dynamical data allocations and error reports. Readers should read the main function first.
In MCML, ID and 2D physical quantities are stored in ID or 2D arrays, respectively. These arrays are dynamically allocated so that users are allowed to specify the array sizes at runtime without wasting computer memory, an advan- tage that static arrays do not provide.
The following list is generated by command cf low -d3 -n –omit-arguments –omit-symbol-names mcml*.c, which shows the structure of the program (MCML 1.2.2) with the nesting depth limited to 3:
1 main() < i n t () at MCMLMAIN.C:186>: 2 ShowVersion() <void () a t MCMLI0.C:71>: 3 C t rPu ts ( ) <void () at MCMLI0.C:48>: 4 p r i n t f ( ) 5 pu ts ( ) 6 GetFnameFromArgvO <void () at MCMLMAIN.C:130>:
APPENDIX 3A. SUMMARY OF MCML 5 9
7 strcpy() 8 GetFileO <FILE * () at MCMLIO.C:111>: 9 printf() 10 scant() 11 strlen() 12 exit() 13 fopen() 14 CheckParm() <void () at MCMLIO.C:531>: 15 ReadNumRuns() <short () at MCMLIO.C:222>: 16 printf() 17 ReadParm() <void () at MCMLIO.C:442>: 18 FnameTaken() <Boolean () at MCMLIO.C:504>: 19 sprintf() 20 free() 21 nrerrorO <void () at MCMLNR.C:17>: 22 FreeFnameListO <void () at MCMLI0.C:517>: 23 rewind() 24 ReadNumRunsO <short () at MCMLIO.C:222>: 25 FindDataLine() <char * () at MCMLIO.C:201>: 26 strcpy() 27 nrerror() <void () at MCMLNR.C:17>: 28 sscanf() 29 ReadParm() <void () at MCMLIO.C:442>: 30 ReadFnameFormatO <void () at MCMLIO.C:242>: 31 ReadNumPhotons() <void () at MCMLI0.C:260>: 32 ReadDzDr() <void () at MCMLIO.C:277>: 33 ReadNzNrNa() <void () at MCMLIO.C:293>: 34 ReadNumLayers() <void () at MCMLI0.C:316>: 35 ReadLayerSpecs() <void () at MCMLIO.C:390>: 36 CriticalAngleO <void () at MCMLIO.C:421>: 37 DoOneRun() <void () at MCMLMAIN.C:145>: 38 InitProfileO 39 cecho2file() 40 InitOutputDataO <void () at MCMLIO.C:562>: 41 Rspecular() <double () at MCMLG0.C:116>: 42 PunchTimeO <time_t () at MCMLMAIN.C:59>: 43 printf() 44 PredictDoneTimeO <void () at MCMLMAIN.C:94>: 45 LaunchPhoton() <void () at MCMLGO.C:140>: 46 HopDropSpin() <void () at MCMLGO.C:734>: 47 exit() 48 ReportResultO <void () at MCMLMAIN.C:115>: 49 FreeData() <void () at MCMLIO.C:598>: 50 fcloseO
6 0 MONTE CARLO MODELING OF PHOTON TRANSPORT IN BIOLOGICAL TISSUE
APPENDIX 3B. PROBABILITY DENSITY FUNCTION
The probability density function (PDF), expressed as p(x), is a function that gives the probability that random variable x assumes a value between jti and X2 as follows:
P{x] < x < xi) = / p(x) dx. (3.79)
The PDF has the following properties:
“+00
/ ” p(x) dx = 1 (3.80)
and p(x) > 0 for x e (-oo, +oo). (3.81)
The relationship between the PDF and the cumulative distribution function (CDF) P(x) is
P(x) = I p{x’) dx J—oo
(3.82)
or
PROBLEMS
p(x) = —P(x). (3.83) ax
3.1 Derive Eqs. (3.45) and (3.46).
3.2 Derive Eq. (3.30). Note that the formula depends on the choice of the local moving coordinate system. The Monte Carlo simulation, however, leads to the same result.
3.3 Find the effective specular reflectance in normal incidence from a water layer (n — 1.33) placed between air and biological tissue (n = 1.37).
3.4 Show that the number of scattering events occurring in path length s follows the Poisson distribution. Neglect absorption.
3.5 Show that the mean and standard deviation of the free path between scat- tering events are both equal to 1/μΛ. Neglect absorption.
3.6 An alternative to Eq. (3.28) when g φ 0 is
cos Θ = — 2g
14V- \+g-2g%)
Why?
PROBLEMS 61
An absorbing slab of thickness d has absorption coefficient ia with neg- ligible scattering. A collimated laser beam of intensity Iin is normally incident on the slab. (a) Assume that the mismatch between the refractive indices of the ambi-
ent media and the absorbing medium is negligible. Calculate the trans- mitted light intensity /out.
(b) Assume that the mismatch between the refractive indices of the ambi- ent media and the absorbing medium causes a specular reflection R on each surface of the slab. Recalculate the transmitted light intensity /out·
A pencil beam is incident on a semiinfinite reference scattering medium. The spatially and temporally resolved diffuse reflectance Ro(r, t) is known for the following optical parameters: absorption coefficient μαο, scattering coefficient μ5ο, and scattering anisotropy go, where r denotes the radial distance between the observation point and the point of incidence, and t denotes time. The speed of light in the medium is c. (a) Write the expression for the new diffuse reflectance /?i(r, /) assuming
that the absorption coefficient is changed from μαο to μα but the other optical parameters are unchanged.
(b) Write the expression for the new diffuse reflectance Ri(r, t) where both the absorption coefficient and the scattering coefficient are scaled by the same factor C and the other optical parameter is unchanged: μ«2 = Ciao and xs2 = Cvs0-
Consider a Gaussian beam with a radius of R, a total energy of 1 mJ, and a radial energy density distribution of S(r) = (2/nR2) exp(—2r2/R2). Derive the sampling expression for the radius r based on the random number ξ that is uniformly distributed between 0 and 1.
Implement a Monte Carlo simulation of photon transport in a semiinfi- nite scattering medium in C/C++ or another programming language. The Henyey-Greenstein function is assumed to be the phase function. Inputs include nrei, μα, μ*, and g. Outputs from the program include the total reflectance R and the depth-resolved fluence. (a) Use the results in Table 3.2 to verify your program. (b) Calculate for nK — 1.37, μα = 0 . 1 cm- 1 , μ5 = 100 cm- 1 , and g =
0.9. (c) Reproduce Figure 3.5.
Use the original Monte Carlo code written for Problem 3.10 and a mod- ified version to verify that the following two algorithms are statistically identical. Compare the total diffuse reflectances and the total absorptions using a statistical test. Use nre — 1, μα = 0.1 cm- 1 , μ5 = 100 cm- 1 , and 8 – 0.9. (a) Sample the step size with s — — 1η(ξ)/μ, and calculate the weight loss
at each step by AW = ^-W.
6 2 MONTE CARLO MODELING OF PHOTON TRANSPORT IN BIOLOGICAL TISSUE
(b) Sample the step size with s = —1η(ξ)/μν and calculate the weight loss at each step by AW = W[ – exp(—μα5)].
READING
Giovanelli RG (1955): Reflection by semi-infinite diffusers, OpticaActa 2: 153-162. (See Section 3.6, above.)
Lux I and Koblinger L (1991): Monte Carlo Particle Transport Methods: Neutron and Photon Calculations, CRC Press, Boca Raton. (See Section 3.2, above.)
Marquez G, Wang LHV, Lin SP, Schwartz JA, and Thomsen SL (1998): Anisotropy in the absorption and scattering spectra of chicken breast tissue, Appl. Opt. 37(4): 798-804. (See Section 3.2, above.)
Prahl SA, Keijzer M, Jacques SL, and Welch AJ (1989): A Monte Carlo model of light propagation in tissue, in Dosimetry of Laser Radiation in Medicine and Biology, Müller GJ and Sliney DH, eds., SPIE Press, IS5: 102-111. (See Section 3.6, above.)
van de Hülst HC (1980): Multiple Light Scattering: Tables, Formulas, and Applications, Academic Press, New York. (See Section 3.6, above.)
Wang LHV, Jacques SL, and Zheng LQ (1995): MCML—Monte Carlo modeling of light transport in multi-layered tissues, Comput. Meth. Prog. Biomed. 47: 131-146. (All sections in this chapter.)
FURTHER READING
Ahrens JH and Dieter U (1972): Computer methods for sampling for the exponential and normal distributions, Commun. ACM 15: 873-882.
Baranoski GVG, Krishnaswamy A, and Kimmel B (2004): An investigation on the use of data-driven scattering profiles in Monte Carlo simulations of ultraviolet light prop- agation in skin tissues, Phys. Med. Biol. 49(20): 4799-4809.
Bartel S and Hielscher AH (2000): Monte Carlo simulations of the diffuse backscattering Mueller matrix for highly scattering media, Appl. Opt. 39(10): 1580-1588.
Boas DA, Culver JP, Stott JJ, and Dunn AK (2002): Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head, Opt. Express 10(3): 159-170.
Carter LL and Cashwell ED (1975): Particle-Transport Simulation with the Monte Carlo Method USERDA Technical Information Center, Oak Ridge, TN.
Cashwell ED and Everett CJ (1959): A Practical Manual on the Monte Carlo Method for Random Walk Problems, Pergamon Press, New York.
Churmakov DY, Meglinski IV, Piletsky SA, and Greenhalgh DA (2003): Analysis of skin tissues spatial fluorescence distribution by the Monte Carlo simulation, J. Phys. D—Appl. Phys. 36(14): 1722-1728.
Cote D and Vitkin IA (2005): Robust concentration determination of optically active molecules in turbid media with validated three-dimensional polarization sensitive Monte Carlo calculations, Opt. Express 13(1): 148-163.
Ding L, Splinter R, and Knisley SB (2001): Quantifying spatial localization of optical map- ping using Monte Carlo simulations, IEEE Trans. Biomed. Eng. 48(10): 1098-1107.
FURTHER READING 6 3
Flock ST, Patterson MS, Wilson BC, and Wyman DR (1989): Monte-Carlo modeling of light-propagation in highly scattering tissues. 1. Model predictions and comparison with diffusion-theory, IEEE Trans. Biomed. Eng. 36(12): 1162-1168.
Flock ST, Wilson BC, and Patterson MS (1989): Monte-Carlo modeling of light- propagation in highly scattering tissues. 2, Comparison with measurements in phan- toms, IEEE Trans. Biomed. Eng. 36(12): 1169-1173.
Gardner CM and Welch AJ (1994): Monte-Carlo simulation of light transport in tis- sue—unscattered absorption events, Appl. Opt. 33(13): 2743-2745.
Hendricks JS and Booth TE (1985): MCNP variance reduction overview, Lecture Notes Phys. (AXV64) 240: 83-92.
Henniger J, Minet O, Dang HT, and Beuthan J (2003): Monte Carlo simulations in com- plex geometries: Modeling laser light transport in real anatomy of rheumatoid arthritis, Laser Phys. 13(5): 796-803.
Igarashi M, Gono K, Obi T, Yamaguchi M, and Ohyama N (2004): Monte Carlo simula- tions of reflected spectra derived from tissue phantom with double-peak particle size distribution, Opt. Rev. 11(2): 61-67.
Jacques SL, Later CA, and Prahl SA (1987): Angular dependence of HeNe laser light scattering by human dermis, Lasers Life Sei. 1: 309-333.
Jacques SL (1989): Time resolved propagation of ultrashort laser-pulses within turbid tissues, Appl. Opt. 28(12): 2223-2229.
Jacques SL (1989): Time-resolved reflectance spectroscopy in turbid tissues, IEEE Trans. Biomed. Eng. 36(12): 1155-1161.
Jaillon F and Saint-Jalmes H (2003): Description and time reduction of a Monte Carlo code to simulate propagation of polarized light through scattering media, Appl. Opt. 42(16): 3290-3296.
Kahn H (1950): Random sampling Monte Carlo techniques in neutron attenuation prob- lems I, Nucleonics 6: 27-37.
Kalos MH and Whitlock PA (1986): Monte Carlo Methods, Wiley, New York. Keijzer M, Jacques SL, Prahl SA, and Welch AJ (1989): Light distributions in artery
tissue—monte-carlo simulations for finite-diameter laser-beams, Lasers Surg. Med. 9(2): 148-154.
Keijzer M, Pickering JW, and Gemert MJCv (1991): Laser beam diameter for port wine stain treatment, Lasers Surg. Med. 11: 601-605.
Kienle A and Patterson MS (1996): Determination of the optical properties of turbid media from a single Monte Carlo simulation, Phys. Med. Biol. 41(10): 2221-2227.
Li H, Tian J, Zhu FP, Cong WX, Wang LHV, Hoffman EA, and Wang G (2004): A mouse optical simulation environment (MOSE) to investigate bioluminescent phenomena in the living mouse with the Monte Carlo method, Acad. Radiol. 11(9): 1029-1038.
Lu Q, Gan XS, Gu M, and Luo QM (2004): Monte Carlo modeling of optical coherence tomography imaging through turbid media, Appl. Opt. 43(8): 1628-1637.
MacLaren MD, Marsaglia G, and Bray AT (1964): A fast procedure for generating expo- nential random variables, Commun. ACM 7: 298-300.
Marsaglia G (1961): Generating exponential random variables, Ann. Math. Stat. 32: 899-900.
McShane MJ, Rastegar S, Pishko M, and Cote GL (2000): Monte carlo modeling for implantable fluorescent analyte sensors, IEEE Trans. Biomed. Eng. 47(5): 624-632.
Meglinski IV (2001): Monte Carlo simulation of reflection spectra of random multilayer media strongly scattering and absorbing light, Quantum Electron. 31(12): 1101 -1107.
6 4 MONTE CARLO MODELING OF PHOTON TRANSPORT IN BIOLOGICAL TISSUE
Meglinski IV and Matcher SD (2001): Analysis of the spatial distribution of detector sen- sitivity in a multilayer randomly inhomogeneous medium with strong light scattering and absorption by the Monte Carlo method, Opt. Spectrosc. 91(4): 654-659.
Mroczka J and Szczepanowski R (2005): Modeling of light transmittance measurement in a finite layer of whole blood—a collimated transmittance problem in Monte Carlo simulation and diffusion model, Optica Applicata 35(2): 311-331.
Nishidate I, Aizu Y, and Mishina H (2004): Estimation of melanin and hemoglobin in skin tissue using multiple regression analysis aided by Monte Carlo simulation, J. Biomed. Opt. 9(4): 700-710.
Patwardhan SV, Dhawan AP, and Relue PA (2005): Monte Carlo simulation of light-tissue interaction: Three-dimensional simulation for trans-illumination-based imaging of skin lesions, IEEE Trans. Biomed. Eng. 52(7): 1227-1236.
Plauger PJ and Brodie J (1989): Standard C, Microsoft Press, Redmond, WA. Plauger PJ, Brodie J, and Plauger PJ (1992): ANSI and ISO Standard C : Programmer’s
Reference, Microsoft Press, Redmond, WA. Qian ZY, Victor SS, Gu YQ, Giller CA, and Liu HL (2003): “Look-ahead distance” of a
fiber probe used to assist neurosurgery: Phantom and Monte Carlo study, Opt. Express 11(16): 1844-1855.
Ramella-Roman JC, Prahl SA, and Jacques SL (2005): Three Monte Carlo programs of polarized light transport into scattering media: Part I, Opt. Express 13(12): 4420-4438.
Sharma SK and Banerjee S (2003): Role of approximate phase functions in Monte Carlo simulation of light propagation in tissues, J. Opt. A 5(3): 294-302.
Swartling J, Pifferi A, Enejder AMK, and Andersson-Engels S (2003): Accelerated Monte Carlo models to simulate fluorescence spectra from layered tissues, J. Opt. Soc. Am. A 20(4): 714-727.
Tycho A, Jorgensen TM, Yura HT, and Andersen PE (2002): Derivation of a Monte Carlo method for modeling heterodyne detection in optical coherence tomography systems, Appl Opt. 41(31): 6676-6691.
van de Hülst HC (1980): Multiple Light Scattering: Tables, Formulas, and Applications, Academic Press, New York.
Wang LHV and Jacques SL (1993): Hybrid model of Monte Carlo simulation diffusion theory for light reflectance by turbid media, J. Opt. Soc. Am. A 10: 1746-1752.
Wang LHV and Jacques SL (1994): Optimized radial and angular positions in Monte Carlo modeling, Med. Phys. 21: 1081-1083.
Wang LHV, Nordquist RE, and Chen WR (1997): Optimal beam size for light delivery to absorption-enhanced tumors buried in biological tissues and effect of multiple-beam delivery: A Monte Carlo study, Appl. Opt. 36(31): 8286-8291.
Wang LHV (2001): Mechanisms of ultrasonic modulation of multiply scattered coherent light: A Monte Carlo model, Opt. Lett. 26(15): 1191-1193.
Wang RKK (2002): Signal degradation by coherence tomography multiple scattering in optical of dense tissue: A Monte Carlo study towards optical clearing of biotissues, Phys. Med. Biol. 47(13): 2281-2299.
Wang XD, Yao G, and Wang LHV (2002): Monte Carlo model and single-scattering approximation of the propagation of polarized light in turbid media containing glucose, Appl. Opt. 41(4): 792-801/
Wang XD, Wang LHV, Sun CW, and Yang CC (2003): Polarized light propagation through scattering media: Time-resolved Monte Carlo simulations and experiments, J. Biomed. Opt. 8(4): 608-617.
FURTHER READING 6 5
Wang XY, Zhang CP, Zhang LS, Qi SW, Xu T, and Tian JG (2003): Reconstruction of optical coherence tomography image based on Monte Carlo method, J. Infrared Millim. Waves 22(1): 68-70.
Wilson BC and Adam GA (1983): Monte Carlo model for the absorption and flux distri- butions of light in tissue, Med. Phys. 10: 824-830.
Wilson BC and Jacques SL (1990): Optical reflectance and transmittance of tissues: Prin- ciples and applications, IEEE J. Quantum Electron. 26: 2186-2199.
Wong BT and Menguc MP (2002): Comparison of Monte Carlo techniques to predict the propagation of a collimated beam in participating media, Num. Heat Transfer B—Fundamentals 42(2): 119-140.
Xiong GL, Xue P, Wu JG, Miao Q, Wang R, and Ji L (2005): Particle-fixed Monte Carlo model for optical coherence tomography, Opt. Express 13(6): 2182-2195.
Yadavalli VK, Russell RJ, Pishko MV, McShane MJ, and Cote GL (2005): A Monte Carlo simulation of photon propagation in fluorescent poly(ethylene glycol) hydrogel microsensors, Sensors Actuators B—Chemical 105(2): 365-377.
Yang Y, Soyemi OO, Landry MR, and Soller BR (2005): Influence of a fat layer on the near infrared spectra of human muscle: Quantitative analysis based on two-layered Monte Carlo simulations and phantom experiments, Opt. Express 13(5): 1570-1579.
Yao G and Wang LHV (1999): Monte Carlo simulation of an optical coherence tomogra- phy signal in homogeneous turbid media, Phys. Med. Biol. 44(9): 2307-2320.
Yao G and Haidekker MA (2005): Transillumination optical tomography of tissue- engineered blood vessels: A Monte Carlo simulation, Appl. Opt. 44(20): 4265-4271.
CHAPTER 4
Convolution for Broadbeam Responses
4.1. INTRODUCTION
The Monte Carlo program MCML, introduced in Chapter 3, computes responses to a pencil beam normally incident on a multilayered scattering medium. These responses are referred to as Green ‘s functions or impulse responses. When a collimated photon beam is of finite width, the Monte Carlo method is still able to compute the responses by distributing the incident positions over the cross section of the beam. Each broad beam, however, requires a new time-consuming Monte Carlo simulation, even if the other parameters are unchanged. Convolu- tion of Green’s functions for the same multilayered scattering medium, however, can efficiently compute the responses to a broad beam. Such convolution was implemented in a program named CONV (Appendix 4A). Like MCML, CONV is written in ANSI Standard C and hence can be executed on various computer platforms. Although convolution is applicable to collimated beams of any inten- sity distribution, only Gaussian and top-hat (flat-top) beams are considered in CONV version 1.
4.2. GENERAL FORMULATION OF CONVOLUTION
Convolution is applicable to a system that is stationary (time-invariant), lin- ear, and translation-invariant. The system here consists of horizontal layers of homogeneous scattering media that have stationary properties (see Figure 3.1 in Chapter 3). The input to the system is a collimated photon beam perpendicu- larly incident on the surface of the scattering medium. The responses can be any observable physical quantities, such as specific absorption, fluence, reflectance, or transmittance. The linearity implies that (1) the responses increase by the same factor if the input intensity increases by a constant factor and (2) any response to two photon beams together is the sum of the two responses to each photon beam alone. The translation invariance in space here means that if the photon beam is
Biomedical Optics: Principles and Imaging, by Lihong V. Wang and Hsin-i Wu Copyright © 2007 John Wiley & Sons, Inc.
67
68 CONVOLUTION FOR BROADBEAM RESPONSES
shifted in any horizontal direction by any distance, the responses are also shifted in the same direction by the same distance. The translation invariance in time indicates that if the photon beam is delayed by a given time, the responses are also shifted by the same delay. Therefore, responses to spatially and temporally broad beams can be computed using the convolution of the impulse responses; only spatial convolution is described in this chapter, however.
Impulse responses to a normally incident pencil beam are first computed using MCML, where a Cartesian coordinate system is set up as described in Chapter 3. The origin of the coordinate system is the incident point of the pencil beam on the surface of the scattering medium, and the z axis is along the pencil beam; hence, the xy plane is on the surface of the scattering medium.
We denote a particular response to a collimated broad photon beam as C(x, y, z) and denote the corresponding impulse response as G(x, y, z). If the broad collimated light source has intensity profile S(x, y), the response to this broad beam can be obtained through the following convolution
/
OO pOO
I G(x – xf, y – yf, z)S(x’, y) dx dy, -co J—oo
(4.1)
which can be reformulated with a change of variables x” = x — x’ and y” = y-y’:
/
CO /»OO
/ G(x’y”,z)S(x~x”,y-y”)dx”dy”. (4.2) -CO «/— OO
Because the multilayered structure has planar symmetry and the photon beam is perpendicularly incident on the surface of the scattering medium, G(x, y,z) possesses cylindrical symmetry. Consequently, the Green function in Eq. (4.1) depends only on the distance ros between the source point (xf, yf) and the obser- vation point (JC, y), rather than on their absolute locations:
ros = > / ( * – * ‘ ) 2 + ( y – / ) 2 . (4.3)
If S(x’, y’) also has cylindrical symmetry about the origin, it becomes a function of only the radius r’
r’ = vV2 + y’1. (4.4)
On the basis of these symmetries, we reformulate Eqs. (4.1) and (4.2) to
C(JC, y,z)= ί ί G (y(x-x’)2 + (y-y’)2, z) S (V*’2 + / 2 ) dx’dy
( 4
/
OO pOO / /
/ G Ux”2 + y”2, z) S U(x – x”)2 + (y- y”)2) dx”dy”.
(4.5)
(4.6)
CONVOLUTION OVER A GAUSSIAN BEAM 6 9
Because C(x, y, z) has the same cylindrical symmetry, Eqs. (4.5) and (4.6) can be rewritten in cylindrical coordinates (r, φ):
C(r, z) = I S(r’)r’ ί G (y/r2 + r’2 – Irr’costy, z) dtf dr’, (4.7)
C(r, z)= I G{r”, z)r” j S i^Jr2 + r”2 – 2rr” cos φ”) </φ” dr”.
(4.8) In Eq. (4.8), the integration over φ” is independent of z and hence needs to be computed only once for all z values. In some cases, the integration over φ” can be solved analytically; thus, the 2D integral in Eq. (4.8) is reduced to a compu- tationally more efficient ID integral. Therefore, Eq. (4.8) is more advantageous computationally than Eq. (4.7).
Example 4.1. Derive Eq. (4.7) from Eq. (4.5).
The differential area element can be changed from dxdy to rdrd§, and the corresponding limits of the integrations are from 0 to -f oo for r and from 0 to 2π for φ. In the polar coordinates aligned with (JC, y), (JC, y) is represented by (r, 0) and (x yf) by (rf, φ’). We convert Eq. (4.3) into ros = yjr2 -f r’2 — 2rrf cos φχ. Thus, Eq. (4.7) can be obtained from Eq. (4.5).
4.3. CONVOLUTION OVER A GAUSSIAN BEAM
For a Gaussian beam, the convolution can be further simplified. The intensity profile of the beam is given by
S(r’) = S0 exp <ff (4-9) Here, R denotes the l/e2 radius of the beam; So denotes the intensity at the center of the beam, and So is related to the total power PQ by
So = 2P0 π/?2 ‘
(4-10)
Substituting Eq. (4.9) into Eq. (4.8), we obtain
C(r,z) = S(r)jC°G(r”,z)exp-2(‘^
D f-2* /4rr”cos( | )” „1 „ ,, (4.11)
7 0 CONVOLUTION FOR BROADBEAM RESPONSES
The inner integral in the square brackets resembles the integral representation of the zeroth-order modified Bessel function, which is defined by
/oW 1 ί2π
2π Jo exp(;c sin φ) d<>
or
1 [2π
2π Jo εχρ0^θ8φ)<Ζφ.
(4.12)
(4.13)
By using Eq. (4.13), we can rewrite Eq. (4.11) as
ί°° {r”Y / 4 r r ” C(r,z) = 2nS(r)j G(r”,z)exp – 2 i – J 70 ί – ^ – J r”dr”. (4.14)
Example 4.2. Derive Eq. (4.10).
The total power
P0 = ί S(r’)2n/dr’ = 2nS0 f exp – 2 ( – )
– 5Ό exp (4.15)
which leads to Eq. (4.10).
Example 4.3. Show that Eqs. (4.12) and (4.13) are equivalent.
Letting φ = φ’ + (π/2) and splitting the integral in Eq. (4.12) into two parts, we obtain
/»2π /·3π/2 I exp(xsin<))d$ — / exp(Jccosφ/)ί/φ/
JO J-n/ ir/2
/
0 /*3π/2
exp(x cos φ’) d φ/ + I exp(;c cos φ’) d(>’ -π/2 JO
1/2
Letting φ’ = φ” + 2π in the first integral on the right-hand side gives
Γ·2π
/
0 /»Z71
exp(jt ^8φ’)οίφ’ = I exp(jc cos φ^) άφ”. -n/2 hn/2
(4.16)
(4.17)
CONVOLUTION OVER A TOP-HAT BEAM 71
Substituting Eq. (4.17) into Eq. (4.16) and merging the two integrals on the right- hand side of Eq. (4.16), we obtain
/»2π /·2π /»3π/2
/ expU sin φ) dφ = / exp(x cos φ”) d<>” + / exp(jt cos φ’) d§’ JO J3n/2 JO
ρ2π
=■ I βχρί^οοβφ) */φ. (4.18) Jo
Since both φ’ and φ” are dummy variables, they are both replaced with φ above.
4 A CONVOLUTION OVER A TOP-HAT BEAM
For a top-hat beam of radius R, the source function becomes
s<H° ϊ ‘»£· <4i9) where So denotes the intensity inside the beam. We have
Po So = ~ , (4.20)
nRz
where Po denotes the total power of the beam. Substituting Eq. (4.19) into Eq. (4.8), we obtain
C{r, z) = 2nSQ / G(r”, Ζ)/Φ(Γ, r”)r”‘ dr (4.21) Jo
where
1 if R > r + r”
^ c o s – ^ ^ + ^ V ^ 2 ) if V-r” < R<r + r”. (4.22) /φ(Γ,Γ”) =
0 if R < r – r”
From Eq. (4.22), the limits of integration in Eq. (4.21) can be changed to a finite range
rr+R it j / J i > ‘ j ‘ f C(r, z) = 2nS0 [ G{r” Ζ)/Φ(Γ, r”)r” dr”, (4.23)
Ja
where
a = max(0, r – R). (4.24)
Function max( ) takes on the greater of the two arguments.
7 2 CONVOLUTION FOR BROADBEAM RESPONSES
If R tends to infinity, Eq. (4.23) becomes /»OO
C(r, z) = 2nSo / G{r z)r” dr”. (4.25) Jo
This equation implies that in MCML, if the absorbed weights in all r grid ele- ments are summed and then divided by the total number of tracked photon packets, the result represents the specific absorption as a function of z for an infinitely wide beam of unit intensity.
4.5. NUMERICAL SOLUTION TO CONVOLUTION
As described in Chapter 3, a grid system is used in the Monte Carlo model. A 2D homogeneous grid system is set up in the r and z directions. The grid element sizes are Ar and Az in the r and z directions, respectively; the total numbers of grid elements in the r and z directions are Nr and Nz, respectively.
When the photon beam is Gaussian or top-hat, the 2D convolution becomes ID. Because the Monte Carlo simulation assigns physical quantities to discrete grid elements, an appropriate integration algorithm is based on the extended trapezoidal rule. This algorithm is ideal for a nonsmooth integrand that is linearly interpolated between available data points; it is implemented in C as a function named t r a p z d ( ) , which is called by another function named q t r a p ( ) .
Another method of integration is to evaluate the integrand at the original grid points. This approach, however, does not offer any control over the integration accuracy. For a top-hat beam, for example, Nr is 50, and R is about 5Ar. If C(0, z) is computed from Eq. (4.23), the integration interval [0, R] covers only 5Ar. Thus, only five function evaluations contribute to the integration and may yield unacceptable accuracy. By contrast, the extended trapezoidal rule continues to perform function evaluations until a user-specified accuracy is reached.
The sequence of integrand evaluations in the extended trapezoidal rule is shown in Figure 4.1a. Subsequent calls to t r apzd( ) incorporate the previous evaluations and evaluate the integrand only at the new points. To integrate f(x) over interval [a,b], we evaluate f(a) and f(b) in the first step as noted by 1 and 2 in Figure 4.1a. To refine the grid, we evaluate / ( | (a + b)) in the second step as noted by 3. This process is repeated until the integral evaluation reaches a specified accuracy. The bottom line shows all function evaluations after four calls.
4.5.1. Interpolation and Extrapolation of Physical Quantities
The physical quantities under discussion have been computed using MCML over a grid system. As discussed in Chapter 3, the optimal r coordinate is
r(ir) = 2 / l2(/r + ±)
Ar, (4.26)
NUMERICAL SOLUTION TO CONVOLUTION
(a)
O
Original data Interpolation Extrapolation
‘ I ‘ I ‘ I ‘ I ‘ I ” I ‘ I 1 2 3 4 5 6 7
(b) r/Ar N -0.5
Figure 4.1. (a) Sequence of integrand evaluations in the extended trapezoidal rule of integration; (b) interpolation and extrapolation of the physical quantities. In this example, Nr = 8. Symbols a and b denote the integration limits, and i denotes the iteration index. Arrows point to where the integrand is evaluated. Solid circles represent the original data points. Dashed and solid lines represent linear interpolation and extrapolation, respectively.
where ir is the index of the grid element (0 < ir < Nr — 1). For ir = 0, r(0) is | Ar instead of ^Ar. The offset between the optimized and the centered coordi- nates in each grid element decreases as ir increases.
In qt rap (), the integrand—of which G is only a part—is evaluated at points that may not fall on the original grid as illustrated in Figure 4.1. Linear interpo- lations are used for those points that fall between two original grid points, and linear extrapolations are used for those points that are located beyond the origi- nal grid system (Figure 4.1b). Extrapolation is extended only up to (Nr – 0.5) Ar because further extrapolation is unreliable. In MCML, the last cells in the r direc- tion are used to collect contributions from photon packets that do not fit into the
74 CONVOLUTION FOR BROADBEAM RESPONSES
grid system and thus do not represent the true local physical quantities. There- fore, the upper limit for extrapolation is (Nr — 0.5)Ar instead of (Nr + 0.5)Ar, and the physical quantity beyond (Nr — 0.5) Ar is set to zero. We denote
rmax = (Nr-0.5)Ar. (4.27)
4.5.2. Integrand Evaluation for a Gaussian Beam
Although the integration in Eq. (4.14) must converge for physical reasons, it may cause overflow in a computer because the modified Bessel function increases rapidly as the argument increases. Therefore, a proper reformulation of the inte- grand is necessary.
We note that the modified Bessel function has the following asymptotic approx- imation for large x values:
expOc) / o W « – ^ . (4.28)
y/2nx
We define the following new function on the basis of IQ
Ioe(x) = Io(x)exp(-x) (4.29)
or
/o(x) = /o*U) expU), (4.30)
where IQ€ is always well bounded. Substituting Eqs. (4.9) and (4.30) into Eq. (4.14), we obtain
poo
C(r, z) = 2π5ο / G(r”, z) exp Jo
Because both the exponential and the I$e terms are well bounded, the integrand can be computed without overflow. Since Eq. (4.28) is not actually used in the computation, Eq. (4.31) does not carry any asymptotic approximation.
Computation speed is another issue. The evaluation of exp( )Ioe() in Eq. (4.31) is a major part of the computation, which can take up to 90% of the total time. For multi-dimensional physical quantities (e.g., the fluence as a function of r and z), the convolution may repeatedly evaluate exp( )Ioe( ) at the same r coordi- nate as the integration is computed for different z coordinates. Therefore, if the values of exp( )/o*( ) are stored during the convolution for one z coordinate, computation time can be reduced. Because the number of function evaluations is unknown in advance, the function evaluations should be saved with dynamic data allocation. Since the evaluation sequence in q t r a p O resembles a binary tree as shown in Figure 4.1a, a binary tree can be used to store the function evaluations. Although the first two nodes are out of balance, the subtree below node 3 is perfectly balanced.
^ Ι Χ Ϊ Κ «->
NUMERICAL SOLUTION TO CONVOLUTION 7 5
4.5.3. Integration Limits for a Gaussian Beam
Although the upper integration limit for a Gaussian beam in Eq. (4.31) is infinity, it can be converted to a finite value by a change of variables; then, the integration can be computed by a function called midexp(). This approach, however, is not computationally efficient. Therefore, q t r a p ( ) is preferred; the upper integration limit, however, must be reduced to a finite value. To this end, the integrand is nonzero if
| r ” – r | <KR (4.32)
or
r-KR<r” <r + KR, (4.33)
where A’ is a constant that can be set in CONV. For example, if K is 4, the exponential term in Eq. (4.31) is exp(—32) ~ 10~~14.
The computation of G covers only interval [0, rm a x] , where rmax is as given by Eq. (4.27). Combining this limit and Eq. (4.33), we rewrite Eq. (4.31) as
fb Γ / r ” – r 2 1 /4rr” C(r,z) = 2nS0j G(r”,z)exp I —2 f—^—1 l0ei ~^- r”dr”% (4.34)
where
a = max(0, r – KR), (4.35)
fc = min(rmax,r + KR). (4.36)
Functions max( ) and min( ) take on the greater and the lesser of the two arguments, respectively.
4.5.4. Integration for a Top-Hat Beam
The integrand for a top-hat beam can be evaluated more easily than that for a Gaussian beam, because the integration limits are finite and the integrand causes no overflow. However, evaluation of /ψ in Eq. (4.23) is also time-consuming. As in the integrand evaluation for Gaussian beams, the evaluated /φ values for q t r a p ( ) are stored in a binary tree for computational efficiency.
Since the physical quantities are computed only in interval [0, rm a x] , Eq. (4.23) can be expressed as
C(r, z) = 2JZSO f G{r Ζ)/Φ(Γ, r”)r” dr”, (4.37) Ja
76 CONVOLUTION FOR BROADBEAM RESPONSES
where
a = max(0, r – /?), (4.38)
fc = min(rmax,r + /?). (4.39)
4.5.5. First Interactions
In MCML, absorption from the first photon-tissue interactions is recorded sep- arately. The first interactions always occur on the z axis and hence contribute to the specific absorption or related physical quantities as a delta function. The total impulse response can be expressed in two parts
G(r, z) = G,(0, z)^- + G2(r, z), (4.40)
where the first term results from the first interactions and the second, from sub- sequent interactions.
For a Gaussian beam, substituting Eq. (4.40) into Eq. (4.34) yields
C(r,z) = G{(0,z)S(r) + 2nS0 f G2{r”,z) Ja
– 2 (̂ Χ )̂’*-· x exp For a top-hat beam, substituting Eq. (4.40) into Eq. (4.37) yields
C(r, z) = Gi(0, z)S(r) + 2nS0 f G2(r Ζ)/Φ(Γ, r”)r” dr”. (4.42) Ja
The numerical results obtained with and without separately recording the first interactions are compared in the next section.
4.5.6. Truncation Error in Convolution
As shown in Eqs. (4.36) and (4.39), the upper integration limits may be bounded by rmax. For a top-hat beam, if
r < W – R, (4.43)
the limited grid coverage in the r direction does not affect the convolution; otherwise, it truncates the convolution and leads to error in the convolution for r > ”max — R- Thus, to convolve reliably for physical quantities at r in response to a top-hat beam, we must ensure that rmax in the Monte Carlo simulation is large enough that Eq. (4.43) holds.
APPENDIX 4A. SUMMARY OF CONV 7 7
For a Gaussian beam, no simple formula similar to Eq. (4.43) exists because a Gaussian beam theoretically extends to infinity. At r ^> R, a Gaussian beam and a top-hat beam of the same R and So have comparable convolution results. Therefore, Eq. (4.43) can be used approximately for Gaussian beams as well.
4.6. COMPUTATIONAL EXAMPLES
In this section, the error that is caused by not recording the first photon-matter interactions separately is illustrated, and a numerical example of convolution is presented. The impulse responses are computed by MCML for a scattering medium described in Table 4.1. The grid element sizes in the r and z directions are both 0.01 cm. The numbers of grid elements in the r and z directions are 50 and 40, respectively. One million photon packets are tracked.
The impulse fluence near the surface of the scattering medium (z = 0.005 cm) is shown in Figure 4.2a, where the first interactions are recorded separately. If they were recorded into the first r grid element instead, it would augment the fluence in the first grid element by 1.95 x 103 cm- 2 , which is significantly greater than the current value of 1.34 x 103 cm- 2 . For comparison, the impulse response is convolved over a top-hat beam of 1-nJ energy and 0.01-cm radius both with, and without, recording the first interactions separately (Figure 4.2b). The convolved results differ at r = 0.015 cm by as much as 120%.
The impulse response is also convolved over a Gaussian beam (1-nJ total energy, 0.1-cm radius), where the convolution error is set to 0.01. The contour lines of the fluence distribution before and after the convolution are shown in Figure 4.3.
APPENDIX 4A. SUMMARY OF CONV
The entire source code of CONV can be found on the Web at ftp://ftp.wuey. com/public/sci_tech_med/biomedical_optics. The program is divided into several files. Header file conv.h defines the data structures and some constants. File convmain.c contains primarily the main function. File convi .c handles data reading. File convo.c handles data writing. File convconv.c implements the
TABLE 4.1. Optical Properties and Structure of a Three-Layered Scattering Medium0
Layer
1 2 3
n
1.37 1.37 1.37
μ« (cm”1)
1.0 1.0 2.0
μ* f (cm” 1)
100.0 10.0 10.0
8
0.9 0 0.7
Thickness (cm)
0.1 0.1 0.2
aRefractive indices for top and bottom ambient media are both 1.0.
78 CONVOLUTION FOR BROADBEAM RESPONSES
0.2 0.3 Radius r (cm)
Figure 4.2. (a) Relative fluence at z = 0.005 cm in response to a pencil beam computed by MCML; (b) fluence at z = 0.005 cm in response to a top-hat beam computed by CONV. Circles and crosses represent data with and without, respectively, the first inter- actions scored separately.
actual convolution. File conviso .c handles calculation of contours. File con- vnr . c contains several functions for dynamical data allocations and error reports. Readers should read the main function first.
The following list is generated by command cf low -d4 -n – -omit-arguments – -omit – symbol- names conv* . c, which shows the structure of the program with
APPENDIX 4A. SUMMARY OF CONV 7 9
(a) 0.2 0.3 Radius r (cm)
0.2 0.3 Radius r (cm)
Figure 4.3. (a) Relative fluence distribution in response to a pencil beam computed by MCML; (b) fluence distribution in response to a Gaussian beam computed by CONV.
the nesting depth limited to 4:
1 main() <int () at C0NVMAIN.C:143>: 2 ShowVersionO <void () at C0NV0.C:37>: 3 pu t s ( ) 4 Cen te rS t r ( ) <char * () a t C0NV0.C:11>: 5 s t r len ( ) 6 strcpy()
8 0 CONVOLUTION FOR BROADBEAM RESPONSES
7 s t r c a t ( ) 8 p r i n t f ( ) 9 ge ts ( )
10 s t r l e n ( ) 11 BranchMainCmd() <void () at CONVMAIN.C:122>: 12 strlen() 13 BranchMaihCmdl() <void () at C0NVMAIN.C:63>: 14 toupper() 15 ReadMcoFileO <void () at CONVI.C:568>: 16 LaserBeam() <void () at C0NVC0NV.C:92>: 17 ConvResolution() <void () at C0NVC0NV.C:127>: 18 ConvError() <void () at CONVCONV.C:156>: 19 ShowMainMenu() <void () at CONVMAIN.C:25>: 20 QuitProgram() <void () at C0NVMAIN.C:44>: 21 puts() 22 BranchMainCmd2() <void () at C0NVMAIN.C:92>: 23 toupper() 24 0u tpu t0 r i gDa ta ( ) <void () at C0NV0.C:784>: 25 OutputConvData() <void () at CONVCONV.C:1017>: 26 ContourOrigData() <void () at C0NV0.C:893>: 27 ContourConvData() <void () at CONVCONV.C:1137>: 28 ScanOrigData() <void () a t CONVO.C:1211>: 29 ScanConvData() <void () at CONVCONV.C:1506>: 30 pu ts ( ) 31 pu ts ( )
PROBLEMS
4.1 Derive Eq. (4.11).
4.2 Derive Eq. (4.21).
4.3 Derive Eq. (4.23).
4.4 Derive Eq. (4.25).
4.5 Derive Eq. (4.31).
4.6 Derive Eqs. (4.41) and (4.42).
4.7 Write a computer program that can convolve the impulse responses over a flat beam.
4.8 Take the Fourier transformation of Eqs. (4.1) and (4.2) with respect to x and y.
4.9 Write the time-domain counterparts of Eqs. (4.1) and (4.2) for an arbi- trary pulse profile S(t) of the incident pencil beam. In this case, an
FURTHER READING 81
impulse response G(x, y, z, /) is defined as a time-resolved response to a temporally infinitely short δ(ί) photon beam. Assuming the response to S(t) to be experimentally measured, explain how to recover the impulse response G(JC, y, z, t) through deconvolution.
4.10 Take the Fourier transformation of the time-domain counterparts of Eqs. (4.1) and (4.2) with respect to t.
4.11 Assume that the incident photon beam is finite in both time and space and can be represented by S(x, y,z,t). Write the convolution over this beam.
Some Matlab Task With Biomedical Science
"Looking for a Similar Assignment? Get Expert Help at an Amazing Discount!"
