Some Matlab Task With Biomedical Science

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.

BIOMEDICAL OPTICS

Save your time - order a paper!

Get your paper written from scratch within the tight deadline. Our service is a reliable solution to all your troubles. Place an order on any task and we will take care of it. You won’t have to worry about the quality and deadlines

Order Paper Now

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

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

http://www.copyright.com

http://www.wiley.com/go/permission

http://www.wiley.com

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 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.

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

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

—

— — —

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)

0.8

0.6

0.4

0.2

–

—r

J – ** S

__ —

/ /

™ “” ” Lol’ ESF 1

, 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)

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

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

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

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

c 3

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)

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.)

Some Matlab Task With Biomedical Science

Save your time - order a paper!

"Looking for a Similar Assignment? Get Expert Help at an Amazing Discount!"

About Us

Writing Services

Contact Us

Talk to us