Photospheric Stereoscopy: Direct Estimation of Solar Surface-Height Variations

Abstract

The orbit of the Solar Orbiter mission carries it and the Polarimetric and Helioseismic Imager (PHI), which is onboard, away from the Sun–Earth line, opening up the first ever possibility of doing stereoscopy of solar photospheric structures. We present a method for a stereoscopic analysis of the height variations in the solar photosphere. This method enables the estimation of relevant quantities, such as the Wilson depression of sunspots and pores. We demonstrate the feasibility of the method using simulated Stokes-$I$ continuum observations of an MHD simulation of the solar-surface layers. Our method estimates the large-scale variations in the solar surface by shifting and correlating two virtual images, mapped from the same surface feature observed from two different vantage points. The resulting vector is then introduced as an initial height estimate in the least-squares Broyden–Fletcher–Goldfarb–Shanno (BFGS) optimization algorithm to reproduce smaller scale structures. The height estimates from the simulated images reproduce well the overall height variations of the MHD simulation. We studied which viewing angles give the best results and found the optimal separation of the view points to be between $10^{\circ }$ and $40^{\circ }$; but neither viewing direction should be inclined by more than $30^{\circ }$ from the vertical to the solar surface. The method yields reliable results if the data have a signal-to-noise ratio of 50 or higher. The influence of the spatial resolution of the observed images is considered and discussed.

1 Introduction

Stereoscopy can provide 3D information from two 2D images if these show the same scene from two different viewing angles. It has been applied to determine depths and distances of many different objects and on a wide range of scales, from microscopy to stellar parallax measurements. Stereoscopic analysis of solar features concentrated mainly on bright loops, filaments, and other coronal objects, such as coronal mass ejections.

The first studies that applied stereoscopy to solar observations used solar rotation to obtain a pair of independent images. This allowed studying the stationary 3D density structure of microwave sources (Aschwanden et al., 1995) and of coronal streamers (Jackson and Hick, 2002). Stereoscopy of non-stationary objects requires coordinated observations from two or more vantage points. The Solar TErrestrial RElations Observatory (STEREO: Kaiser et al., 2008) spacecraft launched on 26 October 2006 constituted the first space mission dedicated and designed to retrieve stereoscopic solar information. Its twin spacecraft orbit the Sun independently at about 1 AU. Most of the stereoscopic studies dedicated to the solar corona and the heliosphere since then used STEREO data (see reviews by Aschwanden et al., 2008; Aschwanden, 2011). They combined the views from the two STEREO spacecraft, sometimes also combined with the view along the Sun–Earth line. Because STEREO observes the Sun and its corona only through EUV imagers and coronagraphs, the photosphere and hence the solar surface could not be accessed by these observations.

The launch of Solar Orbiter in 2020 (Müller et al., 2020), which orbits the Sun on a highly elliptic orbit, provides a new, independent view and the opportunity to apply stereoscopy in combination with observations from another spacecraft. The observations with Solar Orbiter’s Polarimetric and Helioseismic Imager (SO/PHI: Solanki et al., 2020) can be used to extend stereoscopic techniques to also assess photospheric altitude variations. For this purpose, SO/PHI’s high resolution solar-surface observations need to be supplemented by similar observations from, e.g. a near-Earth spacecraft in the same spectral window, or at least covering the same height range.

Stereoscopy is a promising tool when applied to images with clear, identifiable edges or marks. On the solar surface, these marks are most often lacking, and the main difficulty of the application of a stereoscopic analysis to photospheric images is their diffuse brightness variation. For this reason, we present here a test of a compound stereoscopic method that we have developed to measure the photospheric altitude variation. The test uses data for which the true height variation is known so that we can compare our results with ground truth. They consist of synthetic observations calculated from sophisticated magneto-hydrodynamic (MHD) simulations of the solar photosphere. The comparison then allows us to verify if the height estimates that we obtain by our stereoscopic method are reliable.

The most prominent height variation of the solar surface (henceforth assumed to be located at an optical depth [$\tau$] $=1$) is the Wilson depression of sunspots and active regions. It was first observed and described qualitatively by A. Wilson in 1769. This depression below the altitude level of the quiet Sun’s surface is due to the strong magnetic field emanating from the sunspots. The stronger magnetic pressure and curvature forces balance lower gas pressure, temperature, and opacity within a sunspot. The height variation of the Wilson depression is therefore relevant for understanding the magnetic and geometric structure of sunspots and sets constraints on the strength and geometry of the magnetic field.

Quantitative estimates of the depth of the depression by Gokhale and Zwaan (1972) and Prokakis (1974) were based on the varying relative sizes of umbrae and penumbrae parallel and normal to the apparent radial direction of the disk as the sunspot approaches the limb (Wilson effect). From 131 spot observations, Prokakis (1974) obtained depth values in the range between 690 and 2100 km. They noticed that larger spots have deeper depressions than smaller ones. However, the method has the disadvantage that the geometrical altitude of $\tau =1$ depends on the viewing angle so that different altitude surfaces are observed as the sunspot rotates across the solar disk (Solanki, Walther, and Livingston, 1993). In addition, sunspots evolve over time as they rotate across the solar disc center, so that the relative size of umbra to penumbra may change (see Solanki, 2003).

Indirect methods assuming balanced horizontal forces across the sunspot have been used to estimate the Wilson depression. Martínez Pillet and Vázquez (1993) found a linear relation between the enhancement of the magnetic pressure and the reduction of the continuum intensity inside the umbra. For reasonable magnetic tension forces, the authors concluded that a typical depth is of the order of 600 km. Solanki, Walther, and Livingston (1993) investigated the infrared emission from a sunspot and found a linear relation between the vertical magnetic field strength and the gas temperature across cuts through a sunspot. From a quantitative model of the vertical magnetic and thermal force balance they find a depth of about 400 km for the sunspot investigated. However, a somewhat uncertain quantitative estimate of the magnetic curvature force leaves some uncertainty in this depth estimate.

Spectropolarimetric observations of a small penumbral patch were inverted by Puschmann, Ruiz Cobo, and Martínez Pillet (2010) and fitted to a 3D geometrical stationary MHD model of the atmosphere above the patch. They found that in this region, the penumbra surface lowered by about 300 km over a horizontal distance of about 4 Mm. Löptien et al. (2020) extended the work of Puschmann, Ruiz Cobo, and Martínez Pillet (2010) to an entire sunspot. They used two different constraints when fitting the observations, either a vanishing magnetic field divergence or a magnetic force in balance with the gas pressure. For the sunspot that they modeled, these two constraints gave a depression of the umbra of 500 and 700 km, respectively. Hence there is still quite some uncertainty in the depth estimates of the Wilson depression related to uncertainties inherent in the applied techniques.

In general, it was found that larger spots have a stronger magnetic field and a lower temperature in the umbra (see, e.g., Kopp and Rabin, 1992; Schad, 2014; Watson, Penn, and Livingston, 2014; Rezaei et al., 2015). It may, therefore, be expected to also see some dependence of the depth on the spot size, which is indicated by the results of Löptien et al. (2020). However, this relation is still not well known quantitatively due to the uncertainties in the depth estimates available so far. In this article, we propose a new, more direct method for the depth estimation.

The plan of this article is as follows: we first present a brief overview on how the synthetic data were generated. In the following section, we outline our analysis method and explain its assumptions and stereoscopic basics. The method is based on the stereoscopic comparison of high-resolution surface images in identical spectral ranges from different view directions. Since the Sun does not have a solid but a radiative surface layer, we also define here more precisely what we mean by the height of the solar surface. In the following sections, we test our analysis method and then give an overview of the results, exploring the performance of the method with different parameters such as observing geometry, noise, and resolution. In the final sections, we discuss the advantages, possibilities, and limitations of our approach, and we give a brief summary and outlook on the future work that our method will enable.

2 Test Data

Three-dimensional radiation MHD simulations of a typical photospheric magnetic structure (see Riethmüller et al., 2017) were carried out with the MPS/University of Chicago Radiative MHD (MURaM) numerical simulation code (Vögler et al., 2005). The simulation results after a statistical equilibrium is reached are considered to be relatively close to real observational data, since simulations carried out with MURaM have been shown to agree well with observational data (e.g. Schüssler et al., 2003; Keller et al., 2004; Hirzberger et al., 2010; Riethmüller et al., 2014). The output from such a simulation will therefore be used here as test data for our surface-height reconstructions. The test image extracted from the MHD simulations contains small unipolar flux concentrations (pores and sunspots) surrounded by small-scale magnetic structures embedded in a larger quiescent area (Figure 1). We use an $812\times 812$ pixel slice of the MURaM numerical simulations performed over a photospheric height range of 700 km. The grid resolution in the rectangular simulation box simulation was (41.67, 41.67, 15.89) km in $\xi$, $y$, and $z$, respectively. Here we use $x$ exclusively to refer to image pixels, whereas $\xi$ is used for the simulation surface coordinates.

Figure 1

Test data for views corresponding to $\theta =0$° (i.e. from directly overhead), $\theta =30$°, and $\theta =60$°. $\theta$ is defined as the angle between the line-of-sight and the surface normal.

To produce synthetic observations at different viewing angles $\theta$ (i.e. the angle between the line-of-sight and the surface normal), a standard procedure is applied to the forward calculation mode of the Stokes-Profiles INversion-O-Routines (SPINOR) radiative-transfer inversion code (Solanki, 1987; Frutiger et al., 2000) to compute synthetic Stokes spectra around the Fe i 6173 Å line as observed by SO/PHI and the Helioseismic and Magnetic Imager onboard the Solar Dynamics Observatory (SDO/HMI: Scherrer et al., 2012). The applied spectral synthesis procedures have already been described by Valori et al. (2022).

The resulting spectra are used as input to the SO/PHI Software siMulator SOPHISM (see Blanco Rodríguez et al., 2018) to produce SO/PHI-HRT-like observations. These data contain the synthesized contributions of each cell along each ray through the MURaM simulation cube. These contributions can be described in terms of response functions (Landi Degl’Innocenti and Landi Degl’Innocenti, 1977; Ruiz Cobo and del Toro Iniesta, 1994), which we use as weighting functions in the line of sight (LoS) integration with optical depth $[\tau ]$. The response functions are computed from the atmospheric model that SPINOR creates for an inversion of the synthetic observations. To obtain the formation heights of the emitted rays, the response functions must be transformed so that they depend on the geometrical height rather than on the optical depth. This is achieved with the MapTau subroutine of SPINOR, which emulates the radiative transport through the MURaM simulation for each viewing angle $[\theta]$. The output is a height coordinate along the given inclined LoS, which is directly compatible with the optical-depth scale of the magnetic response functions. Rebinning this dataset to match the spatial resolution of the observation-based response functions and the subsequent weighted integration over the LoS for each ray path results in a reference height, i.e. which is the best possible computation of the “true” formation height of the observed light. The filtergrams obtained in this manner at different values of $\theta$ are used as input data to test our stereoscopic surface-height reconstruction.

In the present article, we use only data synthesized in the continuum region near the Fe i 6173 Å line (300 mÅ outside the mean line center as observed with SO/PHI). We are aware that both real and synthetic observations obtained at different inclination angles $[\theta ]$ are formed at different geometrical heights and any direct comparison is prone to intrinsic errors resulting from these formation-height differences. Therefore we try to keep the differences in the inclination angle either small (at least for the reference cases shown in Sections 4.1 and 4.2) or symmetric around $\theta =0$°. We are also aware that any comparison of real SO/PHI observations with other space-borne photospheric imagers, e.g. SDO/HMI (Scherrer et al., 2012) or Hinode/SOT (Tsuneta et al., 2008), is affected by different continuum regions and spectral windows observed by these instruments, which introduces an additional (small) error, which is not reflected in the results presented here.

3 Method

The stereoscopic comparison of images of the same scene taken from two different view points requires a precise geometric calibration of the two observing telescopes (camera calibration; see Faugeras (1993), Chapters 2, 3, and 6 and Trucco and Verri (1998), Chapters 2, 6, and 7 for details of geometrical stereoscopy and epipolar geometry). For real spacecraft data, the necessary parameters are the view direction of the telescope optical axis and the spacecraft roll angle about this axis, both in a WCS coordinate system (e.g. Thompson, 2006). Finally, the camera model of the telescope must be known in terms of internal coordinates, which then provide the relation between a view direction relative to the optical axis and a detector pixel. The comparison of the images from the two spacecraft is then performed along epipolar lines (Trucco and Verri, 1998). The height vectors that we obtain in the end then describe the height variation along epipolar surface profiles. These profiles are defined by the intersection of the associated epipolar plane with the solar surface (see Figure 2). Finally, entire 2D surface areas can be scanned by successively processing stepwise nearby epipolar lines.

Figure 2

Three nearby epipolar planes for spacecraft at positions $\boldsymbol {r}_{1}$ and $\boldsymbol {r}_{2}$ and the epipolar plane intersections with the solar surface. Each of these intersections constitutes an epipolar surface profile on which our stereoscopic analysis is independently applied.

The conditions for application of stereoscopy depend strongly on the difference between the viewing angles from the two spacecraft at the points on the surface profile. These viewing angles $[\theta ]$ from the two spacecraft are defined as the angle between the local surface normal and the respective LoS from the spacecraft. These angles vary along the surface profile, roughly by one degree per degree on the solar surface.

In this article, we use the simulated test data to verify that our method gives reliable results. For the test data used here, the assumed spacecraft configuration and the surface geometry are very much simplified. The simulated surface area is planar and contains the $\xi$–$y$ axes. We call $\xi$ the length coordinate along the surface profile line. The virtual observing spacecraft are therefore in the $\xi$–$z$ plane in the direction of the assumed view angles $\theta$ and at infinite distance.

The origin and the scaling of this surface coordinate is quite arbitrary. Any surface profile line maps to the epipolar line in the images. The images are oriented so that the epipolar lines are oriented along their $x$-pixel axis. Surface profiles on different epipolar planes are parallel and can then be labeled by the $y$-coordinate.

To determine the height variation, we need to identify a unique physical reference point $\boldsymbol {p}^{\mathrm{ref}}$ on the surface-profile line in both images. For real observational data, this point can then be triangulated and its absolute height $h^{\mathrm{ref}}$ determined (see appendix in Zhukov et al., 2021). The precision of this absolute reference height then depends on the precision with which the spacecraft parameters and the view geometry are known. For the test data, the absolute height $h^{\mathrm{ref}}$ refers to the level $z=0$ assumed when the images were created by radiative-transfer integrations from the simulated MHD box. In this article, our main emphasis is on the relative variation of the height $h(\xi )$ that we can obtain for the solar topography.

The procedure that we present below relies on a linearization relative to $\boldsymbol {p}^{\mathrm{ref}}$. This enables the algorithm to perform reasonably rapidly, and the height variations $h(\xi )$ obtained relative to $h^{\mathrm{ref}}$ are much more precise than the absolute height. The drawback is that the length of the surface-profile line has to be limited to a distance from $\boldsymbol {p}^{\mathrm{ref}}$ where the linearization error remains less than errors due to other sources, e.g. data noise and uncompensated image distortion. For the viewing geometry of our test data, the linearization holds exactly.

The next processing step toward the stereoscopic reconstruction is to map the image-intensity data of the two spacecraft to virtual-intensity profiles (VIPs) on the surface-profile line. The goal is to obtain for each image a VIP on the surface at altitude $h^{\mathrm{ref}}$ that would produce, within digitization errors, the observed image intensity.

A basic assumption that we make in the following is that the vertical extent of the response function to temperature $[R_{T}(\xi ,z)]$ with depth $z$ in the MHD cube, at position $\xi$ along the profile can be replaced by

$$ R_{T}(\xi ,z) \rightarrow \delta (z-h(\xi ))\int R_{T}(\xi ,z)\, \mathrm{d}z ,$$

(1)

i.e. by a height-integrated effective response concentrated at height $h(\xi )$. Therefore the height $h(\xi )$ that we reconstruct aims to capture the vertical barycenter of the true $R(\xi ,z)$. Formally, we can consider Equation 1 as the definition of the height variation $h(\xi )$ that we want to obtain. Note that the MapTau routine computes the response functions relative to the depth $z=0$, which denotes the horizontal average of the optical-depth layer $\tau =1$ for vertical emergence (see also Martínez Pillet and Vázquez, 1993; Mathew et al., 2004; Löptien et al., 2018). Because of this arbitrary definition, the $z$-dimension and, thereby, also the height $h$ denote only relative quantities and must not be considered as absolute height levels in the photosphere.

Practically, this assumption can be justified because we expect the horizontal variation of $h(\xi )$ to well exceed the width of the true response function. The problems of this assumption, which ensues when we compare views with different view angles $\theta$, will be discussed in Section 4.3.

By construction, the reference point $\boldsymbol {p}^{\mathrm{ref}}$ must lie on the surface-profile line at some coordinate value $\xi ^{\mathrm{ref}}$. From the triangulation of $\boldsymbol {p}^{\mathrm{ref}}$ we find its respective pixel coordinate $x^{\mathrm{ref}}_{k}$ on the epipolar line of image $k=1$ or 2. For our test data, the epipolar line is exactly along the image $x$-axis, but in general it may be obliquely oriented in the image. In this case, $x_{k}$ is the length coordinate along the epipolar line scaled to one pixel width per $x_{k}$-unit with arbitrary coordinate origin. With these definitions, we have a simple linearized mapping between the pixel coordinate in image $k$ and the surface coordinate (see Figure 3)

$$ x_{k}(\xi ) = x^{\mathrm{ref}}_{k} + \frac{1}{{\Delta \xi }_{k}}((\xi -\xi ^{ \mathrm{ref}})\cos \theta _{k}-h(\xi )\sin \theta _{k}), $$

(2)

where $\xi$ and $h$ are scaled to the same units, and $\Delta \xi _{k}$ is the surface element of one pixel of image $k$ for a vertical view. The VIP along the surface-profile line then is

$$ {J}_{k}(\xi \mid h) = I_{k}(x_{k}(\xi )) \left | \frac{\mathrm{d}x_{k}}{\mathrm{d}\xi}\right |{\Delta \xi }_{k}, $$

(3)

where $I_{k}(x)$ is the intensity along the epipolar line in image $k$. We mention the dependency on the function $h$ explicitly as conditional input for ${J}_{k}$ because it is required in the mapping Equation 2 and the Jacobian derivative $|\mathrm{d}x_{k}/\mathrm{d}\xi |$. The factor ${\Delta \xi }_{k}$ in Equations 2 and 3 compensates different resolutions of different images $k$, whereas the Jacobian of the mapping of Equation 2 takes account of the modified intensity seen on the inclined true surface $h(\xi )$. Obviously, ${J}_{k}(\xi \mid h) = I_{k}(x_{k}( \xi ))$ is independent of $h$ for a vertical view $\theta _{k}=0$. Note also that Equation 2 has a unique inverse only as long as $\tan \theta _{k}<\mathrm{d}h/\mathrm{d}\xi$. If this limit is exceeded, then the local slope of $h(\xi )$ is steeper than the inclination of the view direction causing part of the surface to remain hidden from this view.

Figure 3

Illustration of Equation 2 for one image with inclined view direction.

Now the VIPs ${J}_{1}(\xi )$ and ${J}_{2}(\xi )$ are directly comparable, and if our assumption about the response function holds and $h(\xi )$ was chosen correctly, then they should be identical. Therefore their difference tells us if and where $h(\xi )$ was chosen incorrectly. In particular, an incorrect height in some region would result in a local shift (disparity) of identical structures in ${J}_{1}(\xi )$ and ${J}_{2}(\xi )$. We are therefore only sensitive to $h(\xi )$, where the image intensity along the epipolar line has measurable small-scale variations.

To get rid of different intensity offsets and calibrations, we do not compare ${J}_{1}(\xi \mid h)$ and ${J}_{2}(\xi \mid h)$ directly, but a variant which is windowed, centered, and scaled relative to a point ${\xi _{\mathrm{cent}}}$ on the surface profile. The transformation is conveniently performed in three successive steps (we use <…> to indicate a windowed integration with window $w$)

$$<J_k>({\xi }_{\mathrm{cent}},{\xi }_{\mathrm{width}}\mid h)=\frac{1}{{\xi }_{\mathrm{width}}}\int J_k({\xi }^{\prime }\mid h)w(\frac{{\xi }^{\prime }-{\xi }_{\mathrm{cent}}}{{\xi }_{\mathrm{width}}})\mathrm{d}{\xi }^{\prime }\text{(windowing),}$$

(4)

$$C_k(\xi ,{\xi }_{\mathrm{cent}},{\xi }_{\mathrm{width}}\mid h)=J_k(\xi \mid h)-<J_k>({\xi }_{\mathrm{cent}},{\xi }_{\mathrm{width}}\mid h)\text{(centering),}$$

(5)

$$S_k(\xi ,{\xi }_{\mathrm{cent}},{\xi }_{\mathrm{width}}\mid h)=\frac{C_k(\xi ,{\xi }_{\mathrm{cent}},{\xi }_{\mathrm{width}}\mid h)}{\sqrt{<C_k^2>({\xi }_{\mathrm{cent}},{\xi }_{\mathrm{width}}\mid h)}}\text{(scaling).}$$

(6)

The window function $w$ is symmetric, has support $[-1,1]$, and is normalized to unity so that ${\!< \! {C}_{k} \! >\!}\;=\;{\!< \! {S}_{k} \! >\!}\;=0$ and ${\!< \! {S}_{k}^{2} \! >\!}\;=1$. Then $w(\xi - {\xi _{\mathrm{cent}}})/{\xi _{\mathrm{width}}}$ is symmetric about ${\xi _{\mathrm{cent}}}$ with support $[{\xi _{\mathrm{cent}}}-{\xi _{\mathrm{width}}},{\xi _{\mathrm{cent}}}+{\xi _{\mathrm{width}}}]$. The window that we use here is smoothly tapered at both ends

$$ w(s)= \textstyle\begin{cases} 1+\rho & \text{if}\;|s|\le \rho , \\ \frac{1+\rho}{2}\left (1+\cos (\pi \frac{|s|-\rho}{1-\rho})\right ) & \text{if}\;\rho < |s| < 1, \\ 0 & \text{if}\;\rho < |s| \ge 1. \end{cases} $$

The parameter $\rho$ determines the relative length of the central constant window. We use $\rho =$0.25 throughout.

Here we propose two different approaches to determine the height $h(\xi )$ from a comparison of ${S}_{1}(\xi \mid h)$ and ${S}_{2}(\xi \mid h)$.

3.1 Correlation Method

Here we try to solve locally for a height $h({\xi _{\mathrm{cent}}})$ at the center of the window by minimizing a cost function $f_{\mathrm{cost}}$ by

$$\begin{gathered} h_{\mathrm{corr}}(\xi ,{\xi _{\mathrm{width}}}) ={\underset {h_{\mathrm{cent}}}{\operatorname{\mathrm{argmin}}}\;} f_{ \mathrm{cost}}(\xi ,{\xi _{\mathrm{width}}}\mid h=h_{\mathrm{cent}}), \\ \text{where}\qquad f_{\mathrm{cost}}({\xi _{\mathrm{cent}}},{\xi _{\mathrm{width}}}\mid h=h_{ \mathrm{cent}}), \\ =\;{\!< \! ({S}_{1}-{S}_{2})^{2} \! >\!}({\xi _{\mathrm{cent}}},{\xi _{\mathrm{width}}}\mid h=h_{ \mathrm{cent}}) \\ = 2 - 2{\!< \! {S}_{1}{S}_{2} \! >\!}({\xi _{\mathrm{cent}}},{\xi _{\mathrm{width}}}\mid h=h_{ \mathrm{cent}}). \end{gathered}$$

(7)

Here ${\!< \! {S}_{1}{S}_{2} \! >\!}$ is the weighted correlation coefficient with weight $w$. Instead of the full function $h(\xi )$, we here insert a constant height $h_{\mathrm{cent}}$ inside each window, which is considered representative for the height at the window centre ${\xi _{\mathrm{cent}}}$.

We profit from this simplification because Equations 2 and 3 change to $x_{k}(\xi )=\left .x_{k}(\xi )\right |_{h=0} - h_{\mathrm{cent}} \sin \theta _{k}/{\Delta \xi }_{k}$ and ${J}_{k}(\xi \mid h) = I_{k}(x_{k}(\xi ))$, respectively. Therefore the entire image data inside a window are just shifted in $\xi$ by $h_{\mathrm{cent}}\sin \theta _{k}$. Since $\sin \theta _{1}\ne \sin \theta _{2}$, the surface-profile observations $I_{1}$ and $I_{2}$ are shifted differently until the correlation maximizes for the correct value of $h_{\mathrm{cent}}$.

The fact that we replace $h(\xi )$ by a constant $h_{\mathrm{cent}}$ causes some inconsistency. The level of this inconsistency is larger when the window is wider, meaning a larger ${\xi _{\mathrm{width}}}$. We therefore use Equation 7 initially as a preliminary estimate for the correct $h({\xi _{\mathrm{cent}}})$ for a window size ${\xi _{\mathrm{width}}}$ sufficiently large so that we find a stable minimum with Equation 7. We then successively iterate ${\xi _{\mathrm{width}}}$ downward and compare the computed minima of $h_{\mathrm{corr}}$ with predictions from the previous, wider window sizes. For smaller ${\xi _{\mathrm{width}}}$, the solution of Equation 7 becomes increasingly prone to errors due to accidental local intensity variations and noise. The iteration is stopped, and the latest $h_{\mathrm{corr}}$ accepted as $h({\xi _{\mathrm{cent}}})$ if the minimum ${\xi _{\mathrm{width}}}$ is reached. If the discrepancy between predicted and computed $h_{\mathrm{corr}}$ becomes too large, then the iteration is stopped, and the previous $h_{\mathrm{corr}}$ is accepted as $h({\xi _{\mathrm{cent}}})$.

3.2 Optimization Method

With this approach, we avoid the inconsistency of the correlation method, and we

$$\begin{gathered} \text{minimize}\;\;f_{\mathrm{cost}}({\xi _{\mathrm{cent}}},{\xi _{\mathrm{width}}}\mid h) \quad \text{with respect to}\quad h(\xi ), \\ \text{where}\;\; f_{\mathrm{cost}}({\xi _{\mathrm{cent}}},{\xi _{\mathrm{width}}}\mid h) =\; {\!< \! ({S}_{1}-{S}_{2})^{2} \! >\!}({\xi _{\mathrm{cent}}},{\xi _{\mathrm{width}}}\mid h) +\mu \int \left |\frac{\mathrm{d}^{2}h}{\mathrm{d}\xi ^{2}}\right |^{2} \mathrm{d}\xi . \end{gathered}$$

(8)

The window $w$ is chosen here as a box function with ${\xi _{\mathrm{cent}}}$ in the center of the surface profile and ${\xi _{\mathrm{width}}}$ so that the support of $w((\xi -{\xi _{\mathrm{cent}}})/{\xi _{\mathrm{width}}})$ covers the entire profile.

The price that we have to pay is that we have to determine all height values $h(\xi )$ simultaneously by solving a non-convex optimization problem. Depending on the input data, there may be several minima, e.g. consider periodic input data $I_{k}(x)$ with period $\Delta x$. Then several solutions shifted by $\Delta h=\Delta x\,{\Delta \xi }_{k}/\sin \theta _{k}$ are possible. It is therefore necessary to provide a good initial guess $h(\xi )$ at the start of the iterative solution of Equation 8.

We use the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm (see Nocedal and Wright, 2006), which is based on a gradient search. The required gradient $\nabla _{h} f_{\mathrm{cost}}({\xi _{\mathrm{cent}}}, {\xi _{\mathrm{width}}}\mid h)$ can be calculated analytically tracing the $h$-dependence by the chain rule back to Equations 2 and 3.

For the practical calculations, an equidistant grid for $\xi$ is assumed, which is mapped to the images by Equation 2. A change of a single height element $h(\xi _{i})\rightarrow h(\xi _{i})+\Delta h$ shifts the map of the grid boundaries of cell $i$ in the image and thereby affects the values of ${J}_{k}(\xi _{j})$ only for $j=i-1,i,i+1$. However, these modifications proportional to $\Delta h$ become small in regions where the image data $I_{k}(x_{k}(\xi ))$ has little structure. Here the local gradient $\nabla _{h} f_{\mathrm{cost}}$ becomes small and the minimization procedure becomes insensitive to changes in $h(\xi )$ so that data noise may have a larger influence on $f_{\mathrm{cost}}$. To keep the noise dependence small, the regularization term in Equation 8 proportional to the regularization parameter $\mu$ was added. Depending on $\mu$, this term will smooth $h(\xi )$ and lead to some loss of resolution, but it will stabilize the solution where Equation 8 is too insensitive.

4 Results

In this section, we present the results of applying the stereoscopic analysis described in Section 3 to the test data produced in Section 2. The obtained height variations are then compared to the height variation of the response function $R_{T}(\xi ,z)$ for the vertical view.

The methods described in Section 3 include a number of parameters that may influence the results. We attempt to first determine which set of parameters yields the best results. In both correlation and optimization approaches, a free parameter is the surface grid spacing [$\Delta \xi$] to which the virtual images from the left and right view points are discretized. Another parameter is the width $\xi _{\mathrm{width}}$ of the window function $w$ used for the scaled intensity profiles $S_{k}$. We anticipate that for the optimization method, $\xi _{\mathrm{width}}$ has a minor influence on the results, whereas the regularization coefficient $\mu$ controls the main smoothing of $h(\xi )$.

We present here a few of the test results that we produced. We selected them to explain our choice of optimal parameters. The influence of these parameters on the results is here shown for a standard case, with viewing configuration of $\theta =\pm 20$° to the vertical view and images without noise. After finding the optimal set of parameters for the standard case mentioned, we investigate the parameters given by different observational circumstances such as different viewing angles, different sizes of the image pixels, and different noise levels.

To quantitatively compare our results from the height-variation estimates along a surface profile with the vertical variation of the response function, we characterize the latter by its moments relative to an arbitrary height reference defined as

$$ M_{p}(\xi )=\int (z-z_{\mathrm{ref}})^{p} R_{T}(\xi ,z)\,\mathrm{d}z. $$

The response function $R_{T}(\xi ,z)$ used here was calculated for a vertical view direction. We can now introduce a local, i.e. $\xi$-dependent, parameter $\mathrm{Err}$ that measures the distance between $h(\xi )$ and the vertical mean $h_{\mathrm{mean}}$ of the response function relative to its height extent $h_{\mathrm{width}}$:

$$\begin{gathered} \mathrm{Err}(\xi )= \frac{|h(\xi )-h_{\mathrm{mean}}(\xi )|}{h_{\mathrm{width}}(\xi )}, \end{gathered}$$

(9)

$$\begin{gathered} \text{where}\qquad h_{\mathrm{mean}}(\xi )= \frac{M_{1}(\xi )}{M_{0}(\xi )}\qquad \text{and}\qquad h_{ \mathrm{width}}(\xi )=\sqrt{\frac{M_{2}(\xi )}{M_{0}(\xi )}-h_{ \mathrm{mean}}^{2}(\xi )}. \end{gathered}$$

(10)

An example of a typical response function $R_{T}(\xi ,z)$ is given in Figure 4. The figure shows, for a single point in the quiet Sun, an example of the response-function distribution as a function of the optical depth (left) and of geometrical height (right). Note that $R_{T}(\xi ,z)$ is often strongly skewed so that the maximum of $R_{T}(\xi ,z)$ typically lies below $h_{\mathrm{mean}}$ by about $h_{\mathrm{width}}$. In addition, we define a global measure to assess the quality of our results along an entire $\xi$-profile as

$$ \epsilon =\left ( \frac{ \int \mathrm{Err}^{2}(\xi )\mathrm{d}\xi}{\int \mathrm{d}\xi}\right ) ^{\frac{1}{2}}. $$

(11)

Figure 4

Example of the response function as a function of the optical depth $\tau$ (left), and as a function of the height $z$ (right) at some position $\xi$ along the profile. In the right panel, we have marked $h_{\mathrm{mean}}$ (red vertical line) and the height extent $\pm h_{\mathrm{width}}$ (black horizontal range).

4.1 Test of the Correlation Method

Here we test the influence of the critical parameters in this method, which are the discretization of $J_{k}(x_{k}(\xi ))$ (see Equation 3) and the minimum window size $\xi _{\mathrm{width}}$ used in Equation 4. For this test, we choose as the standard surface grid spacing $\Delta \xi =\mathrm{d}x$, which corresponds to the pixel size $\mathrm{d}x$ (in km) at a vertical view direction $\theta =0$° mapped to the surface. For later reference, we denote this spacing by $\Delta \xi _{\mathrm{or}}$, where the subscript “$\mathrm{or}$” stands for original grid. Note that for an oblique view, this grid requires some interpolation when $J_{k}(x_{k}(\xi ))$ is convolved with the window function to obtain $< J_{k}>$ in Equation 4.

We display the reconstructed height $h(\xi )$ along a surface profile through the simulated active region in Figure 5. The height was calculated for two allowed window-size ranges of $\xi _{\mathrm{width}} = 15 \text{ to } 20 \Delta \xi _{\mathrm{or}}$ (left diagram) and $\xi _{\mathrm{width}} = 30 \text{ to } 60 \Delta \xi _{\mathrm{or}}$ (right diagram). We started the correlation-height estimate with the maximum allowed $\xi _{\mathrm{width}}$ and then reduced $\xi _{\mathrm{width}}$ until the size for which the height could not be determined anymore because of a lack of structure inside of the window or until we reached the smallest $\xi _{\mathrm{width}}$ of the allowed range.

Figure 5

Results for the correlation method using a $\xi _{\mathrm{width}}$ range of 15 to 20 $\Delta \xi _{\mathrm{or}}$ (left) and a range of 30 to 60 $\Delta \xi _{\mathrm{or}}$ (right), where $\Delta \xi _{\mathrm{or}}=41.7\text{ km}$, for the view angles $\theta =-20$° and $\theta =+20$°. For more detail, see Section 4.1.

The accepted height $h(\xi )$ is then successfully obtained from the smallest window size inside the allowed range and shown by the red curves in Figure 5. Because of numerical errors, it is to be expected that at a few points $\xi$ the correlation procedure will return unrealistic, localized outliers. We try to eliminate or reduce these outliers by smoothing $h(\xi )$ using a running average of eight data points, which then yields the height estimate of the green curves.

To compare our reconstructed height $h(\xi )$ with the height from the simulation used to generate the test images (to which we will refer to as the true height), we indicate a range of heights associated with $R_{T}(\xi ,z)$ that denotes the height variation around the response-function maximum, displayed in grayscale. The lightest gray range represents $R_{T}/R_{T\mathrm{max}}\approx 0.5$, the intermediate range represents $R_{T}/R_{T\mathrm{max}}\approx 0.7$, and the darkest represents $R_{T}/R_{T\mathrm{max}}\approx 0.9$, where $R_{T\mathrm{max}}$ is the maximum value of the response function at a fixed ${\xi _{\mathrm{cent}}}$.

The mean error $\epsilon$ from Equation 11 obtained for the two ranges of $\xi _{\mathrm{width}}$ are $\epsilon = 1.44$ and $\epsilon = 2.26$ for $\xi _{\mathrm{width}} = 15 \text{ to } 20 \Delta \xi _{\mathrm{or}}$ and $\xi _{\mathrm{width}} = 30 \text{ to } 60 \Delta \xi _{\mathrm{or}}$, respectively, whereas the minimal cost functions achieved for them are 0.11 and 0.09. Although the agreement of our reconstructed height $h(\xi )$ with the true height is better for the shorter window-width ranges, a smaller cost function was obtained for the wider window-width range. The minimal cost function value, which measures the agreement between the virtual images $< I_{1}>$ and $< I_{2}>$, is the only quality number that we can obtain when using real data when a true height is not known, but it must be used with some care.

This comparison also shows another feature, which is also present in the later examples that we investigated. In the quiet Sun, the calculated height $h(\xi )$ seems to lie systematically below the height of the maximum response function. Since the response function in the quiet Sun is skewed with an upward tail, the vertical barycenter $h_{\mathrm{mean}}$ of the response function lies even higher than the height of the maximum response function. This discrepancy is even larger for more inclined view directions. Surprisingly, this discrepancy is smaller in the deepest parts of the active region (in the umbra). The relative error (Equation 10) is shown in all the diagrams in the lower panels. As $\mathrm{Err}$ measures the mean error relative to the width of the response function, any error less than $\mathrm{Err} = 1$ can be considered a good agreement. $\mathrm{Err} = 1$ is marked in the figures as a dotted line. We can see that the agreement is best in the deepest parts of the active region, and decreases toward the quiet Sun and especially in the penumbral region, where steep height gradients occur.

We also investigated the impact of the grid spacing $\Delta \xi$ at which $J_{1}$ and $J_{2}$ were discretized. We found that the standard spacing $\Delta \xi = \Delta \xi _{\mathrm{or}}$ returned the best result for the correlation method. A finer grid oversamples the data, increasing the computation time and increasing the number of outliers, whereas a coarser grid yields more stable results but at the cost of a reduced spatial resolution. The effect on the variation of $\Delta \xi$ is discussed further in the next sections, where the optimization method is tested.

4.2 Test of the Optimization Method

The optimization method is methodologically more consistent, but it involves a non-convex optimization problem. It is also non-linear in the sense that the matrices to be inverted at each iteration step depend on $J_{1}$ and $J_{2}$ and therefore on the estimate of $h(\xi )$ at each step.

To reduce the computational burden, we use the smoother result from the correlation method (green curves in Figure 6) as the initial $h(\xi )$ of the optimization iterations. To test the performance of this method, we performed our calculations for different grid spacings $\Delta \xi$ and also for different values of the regularization parameter $\mu$.

Figure 6

The optimization of the height vector for the views of $\theta =-20$° and $\theta =+20$° using a grid spacing of $\Delta \xi _{\mathrm{or}}=41.7$ km (left), $2\Delta \xi _{\mathrm{or}}=83.4$ km (right) and different values of $\mu$.

The errors $\epsilon$ with respect to the true height and the final cost function values are summarized in Table 1. The left section of the table shows the results obtained with the standard surface grid spacing $\Delta \xi _{\mathrm{or}}$, whereas the right section corresponds to a coarser grid spacing of $2\Delta \xi _{\mathrm{or}}$.

Table 1 Error and cost function values for the results of Figure 5a with $\xi _{\mathrm{width}}$ 15 to 20 $\Delta \xi _{\mathrm{or}}$ and its optimization using different $\mathrm{\mu}$-values for the regularization parameter for the standard grid spacing $\Delta \xi _{\mathrm{or}}$ and a spacing of $2\Delta \xi _{\mathrm{or}}$ for the view angles $\theta =-20$° and $\theta =+20$°. The row labeled “Correlation” refers to the initial $h(\xi )$ derived from the correlation method, and the results obtained after the optimization iterations with different regularization parameters are listed in the rows below. The error number is given for the full profile, for the quiet Sun, and for the spot; denoted by the subscripts $\mathrm{full}$, $\mathrm{qs}$, and $\mathrm{s}$, respectively. The lowest values of $\epsilon$ and cost function are underlined.

The optimization method reduces the cost function from its initial value obtained with the correlation method. This is expected, since the optimization is designed to achieve this. However, there is a discrepancy between the regularization parameter, which yields the smallest cost function, and the one that yields the smallest $\epsilon$. In both cases, the $\mu$ for which $\epsilon$ is smallest is larger than the one for which the cost function is.

Some of the results for the reconstructed height $h(\xi )$ after the optimization are shown in Figure 6. The figures show the results in the same manner as Figure 5 and along the same surface profile. As for Table 1, the left diagram shows the results obtained with the standard grid spacing, whereas the right one with the coarse grid spacing of $\Delta \xi =2\Delta \xi _{\mathrm{or}}$.

In both cases the curves correspond to the resulting height vectors using $\mu =10^{-6},10^{-5}$, and $10^{-4}$ as the regularization parameter, displayed in red, green, and blue, respectively. The results in Figure 6 show the effect of the regularization parameter on the results: larger $\mu$-values decrease the number of outliers but also can yield an overly smoothed height vector. Since very small $\mu$ values do not make significant changes to the reconstructed height $h(\xi )$ and larger values tend to smooth out many of the small-scale structures, we choose $\mu =10^{-5}$ as the standard parameter for the further tests, along with the standard grid spacing $\Delta \xi _{\mathrm{or}}$.

As was mentioned in Section 4.1, a smaller cost function does not necessarily mean a better agreement with the true height but rather a pair of virtual images that are very similar to each other. A height vector with many outliers can produce very similar virtual images, even if its agreement with the true height is not the best. A smaller $\mu$ produces more details in the height variation. However, in the quiet-Sun region, these details do not always agree with the fine-scale variations of the true height.

4.3 Test of the Effect of Viewing Angles

Once we have determined the set of parameters that yield the best results, with the smallest number of outliers but displaying enough detail (i.e. not over-smoothed), we use them as the standard set of parameters used to evaluate the effect of varying different viewing parameters on the results. In this section, we test the effect of changes in the view angles. The range for the window size $\xi _{\mathrm{width}}$ is set to 15 to 20 $\Delta \xi _{\mathrm{or}}$, the spacing of the surface grid is $\Delta \xi =\Delta \xi _{\mathrm{or}}$, and the regularization parameter is $\mu =10^{-5}$. While these parameters are fixed for all the tests in this section, we present the reconstructed height $h(\xi )$ for different sets of view geometry and study the effects that these have on the results.

The effect of increasing observing angles is displayed in Figure 7. The tests for the correlation and the optimization method were performed for a vertical view at $\theta =0$° and an inclination view of $\theta =-30$° and $\theta =+30$°, respectively. On the respective near-side slope of the sunspot penumbra, we find a large disagreement between our reconstructed height and the true height. The reason for this large error is due to the fact that the $\pm 30$° view direction comes close to the slope of the respective height profile. The slopes on both sides of the sunspot reach atan$\,|\mathrm{d}h/\mathrm{d}\xi | \le 0.5$, which corresponds to an angle as small as 65° from the vertical. Structures on the near-side slope appear therefore compressed in the image, and a correlation with the vertical view is prone to errors. The effect is the same in both viewing directions displayed in Figure 7, but it appears on opposite sides of the depression. Due to the large error in both cases, the mean error of the reconstructed height increased to $\epsilon =3.23$ for the view angles of $\theta =-30$° and 0° and to $\epsilon =2.24$ for the view angles of $\theta = 0$° and 30°.

Figure 7

Results for the reconstruction of the height vector for view angles of $\theta =0$° and −30° (left) and $\theta =+30$° and 0° (right), with a grid spacing of $\Delta \xi _{\mathrm{or}}=41.7$ km and a regularization parameter of $\mu =10^{-5}$.

In Figure 8, we display another example of an asymmetrical view, with the angles −10° and $+20$°, so that the total separation between the view points is again 30° as in Figure 7, but neither of the view directions is vertical. The reconstructed height after the optimization has a mean error of $\epsilon =1.37$, which is similar to the error of the curves of Figure 7. The general agreement of the reconstructed height with the true height is good, but again some of the smaller scale structures have been smoothed after the optimization.

Figure 8

The reconstruction of the height vector for the views of $\theta =-10$° and $+20$° with a grid spacing of $\Delta \xi _{\mathrm{or}}=41.7$ km and a regularization parameter of $\mu =10^{-5}$.

So far we have only considered cases with mutual separation angles of 30°. To further investigate the effect of the view geometry, we reconstruct the height for symmetrical view angles with different separations of $\theta =\pm 10$° and $\theta =\pm 20$°. The results are displayed in Figure 9 as the red and green curves, respectively. We present here only the optimization results for both view angle pairs.

Figure 9

Comparison of the reconstructed height vector for $\theta =\pm 10$° (red) and $\pm 20$° (green) with a grid spacing of $\Delta \xi _{\mathrm{or}}=41.7$ km and a regularization parameter of $\mu =10^{-5}$.

Figure 9 shows that the reconstructed height profile agrees better for symmetric observations. The mean error for $\theta =\pm 10$° is reduced to $\epsilon =1.38$, and for $\theta =\pm 20$° the error is $\epsilon =1.43$. This overall error is better for symmetrical views than for asymmetrical view angles, and the error is also smaller for smaller separation angles. However, for the pixel size of the simulated images of $\mathrm{d}x=\Delta \xi _{\mathrm{or}}=41.67$ km and a viewing angle $\theta$, the depth resolution is limited to about $\Delta h= \mathrm{d}x/\sin \theta =122$ km for $\theta =20$° and 240 km for $\theta =10$°. The reconstructed height for $\theta =\pm 10$° has generally a good agreement with the true height, but the height curve is very smooth, and some of the smaller scale details have been lost due to the lack of resolution and probably over-regularization with our standard value for $\mu$. A large $\mu$ during the optimization may over-smooth the reconstructed height.

The reconstructed height profiles for different sets of view angles give us some information on the effect of the view geometry on the results: the reconstructed height vectors have a smaller $\epsilon$ for symmetrical views and for smaller separation angles down to about $\pm 10$°. In addition, errors might arise on the penumbra of the active regions if one of the view angles is too large so that the view direction comes too close to being parallel to the near-side slope. If the total separation of the view directions is small enough, not exceeding a separation of 40°, then the results from the correlation method seem to be reliable.

4.4 Analysis of a Second Profile of the Active Region

In this section, we analyze another intensity profile, from the same simulated active region but for a different epipolar line profile with different $y$-coordinate in the images. This second profile exhibits more structure and is therefore more complicated. We will show the reconstructed height for the standard view geometry of $\theta =\pm 20$° and for the standard and coarser grid spacing $\Delta \xi$. We also test the results of observing with asymmetrical views of $\theta = 0$° and $\pm 30$° using the standard set of parameters from Section 4.2.

Analogously to the results of Section 4.2, the reconstructed height $h(\xi )$ is better when using the standard grid spacing $\Delta \xi =\Delta \xi _{\mathrm{or}}$ than when using $\Delta \xi =2\Delta \xi _{\mathrm{or}}$. The height curves displayed in Figure 10 show the results of the correlation method (red), the optimization method with the standard (green), and coarse (blue) grid spacing, both with a regularization parameter $\mu =10^{-5}$. The mean error for the correlation results is $\epsilon =2.18$, and those for the optimization are $\epsilon =2.15$ for $\Delta \xi =\Delta \xi _{\mathrm{or}}$ and $\epsilon =2.27$ for $\Delta \xi =2\Delta \xi _{\mathrm{or}}$ (compare with Table 1).

Figure 10

Correlation and optimization results for the reconstruction of the height vector for the views of $\theta =\pm 20$° with different grid spacing of $\Delta \xi _{\mathrm{or}}=41.7$ km and $2\Delta \xi _{\mathrm{or}}=83.4$ km and a regularization parameter of $\mu =10^{-5}$.

On the left penumbral slope of this profile of the active region, there is a structure where the height profile includes a local sharp peak. This structure is not well reproduced in the analysis from $\theta =\pm 20$°, even after optimization of the height vector. However, it can be reproduced with a different viewing geometry as demonstrated in Figure 11a. The curves in Figure 11 correspond to the reconstructed height for a vertical view of $\theta =0$°, $\theta =-30$° (left), and $\theta =+30$° (right). Just as in Figures 7a and 7b, there is a disagreement with the true height at each near-side slope of the sunspot. Nevertheless, in the left panel of Figure 11, the local peak inside the active region is well reproduced because it lies on the far-side slope of the sunspot. However, on the right, near-side slope the reconstruction again fails due to the effect discussed in the previous section, where the view direction comes too close to the slope tangent. Especially after the optimization procedure, shown with the green curve, this structure on the left slope is very well reproduced, and the disagreement with the true height decreases considerably after the optimization in both diagrams of Figure 11.

Figure 11

Results for the reconstruction of the height vector for the views of $\theta =-30$° and 0° (left) and $\theta =30$° and 0° (right) with a grid spacing of $\Delta \xi _{\mathrm{or}}=41.7$ km and a regularization parameter of $\mu =10^{-5}$.

For the simpler standard profile presented in Sections 4.1 and 4.2, the correlation results were generally good, and the optimization method did not change the agreement with the true height. Because of the complex structure inside of the second profile presented in Figure 11, the correlation procedure does not yield such good results. The optimization provides an appreciable improvement here. We therefore suggest that, depending on the level of structure within the active regions, we can choose to use the correlation or the optimization method to achieve the best results. It is also confirmed that the results depend strongly on the view-angle geometry and become erroneous when the angle between one of the view directions and the local surface becomes too small.

4.5 Testing the Effect of Noise in the Reconstructed Height

To study the performance of our method with noisy images, we ran our calculations while increasingly adding noise to the continuum intensity images. We here again use the standard parameters for the calculation and view angles $\theta =\pm 20$°. We determined $\epsilon$ and the cost function of the reconstructed height vector for noise levels that were successively increased by a factor of ten. For each level of the image signal-to-noise ratio (SNR), we determined the mean and standard deviation of both $\epsilon$ and the cost function.

Figure 12 shows the mean $\epsilon$ and cost function as functions of the SNR, where for each point the error bars represent the respective standard deviation. As expected, the mean error and the cost function decrease with the noise level. Above about SNR $\approx 10^{2}$, the decrease is only marginal, and the error and cost functions saturate to the level intrinsic to our method. Conversely, noise starts to significantly influence our height estimates only for an SNR smaller about $10^{2}$.

Figure 12

Curves of $\epsilon$ (left) and the cost function $f_{\mathrm{cost}}$ (right) as functions of the signal-to-noise ratio.

In Figure 13, we show the reconstructed height vector for some of the noise levels analyzed. The red curve represents the reconstructed height when the profiles have no added noise. For the green curve, the profile had an SNR of 500 and for the blue one, an SNR of 50. These two values were selected because SNR = 500 can be expected for the real observations (see Solanki et al., 2020). For SNR = 50, we still obtain reliable results on average, but in regions of steep slopes our method may produce locally larger errors. The number of these local failures increases with decreasing SNR.

Figure 13

Estimated height vector for the intensity profiles of $\theta =-20$° and $\theta =+20$° with different levels of added noise.

4.6 Full 2D Height Maps

So far we displayed the results for individual profiles through the active region to better visualize the performance of our method. By stacking the reconstructed height profiles in the $y$-direction we obtain 2D maps of the reconstructed height. Each profile is calculated individually but with the same parameters: the view angles are $\theta =\pm 20$°, the surface grid size is $\Delta \xi =\Delta \xi _{\mathrm{or}}$, the correlation window has a range of $\xi _{\mathrm{width}}= 15\text{ to }20\Delta \xi$, and a regularization parameter of $\mu =10^{-5}$ is used for the optimization method. As in the previous sections, we obtained a reconstructed height $h(\xi )$ from the correlation method and then smoothed and used it as the initial vector for the optimization. Each horizontal line of the 2D height maps corresponds to the reconstructed height of the respective profile after optimization. No smoothing was applied in the $y$-direction.

To understand the effect of different resolutions of the images on the results, we successively degraded the horizontal ($x$-direction) resolution of the images by a factor of two each time to produce four different image sets with image pixel sizes $\mathrm{d}x=\Delta \xi _{\mathrm{or}}=41.7$ km and its multiples, i.e. 83.4 km, 166.8 km, and 291.9 km. The pixel size in $y$ was maintained. These image sets with the same pixel resolution for both views were used to reconstruct the 2D height maps shown in Figures 14a, b, c, and d.

Figure 14

2D height maps reconstructed from image pairs with view directions $\theta =-20$° and $\theta =+20$° with image pixel sizes of 41.7 km, 83.4 km, 166.8 km, and 291.9 km in the $x$-direction in Panels a, b, c, and d, respectively. For the $y$-direction, the pixel size is always 41.7 km.

The height maps on Figure 14 show the trade-off between resolution and quality. They have to be compared with the map of the barycenter height of the response function $R_{T}(\xi ,z)$ for a vertical view direction of $\theta =0$°, computed directly from the MHD simulation presented in Figure 15. The maps reconstructed at higher spatial resolution contain more visible errors, which mainly appear as bright red or dark blue regions, whereas the lower-resolution maps have fewer errors but smaller scale structures, especially inside the pore, are not reproduced. For real data from SO/PHI and another spacecraft observing from about 1 AU distance, we expect images with the low resolution of Figure 14d at best. It must be mentioned that in Figure 14 the resolution is decreased only in the horizontal direction of the images. In the vertical direction, all the images are composed of 340 reconstructed height vectors with a resolution in the vertical direction equal to $\Delta \xi _{\mathrm{or}}=41.7$ km.

Figure 15

Barycenter height map obtained from the averaged heights associated with the response function $R_{T}(\xi ,z)$ in a range within $R/R_{\mathrm{max}}\approx 0.5$ of the full resolution 2D synthetic MHD data.

4.7 Estimating a Height Vector for Two Views with Different Resolution

Lastly, we evaluate the results obtained when the image resolution is different for the two views. Since real photospheric observations for stereoscopy must be conducted by combining data obtained by SO/PHI with, e.g., SDO/HMI, it is to be expected that the two images to estimate the height with our stereoscopic method will have different resolutions.

We display the respective result for the reconstructed height along one epipolar profile in Figure 16. We used a view direction of $\theta =-20$° for the higher-resolution image with a pixel resolution of 166.7 km and $\theta =+20$° for the lower-resolution image with a pixel resolution of 291.9 km. We reconstruct the height for these images either on a surface grid spacing $\Delta \xi =4\Delta \xi _{\mathrm{or}}=166.7$ km (left) or $\Delta \xi =7\Delta \xi _{\mathrm{or}}=291.9$ km (right). Just as in the previous section, we used the standard set of correlation and optimization parameters.

Figure 16

The optimization method of the height vector from images with different pixel resolution. The views of $\theta =-20$° and $\theta =+20$° have pixel sizes of 166.8 km and 291.9 km, respectively. The surface grid spacings are $\Delta \xi =166.7\text{ km}$ (left) and $\Delta \xi =291.9$ km (right).

When estimating the height vector for intensity profiles with different resolutions, we choose the surface grid to have the same spacing of either one of the intensity profiles. The left panel in Figure 16 shows the reconstructed height obtained by mapping both images to the higher-resolution grid and therefore increasing the size of the low-resolution image, which gave an error number of $\epsilon =4.59$. The right panel shows the reconstructed height when mapping both images to the lower-resolution grid, therefore decreasing the size and resolution of the higher-resolution image profile. This yields a slightly increased error number of $\epsilon =4.71$. The obtained results show that the quality of the reconstructed height is limited by the resolution of the lower-resolution image. Due to the larger pixel sizes used here, both reconstructed height vectors have higher errors than in the previous sections. So the resolution of the data limits the quality that can be achieved, and in addition the lower-resolution image determines the quality or resolution that can be achieved in the reconstructed height vectors. Note that the error roughly tripled compared to the standard case in Figure 6 (left panel), whereas the pixel sizes used here were larger by a factor of four and seven. In comparison with earlier results, the results reproduce the depth of the sunspot relatively well, whereas the feature is estimated to be much larger than it is, and hence the depression of the solar surface is estimated to cover a larger area.

Lastly, we also include 2D height maps for the results from different-resolution images in Figure 17. The geometry and resolution are the same as for Figure 16, but now the heights obtained over the entire image are presented. The resolution of both images in this case is very low, limiting the amount of details that can be reconstructed, but at the same time giving results with only a few outliers. This is another example of the quality/resolution trade-off that is always present in the results of this method.

Figure 17

2D height obtained for the combination $\theta = -20$° with a pixel size of 166.8 km and $\theta = +20$° with 291.9 km. The grid spacing in (a) is of $\Delta \xi =4\Delta \xi _{\mathrm{or}}=166.7$ km, whereas for (b), it is $\Delta \xi =7\Delta \xi _{\mathrm{or}}=291.7$ km.

5 Summary and Discussion

We have developed a stereoscopic method that allows estimating height variations of the solar surface using observations from two different view points. With synthetic test data (obtained from an MHD simulation), this method is capable of estimating with good precision the Wilson depression of a sunspot or large pore and the height structure within such features. The simulated Wilson depression of the active region in our test data was of around 600 km, which agreed with our reconstructed height profiles.

The stereoscopy method yields robust results for view points with a separation ranging approximately from 10° to 40°. For smaller separation angles, the geometrical resolution error of the height vector increases as $1/\sin (|\theta _{1}-\theta _{2}|/2)$ for angles $\theta _{k}$ of view direction $k$ with respect to the surface normal. For larger separation angles, we eventually encounter problems if one of the angles $\pi /2-\theta _{k}$ approaches the angle of the sloping solar surface at the edges of the sunspot.

Our tests with synthetic data also showed stable results when the images were contaminated by a limited amount of noise. The Wilson depression is well estimated for all intensity profiles with a signal-to-noise ratio above about 50. If the noise increases beyond this level, then the smaller-scale structures are reproduced with less precision; only larger-scale structures are still reproduced.

An advantage for our method is that it requires only two surface images to perform the calculations and yields height-variation estimates, as opposed to previous methods, which had to make a series of assumptions to estimate the Wilson depression indirectly (see, e.g., Löptien et al., 2020). Here only the intensity from the two profiles taken from different vantage points in the same spectral window is necessary to perform our calculations.

Our method is limited by the spatial resolution of the data. The higher the resolution, the smaller the scales of the structures that can be reproduced. For the resolution of the images employed in this article, we can expect to reproduce the Wilson depression with good accuracy. For reproducing the height variations of smaller-scale structures within a sunspot and/or in the quiet-Sun region, we require images with higher resolution.

The finite thickness of the photospheric opacity layer sets a limit on the resolution we can achieve. We assume here that this layer is infinitely thin, so that it can be characterized by a single height $h(\xi )$ rather than by a height distribution. Moreover, for real observations, the barycenter of the height distribution varies with the view angle, so that the height $h(\xi )$ we solve for does not compare straight-forwardly with the true distribution of the formation height of the observed radiation from a vertical viewpoint. The larger the angle between view direction and surface normal, the more the formation of the observed radiation is shifted to higher altitudes.

We find in some cases (see Table 1, Figure 7, left) that the agreement between the calculated height and the true height distribution is not always improved by the optimization method as compared to the correlation results. This is because we optimize over the cost function and not over the error number $\epsilon$. Even though we observe the same structure, the different viewing points observe different smaller scale substructures due to changes in photospheric opacity (Carlsson et al., 2004). The optimization procedure copes with these differences in the two intensity profiles with a resulting height vector that yields a more similar mapping from both views, although this does not necessarily produce a better agreement with the true height distribution.

Another effect associated with the photospheric opacity can be observed. A small offset sometimes occurs between the calculated and reference heights. This offset applies to the entire epipolar profile, but it is not the same for different viewing conditions. The offset increases as the viewing angles grow. This may be caused by differences in the formation height (Carlsson et al., 2004) of the observed radiation. The quiet-Sun height ($h=0$) at $\tau =1$ varies with viewing angle because the distance to the solar surface grows and hence the height for which $\tau =1$ does too. We estimate here height variations around the quiet-Sun level and compare them to the reference height from a vertical view. It can be seen throughout the results that regardless of the offset, the quiet-Sun height estimates fall systematically below the level of the reference height, whereas the depth of the umbra displays a better agreement. Therefore our estimates of the Wilson depression, i.e. the height difference between the quiet Sun and the umbra, is a few tens of kilometers smaller than in the reference data. We believe that this reduction of the Wilson depression is a combined effect of the finite width of the applied correlation functions, which is averaging the depth profile within the correlation window, the flattening of the $\tau =1$ layer when the photosphere is observed from an inclined viewpoint, and the asymmetry of that layer when observed from different directions, which might increase the disparity of quiet-Sun structures.

6 Outlook

The synthetic images based on an MHD model used in this article were calculated with virtual spacecraft positions in a simplified epipolar plane geometry with a planar solar surface and epipolar lines parallel to the image $x$-axes. We aim, however, to apply our method to combinations of observations made by PHI and, for example, HMI, which will have more complicated mutual positions and distances form the Sun. As mentioned before, some preprocessing, e.g. rectification of the images, will be necessary to cast the images from arbitrary observer positions into a shape comparable to the synthetic test data used here. Adapting our method so that it can be applied to actual observed images will be the subject to a future article.

An improvement for this method could be to include a dependency of the surface radiance with the angle between view direction and local surface normal. This extension of our model can perhaps help correct the offset that we see in the results presented here.

An extension and application of our method is to stereoscopically determine and compare the height at different wavelengths. This is particularly interesting when considering the core of a spectral line in addition to the continuum. This will allow comparing the line-formation height in different solar features and may allow constraining the temperature gradient.