Multi sky-view 3D aerosol distribution recovery

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 Optical Society of America OCIS codes: (100.3190) Inverse problems; (100.3200) Inverse scattering; (280.1100) Aerosol detection; (280.1310) Atmospheric scattering; (280.4991) Passive remote sensing. 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). 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). 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). 6. M. I. Mishchenko and I. V. Geogdzhayev, “Satellite remote sensing reveals regional tropospheric aerosol trends,” Opt. Express 15, 7423–38 (2007). 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). 8. E. Namer, S. Shwartz, and Y. Y. Schechner, “Skyless polarimetric calibration and visibility enhancement,” Opt. Express 17, 472–93 (2009). 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). 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). #195261 $15.00 USD Received 5 Aug 2013; revised 10 Oct 2013; accepted 10 Oct 2013; published 22 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.025820 | OPTICS EXPRESS 25820 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). 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:10 (2012). 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). 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). 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. 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). 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 Berlin Heidelberg, 2009), pp. 267–293. 23. H. Iwabuchi, “Efficient Monte Carlo Methods for Radiative Transfer Modeling,” Journal of the Atmospheric Sciences 63, 2324–2339 (2006). 24. A. Marshak and A. Davis, 3D Radiative Transfer in Cloudy Atmospheres, Physics of Earth and Space Environments (Springer, 2005). 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). 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).


Introduction
Optical lightfield and integral imaging [1][2][3] 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,5] and monitoring air quality. Mapping aerosol density is significant to aviation safety, which needs real time assessment of conditions and visibility around flight paths.
In remote sensing [6], imaging through air is associated with atmospheric correction, based on aerosol retrieval. In this discipline, the atmosphere is often assumed to be plane-parallel, using one-dimensional vertical radiative transfer (RT). Consequently state-of-the-art aerosol retrieval is done in distinct large lateral blocks [7] with limited height resolution [5]. We, however, seek 3D recovery. Atmospheric correction is related to dehazing [8]. In this paper, however, the medium itself, at all relevant altitudes is the object of interest.
We rely on sky imaging from multiple directions and locations. Such projection is similar to tomography in other scientific domains [9,10]. However, the situation here is distinct. Most tomography setups have a controlled and/or multidirectional radiation source [11,12]. In contrast, our source is the uncontrolled, unidirectional Sun. Moreover, typical tomography relies on a linear model [13]: 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], 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]. 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.

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 = σ aerosol n. 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, dτ = (β aerosol + β air )dl. The optical depth aggregates in extended propagation: where τ air = β air dl. Through an attenuating atmosphere, the transmittance exponentially decays with the optical depth: 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 The angular distribution of scattering by an aerosol is determined by a phase function P aerosol , 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 The phase function is often approximated by a parametric expression. Specifically, the Henyey-Greenstein [16] function is governed by an anisotropy parameter g: Scattering by air molecules follows the Rayleigh model. The phase function is 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] falloff height H air = 8 km. Thus, the coefficients for extinction and scattering by air molecules can be modeled by [17], 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] 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 Thus a photon propagates to a random optical depth If the atmosphere is uniform, then from Eq. (1) the photon reaches a random distance 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. Here, an aerosol density unit is 10 6 particles/m 3 . The atmopsheric domain has area of 50 × 50km 2 , extending from the ground up to altitude of 10km. The domain is divided into a rectilinear grid having N voxels voxels.

3D recovery approach and its assessment
A set of ground-based wide-angle cameras looks upwards ( Fig. 1), performing integral imaging [19] 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 i c (Θ Θ Θ). 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).

Single-scattering forward model
As illustrated in Figs. 1 and 2, the sky is discretized into a grid of N voxels 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 ( Step ii involves the scattering coefficient at voxel k (Sec. 4.2).

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 l SR (q|∆z, Φ SR , k), where ∆z is a voxel's vertical geometric thickness, and Φ SR is the Sun angle. Define a N voxels × N voxels matrix D D D Sun→voxel , whose element (k, q) represents l SR (q): Matrix D D D Sun→voxel is sparse. Analogously, for camera c, denote by [LOS c , k] a LOS bounded between the camera and the center of voxel k (see Fig. 1). The LOS zenith angle is Φ LOS c . Suppose this LOS intersects voxel q. The geometric length of this intersection line segment is l LOS c (m|∆z, Φ LOS c , k). As in Eq. (11), define a N voxels × N voxels sparse matrix whose element (k, m) is 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 Let n n 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: 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:

Scattering
Per camera c and voxel k, the lines [SR, k] and [LOS c , k] intersect at a fixed angle Φ scatter c (k). Column-stacking all angles yields the vector representation of all scattering angles in the domain, Φ Φ Φ scatter c per camera c. Using Eqs. (4) and (6), the angular scattering coefficients across the domain are expressed in vector form bỹ Here the operator denotes the Hadamard (element-wise) product.

Image capture
Compounding the attenuation of irradiance along both [SR, k] and [LOS c , k], and scattering by voxel k towards the camera (Eqs. (2), (4) and (15)) the voxel contributes radiant power where L TOA is the radiance at the TOA. A column-stack vector of all voxel contributions is A camera sensor comprises N pix 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 p c (k) radiating from all voxels. This sum is formulated by a matrix operation Π Π Π c over p p p c : where i i i c is the image, column-stacked to a vector N pix long. Combining Eqs. (15), (16), (18) and (19), the captured image is thus 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.
Geometry: The sun is at zenith angle Φ SR = 45 o . Its L TOA is obtained from [20,21], with respective red-green-blue ratios of L TOA R : L TOA G : L TOA B = 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]. This is a spherical non-absorbing sea-salt/organic particle, whose anisotropy parameter per color channel is [g R , g G , 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] and set the first function as where H aerosol = 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 The true aerosol distribution is then n n n true = S { f f f 1 f f f 2 }, where S 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.

Multiple-scattering forward model by Monte Carlo
The MC simulations rely on the same grid as Sec. 4. We use backward MC: from a detector pixel, photons are 'launched' via the camera lens; then the photons are traced back through the atmosphere. Photons that happen to be back-traced to the Sun are counted, per pixel. In this model, RT takes the following steps, per pixel (Fig. 4a): (i) Launch a photon-packet from camera c to direction Θ Θ Θ. This is the initial ray, denoted R 0 . The packet has intensity I 0 . Per iteration i: (ii) Find the distance on ray R i to which the photon-packet propagates. To do this, Eq. (9) yields τ random . Then, Eq. (10) generalizes to a non-uniform atmosphere: using Eq. (1), numerically seek l random s.t. Distance l random along R i yields 3D position X i . If X i is outside the atmospheric grid, the packet is terminated. If, in addition, R i SR , the packet is counted as contributing to the image pixel.
(iii) The type of particle (air molecule or aerosol) that the photon-packet interacts with at point X i is determined randomly. The relative probabilities of this random step are set by the respective extinction coefficients (β air vs. β aerosol ) at the voxel containing X i . (iv) If the particle is an aerosol, the photon-packet intensity is attenuated to I i+1 = ϖ aerosol I i . If I i+1 is lower than a threshold, the packet is stochastically terminated, following [23].
(v) The photon-packet is scattered to a new direction, determined by inverse transform sampling, according to the phase function of the particle (Eqs. 5 or 6). This yields a new ray, denoted R i+1 , emanating from X i . Local estimation [24] derives the amount of light back-traced to the sun at each scattering. The next iteration of propagation (denoted ii above) proceeds. 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 10 000 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.   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.

Inverse problem
When the type of aerosol above an area is known [22], estimation of the density distribution n n n is needed. The data are N views measured photographs {i i i measured c } N views c=1 . Recovery is phrased as optimization of a cost function that fits the image-formation model to the data. Let C be the set of all distributions complying with some constraints. Particularly, n n n is non-negative, and its spatial support is bounded between the ground and the TOA. The optimization problem iŝ n n n = arg min Here M 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 n 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), ∇ 2 n n n 2 2 . 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)/H smooth . Overall we use Ψ(n n n) = WLn n n 2 where L 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 n n is Here the matrix J J J i i i c (n n n) is the Jacobian of the vector i i i c with respect to n n 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., ∂ i c (Θ Θ Θ)/∂ n(k).
We now detail the derivation of the Jacobian J J J i i i c (n n n). First, we provide some results relating to differentiation. Let a a a(n n n), u u u(n n n) be vector functions: each outputs a vector of length r. Let C C C be a r × N voxels matrix, where N voxels is the length of n n n. Then, a a a u u u) ∂ n n n = ∂ a a a ∂ n n n D {u u u} + ∂ u u u ∂ n n n D {a a a} C C C.
Here D {v v v} is a conversion of a general vector v v v into a diagonal matrix, whose main diagonal elements correspond to the elements of v v v. Now, let a a a(n n n) = exp(−C C Cn n n), where the exponent is element-wise (not raising an operator to some power). Then where τ τ τ air c is independent of n n n. Using Eq. (28), ∂ n n n = ϖ aerosol σ aerosol D P aerosol where

Recovery simulations
To demonstrate recovery, we performed simulations, and tested effects of the following varia- 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] on a computer cluster. Each computer core was dedicated to rendering a modeled image. The optimization was initialized by n n n = 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. 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 = ( n n n 1 − n n n true 1 )/ n n n true 1 . To also sense local errors, we use ε = n n n − n n n true 1 / n n n true 1 . These results are listed in Table 1.
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 singlescattering 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.

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 N views < 36. All rendered images used MC. Figure 9 plots the performance. Diminishing returns are obtained beyond N views ≈ 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.
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., [g R , g G , g B ] = [0, 0, 0]. In the other variation, the aerosol is partly absorbing: its single scattering albedo per color channel is  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.
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 N views = 36. Still, Fig. 10(c) shows that the spatial distribution outlines are recovered similarly to the other tests.

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] should be set to optimally recover a wide range of aerosol distributions.