Fitting Local, Low-Dimensional Parameterizations of Optical Turbulence Modeled from Optimal Transport Velocity Vectors

Fitting Local, Low-Dimensional Parameterizations of Optical Turbulence Modeled from Optimal Transport Velocity Vectors

Journal Pre-proof

Fitting Local, Low-Dimensional Parameterizations of Optical Turbulence Modeled from Optimal Transport Velocity Vectors Tegan H. Emerson, Jonathan M. Nichols PII: DOI: Reference:

S0167-8655(19)30302-2 https://doi.org/10.1016/j.patrec.2019.10.023 PATREC 7671

To appear in:

Pattern Recognition Letters

Received date: Revised date: Accepted date:

25 July 2018 8 September 2019 20 October 2019

Please cite this article as: Tegan H. Emerson, Jonathan M. Nichols, Fitting Local, Low-Dimensional Parameterizations of Optical Turbulence Modeled from Optimal Transport Velocity Vectors, Pattern Recognition Letters (2019), doi: https://doi.org/10.1016/j.patrec.2019.10.023

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.

Highlights • Proposes physics-based model for image sequences degraded by atmospheric turbulence. • Shows the model to be low-dimensional over spatially localized regions of the image. • Provides algorithm for identifying local image models and method for combining to form global image model. • Demonstrates the ability to recover clean image from surrogate global image model.

1

Fitting Local, Low-Dimensional Parameterizations of Optical Turbulence Modeled from Optimal Transport Velocity Vectors Tegan H. Emersona,∗, Jonathan M. Nicholsb a Pacific

Northwest National Laboratory, Data Science and Analytics Group, Seattle, WA Research Laboratory, Optical Sciences Division, Washington, DC.

b Naval

Abstract This work exploits a connection between optimal transport theory and the physics of image propagation to yield a locally low-dimensional model of turbulence-corrupted imagery. Optimal transport produces an invertible, pixel-wise linear trajectories to approximate the globally nonlinear turbulence between a clean and turbulence corrupted image pair. We use the low-dimensional model to fit subsets of the optimal transport vector fields and stitch the local models into a surrogate for the global map to be used for image cleaning. Experiments are performed on laboratory generated data of beam propagation using different values of the Fried parameter (a scale measuring turbulence coherence) as well as a toy data set. The results suggest this is a fruitful direction, and first step, towards using multiple realizations of turbulence corrupted images to learn a blind surrogate for the optimal transport vector field for image cleaning. Keywords: Turbulence, Optimal Transport, Velocity Fields, Geometric Optics 1. Introduction Understanding, and subsequently mitigating, the influence of a turbulent atmosphere on a propagating electric field remains a challenge to those working on radar [1] and imaging [2] applications. Despite early theoretical advances (see e.g., [3, 4]), researchers continue to improve propagation models [5, 6], the computational tools required to solve those models [7, 8], and the techniques to mitigate the effects of turbulence on acquired signals [9, 10]. In prior work, we formulate a transport-based propagation model for imagery, and then show how one can quickly and efficiently estimate that model using only a pair of turbulence-corrupted images [11]. The model predicts the spatially localized image “dancing” observed in video sequences as well as the blurring effects of turbulence [2]. Importantly, the solution of this model is also shown to approximately solve the equations governing the problem physics. Thus, our estimated transport models provide us a forward model consistent with the geometric effects of optical turbulence. By contrast, existing methods are phenomenological in nature, correcting for the effects of turbulence using global registration and/or deblurring operators [10, 12]. The Optimal Transport (OT) model is a global, inevitable coordinate transformation that captures the geometric optics effects of turbulence between a clean and ∗ Corresponding

author. Email addresses: [email protected] (Tegan H. Emerson), [email protected] (Jonathan M. Nichols) Preprint submitted to Pattern Recognition Letters

corrupted image pair. The coordinate transformation approximates the highly nonlinear global behavior of turbulence with pixel-wise linear trajectories described by a velocity vector field. Current model estimates require both clean and turbulence-corrupted images; one can invert the resulting vector field to obtain the former from the latter. However, to perform the more practically useful blind transport one would require a means of estimating the vector field using only corrupted imagery. As with all estimation problems, the quality of the estimates is improved by maximizing the number of data while minimizing the number of parameters to be identified. We view development of a low dimensional representation of the vector field as critical to this end goal. This paper therefore develops a low-dimensional model of the OT vector fields associated with turbulencecorrupted imagery. We show that locally, a Discrete Space Dynamical Systems (DSDS) model captures the coordinate transformation of a single ray of light (i.e., geometric optics [13]). The coherence of the electric field provides a natural length scale over which the model can be applied. Each local DSDS is shown to be governed by a single inhomogeneous, first-order, linear Ordinary Differential Equation (ODE) and therefore admits only a few possible solutions described by a few parameters. We then show how one can stitch these local, low-dimensional models into a global transport map appropriate for image inversion/turbulence correction. We begin with a brief introduction to OT and imaging in Section 2.1. Next, in Section 2.2 we describe turbulence as a complex optical system, the effects of which can be locally described using the ray transfer matrices of geometric October 21, 2019

optics. A mathematical model tying the OT vector fields and ray transfer matrices to a low-dimensional parameterized system is presented in Section 3. Some preliminary experiments and their results are reported Section 4. Given the resolution and noise performance of modern imaging systems, atmospheric turbulence is now often the limiting factor in both long range imaging and free-space optical communications [14, 15, 16]. Overcoming atmospheric effects is therefore essential to improving performance in either application. By incorporating the physics of corruption in our low-dimensional image model, we hope to ultimately mitigate its influence on bit error rates (for communications) and/or the ability to recognize objects in a scene (long range imaging). We have therefore included a discussion of the advantages, disadvantages, and implications of the proposed model in the context of turbulence mitigation in Section 5. Future directions and conclusions are presented in Section 6.

Such a solution can be realized by finding the map that minimizes Z f ∗ (~x0 , z) = kf ∗ (~x0 , Z) − ~x0 k2 ρ(~x0 )d~x (3) f

subject to the continuity constraint, Eqn. (2a). The coordinate transformation f ∗ (~x0 , Z) can be viewed as a set of instructions, governing the movement of each parcel of intensity comprising an image. The mapping is exact in a homogeneous atmosphere and is providing a linear approximation to intensity movement in the event that the atmosphere is turbulent. A benefit of this formulation is that solutions to (3) can be numerically obtained using recent computational advances from the field of optimal transport and a single pair of clean/corrupted images as input [18]. Thus, we have a practical means of estimating an invertible model that approximately solves Equation (1) and can be easily applied to observed imagery. An example of the output optimal transport velocity vector field can be seen in Figure 1.

2.1. Optimal Transport and Imaging As a starting point, we take the homogeneous wave equation governing the propagation of an electric field through a medium with real-valued refractive index n(~x, z)2 = 1 + 2η(~x, z). Assuming a wave model who’s propagation is dominated by a transverse planar form in the direction of propagation, ρ(~x, z) exp[i(ωt + k0 z + φ(~x, z)], the complex wave equation becomes two, realvalued Partial Differential Equations (PDEs) governing the magnitude ρ(~x, z) and phase φ(~x, z) of the field. The PDEs are given by

2.2. Discrete Space Dynamical System Over the last 50 years it has been observed that coherent light propagation through a turbulent atmosphere is comparable to light propagation through an optical system/component [3, 19]. Loosely we can think of one realization of imaging through atmospheric turbulence as the superposition of multiple coherent beams propagating through an optical component with spatially varying index of refraction. Furthermore, this global, spatially varying index of refraction optical component can locally be approximated by thick lenses. Let us begin by assuming that we can assign a single coherent beam of light, i.e. a ray, to each pixel in the image. Simultaneously consider the magnitude of that ray to be a packet of light/intensity. As that ray moves in the direction of propagation it encounters the spatial variations in the index of refraction and deviates from its turbulencefree straight trajectory. This is the geometric optics interpretation of atmospheric turbulence acting as an optical system. In our approach we assume that the ray encounters N discrete changes in the index of refraction. After each change the new location, in the transverse plane, can be described by a ray transfer matrix. We use Bn−1 to denote

(1a) (1b)

where v(~x, z) ≡ ∇X φ(~x, z)/k0 is the spatial gradient of the phase in the transverse plane, denoted by coordinates ~x = (x, y). Thus, we have that evolution of the electric field defining the image is being “driven” by the perturbations to the refractive index while the image intensity obeys the continuity equation (see [11, 17] for a complete derivation and discussion of the model given by Eqn. 1). In the absence of refractive index variations, the solution to (1) can be defined by an invertible transport map, f ∗ (~x0 , z) that evolves the initial intensity locations, ~x0 forward along the axis of propagation. In this Lagrangian view of the problem, the magnitude and phase at propagation distance Z are given by ρ(~xZ ) = |Jf ∗ (~x0 , Z)|−1 ρ(~x0 ) k0 ∗ (f (~x0 , Z) − ~x0 ) ∇X φ(~xZ ) = Z

𝜌(𝑥⃗/ )

Figure 1: (left) Turbulence-free image of a coherent beam, (middle) the optimal transport velocity vector field, and (right) the turbulence-corrupted image of the beam.

2. Theory and Calculations


∂ρ(~x, z) + ∇X · ρ(~x, z)v(~x, z) = 0 ∂z ∂v(~x, z) + (v(~x, z) · ∇X )v(~x, z) = ∇X η(~x, z), ∂z

𝑓 𝑥⃗% , 𝑧 = 𝛻, 𝜙(𝑥⃗. )

𝜌(𝑥⃗% )

(2a) (2b) 3

the ray transfer matrix for the nth change in the index of refraction. Consequently, the location of the pixel intensity can be modeled by a discrete-space dynamical system given by       xn xn−1 xn−1 = + Bn−1 . (4) yn yn−1 yn−1

to the transverse coordinates of the region’s optical axis to yield # "     ∂φ 1 x Ox ∂x = F + (10) ∂φ y Oy Z ∂y It was shown in [11] that under the paraxial assumption #  "  ∂φ dx ∂x dz = . (11) ∂φ dy

After performing a few iterations of the dynamical system one can see that the dynamical system can be written in terms of the initial coordinates and some 2 × 2 matrix, Fn−1 that is a function of all the prior Bi . Generically we have that       xn x0 x0 = + Fn−1 (5) yn y0 y0

∂y

Finally, we can write Equation 10 as  dx      1 x(z) Ox (z) dz = F + . dy y(z) Oy (z) Z dz

Within the imaging context, all of the intermediate positions of the pixel intensity are unknown and we would only observe the final coordinates in the transverse plane ~xZ for the intensity that originated at ~x0 . Thus, we approximate the real non-linear trajectory of the pixel intensity with a single linear trajectory described by one ray transfer matrix F~x0 . Another interpretation of this model is that the complex optical system comprised of multiple simple optical elements (thin lenses) can be described by a single thick lens. The final coordinates for a packet of light originating at ~x0 are equivalently given by the optimal transport coordinate map,   xN = f (~x0 , Z). (7) yN

3. Materials and Methods: Parameterization Model

x0 y0



.

A Low-dimensional

In the context of image processing the Fried parameter provides a number of pixels over which coherent structure may be expected. This critical number, henceforth to be denoted r0 , provides an upper limit on the diameter of a circle where the velocity vectors within the grid can be assumed to be associated to the same ODE (i.e. the optical turbulence can be can be modeled as the same lens due to coherence). Thus, each incident ray of light can locally (within r0 ) be modeled as a different initial condition to a single linear ODE (12). For an autonomous, first-order, linear, homogeneous ODE of the form ~x(z)0 = A~x(z) the solutions are well understood and can be described by exactly one of ten possible forms; these are shown in Figure 2 as phase portraits [21]. Each phase portrait class is tied to properties of the linear operator A governing the ODE. Specifically, the phase portrait can be determined based on the trace and determinant of the linear operator. Solving an autonomous inhomogeneous first-order linear ODE, (~x(z)0 = A~x(z) +~s(z) with ~s(z) 6= 0), first requires solving the homogeneous case. The general solution can then be identified based on the form of ~s(z) for several function forms. Specifically, we are interested in the case where ~s(z) is a constant vector. From the optimal transport model we obtain a vector field (∂φ/∂x, ∂φ/∂y). We have one matrix containing all the partial derivatives for x and one matrix for the partial derivatives with respect to y. As we showed in Section 2.2, the optimal transport velocity vector field can be associated to a linear, first-order inhomogeneous ODE. Over a

We can further relate Equations (6) and (7) through 

(12)

Equation 12 is a first-order, linear, autonomous, inhomogeneous ODE. Thus, under the paraxial assumption the evolution of the coordinates of a parcel of light intensity can be described by a DSDS which in turn can be approximated by a constant velocity path equivalent to the constant velocity solution from optimal transport. The constant velocity solution admits a first-order linear ODE which can be characterized by one of ten behaviors as described subsequently in Section 3.

Let the coordinate displacement of the incident ray over it’s entire path be given by       xN x0 x0 − = F~x0 . (6) yN y0 y0

1 f (~x0 , Z) − ~x0 = F~x0 Z Z

dz

(8)

As a result of [11]) we can write Equation (8) in terms of the velocity vector field, or equivalently the gradient of the electric field phase by   1 x0 v(~x0 , Z) = ∇X φ(~x, z) = F~x0 (9) y0 Z where ∇X denotes the gradient in the transverse plane. Each incident ray’s path could be characterized by a different matrix. However, over a region defined by the Fried parameter [20] the paths of all rays initially within the region should be related, i.e. the rays over a region in the transverse plane pass through a single optical system with the same optical axis. Consequently, we can drop the subscripts in Equation (9) by adding a constant corresponding 4

ˆ = VX ˜ † . We employ this pseudoinverse can be found as M based solution in the experiments that follow. This solution is a maximum likelihood estimate for M in the case that the noise values are independent, jointly Gaussian distributed. Further, note that so long as the patch contains more than six pixels this problem will be over-determined.

small spatial region of size r0 , the optimal transport velocity vectors will take on the form of one of the ten phase portrait types. Thus, we seek to exploit this relationship to correct for the effects of turbulence in imaging. Other work has leveraged the finite number of phase portraits corresponding to a first order, linear, inhomogeneous ODE [22, 23]. In [22] they analyze flow fields and flow-like textures that can arise in a variety of contexts. The authors state that optical flow can admit the types of flow fields considered but no formal model is presented related to turbulence in imaging. Moreover, the proposed model is locally fit and used to describe only local features. Alternatively, [23] is interested in both local and global modeling and approaches the problem by identifying local singular pattern regions and linearly combining these local regions through weighted linear combination to approximate a global fit. A key difference between the two works mentioned and our proposed approach is that the referenced works are exclusively interested in modeling a known vector field. As such, both papers describe metrics for measuring the goodness of the vector field fit. These metrics may be incorporated in future work to determine “good” training pairs for learning blind transport patch parameters. In future work the authors intend to consider using singular patterns and weighted linear combinations as alternative stitching algorithms and evaluate the viability for blind constructions. All of the velocity vectors from inside a circle with diameter r0 , the Fried parameter, can be used to approximate a linear operator and constant who’s corresponding phase portrait matches the velocity vector field. A linear operator is defined by four parameters and the optical axis shift is defined by two parameters. As a consequence if r0 admits a grid containing at least six velocity vectors/pixel the estimation problem will be well defined. The smaller the grid the better the fit expected. However, the smaller the grid the more stitching together of the patches there will be. A starting place for fitting the field over a patch is to solve a least squares problem for each patch. Assume r0 admits a grid of size m × m, i.e. m2 velocity vectors. Let ~vi be the ith velocity vector at the corresponding ~xi position in the grid. Define V = [~v1 , ~v2 , . . . , ~vm2 ] and X = [~x1 , ~x2 , . . . , ~xm2 ]. We seek to find F and O such that V = FX + O. This can be written as a least squares optimization problem as ˜ 2 min ||V − MX|| F

M = [F, O]

˜ = and X



We present the results of two experiments. First, we create a synthetic data set where turbulence is modeled as a horizontal shift and the resulting optimal transport map matches a human intuition of moving between the images. In this case local and global behavior are the same and no stitching is necessary. Alternatively, we consider optical fiber mode transmission data. Turbulence is modeled by passing the modes through phase screens emulating the the Kolmogorov model. For both experiments we start by computing the optimal transport map between the clean and simulated turbulence effected image. This optimal transport map appears as a large vector field. We then subdivide the field into patches identified by the grid size and the overlap/stride. Each patch of the full field is a V and the corresponding low-dimensional model parameters are learned by solving Equation 13. The set of learned parameters in turn produces a local vector field (see Figure 2). These learned fields are then stitched together to form a global surrogate for the original optimal transport map. Stitching is performed by allowing patch overlap and averaging velocity vector estimates over regions of overlap. 4.1. Toy Problem The turbulence free image in this toy example is a Gaussian column of intensity set in a near zero background. Turbulence is modeled as a basic horizontal shift to the right by 5 pixels. A 20×20 resolution is used and the patch size is 20 pixels square. First, the optimal transport map is computed between the clean and the corrupted images. Next, Equation 14 is solved using right multiplication by ˜ Finally, the corrupted image is the pseudoinverse of X. cleaned using the fitted field. Figure 3 shows the results of fitting the optimal transport velocity vectors to Equation 10 by solving Equation 14. The upper left image is the clean image while the lower left is the turbulence corrupted image. In the center column contains the cleaned images resulting from the velocity vector fields resulting from optimal transport (bottom) and the DSDS (top) to their right. Turbulence in this example was modeled by a horizontal shift and the optimal transport velocity vectors are all horizontal arrows to the right (the direction of the shift). The fitted vector field also captures this horizontal shift although it minimally disagrees with optimal transport field near the left edge.

(13)

M

where

4. Experiments and Results

X 12×m2



.

(14)

ˆ be the approximation to the velocity vector field. Let M ˜ † be the pseudoinverse of X. ˜ The solution Further, let X to the least squares optimization problem in Equation 14 5

𝐷 = 𝑇 $ /4

𝐷 = 𝑇 $ /4

Nodal Sink

Spiral Source Nodal Source

Center

T

(Trace of A )

Spiral Sink

Stable Saddle-Node

Unstable Saddle-Node

(0,0)

Saddle

D (Determinant of A)

Figure 2: The classification of the phase portrait of a first-order, linear homogeneous ODE based on the trace (T) and determinant (D) of the linear operator A. Examples of each of the canonical phase portraits for first-order, linear, homogeneous ODEs produced by linear operators of with each possible combination of determinant and trace relationship are also shown.

Turbulence Free

DSDS Cleaned

DSDS Vector Field

Turbulence Corrupted

Optimal Transport Cleaned

Optimal Transport Velocity Vectors

(far right column) for Fried parameters r0 ≈ 0.38cm and ◦ 0 ≈ 0.12cm. The reconstructions generated from stitched local fields of size 32, 16, 8, 4, and 2 are shown between each clean-corrupted image pair (from left to right). The corresponding global stitched vector fields are provided in the row below the images. Based on the Fried parameter, aperture size, and resolution of the images we determine that the maximal diameter for turbulence coherence for the experiment with r0 ≈ 0.38cm is 50 pixels. For the experiment with r0 ≈ 0.12cm the maximal diameter is 17 pixels. Inscribing a square grid inside a circular diameter of the indicated pixels identifies rectangular grids of size 35 pixels and 12 pixels, respectively, over which coherence in the turbulence is expected. The grid size bounds for coherence provide upper limits for which vector field fitting should begin to capture the patterns and structure of the turbulence. Any grid size lower than this bound is valid and the finer the resolution the better the expected fit. One can see that by a square grid of size 8 the global structures of the turbulence map/OT vector field are being captured by the stitched fields. As the size of the grid decreases we see the intensity returning to the correct regions (using the turbulence free image as reference) of the image. In future work we seek to explore more sophisticated stitching and overlap schemes to improve the quality of reconstructions for grid sizes nearer to the Fried parameter determined upper limit.

Figure 3: Reconstructions and vector fields generated for a toy problem. The clean image is a blurred column of intensity and the turbulence corrupted image is a horizontal shift of the column. The OT vector field is shown in red and the DSDS fitted vector field is shown in blue.

4.2. Free-Space Optical Communications Consider next an optical communications example whereby different modes of a coherent Laguerre-Gauss beam are used to encode bits of information [24]). The challenge is to identify the clean image given only the turbulence-corrupted image [15]. Experimentally we generate several such clean/corrupted image pairs by passing the beam through a Kolmogorov phase screen (see e.g., [25, 26]). For each clean/corrupted image pair we estimated the vector field that solves Equation 1 using the variational approach described in [27]. For each OT vector field we then fit local regions of the velocity field using the pseudoinverse method proposed in Section 2.2. These local regions are stitched together using vector averaging over a sliding window. We consider local patches of size 32 × 32, 16 × 16, 8 × 8, 4 × 4, and 2 × 2 with overlaps/stride of size 15, 7, 3, 1, and 0 respectively. The quality of the the reconstructions derived from the inverted vector fields are also shown. Figure 4 shows the turbulence corrupted beam profiles (far left column) and the turbulence free beam profiles

5. Discussion One of the chief motivations for developing lowdimensional models is to facilitate their identification from observed data. Indeed, the ability to successfully identify a model is predicated largely on maximizing observations while minimizing the number of parameters to be identified. The local, low-dimensional models proposed in this work provide one reasonable path toward this goal for com6

Turbulence Corrupted

32 x 32 Grid Size, Overlap=15

16 x 16 Grid Size, Overlap=7

8 x 8 Grid Size, Overlap=3

4 x 4 Grid Size, Overlap=1

2 x 2 Grid Size, Overlap=0

Turbulence Free

𝑟" ≈ 0.38𝑐𝑚 Optimal Transport Velocity Vectors Corrupted

Grid Size=32

Grid Size=16

Grid Size=8

Grid Size=4

Grid Size=2

Clean Image

𝑟" ≈ 0.12𝑐𝑚 Optimal Transport Velocity Vectors

Figure 4: Experimental results for turbulence with Fried parameters r0 ≈ 0.38cm and r0 ≈ 0.12. The cleaned images resulting from the stitched, fitted vector fields. The far left column shows the turbulence corrupted image and the far right column shows the turbulence free mode. Stitched fields producing the cleaned images are located below the reconstructed image. In red are the optimal transport derived vector fields.

plex (high-dimensional) image data that have been influenced by atmospheric dynamics. The intended application is the removal of turbulenceinduced artifacts from imagery. Accomplishing this “blindly” requires one to identify an image model in accordance with a chosen measure of image quality. To this end, blind de-convolution is a frequently used technique (see e.g., [10]) which couples a linear image model with some measure of image sharpness. Such an approach is phenomenological, however, in that it models qualitative features of the corrupted image rather than the physics that produced the image. In this sense, the models we propose offer an advantage. They retain the low-dimensionality, invertibility, and ease of computation afforded by the linear model while matching the physics of image distortion (see section 2.1 and reference [17]). In fact, for the application to free-space optical communications we have shown previously that this approach to modeling yields a greatly simplified classifier over one based on machine learning [15]. The chief drawbacks to this local model are 1) the need to select a local length scale over which to construct the model and 2) the need to “stitch” together multiple local patches in order to construct the global model. As

discussed in the manuscript, length scale can be appropriately selected as the Fried parameter associated with the turbulence. The stitching problem, however, will require an additional computational step whereby local model parameters are adjusted to match the OT vector fields of nearby patches. Our current research is focused on precisely this problem. We hope to utilize theory related to segmenting discrete vector field and Hodge decomposition using groundwork laid in [28], and by others, towards this end. 6. Conclusion We have posed a local, low-dimensional model for turbulence-corrupted imagery and shown that it matches the physics of the process that produced the imagery. We then used experimentally obtained, clean-corrupted image pairs to estimate OT velocity fields. Inversion of these fields is shown to recover the clean image from the corrupted image suggesting that the influence of turbulence can be well approximated via OT maps applied to an image. Generating the OT velocity vectors requires knowledge of both the turbulence-free image and turbulencecorrupted image pair. Given a local, low-dimensional pa7

rameterization of geometric optics we can stitch together local models into a global vector field that closely matches the optimal transport map and resultantly the geometric optics effects of turbulence. In future work the authors seek to develop algorithmic frameworks that utilize multiple turbulence corrupted images and optimal transport maps, in the absence of a clean image, together with machine learning to learn the local models that can be stitched together to recover a clean image. With sufficient training data, the authors believe it will be possible to use machine and deep learning techniques to learn the local, low-dimensional parameterizations of OT velocity fields without clean reference images. More sophisticated computer vision methods will be used to stitch the learned local models into a surrogate for the global OT map. If such an approach is fruitful it would enable “blind transport.” Blind transport would combine aspects of fluid mechanics, OT, dynamical systems, machine learning, and computer vision into a single algorithm capable of correcting for image corruption due to atmospheric turbulence.

[16] R. Arora, G. Singh, Key challenges and resolutions for atmospheric turbulence effect removal, International Journal of Grid and Distributed Computing 11 (11) (2018) 39–52. [17] J. M. Nichols, T. H. Emerson, G. K. Rohde, A transport model for broadening of a linearly polarized, coherent beam due to inhomogeneities in a turbulent atmosphere, Journal of Modern Optics 66 (8) (2019) 835–849. [18] S. Kolouri, A. B. Tosun, J. A. Ozolek, G. K. Rohde, A continuous linear optimal transport approach for pattern analysis, Pattern Recognition 51 (2016) 453–462. [19] M. Charnotskii, Propagation of polarized waves in inhomogeneous media, Journal of the Optical Society of America – A 33 (7) (2016) 1385–1394. [20] D. L. Fried, Statistics of a geometric representation of wavefront distortion, J. Opt. Soc. Am. 55 (11) (1965) 1427–1435. [21] S. H. Strogatz, Nonlinear Dynamics and Chaos: with applications to physics, biology, chemistry, and engineering, Addison Wesley, Reading, MA, 1996. [22] A. Rao, R. Jain, Computerized flow field analysis: Oriented texture fields, IEEE Transactions on Pattern Analysis and Machine Intelligence 14 (7) (1992) 693–709. [23] W. Liu, E. Ribeiro, Detecting singular patterns in 2d vector fields using weighted laurent polynomial, Pattern Recognition 45 (11) (2012) 3912–3925. [24] T. Doster, A. Watnik, Laguerre–gauss and bessel–gauss beams propagation through turbulence: analysis of channel efficiency, Applied optics 55 (36) (2016) 10239–10246. [25] R. Lane, A. Glindemann, J. Dainty, et al., Simulation of a kolmogorov phase screen, Waves in random media 2 (3) (1992) 209–224. [26] T. Doster, A. Watnik, Machine learning approach to oam beam demultiplexing via convolutional neural networks, Appl Opt 56 (12) (2017) 3386–3396. [27] E. Haber, T. Rehman, A. Tannenbaum, An efficient numerical method for the solution of the l 2 optimal mass transfer problem, SIAM Journal on Scientific Computing 32 (1) (2010) 197–211. [28] H. Li, W. Chen, I.-F. Shen, Segmentation of discrete vector fields, IEEE Transactions on Visualization and Computer Graphics 12 (3) (2006) 289–300.

References [1] F. J. Ryan, Analysis of electromagnetic propagation over variable terrain using the parabolic wave equation, PN, 1991. [2] F. G. Smith, Atmospheric Propagation of Radiation, The infrared and electro-optical systems handbook, SPIE Optical Engineering Press, Bellingham, Washington USA, 1993. [3] V. Tatarskii, The effects of the turbulent atmosphere on wave propagation, 1971. [4] Y. A. Kravtsov, Y. I. Orlov, Geometric Optics of Inhomogeneous Media, Springer-Verlag, Berlin Heidelberg, 1990. [5] M. Charnotskii, Extended Huygens-Fresnel principle and optical waves propagation in turbulence: discussion, Journal of the Optical Society of America – A 32 (7) (2015) 1357–1365. [6] J. Rubinstein, G. Wolansky, Geometrical optics and optimal transport, JOSA A 34 (10) (2017) 1817–1823. [7] J.-D. Benamou, O. Lafitte, R. Sentis, I. Solliec, A geometric optics method for high-frequency electromagnetic fields computations near fold caustics– part i, Journal of Computational and Applied Mathematics 156 (2003) 92–125. [8] H. Liu, S. Osher, R. Tsai, Multi-valued solution and level set methods in computational high frequency wave propagation, Communications in computational physics 1 (5) (2006) 765– 804. [9] Y. Mao, J. Gilles, Non rigid geometric distortions correctionapplication to atmospheric turbulence stabilization, Inverse Problems and Imaging 6 (3) (2012) 531–546. [10] X. Zhu, P. Milanfar, Removing atmospheric turbulence via space-invariant deconvolution, IEEE Transactions on Pattern Analysis and Machine Intelligence 35 (1) (2013) 157–170. [11] J. Nichols, T. Emerson, L. Cattell, S. Park, A. Kanaev, F. Bucholtz, A. Watnik, T. Doster, G. Rohde, Transport-based model for turbulence-corrupted imagery, Applied optics 57 (16) (2018) 4524–4536. [12] R. C. Hardie, D. A. LeMaster, On the simulation and mitigation of anisoplanatic optical turbulence for long range imaging, in: Proceedings of the SPIE, Vol. 102040B, 2017, p. doi: 10.1117/12.2264064. [13] B. E. Saleh, M. C. Teich, B. E. Saleh, Fundamentals of photonics, Vol. 22, Wiley New York, 1991. [14] Photonics Research. [15] De-multiplexing vortex modes in optical communications using transport-based pattern recognition.

8