# Quantitative shadowgraphy and proton radiography for large intensity modulations

## Abstract

Shadowgraphy is a technique widely used to diagnose objects or systems in various fields in physics and engineering. In shadowgraphy, an optical beam is deflected by the object and then the intensity modulation is captured on a screen placed some distance away. However, retrieving quantitative information from the shadowgrams themselves is a challenging task because of the non-linear nature of the process. Here, a novel method to retrieve quantitative information from shadowgrams, based on computational geometry, is presented for the first time. This process can also be applied to proton radiography for electric and magnetic field diagnosis in high-energy-density plasmas and has been benchmarked using a toroidal magnetic field as the object, among others. It is shown that the method can accurately retrieve quantitative parameters with error bars less than 10%, even when caustics are present. The method is also shown to be robust enough to process real experimental results with simple pre- and post-processing techniques. This adds a powerful new tool for research in various fields in engineering and physics for both techniques.

###### pacs:

Valid PACS appear here^{1}

## I Introduction

Shadowgraphy is a technique to visualise modulations in discrete objects (1); (2) and is used extensively in our daily life. For example, when sun rays propagate through a transparent object with non-flat surface, one can readily observe the modulation of the light intensity behind the object on a screen placed a suitable distance away. This occurs because when light rays propagate, different refractive indices in the object cause the rays’ paths to be deflected, resulting in the intensity modulation. With its simplicity, the shadowgraphy technique has become a widely used diagnostic tool in many different fields in physics and engineering. Examples are diagnosing plasma wakefields (3), measuring temperatures in combustion processes (4), and characterization of optical systems (5). A similar technique, that of proton radiography, is also widely employed to diagnose the structure in laser-plasma experiments (6); (8); (7); (9); (10); (11). In proton radiography, instead of using light rays, a proton beam is fired into the plasma. The electric and magnetic fields inside the plasma deflect the protons’ trajectories. Proton beams are both highly laminar and have discrete divergence angles that allow magnification of the object, provided that the screen is placed far enough away from the object. By looking at the intensity modulation of the proton beam on the screen, one can see the structure inside the plasma with resolution. Among the applications of proton radiography are studies of experimental magnetic reconnection phenomena (6); (7), observing solitons (8), laser channeling in plasmas (9); (10), and observing the Weibel instability (11).

One emphasises here that both the shadowgraphy and proton radiography techniques share the same underlying principle. Thus, one can refer to proton radiography as shadowgraphy and vice-versa, without losing generalities. We will do this throughout this paper. By doing so, we show that this new approach provides a powerful new quantitative diagnostic tool for high-energy-density plasma science.

Although shadowgraphy is widely used in plasma science, in many cases it is used as a qualitative analysis tool (9); (10). There have been many efforts in the past to retrieve the quantitative information from shadowgrams, but it has only been possible, so far, in limited cases where the intensity modulations are small. This is done mainly by employing Poisson’s equation solver (1); (13); (14); (15) or by using the diffusion equation (16) for specific cases (5); (17). The equation for small intensity modulation of shadowgraphy was also obtained by Pogany, et al. (13) using phase contrast approach and Fresnel diffraction. The non-linear nature of shadowgraphy makes it a challenging task for large modulation cases. Some experiments also make use of a grid to estimate the deflection of the beam (12); (7). However, the technique depends on the grid resolution and it becomes harder to estimate when the feature to be observed is about the same size as the grid resolution or smaller (9); (10).

In this paper a method to retrieve quantitative information from shadowgraphic images for large intensity modulations, without using a grid, is presented. A coherent beam for optical shadowgraphy is also assumed throughout. By retrieving the quantitative information one can interpret phenomena in much greater detail, and thus provide a greater understanding of the diagnosed system. Section II provides equations underlying the shadowgraphy and proton radiography techniques as well as the basic tools used in the methods section. Then the new method to retrieve the quantitative information is explained, as well as its implementation, in section III. Benchmarking with simulations is presented in section IV and tests on real experimental results in section V. Section VI concludes the paper.

## Ii Theory

### ii.1 Deflectometry

If beams of light or charged particles are fired into deflecting objects along the -axis, they will be deflected by an amount of

(1) |

where is the deflection potential. The deflection potentials for optical shadowgraphy and proton radiography cases, respectively, are (14); (15)

(2a) | ||||

(2b) | ||||

(2c) |

where is the refractive index of the object in light shadowgraphy cases, and respectively are the electric and magnetic potential, , , and are the charge, energy, and mass of the particle in the beam, respectively. It is assumed that the beam propagates in straight lines during the interaction with the object.

With each deflection, beams at position on the object plane are mapped to position on the screen according the equations below,

(3) |

where is the distance between the object and the screen. These equations assume the beams are collimated before the interaction with the objects. For diverging beams using the paraxial approximation, one can simply replace and , where is the distance from the beam source to the object.

From the mapping equations, one can obtain the intensity of the beam on the screen as (15)

(4) |

where is the beam intensity on the screen without deflections. The term is the determinant of the Jacobian matrix of with respect to . The Jacobian in the denominator is what makes shadowgraphy cases non-linear for relatively large or . Moreover, if or is large enough, it can make the determinant of the Jacobian matrix very small, hence it causes very high intensity at some positions on the screen. This is called caustic.

It is assumed that the object does not emit or absorb the beam, so the total flux on the screen without the object (source profile) is the same as the total flux on the screen with the object (target profile). With this assumption, the problem can be restated as the Monge transport problem (18): how are the particles transported from the source profile to the target profile such that the total distance for all particles is minimised? This can be solved using a combination of Lloyd’s algorithm (19), Voronoi and power diagram (20), and optimization (21).

### ii.2 Voronoi and power diagram

Consider a 2D plane with several sites located on the plane. For every point on the plane, there is a site which is closest to the corresponding point. As an example, Figure 1(a) shows a plane with 3 sites and point A. Compared to the other sites, site 1 is the closest to the point A. Therefore point A belongs to site 1.

In the construction of a Voronoi diagram (20), the plane is divided by some regions. All points in a region belong to the site in the same region. Figure 1(b) is an example of a Voronoi diagram. Mathematically, the -th site at occupies a region or cell on the source plane, , where for all ,

(5) |

The equation above applies only for a case where all sites have the same weights. However, in some cases, this does not apply. A site with a larger weight tends to have a larger region compared to sites with smaller weights. A diagram resulting from weighted sites is called as weighted Voronoi diagram or power diagram. A region in power diagram is called as a power cell. In the power diagram with weights , the -th site at occupies a region or power cell on plane where

(6) |

for all . Figure 1(c) shows an example of a power diagram with more weight on site 1. In a power diagram, it is possible for a site to not be located inside its region or even have no region. Setting all weights to be uniform or zero produces the Voronoi diagram.

### ii.3 Lloyd’s algorithm

Lloyd’s algorithm is a method of dividing a bounded plane into several regions with approximately the same area. The algorithm starts by deploying randomly a number of sites on a bounded plane. Then a Voronoi diagram is constructed to divide the plane into several regions. For every region, the algorithm calculates its centroid position. The sites are then moved to the centroid position of its region, and constructing the Voronoi diagram for the new positions. The process is then repeated until any stopping conditions are reached, e.g. maximum number of iterations, minimum displacement, etc. An illustration of the algorithm can be found on Figure 1(d)-(g).

There are some cases where the plane is not uniform. If this is the case, then there are several improvements that can be made. First, the site can be deployed randomly using a simple rejection method (22). Positions on the plane with lower values tend to reject a site with higher probabilities. The rejected sites are deployed to other positions until they are accepted. Second, the centroid can be calculated by adjusting the values on the plane. It is similar to calculating the centre of mass of a 2D object with a non-uniform density.

## Iii Method

In order to retrieve quantitative information of the object from the screen, one needs the beam profiles both with and without the object in position. We refer to the beam profile without the object as the source profile, , and the profile with the object as the target profile, .

Initially, a number of sites are deployed randomly on the source plane profile with a simple rejection method mentioned above. Then, Lloyd’s algorithm is applied on the source plane profile to distribute the sites so that each site has approximately the same flux. This produces a Voronoi diagram, or a power diagram with weights . Once the Lloyd’s algorithm finishes, then the algorithm performs optimization on the weights.

Denote as the -th region on the source plane and as the -th region on the target plane as a function of all sites’ assigned weights, . Note that . Also denote and as the flux of the -th region on the source and target planes, respectively. The objective of the algorithm is to find the weights, such that transporting the flux from the source plane with intensity profile to the target plane produces the same intensity profile as the target profile, , and the total distance travelled by all regions from the source plane to the target plane is minimised. Aurenhammer (21) found that the weights can be found in the minimum of a convex function,

(7) |

where and are the -th site position and the assigned weight, respectively. It is noted that . The gradient of the function is given by

(8) |

so any gradient based optimization methods can be employed. Note that in the optimization process, the sites positions do not change. It is only the assigned weights that are changed. These equations have been employed to design surfaces of transparent objects that produce caustic designs (23).

Once the minimum of equation 7 is reached, the centroid position of each power cell in the power diagram, , is computed. From the -th power cell’s centroid position on the target plane, , and the site’s position on the source plane, , the displacement in the and directions can be obtained by . However, the displacement from source plane to target plane is obtained only at positions where the sites on the source plane are located. To fill in the displacement as a function of every position on the source plane, , sites closest to the 4 corners are first moved to the corners and then the natural neighbour interpolation is used. The sites need to be moved to the corners so that the convex hull of the sites covers all the source plane and thus natural neighbour interpolation can be used. This causes some distortion near the corners, but this can be minimised by having more sites. The result of this method is curl-free at most positions, thus the deflection potential, , can then be obtained by integrating the deflection in the or direction. We call this method ‘the power diagram method’ in the remaining sections of this paper. The complete pseudocode of this algorithm is given in Algorithm 1 where all the bold face variables show the vectors of variables for all sites. The illustration is shown in Figure 2

### iii.1 Implementation

There are a lot of basic computational geometry algorithms employed in the implementation of this method. First, to obtain the power diagram of sites, algorithms that use convex hull and transformation to dual space are employed (24). Voronoi diagram can be obtained by the same algorithm by setting all weights to zero. Bounded Voronoi and power diagrams inside a rectangle are obtained by clipping the diagram with the rectangle using the Sutherland-Hodgman algorithm (25). The Sutherland-Hodgman algorithm is employed for all polygon clippings in the implementation, since all polygons are convex in this case.

To calculate the function in equation 7, one needs to compute the weighted area (i.e. and ), weighted centroid position (i.e. and ), and the weighted moment of inertia (i.e. ) of each cell in the power diagram. In order to simplify the problem, it is reasonable to pixelate the intensity profile and assume the intensity within one pixel is constant. Thus, the above parameters can be computed by splitting the cell into several polygons with uniform density within a pixel, computing the parameters for each polygon, and merging the parameters to give the parameters for the given cell (29). The area, centroid position, and moment of inertia with respect to the origin of a 2D convex polygon with vertices can be shown to be

(9a) | ||||

(9b) | ||||

(9c) | ||||

(9d) | ||||

where is the vertex position of each polygon and they are ordered in the clockwise direction. Note that . The cells’ centroids for Lloyd’s algorithm are also computed by this method.

To obtain faster convergence to the global minimum of the function in equation 7, one can use a quasi-Newton gradient descent algorithm (26). However, using a quasi-Newton algorithm requires memory, where is the number of sites, and it can be very large computationally. Thus, using the limited memory BFGS (L-BFGS) method (27); (28) can save memory while still achieving fast convergence. One can also use a multi-stage approach to minimise equation 7 faster (29). The complete implementation code of the algorithm on this paper can be found at https://github.com/mfkasim91/invert-shadowgraphy.

## Iv Benchmark with simulations

### iv.1 Magnetic field proton radiography

The first test for this method considers the case of a proton beam with energy of propagating in the positive -direction and going through a toroidal magnetic field. The toroidal magnetic field around the centre gives a line-integrated magnetic field on the object plane of

(10) |

where is the maximum value of line-integrated value of the toroidal magnetic field. This basic structure has been found in laser-plasma experiments, such as in magnetic reconnection experiments (6); (7). Even though only magnetic field cases are considered here, it can be expanded into light shadowgraphy and electric field cases using equations 2.

The transverse size of the toroidal magnetic field is assumed to be . The beam is deflected by the magnetic field and captured on the screen away. The distance from the source to the magnetic field is , thus giving magnification of 15. It is assumed that the magnetic field extent in the -direction is very small compared to and . Visualisation of the test case can be seen in Fig. 3.

The beam’s deflected velocity is where is the charge-to-mass ratio of the proton beam. Thus, the deflected angle is . This gives the deflection potential as in equation 2, given .

The value of is varied from to . These values cover the cases from small intensity modulation to the cases where caustics are formed. Caustics start to appear on the screen at . The beam’s intensity modulation is shown in Fig 4. Gaussian noise is added to each image with variance about of the average intensity.

From each image of the intensity modulation, the deflection potentials are retrieved using the method explained in this paper. Then one calculates the magnitude of the line integrated magnetic field, , from the deflection potential. The retrieved value is then compared with the original value to benchmark the method.

The images of the retrieved line-integrated magnetic field are shown in Fig 5(a). Comparison between the peaks of the retrieved values of line-integrated magnetic field and the original values are presented, as well as their relative errors. To see the improved performance of the method described in this paper, the line-integrated magnetic field profiles are retrieved using Poisson’s equation solver and compared. Note that no noise has been added to the intensity images for the Poisson’s equation solver case. These comparison results are shown in Fig 5(b).

It can be seen that the retrieved line integrated magnetic field gives very good agreement with the original profile, even when caustics are formed. The error on the retrieved value increases just when the branches of the caustics are relatively distinguishable. This is because in regions between the caustics branches, the beams are coming from more than one different position on the object plane, while this method assumes that each region on the target plane is formed from one region on the object plane only. However, one can still infer magnitude of the deflection potential in this case within some error.

A slightly higher relative error at small values of is caused by the noise. The intensity modulation at that point is comparable to the noise. As the intensity modulation gets larger, the effect of noise seems to be weaker.

It is observed that this method amplifies low-frequency components of the image and reduces the high-frequency components. It makes the power diagram method somewhat less robust to low frequency noise, but robust to high frequency noise. This can be solved by applying a high frequency filter to the image before it is processed using the power diagram method.

On the other hand, the retrieved value using the Poisson’s equation solver deviates significantly from the original value. The Poisson’s equation solver gives relative error of 10% when while the power diagram method gives the same relative error when . If one wants an accuracy of less than 10%, the power diagram method gives 10 times larger working range than the Poisson’s equation solver’s working range.

As an additional benchmark to the case described above, we also tried retrieving quantitative information for arbitrary magnetic field structures. Fig 6 shows the retrieval results of the magnetic field with various structures with the same setup as the previous case. Each image has a size of pixels with depth of 16 bits. The code was run on a highly parallel computer cluster using 32 cores. It takes around 2-3 hours to process one image.

In Fig 6, one can see that the retrieved magnetic field structures agree very well with the original magnetic field. It is also apparent that the retrieved magnetic field can be different from what it seems in the proton radiography image. Moreover, the size of the structure can also be different, as shown in Fig 6(e) where the structure’s size is actually smaller than it seems in the proton radiography image.

## V Tests with experimental results

An analysis method will only be useful if it can be shown to work on real experimental data and give reasonable results. In this section, the power diagram method is used to analyse experimental data from Sävert, et al. in cases of plasma wakefield shadowgraphy (3). In the experiment, a laser pulse was fired into a plasma to generate an electron density modulation wave associated with a laser-driven wakefield. Another laser pulse with much lower intensity was fired perpendicularly to the wakefield as a probe for the shadowgraphy method. The electron density fluctuations of the wakefield caused local modulations of the refractive index in the plasma. The refractive index modulation in the plasma caused the probe’s path to bend so that some parts of the probe were brighter than others at the detector. The refractive index of a plasma with density profile is

(11) |

where is the vacuum permittivity constant, and are the electron’s charge and mass, respectively, and is the frequency of the light. Using equations 2 and 11 with approximation , the deflection potential for light in a plasma is

(12) |

Thus, the information that can be retrieved from shadowgrams using the power diagram method is .

One of the main challenges in inverting the shadowgrams for real experimental data is the non-uniformity of the probe’s unmodulated intensity, i.e. the intensity profile without deflection from objects. This can be a big problem because the inversion processes from shadowgrams to deflection potentials amplify low frequency components. Even though this can be solved by taking the intensity profile without the object, the data is not usually available or reliable because of shot-to-shot variations. Therefore, a straightforward solution is to apply a high-pass filter to either the shadowgrams and/or the resulting deflection potentials.

The pulse that drives the wakefield has a wavelength of 810 nm, duration of 35 fs, and peak intensity . The probe pulse has the same wavelength, but with shorter duration, 5.9 fs. The shadowgrams for a plasma with density are shown in Fig. 7(a). Using the power diagram method and equation 12, it is possible to infer the line-integrated relative electron density modulation, from the shadowgrams. The inverted results from the shadowgrams are shown in Fig. 7(b), where the grey scale shows the value of the line-integrated relative electron density modulation, , as well as their 1D-cross-sections at the centre of the wakefield in Fig. 7(c). The figures still show the wakefield features with additional information of . It is shown from Fig. 7 that the power diagram method in this paper is robust enough to analyse real experimental results with additional preprocessing and post-processing. It should be noted that we have neglected relativistic effects in the plasma, e.g. the plasma electronsâ mass increase, which may occur during their interaction with the high-intensity driver pulses. To model this, multi-dimensional numerical simulations need to be applied (30)

## Vi Conclusions

We have presented a new method to retrieve quantitative data from shadowgraphic images. In the cases considered in this paper, a beam propagates through an object, gets deflected by it, and is then captured on a screen. The intensity modulation on the screen acts as the input and the deflection potential of the object is regarded as the output of this method. It assumes the beam propagates in straight lines while interacting with the object. Besides shadowgraphy, the method in this paper can also be applied to proton radiography cases.

The method has been benchmarked for a toroidal magnetic field case, which has been found in some laser plasma experiments, and a plasma wakefield shadowgraphy case. In some test cases, the method successfully retrieved the deflection potential profiles with relative error less than 10% for large intensity modulation, even for cases where caustics are present. It is also tested using arbitrary structures of the diagnosed objects and gives very good results in retrieving structures with their quantitative parameters. Moreover, it has been shown that the method is also robust to noise, especially high-frequency noise. This extends the working range of the Poisson’s solver equation by an order of magnitude. It is also shown that the method can be applied to real experimental results, with some additional pre-processing and post-processing. By applying this method, one can infer quantitative information from shadowgraphy images with high accuracy. This opens up a new dimension of research in a wide range of areas in engineering and physics.

###### Acknowledgements.

The authors would like to acknowledge the support from the plasma physics HEC Consortium EPSRC grant number EP/L000237/1, as well as the Central Laser Facility and the Computer Science Department at the Rutherford Appleton Laboratory for the use of SCARF-LEXICON computer cluster. We would also like to thank ARCHER UK National Supercomputing Service for the use of the computing service. The authors would like to thank M. C. Levy for the useful discussions. We also wish to thank the UCLA/IST OSIRIS consortium for the use of OSIRIS. M.F.K. would like to gratefully thank Indonesian Endowment Fund for Education for its support. The authors gratefully acknowledge support from OxCHEDS and P.A.N. for his William Penney Fellowship with AWE plc. A.S. and M.C.K. are grateful for support by DFG (Grants No. TR18 B9, and No. KA 2869/2-1), and by BMBF (Contracts No. 05K10SJ2 and No. 03ZIK052).### Footnotes

- preprint: APS/123-QED

### References

- G. S. Settles, Schlieren and Shadowgraph Techniques: Visualizing Phenomena in Transparent Media (Heidelberg: Springer, 2001).
- P. K. Panigrahi and K. Muralidhar, Schlieren and Shadowgraph Methods in Heat and Mass Transfer (New York: Springer, 2012).
- A. Sävert, S. P. D. Mangles, M. Schnell, E. Siminos, J. M. Cole, M. Leier, M. Reuter, M. B. Schwab, M. MÃ¶ller, K. Poder, O. Jäckel, et al., Phys. Rev. Lett. 115, 055002 (2015).
- R. W. Lewis, R. E. Teets, J. A. Sell, and T. A. Seder, Applied Optics Vol. 26, Issue 17, pp. 3695-3704 (1987).
- H. Canabal, J. Alonso, and E. Bernabeu, Opt. Eng. 40 (11) pp. 2517-2523 (2011).
- P. M. Nilson, L. Willingale, M. C. Kaluza, C. Kamperidis, S. Minardi, M. S. Wei, P. Fernandes, M. Notley, S. Bandyopadhyay, M. Sherlock, et al., Phys. Rev. Lett. 97, 255001 (2006).
- C. K. Li, F. H. Séguin, J. A. Frenje, J. R. Rygg, R. D. Petrasso, R. P. J. Town, O. L. Landen, J. P. Knauer, and V. A. Smalyuk, Phys. Rev. Lett. 99, 055001 (2007).
- M. Borghesi, S. Bulanov, D. H. Campbell, R. J. Clarke, T. Zh. Esirkepov, M. Galimberti, L. A. Gizzi, A. J. MacKinnon, N. M. Naumova, F. Pegoraro, et al., Phys. Rev. Lett. 88, 135002 (2002).
- L. Romagnani, A. Bigongiari, S. Kar, S. V. Bulanov, C. A. Cecchetti, T. Zh. Esirkepov, M. Galimberti, R. Jung, T. V. Liseykina, A. Macchi, et al., Phys. Rev. Lett. 105, 175002 (2010).
- L. Willingale, P. M. Nilson, A. G. R. Thomas, J. Cobble, R. S. Craxton, A. Maksimchuk, P. A. Norreys, T. C. Sangster, R. H. H. Scott, C. Stoeckl, et al., Phys. Rev. Lett. 106, 105002 (2011).
- C. M. Huntington, F. Fiuza, J. S. Ross, A. B. Zylstra, R. P. Drake, D. H. Froula, G. Gregori, N. L. Kugland, C. C. Kuranz, M. C. Levy, et al., Nat. Phys. 11, 173â176 (2015).
- L. Willingale, A. G. R. Thomas, P. M. Nilson, M. C. Kaluza, S. Bandyopadhyay, A. E. Dangor, R. G. Evans, P. Fernandes, M. G. Haines, C. Kamperidis, et al., Phys. Rev. Lett. 105, 095001 (2010).
- A. Pogany, D. Gao, and S. W. Wilkins, Rev. Sci. Instrum. 68 (7) (1997).
- G. Izarra and C. Izarra, Eur. J. Phys. 33, pp. 1821â1842 (2012).
- N. L. Kugland, D. D. Ryutov, C. Plechaty, J. S. Ross, and H.-S. Park, et al., Rev. Sci. Instrum. 83, 101301 (2012).
- C. Graziani, P. Tzeferacos, D. Q. Lamb, and C. Li, arXiv:1603.08617 [physics.plasm-ph] (2016).
- A. Boné, N. Lemos, G. Figueira, and J. M. Dias, J. Phys. D: Appl. Phys. 49, 155204 (2016).
- G. Monge, Histoire de l’Académie Royale des Sciences de Paris, avec les Mémoires de Mathématique et de Physique pour la même année, pp. 666-704 (1781).
- M. McCool and E. Fiume, Proc. of the Graphics Interface, pp. 94-105 (1992).
- F. Aurenhammer, ACM Computing Surveys, Vol. 23, No. 3 (1991).
- F. Aurenhammer, F. Hoffmann, and B. Aronov, Algorithmica, Vol. 20, Issue 1, pp. 61-76 (1998).
- G. Casella, C. P. Robert, and M. T. Wells, Lecture Notes - Monograph Series, Vol. 45, pp. 342-347 (2004).
- Y. Schwartzburg, R. Testuz, A. Tagliasacchi, and M. Pauly, ACM Trans. Graph. 33, 4, Article 74 (2014).
- A. Nocaj and U. Brandes, Computer Graphics Forum, Vol. 31, No. 3, pp. 855-864 (2012).
- I. E. Sutherland and G. W. Hodgman, Communications of the ACM, Vol. 17, Issue 1, pp. 32-42 (1974).
- A. G. Buckley, Mathematical Programming, Vol. 15, Issue 1, pp. 200-210 (1978).
- D. C. Liu and J. Nocedal, Mathematical Programming, Vol. 45, Issue 1, pp. 503-528 (1989).
- M. Schmidt. minFunc: unconstrained differentiable multivariate optimization in Matlab (2005).
- Q. Mérigot, Computer Graphics Forum, Vol. 30, Issue 5, pp. 1583-1592 (2011).
- E. Siminos, S. Skupin, A. Sävert, J. M. Cole, S. P. D. Mangles, and M. C. Kaluza, Plasma Physics and Controlled Fusion, Vol. 58, No. 6 (2015).