OSA's Digital Library

Optics Express

Optics Express

  • Editor: Andrew M. Weiner
  • Vol. 21, Iss. 22 — Nov. 4, 2013
  • pp: 25820–25833
« Show journal navigation

Multi sky-view 3D aerosol distribution recovery

Amit Aides, Yoav Y. Schechner, Vadim Holodovsky, Michael J. Garay, and Anthony B. Davis  »View Author Affiliations


Optics Express, Vol. 21, Issue 22, pp. 25820-25833 (2013)
http://dx.doi.org/10.1364/OE.21.025820


View Full Text Article

Acrobat PDF (5854 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

Aerosols affect climate, health and aviation. Currently, their retrieval assumes a plane-parallel atmosphere and solely vertical radiative transfer. We propose a principle to estimate the aerosol distribution as it really is: a three dimensional (3D) volume. The principle is a type of tomography. The process involves wide angle integral imaging of the sky on a very large scale. The imaging can use an array of cameras in visible light. We formulate an image formation model based on 3D radiative transfer. Model inversion is done using optimization methods, exploiting a closed-form gradient which we derive for the model-fit cost function. The tomography model is distinct, as the radiation source is unidirectional and uncontrolled, while off-axis scattering dominates the images.

© 2013 OSA

1. Introduction

Optical lightfield and integral imaging [1

1. A. Stern and B. Javidi, “Three-dimensional image sensing, visualization, and processing using integral imaging,” Proceedings of the IEEE 94.3 (2006).

3

3. M. Levoy, R. Ng, A. Adams, M. Footer, and M. Horowitz, “Light field microscopy,” ACM Transactions on Graphics (TOG) 25, 924–934 (2006). [CrossRef]

] sample the optical radiance distribution in location and direction. These are mainly used in small-scale setups. This paper deals with a far larger scale, to estimate scatterer 3D spatial distribution in the sky. The atmosphere allows light to pass through in multiple locations and directions. Light is affected by this medium. Therefore, in this paper we lay out a principle to recover this 3D distribution using measured and modeled light-fields (Fig. 1). 3D recovery of this medium has direct implications to various scientific communities that either rely on remotely-sensed imagery, study the atmosphere, or overcome the medium to see beyond. These include meteorology, atmospheric sciences, volcanology, and climatology. Aerosol retrieval is important for understanding climate evolution [4

4. U. Dayan, B. Ziv, T. Shoob, and Y. Enzel, “Suspended dust over southeastern Mediterranean and its relation to atmospheric circulations,” International Journal of Climatology 924, 915–924 (2008). [CrossRef]

, 5

5. O. V. Kalashnikova, M. J. Garay, A. B. Davis, D. J. Diner, and J. V. Martonchik, “Sensitivity of multi-angle photo-polarimetry to vertical layering and mixing of absorbing aerosols: Quantifying measurement uncertainties,” Journal of Quantitative Spectroscopy and Radiative Transfer 112, 2149–2163 (2011). [CrossRef]

] and monitoring air quality. Mapping aerosol density is significant to aviation safety, which needs real time assessment of conditions and visibility around flight paths.

Fig. 1 Integral (lightfield) imaging through a volumetric distribution in the atmosphere, using ground-based cameras.

We rely on sky imaging from multiple directions and locations. Such projection is similar to tomography in other scientific domains [9

9. A. Bluestone, G. Abdoulaev, C. Schmitz, R. Barbour, and A. Hielscher, “Three-dimensional optical tomography of hemodynamics in the human head,” Opt. Express 9, 272–86 (2001). [CrossRef] [PubMed]

, 10

10. I. Ihrke, K. N. Kutulakos, H. P. Lensch, M. Magnor, and W. Heidrich, “State of the art in transparent and specular object reconstruction,” in EUROGRAPHICS 2008 STAR–STATE OF THE ART REPORT , (2008).

]. However, the situation here is distinct. Most tomography setups have a controlled and/or multidirectional radiation source [11

11. D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” Sig.Proc. Magazine 18, 57–75 (2001). [CrossRef]

, 12

12. H. Messer, A. Zinevich, and P. Alpert, “Environmental sensor networks using existing wireless communication systems for rainfall and wind velocity measurements,” IEEE Instrumentation & Measurement Magazine 15, 32–38 (2012). [CrossRef]

]. In contrast, our source is the uncontrolled, unidirectional Sun. Moreover, typical tomography relies on a linear model [13

13. J. Gregson, M. Krimerman, M. B. Hullin, and W. Heidrich, “Stochastic tomography and its applications in 3D imaging of mixing fluids,” ACM Trans. Graph. 31, 52:1–52 (2012). [CrossRef]

]: the pixel value is a linear combination of components along a line of sight (LOS), or a multiplicative combination (linearized by a logarithm). Linear tomography can detect gases [14

14. B. R. Cosofret, D. Konno, A. Faghfouri, H. S. Kindle, C. M. Gittins, M. L. Finson, T. E. Janov, M. J. Levreault, R. K. Miyashiro, and W. J. Marinelli, “Imaging sensor constellation for tomographic chemical cloud mapping,” Appl. Opt. 48, 1837–52 (2009). [CrossRef] [PubMed]

], which absorb UV or IR radiation. However, aerosol attenuation of visible light is typically dominated by scattering, rather than absorption. Since the radiation source is single and fixed in space and time, we cannot rely on direct illumination for tomographic recovery of attenuation fields [15

15. J. A. Aviles, “The Development and Validation of a First Generation X-Ray Scatter Computed Tomography Algorithm for the Reconstruction of Electron Density Breast Images Using Monte Carlo Simulation,” Ph.D. thesis (2011).

]. We may only use sunlight scattered into the LOS. The model is nonlinear, yet tractable. We model passive optical tomographic imaging of 3D atmospheric scatterer distributions in cloudless conditions. Then, we solve this tomography problem, to recover the distribution. Recovery is formulated as an optimization that minimizes a cost function. We derive the gradient of this cost function, to enable efficient optimization.

2. Theoretical background

Extinction: Sun rays (SR) irradiate a small volume that includes particles of a certain type. Each particle has an extinction cross section for interacting with the irradiance. An aerosol extinction cross section is denoted σaerosol. The aerosol density is n. Per unit volume, the extinction coefficient due to aerosols is βaerosol = σaerosoln. In addition, the atmosphere contains air molecules. The extinction coefficient due to the molecules is βair. The volume has infinitesimal length dl. Then, the relative portion of extinct SR irradiance is the unitless differential optical depth, = (βaerosol + βair)dl. The optical depth aggregates in extended propagation:
τ=dτ=(βaerosol+βair)dl=(σaerosoln+βair)dl=τair+σaerosolndl,
(1)
where τair = ∫ βairdl. Through an attenuating atmosphere, the transmittance exponentially decays with the optical depth:
t=exp(τ).
(2)

Scattering: Interaction of a single particle with the irradiance is by absorption and scattering. The weight of scattering (to all directions), relative to the total extinction is given by the unitless single scattering albedo of the particle. For an aerosol, the single scattering albedo is denoted ϖaerosol. The scattering coefficient due to aerosols in the volume is
αaerosol=ϖaerosolβaerosol=ϖaerosolσaerosoln.
(3)
The angular distribution of scattering by an aerosol is determined by a phase function Paerosol, which is normalized: its integral over all solid angles is unit. Part of the light scatters towards a camera’s LOS, as illustrated in Fig. 1. The angle between the SR and LOS is the scattering angle Φscatter. The angular scattering coefficient due to aerosols is
α˜aerosol(Φscatter)=αaerosolPaerosol(Φscatter)=ϖaerosolσaerosolnPaerosol(Φscatter).
(4)
The phase function is often approximated by a parametric expression. Specifically, the Henyey-Greenstein [16

16. W. M. Cornette and J. G. Shanks, “Physically reasonable analytic expression for the single-scattering phase function,” Appl. Opt. 31, 3152–3160 (1992). [CrossRef] [PubMed]

] function is governed by an anisotropy parameter g:
Pgaerosol(Φscatter)38π(1g2)[1+(cosΦscatter)2](2+g2)(1+g22gcosΦscatter)32
(5)

Scattering by air molecules follows the Rayleigh model. The phase function is
Pair(Φscatter)316π(1+cos2Φscatter)
(6)
and the single scattering albedo in visible light is ϖair=1. Air molecular density falls off approximately exponentially with altitude h, with a characteristic [17

17. L. Levi, Applied Optics (John Wiley & Sons, Inc., 1980).

] falloff height Hair = 8 km. Thus, the coefficients for extinction and scattering by air molecules can be modeled by [17

17. L. Levi, Applied Optics (John Wiley & Sons, Inc., 1980).

],
αair(h,λ)=βair(h,λ)1.09×103λ4exp(h/Hair),
(7)
where λ is the wavelength, in microns.

Inverse transform sampling: A tool for RT modeling is the Monte-Carlo (MC) method. MC considers scattering and extinction as random phenomena, that are sampled from their probability distributions. Inverse transform sampling [18

18. L. Devroye, “Sample-based non-uniform random variate generation,” in Proceedings of the 18th conference on Winter simulation , (ACM, 1986), pp. 260–265. [CrossRef]

] is used for sampling random numbers that comply with a specified probability density function. Let u be a random number drawn from a uniform distribution in the interval [0, 1]. The number u can be transformed into a random variable X, whose cumulative distribution function (CDF) is F(X). The transform is X = F−1(u), where F−1 denotes the inverse function of F.

For example, consider a photon propagating in the atmosphere. The photon has high probability of propagating as long as t is high, and the probability diminishes as t → 0. Thus (Eq. (2)) can be viewed as a probability density function, whose CDF is
F(τ)=0τexp(τ)dτ=1exp(τ).
(8)
Thus a photon propagates to a random optical depth
τrandom=F1(u)=ln(1u).
(9)
If the atmosphere is uniform, then from Eq. (1) the photon reaches a random distance
lrandom=(βaerosol+βair)1τrandom=(βaerosol+βair)1ln(1u),
(10)
before interacting with a particle. If the number of photons launched is very high, their average number falls off in consistency with Eqs. (1) and (2). Analogous transforms generate scattering at random angles, according to the phase function.

3. 3D recovery approach and its assessment

A set of ground-based wide-angle cameras looks upwards (Fig. 1), performing integral imaging [19

19. S.-H. Hong, J.-S. Jang, and B. Javidi, “Three-dimensional volumetric object reconstruction using computational integral imaging,” Opt. Express 12, 483–491 (2004). [CrossRef] [PubMed]

] of the sky. Per viewpoint, the view azimuth and elevation angles (relative to the zenith and north) are encapsulated in vector Θ. Any camera c is at a distinct fixed location, capturing raw image ic(Θ). This paper reports a first step: formulating a principle for estimating 3D atmospheric aerosols distributions, using such data. To theoretically study the feasibility and cross validate it, we perform RT using two different forward models. A forward model takes as input a 3D aerosol distribution, and outputs images as if captured at various viewpoints. One model uses MC (Sec. 6): it is stochastic, slow but naturally expresses multiple-scattering. Hence, MC rigorously simulates realistic scenes of arbitrary complexity. The second model uses the single-scattering approximation (Sec. 4). It is less accurate than MC, but solves RT in an analytic, closed form. This form enables simple inversion of the model, to estimate the underlying aerosol distribution (Sec. 7).

4. Single-scattering forward model

As illustrated in Figs. 1 and 2, the sky is discretized into a grid of Nvoxels rectangular cuboid volume elements (voxels), indexed by k, m or q. In a single-scattering model, any SR changes direction only once on the way to a camera. This model is valid in atmospheres that are not very dense (inside fog and clouds this model does not apply). Here, RT has three steps (Fig. 1): ( i) Attenuation of a SR propagating from the top of the atmosphere (TOA) to voxel k. ( ii) Light scattering at voxel k, towards a camera. ( iii) Light attenuation on the LOS from voxel k to the camera. Steps i, iii involve the optical depth along the light path to and from voxel k (analyzed in Sec. 4.1). Step ii involves the scattering coefficient at voxel k (Sec. 4.2).

Fig. 2 Simulated 3D aerosol distributions: [a] Haze blobs and [b] Haze Front. Here, an aerosol density unit is 106 particles/m3. The atmopsheric domain has area of 50 × 50km2, extending from the ground up to altitude of 10km. The domain is divided into a rectilinear grid having Nvoxels voxels.

4.1. Optical depth

Denote a SR line segment between the TOA and the center of voxel k by [SR, k] (See Fig. 1). This SR intersects voxel q. The intersection line segment has geometric length lSR(qz, ΦSR, k), where Δz is a voxel’s vertical geometric thickness, and ΦSR is the Sun angle. Define a Nvoxels × Nvoxels matrix DSun→voxel, whose element (k, q) represents lSR(q):
DSunvoxel(k,q)={lSR(q)ifq[SR,k]0otherwise
(11)
Matrix DSun→voxel is sparse. Analogously, for camera c, denote by [LOSc, k] a LOS bounded between the camera and the center of voxel k (see Fig. 1). The LOS zenith angle is ΦcLOS. Suppose this LOS intersects voxel q. The geometric length of this intersection line segment is lLOSc(mz, ΦcLOS, k). As in Eq. (11), define a Nvoxels × Nvoxels sparse matrix whose element (k, m) is
Dcvoxelcam(k,m)={lLOSc(m)ifm[LOSc,k]0otherwise.
(12)

We focus on situations where, in addition to air molecules, there is a single type of aerosol over a scene. Hence, the three-element vector [σaerosol, ϖaerosol, g] is uniform across the scene, but the density distribution n(k) is variable. As a numerical approximation, assume that within any voxel, the molecular parameters {βair(k), αair(k)} and the aerosol density n(k) are constants, e.g., corresponding to the value at each voxel center. Following Eq. (1), the optical depths between the TOA and voxel k, and from voxel k to camera c are
τSR(k)=q[SR,k]lSR(q)[βair(q)+σaerosoln(q)]τLOSc(k)=m[SR,k]lLOSc(m)[βair(m)+σaerosoln(m)].
(13)
Let n, τSR and βair be column stack vector representations of n(k), τSR(k) and βair(k), respectively. Then, we can write Eq. (13) using matrix notation:
τSR=DSunvoxelβair+σaerosolDSunvoxelnτLOSc=Dcvoxelcamβair+σaerosolDcvoxelcamn.
(14)

The SR and LOS paths joining at each voxel form a single path from the TOA to camera c. Hence, in matrix notation, Eq. (1) yields the total optical depths corresponding to LOSs (of camera c) and SRs that cross all voxels:
τc=τcair+σaerosolDcn,
(15)
where Dc=DSunvoxel+Dcvoxelcam and τcair=Dcβair.

4.2. Scattering

Per camera c and voxel k, the lines [SR, k] and [LOSc, k] intersect at a fixed angle Φcscatter(k). Column-stacking all angles yields the vector representation of all scattering angles in the domain, Φcscatter per camera c. Using Eqs. (4) and (6), the angular scattering coefficients across the domain are expressed in vector form by
α˜cair=βairPair(Φcscatter),α˜caerosol=ϖaerosolσaerosolnPgaerosol(Φcscatter).
(16)
Here the operator ⊙ denotes the Hadamard (element-wise) product.

4.3. Image capture

Compounding the attenuation of irradiance along both [SR, k] and [LOSc, k], and scattering by voxel k towards the camera (Eqs. (2), (4) and (15)) the voxel contributes radiant power
pc(k)=LTOA[α˜cair(k)+α˜caerosol(k)]exp[τc(k)],
(17)
where LTOA is the radiance at the TOA. A column-stack vector of all voxel contributions is
pc=LTOA(α˜cair+α˜caerosol)exp(τc).
(18)

A camera sensor comprises Npix pixels. Each pixel collects light from a narrow cone in the atmosphere. The cone contains or intersects some voxels, while being oblivious to all the rest. Overall, light power captured at the pixel is a weighted sum of the power pc(k) radiating from all voxels. This sum is formulated by a matrix operation Πc over pc:
ic=Πcpc,
(19)
where ic is the image, column-stacked to a vector Npix long. Combining Eqs. (15), (16), (18) and (19), the captured image is thus
ic=LTOAΠc{[α˜cair+ϖaerosolσaerosolPgaerosol(Φcscatter)n]exp[(τcair+σaerosolDcn)]}.
(20)

5. Examples

Figure 3 shows examples of photographs rendered using the single-scattering model. To create such photographs we simulated several scenarios, whose details are as follows.

Fig. 3 Different viewpoints of a sky that includes Haze Blobs. Images are rendered using the single-scattering model. A yellow dot marks the Sun location. Dashed white circles mark zenith angles.

Geometry: The sun is at zenith angle ΦSR = 45°. Its LTOA is obtained from [20

20. M. Charity, “Blackbody color datafile,” (2001), http://www.vendian.org/mncharity/dir3/blackbody/UnstableURLs/bbr\_color.html.

, 21

21. Wikipedia, “Sunlight — Wikipedia, The Free Encyclopedia,” (2012), http://en.wikipedia.org/w/index.php?title=Sunlight\&oldid=502554571.

], with respective red-green-blue ratios of LRTOA:LGTOA:LBTOA=255:236:224.

The atmospheric domain is described in Fig. 2. We tested domain division to various grids, from 10 × 10 × 20 voxels, up to 50 × 50 × 100. The corresponding voxels had width, length and thickness of 5km × 5km × 0.5km and 1km × 1km × 0.1km. In all cases, the vertical resolution was higher than horizontal, via use of voxels that are vertically thin but laterally wide. This expresses the natural higher variability of air and aerosols in the vertical dimension. Each ground-based camera has a hemispherical field of view, and low resolution, 128px × 128px, as cloudless (yet hazy) sky images are rather smooth.

Aerosol: We used particle-type 6 from the aerosol list in [22

22. J. V. Martonchik, R. A. Kahn, and D. J. Diner, “Retrieval of aerosol properties over land using MISR observations,” in Satellite Aerosol Remote Sensing over Land,, A. A. Kokhanovsky and G. Leeuw, eds. (Springer BerlinHeidelberg, 2009), pp. 267–293. [CrossRef]

]. This is a spherical non-absorbing sea-salt/organic particle, whose anisotropy parameter per color channel is [gR, gG, gB] = [0.763, 0.775, 0.786]. Its corresponding extinction cross sections are [σRaerosol,σGaerosol,σBaerosol]=[16.5,16.2,15.9]μm2. At all channels, ϖaerosol = 1. We simulated spatial distributions using a product of two functions. To express the general trend of reduced density with altitude h(k) of voxel k, we follow [17

17. L. Levi, Applied Optics (John Wiley & Sons, Inc., 1980).

] and set the first function as
f1[h(k)]=nsealevelexp[h(k)/Haerosol],
(21)
where Haerosol = 2 km. To express a clustered nature of aerosol distributions (as clouds), we define blobs in the form of 3D ellipsoids. There may be multiple ellipsoids suspended. Then
f2(k)={1ifkanyblobcluster0otherwise.
(22)
The true aerosol distribution is then ntrue = 𝒮 {f1f2}, where 𝒮 is 3D spatial smoothing, obtained by narrow 3D Gaussian filtering. It expresses non-sharp boundaries of typical distributions. The Haze Blobs scene in Fig. 2(a) uses two ellipsoids: one is 32km wide, 2.8km thick and centered at altitude 2.5km; the other is 24km wide, 2.1km thick and centered at altitude 5km. Horizontal widths are much larger than the vertical thickness, in consistency with atmospheric scales. Figure 3 shows photographs of the Haze Blobs scene, as simulated by Eq. (20). For clarity of display, the brightness of the displayed pictures in this paper underwent the same standard contrast enhancement (stretching), including gamma-correction. The recovery algorithm, of course, uses the raw (not brightness enhanced) images.

The Haze Front scene in Fig. 2(b) uses an ellipsoid degenerated to an elliptic cylinder that partly enters the analyzed volume. It is 32km long (only part of which enters the volume), 4km thick and centered at altitude 5km. The front stretches across the domain width.

6. Multiple-scattering forward model by Monte Carlo

Fig. 4 [a] Multiple-scattering forward model by MC. [b] Photograph of the same sky and view point as in Fig. 3(a), rendered using MC.

The quality of MC increases with the number of photons launched. Figure 4(b) shows a photograph of the Haze Blobs scene, derived by the MC process, using 10000 photons per pixel, from the same viewpoint as Fig. 3(a). Similarly, Fig. 5(a) corresponds to the same viewpoint as Fig. 3(d). Figure 6 shows the contribution by successive scattering orders. They are derived from the MC image. Photons contributing to any pixel are accumulated in steps ii and v above. The contributing photons that stem from just a single scattering event yield Fig. 6(a). Similarly, contributing photons that had experienced exactly two or three scattering events yield respectively Fig. 6(b) and Fig. 6(c). First-order scatter yields most of the energy in the MC image. Radiance by higher order scattering is mostly evident at the horizon. A horizontal LOS passes through a long and dense part of the atmosphere, while a vertical LOS cuts short through the dense, lower part. Hence, toward the horizon the probability of high order scattering increases.

Fig. 5 Images simulated by MC, from the same viewpoint, but different aerosol characteristics. (a) Haze blobs, characterized in Sec. 5. (b) Partly absorbing aerosol. (c) Aerosol having an isotropic phase function. (d) High aerosol density. Details of the scenarios used in creating (b) and (c) are given in Sec. 8.
Fig. 6 Separation of the image shown in Fig. 4(b) to contributions by successive orders of scattering: (a) first, (b) second, (c) third and (d) forth scattering order.

Figure 7 plots the image intensity profile along a straight line from left to right, via the image center (zenith). The plot compares photographs simulated using the single-scattering and the MC forward models. Figure 7(a) is extracted from Fig. 3(d) and Fig. 5(a) of the Haze Blobs scene. This plot demonstrates consistency in this scene, for most pixels. Figure 7(b) plots the same profiles, where the aerosol density was increased tenfold (see Sec. 8 and the corresponding MC photograph, Fig. 5(d)). In Fig. 7(b) the plots differ more than in Fig. 7(a), due to high order scattering.

Fig. 7 Cross-sections along the X-axis of photographs rendered by the single-scattering (dotted) and MC (solid line) forward models. (a) cross sections of the Haze Blobs scene (described in Sec. 5). The difference between the models is much more pronounced (b) when the scatterers are ten times denser (see Sec. 8).

7. Inverse problem

When the type of aerosol above an area is known [22

22. J. V. Martonchik, R. A. Kahn, and D. J. Diner, “Retrieval of aerosol properties over land using MISR observations,” in Satellite Aerosol Remote Sensing over Land,, A. A. Kokhanovsky and G. Leeuw, eds. (Springer BerlinHeidelberg, 2009), pp. 267–293. [CrossRef]

], estimation of the density distribution n is needed. The data are Nviews measured photographs {icmeasured}c=1Nviews. Recovery is phrased as optimization of a cost function that fits the image-formation model to the data. Let 𝒞 be the set of all distributions complying with some constraints. Particularly, n is non-negative, and its spatial support is bounded between the ground and the TOA. The optimization problem is
n^=argminn𝒞E(n)
(24)
where
E(n)=c=1Nviews[icmeasuredic(n)]22+ηΨ(n).
(25)
Here masks the area around the sun: in real-world images it is indeed blocked. We mask the horizon, by using a camera whose field of view is a little narrower than hemispherical. This reduces the influence of LOSs that are more strongly affected by high-order scattering (Fig. 6(b)). Consequently, the fit of the imaging model focuses on the bulk image region, where the single scattering model is more valid. In Eq. (25), Ψ(n) is a regularization term that expresses some prior knowledge on the distribution, while η is the regularization weight.

Since aerosol distributions are usually fuzzy, useful regularization is by a smoothness term, which penalizes for energy in the second order spatial derivatives (3D Laplacian), 2n22. Regularization is not required when data is sufficient and reliable. Voxels at high altitude have good coverage by many cameras at multiple directions (Fig. 1), while low altitude voxels are mainly observed by local cameras. Hence, regularization can be weakened with altitude. We accomplish this weakening using a weight w(k) = exp[−h(k)/Hsmooth]. Overall we use
Ψ(n)=Wn22
(26)
where is a matrix representation of the 3D Laplacian operator, and matrix W is diagonal, whose elements are w(k).

We solve Eq. (24) using standard optimization tools. The gradient of E with respect to n is
nE=2c=1Nviews[Jic(n)][icmeasuredic(n)]+2ηWWn.
(27)
Here the matrix Jic(n) is the Jacobian of the vector ic with respect to n. Element (Θ, k) of the Jacobian differentiates the intensity in pixel Θ (in viewpoint c) with variation of the aerosol density at voxel k, i.e., ∂ic(Θ)/∂n(k).

We now detail the derivation of the Jacobian Jic(n). First, we provide some results relating to differentiation. Let a(n), u(n) be vector functions: each outputs a vector of length r. Let C be a r × Nvoxels matrix, where Nvoxels is the length of n. Then,
C(au)n=[an𝔻{u}+un𝔻{a}]C.
(28)
Here 𝔻{v} is a conversion of a general vector v into a diagonal matrix, whose main diagonal elements correspond to the elements of v. Now, let a(n) = exp(−Cn), where the exponent is element-wise (not raising an operator to some power). Then
an=C𝔻{exp(Cn)}.
(29)
Using Eq. (29),
exp[(τcair+σaerosolDcn)]n=σaerosolDc𝔻{exp[(τcair+σaerosolDcn)]},
(30)
where τcair is independent of n. Using Eq. (28),
[ϖaerosolσaerosolPgaerosol(Φcscatter)n]n=ϖaerosolσaerosol𝔻{Pgaerosol(Φcscatter)}.
(31)
Based on Eqs. (28), (30) and (31), we derive the Jacobian of Eq. (20) in close-form,
Jic(n)=LTOAσaerosol(AB)𝔻{exp[(τcair+σaerosolDcn)]}Πc,
(32)
where
A=ϖaerosol𝔻{Pgaerosol(Φcscatter)},B=Dc𝔻{[α˜cair+ϖaerosolσaerosolPgaerosol(Φcscatter)n]}
(33)

8. Recovery simulations

To demonstrate recovery, we performed simulations, and tested effects of the following variations: {1} Multiple or only single scatter in the forward model; {2} Spatial distribution; {3} Density of viewpoints; {4} Aerosol characteristics; {5} Density scale.

In all cases, images were created as described in Secs. 5 and 6. These images served as input to the reconstruction algorithm. The optimization used an L-BFGS-B solver [25

25. C. Zhu, R. H. Byrd, P. Lu, and J. Nocedal, “Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization,” ACM Trans. Math. Softw. 23, 550–560 (1997). [CrossRef]

] on a computer cluster. Each computer core was dedicated to rendering a modeled image. The optimization was initialized by nnn = 0. Satisfactory convergence occurred in several hundred iterations. Depending on the resolution it took between minutes to a couple of hours to complete, in total. The total estimation error can be quantified by the aerosol mass that is over and under-estimated in all voxels, relative to the total aerosol mass in the scene. Using the 1 norm, the total mass relative difference is δmass = (||||1 − ||ntrue||1)/||ntrue||1. To also sense local errors, we use ε = ||ntrue||1/||ntrue||1. These results are listed in Table 1.

Table 1. Relative errors in various simulations. Here 𝒪 denotes order of magnitude.

table-icon
View This Table

Forward models and distributions: Figure 8 shows estimation results, corresponding to the original distributions of Fig. 2. Here, the atmosphere domain was defined over a 20 × 20 × 40 grid. Cameras were placed ∼ 7km apart on a 6 × 6 grid. In Figs. 8(a) and 8(b), the single-scattering model was used both in rendering of the raw images, and in the estimation algorithm. In Figs. 8(c) and 8(d), the MC model (multiple scattering) rendered the raw images, over which our recovery (Sec. 7) was tested. Overall, the distributions, the density order of magnitude and values are consistent. Errors in Figs. 8(a) and 8(b) stem from random noise, which had been induced in the raw images. Errors are larger when MC is used to render the images, since high-order scattering is not accounted for in the algorithm of Sec. 7.

Fig. 8 Recoveries of scenes shown in Fig. 2. Color represents aerosol density, in units of 106 particles/m3. In (a) and (b) images simulated by single-scattering are used as input for the recovery. In (c) and (d), MC simulated images are used as input for the recovery.

Density of viewpoints: Keeping the domain and scene (Haze blobs) fixed, we repeated the reconstruction many times, each using a different subset of the above mentioned 36 simulated cameras. Each test randomly selected viewpoints, for various fixed values of Nviews < 36. All rendered images used MC. Figure 9 plots the performance. Diminishing returns are obtained beyond Nviews ≈ 20, which corresponds to ≈ 11km separation between neighboring viewpoints. These tests point to fundamental questions about spatial-angular sampling: how spatially dense should cameras be, and how dense should a camera’s field of view be sampled (in pixels)? This likely depends on the finesse of atmospheric structure sought, and the level of optical diffusion achieved by various aerosols. The novel tomography domain introduced here thus raises new theoretical questions that will require extensive further research.

Fig. 9 The relative error ε measure decreases with the number of cameras. Bars represent the standard deviation of ε, over our random tests.

Aerosol characteristics: We used the Haze blobs scene, but varied the aerosol characteristics. Compared to the aerosol described in Sec. 5, we tested two variations. In one variation, the phase function was isotropic, i.e., [gR, gG, gB] = [0, 0, 0]. In the other variation, the aerosol is partly absorbing: its single scattering albedo per color channel is [ϖRaerosol,ϖGaerosol,ϖBaerosol]=[0.69,0.74,0.78]. All rendered images used MC (See Figs. 5(a)–5(c)). Recovery results using Nviews = 36 are shown in Figs. 10(a) and 10(b) and in Table 1. The results indicate that if the aerosol’s phase function has smaller anisotropy, recovery of the aerosol’s density is easier. We hypothesize that tomography is better served by scattering isotropy, since then significant radiance from any voxel reaches cameras in a wide range of directions. This may make back-projection (recovery from cameras to voxels) more accurate. On the other hand, if the phase function is too strongly peaked around one direction, most cameras would not receive much scattered radiance from a voxel, undermining tomography. A very broad research is needed to verify and quantify such dependencies.

Fig. 10 Recovered distributions when the aerosol is either party absorbing (a), has isotropic phase function (b) or has high density (c).

Density scale: We used the aerosol Haze blobs distribution and characteristics as described in Sec. 5, but increased the aerosol density tenfold everywhere. This significantly increases the effects of high order scattering, as shown in Fig. 7(b) and Fig. 5(d). Consequently, ε increases (Table 1), when using Nviews = 36. Still, Fig. 10(c) shows that the spatial distribution outlines are recovered similarly to the other tests.

9. Discussion

We describe a way to sense the 3D atmosphere. It can help remote sensing depart from common one-dimensional aerosol distribution models, in favor of detailed 3D RT analysis. The paper describes both a novel data acquisition system concept for this recovery task (ground-based camera network), and a dedicated algorithm for reconstructing aerosol distributions. Using more prior knowledge on the nature of distributions will improve the recovery.

The estimation algorithm uses a single scattering model. This approximation is valid when the optical depth is small. However, as MC simulations show, it is important to generalize the recovery algorithm so that high-order scattering is included in it. Such generalization would be computationally complex: in our tests, the MC forward model was slower by 2–3 orders of magnitude than a single-scattering analytic model. Thus, inverting a 3D multiple-scattering model requires research into numerical schemes for increasing the solution efficiency of both forward and inverse problems. The remote sensing community has devised several approaches to better approximate multiple-scattering problems, mainly in one-dimensional models (e.g., successive orders of scattering and the discrete ordinates method). We may tap into some of these approaches. Deploying an experimental system will be needed: camera specifications [26

26. N. J. Pust, A. R. Dahlberg, M. J. Thomas, and J. a. Shaw, “Comparison of full-sky polarization and radiance observations to radiative transfer simulations which employ AERONET products,” Opt. Express 19, 18602–18613 (2011). [CrossRef] [PubMed]

] should be set to optimally recover a wide range of aerosol distributions.

Acknowledgments

We are grateful to David J. Diner for fruitful discussions and his support at all levels. We thank the reviewers for useful comments. Yoav Schechner is a Landau Fellow - supported by the Taub Foundation. This work is supported by the US-Israel Binational Science Foundation (BSF grant 2012202) and conducted in the Ollendorff Minerva Center. Minerva is funded through the BMBF. Part of this work was performed at the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA and with support by NASA’s Earth Science Technology Office Advanced Information Systems Technology Program.

References and links

1.

A. Stern and B. Javidi, “Three-dimensional image sensing, visualization, and processing using integral imaging,” Proceedings of the IEEE 94.3 (2006).

2.

J. Kim, D. Lanman, Y. Mukaigawa, and R. Raskar, “Descattering transmission via angular filtering,” in Proc. ECCV’ 10, (Springer-Verlag, Berlin, Heidelberg, 2010), pp. 86–99.

3.

M. Levoy, R. Ng, A. Adams, M. Footer, and M. Horowitz, “Light field microscopy,” ACM Transactions on Graphics (TOG) 25, 924–934 (2006). [CrossRef]

4.

U. Dayan, B. Ziv, T. Shoob, and Y. Enzel, “Suspended dust over southeastern Mediterranean and its relation to atmospheric circulations,” International Journal of Climatology 924, 915–924 (2008). [CrossRef]

5.

O. V. Kalashnikova, M. J. Garay, A. B. Davis, D. J. Diner, and J. V. Martonchik, “Sensitivity of multi-angle photo-polarimetry to vertical layering and mixing of absorbing aerosols: Quantifying measurement uncertainties,” Journal of Quantitative Spectroscopy and Radiative Transfer 112, 2149–2163 (2011). [CrossRef]

6.

M. I. Mishchenko and I. V. Geogdzhayev, “Satellite remote sensing reveals regional tropospheric aerosol trends,” Opt. Express 15, 7423–38 (2007). [CrossRef] [PubMed]

7.

J. Martonchik, D. D. J. Diner, J. M. David, R. Kahn, T. P. Ackerman, M. M. Verstraete, B. Pinty, and H. R. Gordon, “Techniques for the Retrieval of Aerosol Properties over Land and Ocean Using Multi-angle Imaging,” IEEE Trans. on Geoscience and Remote Sensing 36, 1212–1227 (1998). [CrossRef]

8.

E. Namer, S. Shwartz, and Y. Y. Schechner, “Skyless polarimetric calibration and visibility enhancement,” Opt. Express 17, 472–93 (2009). [CrossRef] [PubMed]

9.

A. Bluestone, G. Abdoulaev, C. Schmitz, R. Barbour, and A. Hielscher, “Three-dimensional optical tomography of hemodynamics in the human head,” Opt. Express 9, 272–86 (2001). [CrossRef] [PubMed]

10.

I. Ihrke, K. N. Kutulakos, H. P. Lensch, M. Magnor, and W. Heidrich, “State of the art in transparent and specular object reconstruction,” in EUROGRAPHICS 2008 STAR–STATE OF THE ART REPORT , (2008).

11.

D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” Sig.Proc. Magazine 18, 57–75 (2001). [CrossRef]

12.

H. Messer, A. Zinevich, and P. Alpert, “Environmental sensor networks using existing wireless communication systems for rainfall and wind velocity measurements,” IEEE Instrumentation & Measurement Magazine 15, 32–38 (2012). [CrossRef]

13.

J. Gregson, M. Krimerman, M. B. Hullin, and W. Heidrich, “Stochastic tomography and its applications in 3D imaging of mixing fluids,” ACM Trans. Graph. 31, 52:1–52 (2012). [CrossRef]

14.

B. R. Cosofret, D. Konno, A. Faghfouri, H. S. Kindle, C. M. Gittins, M. L. Finson, T. E. Janov, M. J. Levreault, R. K. Miyashiro, and W. J. Marinelli, “Imaging sensor constellation for tomographic chemical cloud mapping,” Appl. Opt. 48, 1837–52 (2009). [CrossRef] [PubMed]

15.

J. A. Aviles, “The Development and Validation of a First Generation X-Ray Scatter Computed Tomography Algorithm for the Reconstruction of Electron Density Breast Images Using Monte Carlo Simulation,” Ph.D. thesis (2011).

16.

W. M. Cornette and J. G. Shanks, “Physically reasonable analytic expression for the single-scattering phase function,” Appl. Opt. 31, 3152–3160 (1992). [CrossRef] [PubMed]

17.

L. Levi, Applied Optics (John Wiley & Sons, Inc., 1980).

18.

L. Devroye, “Sample-based non-uniform random variate generation,” in Proceedings of the 18th conference on Winter simulation , (ACM, 1986), pp. 260–265. [CrossRef]

19.

S.-H. Hong, J.-S. Jang, and B. Javidi, “Three-dimensional volumetric object reconstruction using computational integral imaging,” Opt. Express 12, 483–491 (2004). [CrossRef] [PubMed]

20.

M. Charity, “Blackbody color datafile,” (2001), http://www.vendian.org/mncharity/dir3/blackbody/UnstableURLs/bbr\_color.html.

21.

Wikipedia, “Sunlight — Wikipedia, The Free Encyclopedia,” (2012), http://en.wikipedia.org/w/index.php?title=Sunlight\&oldid=502554571.

22.

J. V. Martonchik, R. A. Kahn, and D. J. Diner, “Retrieval of aerosol properties over land using MISR observations,” in Satellite Aerosol Remote Sensing over Land,, A. A. Kokhanovsky and G. Leeuw, eds. (Springer BerlinHeidelberg, 2009), pp. 267–293. [CrossRef]

23.

H. Iwabuchi, “Efficient Monte Carlo Methods for Radiative Transfer Modeling,” Journal of the Atmospheric Sciences 63, 2324–2339 (2006). [CrossRef]

24.

A. Marshak and A. Davis, 3D Radiative Transfer in Cloudy Atmospheres, Physics of Earth and Space Environments (Springer, 2005). [CrossRef]

25.

C. Zhu, R. H. Byrd, P. Lu, and J. Nocedal, “Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization,” ACM Trans. Math. Softw. 23, 550–560 (1997). [CrossRef]

26.

N. J. Pust, A. R. Dahlberg, M. J. Thomas, and J. a. Shaw, “Comparison of full-sky polarization and radiance observations to radiative transfer simulations which employ AERONET products,” Opt. Express 19, 18602–18613 (2011). [CrossRef] [PubMed]

OCIS Codes
(100.3190) Image processing : Inverse problems
(280.1100) Remote sensing and sensors : Aerosol detection
(280.1310) Remote sensing and sensors : Atmospheric scattering
(280.4991) Remote sensing and sensors : Passive remote sensing
(100.3200) Image processing : Inverse scattering

ToC Category:
Remote Sensing

History
Original Manuscript: August 5, 2013
Revised Manuscript: October 10, 2013
Manuscript Accepted: October 10, 2013
Published: October 22, 2013

Virtual Issues
December 3, 2013 Spotlight on Optics

Citation
Amit Aides, Yoav Y. Schechner, Vadim Holodovsky, Michael J. Garay, and Anthony B. Davis, "Multi sky-view 3D aerosol distribution recovery," Opt. Express 21, 25820-25833 (2013)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-21-22-25820


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. A. Stern and B. Javidi, “Three-dimensional image sensing, visualization, and processing using integral imaging,” Proceedings of the IEEE94.3 (2006).
  2. J. Kim, D. Lanman, Y. Mukaigawa, and R. Raskar, “Descattering transmission via angular filtering,” in Proc. ECCV’ 10, (Springer-Verlag, Berlin, Heidelberg, 2010), pp. 86–99.
  3. M. Levoy, R. Ng, A. Adams, M. Footer, and M. Horowitz, “Light field microscopy,” ACM Transactions on Graphics (TOG)25, 924–934 (2006). [CrossRef]
  4. U. Dayan, B. Ziv, T. Shoob, and Y. Enzel, “Suspended dust over southeastern Mediterranean and its relation to atmospheric circulations,” International Journal of Climatology924, 915–924 (2008). [CrossRef]
  5. O. V. Kalashnikova, M. J. Garay, A. B. Davis, D. J. Diner, and J. V. Martonchik, “Sensitivity of multi-angle photo-polarimetry to vertical layering and mixing of absorbing aerosols: Quantifying measurement uncertainties,” Journal of Quantitative Spectroscopy and Radiative Transfer112, 2149–2163 (2011). [CrossRef]
  6. M. I. Mishchenko and I. V. Geogdzhayev, “Satellite remote sensing reveals regional tropospheric aerosol trends,” Opt. Express15, 7423–38 (2007). [CrossRef] [PubMed]
  7. J. Martonchik, D. D. J. Diner, J. M. David, R. Kahn, T. P. Ackerman, M. M. Verstraete, B. Pinty, and H. R. Gordon, “Techniques for the Retrieval of Aerosol Properties over Land and Ocean Using Multi-angle Imaging,” IEEE Trans. on Geoscience and Remote Sensing36, 1212–1227 (1998). [CrossRef]
  8. E. Namer, S. Shwartz, and Y. Y. Schechner, “Skyless polarimetric calibration and visibility enhancement,” Opt. Express17, 472–93 (2009). [CrossRef] [PubMed]
  9. A. Bluestone, G. Abdoulaev, C. Schmitz, R. Barbour, and A. Hielscher, “Three-dimensional optical tomography of hemodynamics in the human head,” Opt. Express9, 272–86 (2001). [CrossRef] [PubMed]
  10. I. Ihrke, K. N. Kutulakos, H. P. Lensch, M. Magnor, and W. Heidrich, “State of the art in transparent and specular object reconstruction,” in EUROGRAPHICS 2008 STAR–STATE OF THE ART REPORT, (2008).
  11. D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” Sig.Proc. Magazine18, 57–75 (2001). [CrossRef]
  12. H. Messer, A. Zinevich, and P. Alpert, “Environmental sensor networks using existing wireless communication systems for rainfall and wind velocity measurements,” IEEE Instrumentation & Measurement Magazine15, 32–38 (2012). [CrossRef]
  13. J. Gregson, M. Krimerman, M. B. Hullin, and W. Heidrich, “Stochastic tomography and its applications in 3D imaging of mixing fluids,” ACM Trans. Graph.31, 52:1–52 (2012). [CrossRef]
  14. B. R. Cosofret, D. Konno, A. Faghfouri, H. S. Kindle, C. M. Gittins, M. L. Finson, T. E. Janov, M. J. Levreault, R. K. Miyashiro, and W. J. Marinelli, “Imaging sensor constellation for tomographic chemical cloud mapping,” Appl. Opt.48, 1837–52 (2009). [CrossRef] [PubMed]
  15. J. A. Aviles, “The Development and Validation of a First Generation X-Ray Scatter Computed Tomography Algorithm for the Reconstruction of Electron Density Breast Images Using Monte Carlo Simulation,” Ph.D. thesis (2011).
  16. W. M. Cornette and J. G. Shanks, “Physically reasonable analytic expression for the single-scattering phase function,” Appl. Opt.31, 3152–3160 (1992). [CrossRef] [PubMed]
  17. L. Levi, Applied Optics (John Wiley & Sons, Inc., 1980).
  18. L. Devroye, “Sample-based non-uniform random variate generation,” in Proceedings of the 18th conference on Winter simulation, (ACM, 1986), pp. 260–265. [CrossRef]
  19. S.-H. Hong, J.-S. Jang, and B. Javidi, “Three-dimensional volumetric object reconstruction using computational integral imaging,” Opt. Express12, 483–491 (2004). [CrossRef] [PubMed]
  20. M. Charity, “Blackbody color datafile,” (2001), http://www.vendian.org/mncharity/dir3/blackbody/UnstableURLs/bbr\_color.html .
  21. Wikipedia, “Sunlight — Wikipedia, The Free Encyclopedia,” (2012), http://en.wikipedia.org/w/index.php?title=Sunlight\&oldid=502554571 .
  22. J. V. Martonchik, R. A. Kahn, and D. J. Diner, “Retrieval of aerosol properties over land using MISR observations,” in Satellite Aerosol Remote Sensing over Land,, A. A. Kokhanovsky and G. Leeuw, eds. (Springer BerlinHeidelberg, 2009), pp. 267–293. [CrossRef]
  23. H. Iwabuchi, “Efficient Monte Carlo Methods for Radiative Transfer Modeling,” Journal of the Atmospheric Sciences63, 2324–2339 (2006). [CrossRef]
  24. A. Marshak and A. Davis, 3D Radiative Transfer in Cloudy Atmospheres, Physics of Earth and Space Environments (Springer, 2005). [CrossRef]
  25. C. Zhu, R. H. Byrd, P. Lu, and J. Nocedal, “Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization,” ACM Trans. Math. Softw.23, 550–560 (1997). [CrossRef]
  26. N. J. Pust, A. R. Dahlberg, M. J. Thomas, and J. a. Shaw, “Comparison of full-sky polarization and radiance observations to radiative transfer simulations which employ AERONET products,” Opt. Express19, 18602–18613 (2011). [CrossRef] [PubMed]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited