An Efficient Approach for Optical Radiative Transfer Tomography using the Spherical Harmonics Discrete Ordinates Method

0.75 MB·English

by  Aviad Levis

An Efficient Approach for Optical Radiative
Transfer Tomography using the Spherical
5
1
Harmonics Discrete Ordinates Method
0
2

n AviadLevis1,YoavY.Schechner1,AmitAides1,AnthonyB.Davis2
a
J

4
2
1ELECTRICALENGINEERINGDEPARTMENT, TECHNION-ISRAELINSTITUTEOFTECHNOLOGY,

HAIFA32000,ISRAEL
]
h 2 JET PROPULSION LABORATORY, CALIFORNIA INSTITUTE OF TECHNOLOGY, PASADENA, CA
p 91109,USA
-
o
Abstract
a
.
s Thispaperintroducesamethodtopreformopticaltomography,using3Dra-
c diativetransferastheforwardmodel. Weuseaniterativeapproachpredicatedon
i
s theSphericalHarmonicsDiscreteOrdinatesMethod(SHDOM)tosolvetheopti-
y mizationprobleminascalablemanner.Weillustratewithanapplicationinremote
h sensingofacloudyatmosphere.
p
[

1 Introduction
1
v
3 Opticaltomographyisa3Dimagingtechniquethatusesopticalmeasurementsonthe
9 boundaryof a domainto findthe spatialdistribution ofparameters within. It hasnu-
0
merous applications in bio-medical imaging. For a list of applications we refer the
6
readerto[1,2]andreferencestherein. Insharpcontrast,operationalremotesensingof
0
. cloudsintheEarth’satmosphereusingspace-borneopticalsensorsispredicatedonthe
1
classic plane-parallel medium geometry with horizontally uniform properties, hence
0
5 the1Dradiativetransferequation(RTE);thisreducestheproblemtoestimatingasin-
1 glecloudopticalthickness,andshiftsthefocustowardmicrophysicalproperties[3–5].
: Someofthebiasesinretrievedcloudpropertiescausedbythissometimesverycrude
v
i approximationofrealcloudstructurearedocumentedin[6].
X
Solvingtheinverseproblemusingthefull3DRTEasaforwardmodelcanbedif-
r ficultandcomputationallydemanding. Forsomeapplications, itispossibletousean
a
approximatemodel. Indensemedia(meanfreepathsmallcomparedtheoverallsize),
with scatteringdominating absorption, it is possibleto use adiffusion approximation
of the 3D RTE. This results in the inverse problem of Diffuse Optical Tomography
(DOT) [7–9]. When the mean free path is large compared to the outer size of the
medium, the measured radiant energy is dominated by direct and single scattered in-
tensities. The resulting inverse problem is single scattering tomography [2,10,11].
1
I
incident
s
I(x,ω)
I
∂Ω s
ϑ
∂Ω
Ω Ω
w
(a) (b)
Figure1: (a)Domainandboundaryconditions;(b)Aperturefunction. Bluemarksthe
spatialsupportandredmarkstheangularsupport
Otherapproximationsandtheirderivationscanbefoundin[12]. Weproposetosolve
the inverse problem using the RTE, with neither single scattering, nor diffusion ap-
proximations. However,tomaketheinverseproblemtractable,wederiveaniterative
optimizationframeworkthatusesdifferentpartsoftheSphericalHarmonicsDiscrete
OrdinatesMethod(SHDOM)3DRTEsolver[13].
2 Theoretical background
2.1 Radiativetransfer
Ourforwardmodelissteady-state3Dradiativetransfer. Thismodelisusedinpassive
imaging when source gating is sufficiently slow. The 3D RTE governs the transport
oflightthroughamediumwithdistributionofscatteringand/orabsorption. Consider
a domain Ω having boundary ∂Ω whose outward facing normal is ϑ (Fig. 1a). The
transport domain is indicated by position x ∈ R3 and direction of propagation ω ∈
S2 (unit sphere). The radiation (light) field is I (x,ω). The subscript λ indicates
λ
wavelengthdependencyandwillbeomittedfromhereon,forsimplicity.Theboundary
condition(Fig.1a)is
I(x,ω)=I (x,ω) if ω·ϑ<0, x∈∂Ω, (1)
incident
whereω·ϑ<0definesincomingradiation,ϑbeingtheoutgoingnormalto∂Ω. The
transportoflightisformulatedintermsofthefollowingconservationlaw[14,15]:
ω·(cid:79)I(x,ω)=β(x)[J(x,ω)−I(x,ω)] x∈Ω, (2)
where β is the extinction coefficient. Here J(x,ω) is the source function. In the
absence of emission within the medium, it is also known as the in-scattering term,
2
since it accounts of an increase of radiation due to in-scattering. In this case, it is
expressedas
(cid:36) (cid:90)
J(x,ω)= p(ω·ω(cid:48))I(x,ω(cid:48))dω(cid:48), (3)
4π
S2
where(cid:36)isthesinglescatteringalbedoandp(ω·ω(cid:48))1isthephasefunction.Thephase
function describes the probability of a photon traveling in direction ω(cid:48) into scatter to
directionω. Thephasefunctionsatisfiesanormalizationcondition
1 (cid:90)
p(ω·ω(cid:48))dω(cid:48) =1. (4)
4π
S2
Eqs.(2)–(3),andtheboundarycondition(1),defineacompleteradiativetransferfor-
wardmodelforanexternallyilluminatedmedium.
Integrating Eq. (2) along a specific direction ω results in the integral form of the
3DRTE[14,16]
 (cid:90)x0 
I(x,ω)=I0exp− β(r)dr
x
(cid:90)x0  (cid:90)x(cid:48) 
+ J(x(cid:48),ω)β(x(cid:48))exp− β(r)drdx(cid:48), (5)
x x
bearing in mind the expression of J(x,ω) in terms of I(x,ω) in (3). Here x is
0
a point on the boundary and I holds the boundary condition (1). The operation we
0
denoteas
x(cid:48)
(cid:90)
f(r)dr (6)
x
issimplya1D(line)integraloverafieldf(x)alongthesegmentextendingfromxto
x(cid:48). Numerically,thisispreformedbyback-projecting(BP)araythroughthemedium.
TheexpressionforI(x,ω)in(5),foragivenJ(x,ω),isoftenreferredtoasthe
up-windsweepofradiantenergysources. WeadoptedtheSphericalHarmonicsDis-
creteOrdinatesMethod(SHDOM)[13,17]fortheforward3DRT,whichisavailableas
open-sourcecode. InSHDOM,(5)isindeedusedinastraightforwardpost-processing
stage for estimating the observed quantity I(x,ω) after a nontrivial numerical pro-
cedure is performed to obtain J(x,ω) throughout the grid. This staged approach is
exploitedfurtheron. Wenowformalizeourforwardmodelinoperatorlanguage.
2.2 OperatorNotation
We follow the definitions of [2] to formulate the forward model using operator nota-
tions (Table 1 summarizes the notations). Denote by Θ the space of extinction fields
1Inprinciple,thephasefunctiondependsonbothincomingandoutgoingdirections(ω,ω(cid:48)). However,
dependencyisoftenassumedtobesolelyonthescatteringangle,equivalenttoω·ω(cid:48).
3
Operator Symbol Domain Range
Radianceforwardmapping I Θ-Extinctionfields Z -Radiancefields
Measurement M Z -Radiancefields Y -Measurements
Forwardoperator F Θ-Extinctionfields Y -Measurements
Scatterforwardmapping J Θ-Extinctionfields V -Scatterfields
Transformationoperator T V -Scatterfields Z -Radiancefields
Table1: Operatornotationsummary
over the domain Ω. Let β(x) ∈ Θ be an extinction field over the domain Ω. Let S
representasetofradiationsourcesover∂Ω.Foragivensources∈S,denoteZ asthe
s
setofallpossibleradiationfieldsthatsatisfyEqs. (1,2)acrossallpossibleβ ∈Θ. The
setZ isinfinite,sinceforeachextinctionfieldβ(x)thereisacorrespondingradiation
s
field I = I (β). More generally the radiance forward mapping, I : Θ → Z, is a
s s
(cid:83)
mappingfromtheopticalparameterdomain, totheradiancefieldrangeZ = Z .
s s
MeasurementsareanoperatorM : Z → Y,mappingthecontinuousfunctionspace
toavectorspaceY. ConsequentlytheforwardoperatorF isdefinedas
F =MI :Θ→Y. (7)
ThemeasurementoperatorMisdefinedbythedetector’saperturefunction,w ∈Ω×S2
(Fig.1b). Theaperturefunctiondefinesthemannerinwhichthedetectorcollectsradi-
ance,overaspatialandangularsupport. Agivenapertureandsourcepair(w,s)yields
anelement
y =F (β)=M I (β)=(cid:104)w,I (cid:105) , (8)
w,s w,s w s s Ω
foraspecificextinctionfieldβ,where
(cid:90) (cid:90)
(cid:104)·, ·(cid:105) ≡ ··dωdx. (9)
Ω
Ω S2
Foranidealizedsingle-pixeldetectorpositionedatx∗,collectingradiationflowingin
directionω∗,
y =(cid:104)δ(x−x∗)δ(ω−ω∗) , I (cid:105) =I (x∗,ω∗). (10)
w,s s Ω s
WecandefineamatrixY whoseelementsy correspondtodifferentsource-detector
w,s
configurations. AcolumnofY representsmeasurementsbyasingledetector,formul-
tiplesources. Arowrepresentsmeasurementsbymultipledetectors,foroneparticular
lightsource. ThevectoryisthecolumnstackofY.
2.3 Opticaltomography
UsingtheoperatorsdefinedinSec.2.2weexpresstheprocessofopticaltomography.
Tomographicreconstructionisanestimatorofβ thatminimizesadefinedcost
βˆ=argmin{E[y,F(β)]+αΨ(β)}, (11)
β
4
θ(x)
b(x)
1 b2 b b4
b 3
1
x x
x x x x
1 2 3 4
Figure2: Illustrationofa1Dlinearinterpolationkernel
whereE[y,F(β)]isthedatafit(fidelity)functionalandΨ(β)isaregularizationon
theopticalparameters. Hereαisatunableparameter, choseninaccordancewiththe
noiselevel,tobalancethetwoterms.
Exceptforafewspecialgeometricconfigurations[2,18],3Dtomographicrecovery
cannotbedoneanalytically. Thus,thecontinuousfunctionisoftendiscretized(Fig.2)
Ngrid
(cid:88)
β(x)= β b (x), (12)
k k
k=1
where{β }Ngrid arediscreteparameters,b (x)isaninterpolationkernelandN is
k k=1 k grid
thenumberofgridpointsusedfordiscretization. Weapply(12)toEq.(11)toseekan
estimatorforβ =(cid:0)β ,...,β (cid:1)T,where(·)T denotestransposition.
1 Ngrid
We focus here on gradient-based optimization methods. The forward operator F
is a non-linear function of the optical parameters. Thus, one approach [2] estimates
(11)bylinearizingF. Supposethesolutionisβδ = β +δβ,whichisaperturbation
0
ofaninitialguessβ . ItispossibletolinearizeF andsolveδβ usingalinearsetfor
0
equations. However, this approach requires an initial guess very close to the true so-
lution. Anotherapproach[19–22]iterativelyestimatesthegradientwithrespecttoβ.
However, this approach has high computational complexity, as it requires O(N )
grid
simulations of the forward model per iteration (i.e., for a single computation of the
gradient). Forwardmodelsimulationistimeconsumingparticularlyinthepresenceof
multiplescattering. Simulatingtheforwardmodelperelementofthegradientdoesnot
scalewellasN increases.
grid
Our approach is also iterative, however each iteration is computationally simple,
having run-time independent of N . The approach does not directly optimize the
grid
nonlinearF. Instead,I isdecomposedintoaproductoftwooperators. Eachisopti-
mizedintime. Thisconceptisanalogoustoexpectationminimization(EM)optimiza-
tion,whichisusedinDOT[23–25].
5
I (x ,ω∗)
s 0
∂Ω
J (x,ω∗)
s Ω
(cid:96)
w
x∗
(a) (b)
Figure3: (a)Alineintegraloverthein-scatterfieldJ . (b)Integrationoverthefieldf
s
foranarbitraryaperturespatialandangularsupportw
3 Optimization Approach
3.1 Forwardmapdecomposition
In Sec. 2 the forward operator was defined in terms of the radiance field I. We now
show that defining F in terms of the in-scatter field J considerably speeds up the
computation time. We decompose the radiance forward mapping into two operators,
I =TJ,whichwenowintroduce.
Foraparticularradiationsources ∈ S,letJ =J (β)bethein-scatterfieldthat
s s
satisfiesEq.(3),whileI satisfiesEqs. (1,2). DefineV asthesetofpossiblein-scatter
s s
fields for a particular source. The in-scatter forward mapping, J : Θ → V, is a
(cid:83)
mapping from the optical parameter domain to the in-scatter field range V = V .
s s
Wedefinethefollowingtransformationoperator
T :V →Z, I =T (β)J . (13)
s s
Eq. (3)defines a relationbetween a givenlight field I and acorresponding in-scatter
fieldJ.ThetransformationoperatorT,definedbyEq.(5),definestheinverserelation:
for a given in-scatter field J, a corresponding light field I is attained. Consequently,
wedefinetheforwardoperatorintermsofthein-scatterforwardmapping
F =MT (β)J (β):Θ→Y. (14)
Define(cid:96)asalineofsightfromx∗indirection−ω∗(Fig.3a). Foraconvexspatialdo-
main,theintersectionpointbetweenlineofsightandthedomain’sboundaryisunique
anddenotedbyx =∂Ω∩(cid:96)(Fig.3a).
0
OperatorT isdefinedbyEq.(5). Forsources ∈ S, theresponseofanidealized
6
E[y,F(θ)] E[y,F(θ)]
E[y,F (θ)]
1
E[y,F (θ)]
0
E[y,F (θ)]
0
θ θ
θ θ θ θ θ θ θ θ
0 min min 0 0 1 1 0
(a) (b)
Figure4: (a)Gradientdecentoptimization. (b)Surrogatefunctioniterativeoptimiza-
tion. Bothdependontheinitialguessβ
0
singlepixelwithaperturew =δ(x−x∗)δ(ω−ω∗)is
y =M TJ
w,s w s
 (cid:90)x0  (cid:90)x0  (cid:90)x 
=Is(x0,ω∗)exp− β(r)dr+ Js(x,ω∗)β(x)exp− β(r)drdx
x∗ x∗ x∗
(15)
Eq. (15) is simply an accumulation of the scattered radiance along the line of sight
weighted by its corresponding attenuation factor (Fig. 3a). Applying a general mea-
surementoperationisamannerofintegratingovertheaperturefunctionofthedetector
(Fig.3b).
3.2 Iterativesurrogatefunctionoptimization
OptimizingEq.(11)directlyoverβiscomputationallyexpensive.Findingthegradient
direction in each iteration requires O(N ) numerical computations of the forward
grid
model. With our approach, however, we iteratively keep J constant and estimate β
s
solelybasedonT.Insteadofoptimizingthroughafunctionthatisdifficulttocompute
(Fig.4b)weoptimizeβusingasurrogatefunction[26],whichisefficientlycomputed
(Fig. 4a). We define the following iterative optimization process. Define β as an
n
estimate of β in the “n” iteration. The first step consists of computing the in-scatter
field,whichcorrespondstothecurrentestimateβ
n
J =J (β ).
n n
In the second step, keeping the in-scatter field constant, we solve the following opti-
mizationproblemtofindournextestimateofβ
β =argmin{E[y,F (β)]+αΨ(β)}, with (16)
n+1 n
β
E[y,F (β)]=[y−MT (β)J ]T Σ−1 [y−MT (β)J ],
n n meas n
7
(a) (b)
Figure5: (a)Eachgradientelementiscomputedbytracingalltherays. Onlythered
rayscontributetothecomputationoftheredgridpointgradient. (b)Asingletraceof
eachrayisperformed. Whiletracingeachraywecomputethegradientcontributionto
eachrelevantgridpoint
where
F (β)=MT (β)J , (17)
n n
isthenthsurrogatefunction,andΣ isthecovariancematrixofourmeasurements.
meas
Intheremainder,wewillnotmakeuseofregularization(setα=0).
Eq.(17)definesaniterativeoptimizationprocess,
i. Startwithaninitialguessβ .
0
ii. Basedonthecurrentestimateβ ,numericallyfindthein-scatterfieldJ =J (β ).
n n n
iii. Optimize(17)tofindthenextestimateβ .
n+1
iv. Repeatstepsii-iiiuntilconvergence.
3.3 Scalability
For a given measurement vector y ∈ Y, we solve the optimization of (17) with a
gradient-basedmethod. Weusethediscretizationdescribedin(12). Assuminguncor-
relatedmeasurements
Σ =diag(σ2 ), σ2 =(cid:0)σ2,...,σ2 (cid:1)T , (18)
meas meas meas 1 Nmeas
where N = N ×N is the number of measurements, N being the number of
meas w s w
detectors and N the number of sources. In this study based on radiometry, we use
s
σ =0.03×y (a3%noiseduemostlytocalibrationerror).
i i
The task is thus to recover the gridded extinction {β }Ngrid. For this purpose we
k k=1
findthek’thgradientelement. Withoutlossofgeneralityweformulatetheproblemin
termsofasinglesourceandmanydetectors.
∂ N(cid:88)meas 1 (cid:20) ∂ (cid:21)
E[y,F (β)]= [F (β)−y ]M T (β) J , (19)
∂β n σ2 w,n w w ∂β n
k w=1 w k
8
For anideal single-pixel detector positionedat x∗, collectingradiation flowing indi-
rectionω∗weget
(cid:20) (cid:21)
∂
M T (β) J =A +B (20)
w ∂β n w,k w,k
k
where
 (cid:90)x0   (cid:90)x0 
Aw,k =− bk(r)drI(x0,ω∗)exp− β(r)dr (21)
x∗ x∗
and
(cid:90)x0  (cid:90)x   (cid:90)x 
Bw,k = Jn(x,ω∗)bk(x)−β(x) bk(r)drexp− β(r)drdx. (22)
x∗ x∗ x∗
Term A results from the boundary illumination and can be readily computed. Term
B defines a line integral over a field, computed using back-projection of rays from
thedetectorthroughthemedium. Astraightforwardapproachcalculateseachelement
of the gradient vector by summing all the back-projected rays from all the detectors
(Alg.2a). ThecomplexityofthisapproachisO(N ×N )back-projectingop-
meas grid
erations. However,mostofback-projectedraysdonotcontributetoB sincetheinter-
polationkernelb (x)typicallyhasasmallsupportregion(Fig.5a). Instead,foreach
k
point along a ray we define a support region of grid points, according to the support
regionoftheinterpolationkernel(Fig.5b).
A specific ray will contribute to the computation of gradient elements associated
with the grid points in the support region. In the case of a linear interpolation ker-
nel,asupportregioniscomposedoftheeightcornergridpointsthatmakeagridcell.
Sincethesupportregionistypicallysmall,foreachpointalongtherayweonlyneedto
lookatafiniteamountofneighboringgridpoints(Alg.2b). Thecomplexityofthisap-
proachisO(N )back-projectingoperations.Hence,thecomputationalcomplexity
meas
isindependentofthenumberofgridpointswhichmakesthisapproachscalable.
4 Application to remote sensing
4.1 Generalsetup
Weuselargeeddysimulation(LES)[27,28]togeneratethemicrophysicalquantitiesof
arealisticcloudfield(Fig.6). WeuseMiescatteringtheory(implementedintheavail-
ablecode[29])toconvertthemicrophysicalquantities(liquidwatercontentfromLES,
assumed effective droplet radius r and effective variance of the droplet size distri-
eff
butionv )intotherequiredopticalquantities(extinctioncoefficient,phasefunction,
eff
singlescatteringalbedo). Wetaker =10µmandv =0.1forGamma-distributed
eff eff
dropletsizes,andthewavelengthissetto672nm. Finally,weaddtheextinctionand
phasefunctionduetoRayleighscatteringbyairmolecules;aerosolsareignoredinthis
simulation.
9
1: fork =1→Ngriddo 1: Grad=0
2: Gradk =0 2: forw =1→Nmeasdo
3: forw =1→Nmeasdo 3: Err= σ12 [Fn(β)−yw]
4: BPtocomputeAw,k,Bw,k 4: fork ∈Swupportregiondo
5: Err= σ12 [Fn(β)−yw] 5: BPtocomputeAw,k,Bw,k
w
6: Gradk =Gradk + 6: Gradk =Gradk +
+Err(A +B ) +Err(A +B )
w,k w,k w,k w,k
7: endfor 7: endfor
8: endfor 8: endfor
(a)Straightforwardapproach (b)Scalableapproach
Alg. 2: Comparison of the two approaches to compute the gradient of the surrogate
(cid:0) (cid:1)
function. Grad= Grad ,...,Grad
1 Ngrid
Weseektorecovertheextinctionoftheclouddropletsonacartesiangrid, where
theairistakentobeaknownparameterdependentonlyontheheight. Theboundary
conditionsforthedomainare
• Collimatedsolarradiationatthetopoftheatmosphere(TOA),atazenithangle
of60◦.
• Openboundaryforthesidefaces,toapproximateanisolatedcloudobservedata
wavelengthwheretheRayleighskyisalmostblack(Rayleighopticalthickness
oftheatmosphereisonly0.04532at672nm).
• Lambertian reflectance at the surface (Earth) with an albedo of 0.05, to mimic
thedarkoceansurface.
Atλ=672nm,liquidwaterandairarenon-absorbing,sothesinglescatteringalbedo
(cid:36)iseverywhereunity. Inordertofindthein-scatterfieldJ ,weuseSHDOM[13].
n
4.2 SHDOM
Running the forward model is a balance between speed and accuracy. Monte Carlo
methods can handle very complex optical media. However, they compute radiomet-
ric quantities by random sampling the infinite domain of possible light paths. This
introduces stochastic noise at the output, which can be controlled by increasing the
numberofsamples(photonpaths)andvariancereductiontechniques. Whenmanyra-
diometric outputs of the same scene are sought, as in the case of many viewpoints, a
model that solves the RTE directly has a much more favorable runtime [32]. A dis-
crete ordinates representation models the flow of radiant energy in the domain easily
andintuitively[13,17,33,34]. TheSHDOMmodelusesagridforthespatialdepen-
dency,andalinearinterpolationkernelb (x). SphericalHarmonicexpansion[12,35]
k
is used to compute angular integrals. SHDOM solves the forward model 3D RTE in
an integral form for the source function J(x,ω). It can adaptively truncate negligi-
blecoefficientsintheseriesexpansion. Italsohasadaptivespatialgridcapabilityand
multi-gridacceleration.
10