# Reconstructing the Normal and Shape at Specularities in Endoscopy

###### Abstract

Specularities are numerous in endoscopic images. They occur as many white
small elliptic spots, which are generally ruled out as nuisance in image
analysis and computer vision methods. Instead, we propose to use specularities
as cues for 3D perception. Specifically, we propose a new method to
reconstruct, at each specularity, the observed tissue’s normal direction
(i.e., its orientation) and shape (i.e., its curvature) from a single image.
We show results on simulated and real interventional images.

Index Terms—  Specular, shape operator, endoscopy

##  1 Introduction

Specularities on smooth reflecting surfaces depend on the surface’s local
orientation and curvature, on the light source and on the camera position [1,
Chapter 2]. In the general setting, all these contain far too many unknowns to
allow one to exploit specularities effectively in 3D reconstruction. However,
endoscopy is a particular case where the light source and camera form a
collocated rig. Specularities are numerous in endoscopy, where the moist
tissue causes specular reflection. However, they are generally considered as
artifacts in endoscopic 3D reconstruction methods, both in classical methods
such as Shape-from-Shading (SfS) [2] and learning-based methods such as direct
depth estimation [3]. More recently, it has been understood that specularities
could form useful cues for 3D reconstruction. The study of [4] has shown that
training a learning-based reconstruction model from synthetic data which can
eventually deal with real data requires one to generate specularities.
Concretely, we have recently shown that the surface normal could be
reconstructed at each specularity [5, 6], with a two-fold ambiguity. The
reconstruction method uses the specular isophote, which is the curve enclosing
the specularity in the image. Based on the assumption that the surface
underlying the specularity is flat, we have used a mathematical model to show
that the isophote is elliptic. This approach has advantages: the method runs
fast and can process full high image resolution, in contrast to the learning-
based approach which is currently limited to lower image resolution to keep
runtime reasonable (the available datasets such as the EndoSLAM [7] and UCL
[8] ones are at lower resolution). However, it also has strong limitations.
First, the reconstruction is up to a two-fold ambiguity: only one of the two
directions found at each point is the correct normal. Second, the local
planarity assumption is unreasonable, as specularities typically appear near
highly curved surface patches. Is an ambiguous surface normal all one can
reconstruct from a specularity?

We propose to move a step ahead by using a more advanced local surface model.
Specifically, we leave the restrictive flat surface to adopt a curved surface
model. We show that with this model, the specular isophote is still elliptical
in the image, albeit that its eccentricity is not just due to surface
orientation as in past work, but to both surface orientation and curvature. We
show that one can reconstruct the surface normal from just the specularity’s
Brightest Point (BP) and the surface curvature from the elliptical isophote.
Specifically, we show that the principal curvature directions (the directions
of minimum and maximum curvatures on the tangent plane) and the curvature
ratio can be reconstructed. Therefore, the specular isophote contains much
more information than one could think from past work. We propose an automatic
reconstruction method, which, from an endoscopic image and the camera’s
intrinsic parameters, detects the specularities and reconstructs the local
surface orientation and curvature ratio at each specularity. This
reconstruction can be used to implement downstream tasks such as computer-
aided navigation, to bootstrap other reconstruction methods such as SfS and to
train depth estimation networks. We evaluate the method on realistic simulated
colonoscopic data from Blender and surgical laparoscopy images.

##  2 Theory and methods

###  2.1 Mathematical Model

We consider a calibrated perspective camera model with intrinsic parameters in
matrix $\mathbf{K}\in\mathbb{R}^{3\times 3}$ and a point light source near the
camera centre located at the origin. We consider a smooth specular surface and
assume that the brightness reﬂected along the sightline direction $V$ is a
function of $V\cdot P$ , where $P$ is the specular direction, also known as
the direction of a perfect specular reflection. This assumption is commonly
used in reflection models such as Blinn-Phong [9] and the physics-based BSDF
model popularised by Disney [10]. In these models, each specularity has a BP,
approximately located at its centre, where the reflected light energy is the
highest. Let $L$ be the lighting direction; it is known that, from the fact
that the vectors $L$, $V$ and $P$ have unit length [1, 9, 10], the half-angle
direction $H=(L+V)/2$ is collinear with the surface normal $\mathbf{N}$ at the
specularity’s BP. In endoscopy, as we have $L\approx V$, the surface normal is
thus collinear with the sightline $V$.

We model the observed tissue as a collection of 3D surface patches. We model
each patch by three components, as shown in figure 1 (left): 1) its normal
direction $\mathbf{N}$, 2) its tangent plane, with basis vectors
$\mathbf{u}_{1},\mathbf{u}_{2}$ chosen along the patch’s principal directions
(i.e., the directions of minimal and maximal curvature), and 3) its curvature
ratio, which is the ratio of the minimal to the maximal curvature. We
instantiate one patch per specularity. We assume the specularity to conform to
the elliptic model, which is widely applicable in endoscopic images, where
they occur as small isolated white blobs [1, 5].

Fig. 1:  Left: local coordinate system at the BP oriented using the normal and
principal directions. The light and camera are collocated. Elliptical specular
isophotes are shown in cyan. Right: an example of three neighbouring surface
patches from a real endoscopic image, showing the difficulty to determine the
BP directly from image intensities.

###  2.2 Reconstruction Method

Our reconstruction method has three main steps given as Algorithm 1. The first
step, which is as in [6], uses a learning-based specularity detector and
ellipse fitting to retrieve an isophote for each specularity. We perform a
simple ellipticity test based on the fitting residual. The second and third
steps are geometry-based reconstruction given directly below.

####  2.2.1 Sightline-based Normal Reconstruction

We exploit the property that, at the BP, the surface normal is collinear with
the sightline. We approximate the BP by the elliptic isophote’s centre
$(x_{0},y_{0})$ in image pixel coordinates, as shown in figure 1 (right). The
sightline is trivially found using the camera model, giving the surface normal
direction $\mathbf{N}\in S^{2}$ as
$\mathbf{N}\propto\mathbf{K}^{-1}(x_{0},y_{0},1)^{\top}$, where $S^{2}$ is the
set of direction vectors in 3-space. Obviously, the reconstruction accuracy
directly depends on the accuracy of the BP estimate. However, determining the
exact localtion of the BP is challenging because of camera saturation at
specularities. As seen in the right side of figure 1, all the red pixels have
the maximal intensity for the selected surface patch. In contrast, it is much
easier to detect an isophote whose centre and the unknown BP are nearly
coincident.

####  2.2.2 Ellipse-based Shape Reconstruction

We represent a specular isophote ellipse by its standard $3\times 3$ real
symmetric matrix form $C$.
This satisfies $(x,y,1)C(x,y,1)^{\top}=0$ for the image pixel coordinates
$(x,y)$ of an isophote point. Next, we map the conic $C$ to the normalized
image plane111The origin is at the principal point and the distances are
measured in units of focal length. as
$C^{\prime}=\mathbf{K}^{\top}C\mathbf{K}$. Note that $C^{\prime}$, in its
turn, is a $3\times 3$ real symmetric matrix. It has therefore the following
properties: 1) always diagonalizable; 2) it has real eigenvalues; and 3) its
eigenvectors are orthogonal. We then diagonalise $C^{\prime}$ as
$C^{\prime}=\mathbf{R}^{\top}diag(\lambda_{1},\lambda_{2},\lambda_{3})\mathbf{R}$,
where $\\{\lambda_{1},\lambda_{2},\lambda_{3}\\}$ are the eigenvalues sorted
in ascending order. The corresponding eigenvector set
$\\{\vec{v_{1}},\vec{v_{2}},\vec{v_{3}}\\}$ forms an orthonormal basis of
$\mathbb{R}^{3}$. In fact, $\mathbf{R}=(\vec{v_{1}},\vec{v_{2}},\vec{v_{3}})$
is a 3D rotation of the standard camera coordinate system to the so-called
eigenvector-frame. It aligns the principal directions on the tangent plane at
the BP with the ellipse’s major and minor axes, with the normal direction
$\mathbf{N}$ being aligned with the cross product
$\vec{v_{3}}=\vec{v_{1}}\times\vec{v_{2}}$. The principal directions on the
tangent plane are therefore given by
$\mathbf{u}_{1},\mathbf{u}_{2}=\vec{v_{1}},\vec{v_{2}}$. The ratio of
principal curvatures $|\lambda_{1}/\lambda_{2}|$ can be derived from the
ellipse’s eccentricity defined as
$\epsilon=\sqrt{1-(\lambda_{1}/\lambda_{2})^{2}}$.

Algorithm 1 Single-image specular normals and shapes reconstruction

Image $I$, camera intrinsics $\mathbf{K}$

Sets of Brightest Points (BP), normals $\mathbf{N}$, principal directions
$\mathbf{u}_{1},\mathbf{u}_{2}$, curvature ratios $\epsilon$

\- Use the neural network to segment the specular mask $M$ from $I$
$\triangleright$ [6]

for each connected component $c\in M$ do

\- find the specular isophote as the outer boundary of $c$

\- smooth the isophote using a cubic B-spline

\- fit an ellipse to the isophote, with matrix representation $C$ and fitting
residual $r$

if $r\leq t$ then $\triangleright$ ellipticity test with threshold $t$

\- find the BP $(x_{0},y_{0})$ as the ellipse centre

\- compute the normal direction
$\mathbf{N}=\mathbf{K}^{-1}(x_{0},y_{0},1)^{\top}$

\- normalise it as $\mathbf{N}=\mathbf{N}/\|\mathbf{N}\|$

\- transfer the ellipse to the normalised image plane as
$C^{\prime}\leftarrow\mathbf{K}^{\top}C\mathbf{K}$

\- retrieve the ellipse axes by eigen-decomposition as
$C^{\prime}=\mathbf{R}^{\top}diag(\lambda_{1},\lambda_{2},\lambda_{3})\mathbf{R}$,
with eigenvalues $\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$ and eigenvectors
$\mathbf{R}=(\vec{v_{1}},\vec{v_{2}},\vec{v_{3}})$

\- determine the principal directions on the tangent plane as
$\mathbf{u}_{1},\mathbf{u}_{2}=\vec{v_{1}},\vec{v_{2}}$

\- compute the curvature ratio from the ellipse eccentricity $\epsilon$

end if

end for

##  3 Experimental results

###  3.1 Synthetic Colonoscopic Data

####  3.1.1 Model

We use the modelling and rendering software suite Blender to synthesise images
using a wide-angle perspective camera of 10 mm focal length. We use a
spotlight located at the camera centre with a radius of 0.1 mm and projected
to a 120∘ cone. We generate a colon-like deforming surface with free-form
deformation [11]. The generated surface is represented by a triangular mesh
with 502,502 vertices and 1,002,000 faces. We deform it over 250 frames. The
surface is imparted with a strong ‘clearcoat’ to boost specular reflections.
We use the anisotropic principled BSDF model [10], mimicking real-world
shading and specularities.

####  3.1.2 Normal Reconstruction

Our method found 1,500 elliptic specularities. The reconstructed normals using
the proposed sightline based method are compared with the true mesh normals,
at mesh vertices corresponding to the set of detected BPs in the image. The
histogram of errors in normal estimates is shown in figure 4.

####  3.1.3 Shape Reconstruction

We estimate the true principal curvatures $\kappa_{1}$ and $\kappa_{2}$ as the
eigenvalues of the second fundamental form that we compute numerically by the
method [12], from which we obtain the true curvature ratio. We also compute
the true principal directions, as the eigenvectors of the second fundamental
form. A representative example is shown in figure 2 (top right), showing that
the specular ellipse’s axes match the principal directions. We perform a
sanity check by verifying that the reconstructed normal is orthogonal to the
reconstructed principal directions on the tangent plane, confirming that they
together form an orthonormal local coordinate system. The histogram of errors
in principal direction estimates is shown in figure 4.

Fig. 2:  Colon-like synthetic sequence. The available 2D-3D correspondences
(in magenta) are used to compare between our reconstruction of the normal, the
principal curvature ratio and the principal directions for the surface at
specularity’s brightest point and the same quantities estimated from the 3D
model via the method of Rusinkiewicz [12].

Fig. 3: Error histogram: angle between reconstructed and true normals (in
degrees).

Fig. 4: angle between reconstructed normals and $\vec{v_{1}}\times\vec{v_{2}}$
(in degrees).

###  3.2 Real Data

The applicability of our reconstruction approach to real data is evaluated in
2D laparoscopic views of the liver, where a sparse set of specular highlights
is detected and for which we have 3D reference models. These references models
are obtained by 3D-to-2D registration of a virtual patient model segmented
from a preoperative CT scan by experts. The registration is solved by the
method from [13] to a satisfying accuracy. We use three images extracted from
three procedures collected in our hospital under ethical approval
IRB00008526-2019-CE58 issued by CPP Sud-Est VI in Clermont-Ferrand, France. A
representative example is shown in figure 5.

Fig. 5:  Qualitative results for the liver data.

####  3.2.1 Normal Reconstruction

Contrary to [5] in which a small set of isolated specularities is manually
selected to then detect isophotes using the marching squares, we automatically
detect numerous specular isophotes from our neural network [6]. An example
showing elliptical isophotes is shown in figure 5. 193 small elliptical
specularities are kept from the three selected images according to the
ellipticity criterion and used to compare the reconstructed normals with two
other estimates: 1) the normals produced by the pose from circle method [5],
and 2) the normals directly determined from the registered virtual 3D model.
The histograms of figure 6 show a good agreement between the three methods.
The green histogram shows an angular error centred at 2.5∘, which confirms the
efficiency of both the proposed and the registration methods [14]. The
improvement is partly due to the ellipticity, which is not present in [6].
Nevertheless, the magenta histogram shows an angular error centred at 3∘,
which shows a fair agreement between the proposed and the previous methods.

Fig. 6:  Evaluation of normal reconstruction for the laparoscopic liver
images. The histograms show the absolute angular difference (in degrees).

####  3.2.2 Principal directions principal curvature ratio

First, in order to evaluate the quality of our reconstructed principal
directions $\\{\vec{v_{1}},\vec{v_{2}}\\}$, we have compared them to the
principal directions $\\{\mathbf{u}_{1},\mathbf{u}_{2}\\}$ estimated using
[12] from the triangle mesh. An easy way to do that is to quantify the angle
between their cross products since both of them are supposed to be aligned
with the normal direction. Let
$\theta_{e}=\angle(\mathbf{u}_{1}\times\mathbf{u}_{2},\vec{v_{1}}\times\vec{v_{2}})$
be this angle expressed in degrees. We then consider its absolute value
$\left|\theta_{e}\right|$ to quantify the degree of agreement between the two
methods, leading to the histogram error illustrated in figure 7(a). This
histogram shows a small error which is always less than 3∘ for the same
specularities used in Section 3.2.1 for normal reconstruction evaluation.

Second, we evaluated the differences between the two estimated curvature
ratios (proportional to the eccentricity) as
$d=\left|\left|\lambda_{1}/\lambda_{2}\right|-\left|\kappa_{1}/\kappa_{2}\right|\right|$
where $\kappa_{1}$ and $\kappa_{2}$ are the minimum and maximum curvatures
estimated using [12]. The histogram of figure 7(b) shows that the error was
almost comprised between 0 and 0.2 (and mostly close to zero) which is still
very acceptable since the eccentricity of an ellipse always lies between 0 (a
perfect circle) and 1 (a straight line). The example presented in figure 5
(top right) show that the reconstructed normals and principal direction are of
good quality.

(a) principal directions (differences $\left|\theta_{e}\right|$ in degrees)

(b) principal curvature ratio (differences $d$ in a.u.)

Fig. 7: Evaluation of shape reconstruction: direct comparison between our
method and [12].

We finally show qualitative results for a colonoscopic image extracted from
the EndoMapper database [15] and for which high quality 3D models are not
available.

Fig. 8:  Colonoscopy. From left to right: image $I$, neural mask $M$, the
third and fourth images show some reconstructed principal directions (red) and
normals (blue).

##  4 Conclusion

We have presented a mathematical model based on a second-order model of the
tissue surface at specularities. This model leads to a reconstruction method
which can infer the local surface normal, principal directions and curvature
ratio. The method is simple and fast. It produces normal estimates with
accuracy on par with previous methods, solves the local two-fold ambiguity and
estimates an additional level of local shape information. As it works from a
single image, we plan to use it to boost the NRSfM method in the endoscopic
case.

## References

  * [1] David A Forsyth and Jean Ponce,  Computer vision: a modern approach,  prentice hall professional technical reference, 2002.
  * [2] Gastone Ciuti, Marco Visentini-Scarzanella, Alessio Dore, Arianna Menciassi, Paolo Dario, and Guang-Zhong Yang,  “Intra-operative monocular 3d reconstruction for image-guided navigation in active locomotion capsule endoscopy,”  in 2012 4th IEEE RAS & EMBS International Conference on Biomedical Robotics and Biomechatronics (BioRob). IEEE, 2012, pp. 768–774.
  * [3] Yuying Liu and Siyang Zuo,  “Self-supervised monocular depth estimation for gastrointestinal endoscopy,”  Computer Methods and Programs in Biomedicine, p. 107619, 2023.
  * [4] Rema Daher, O León Barbed, Ana C Murillo, Francisco Vasconcelos, and Danail Stoyanov,  “Cyclesttn: A learning-based temporal model for specular augmentation in endoscopy,”  in International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, 2023, pp. 570–580.
  * [5] Karim Makki and Adrien Bartoli,  “Normal reconstruction from specularity in the endoscopic setting,”  in 2023 IEEE 20th International Symposium on Biomedical Imaging (ISBI). IEEE, 2023, pp. 1–5.
  * [6] Karim Makki, Kilian Chandelon, and Adrien Bartoli,  “Elliptical specularity detection in endoscopy with application to normal reconstruction,”  International Journal of Computer Assisted Radiology and Surgery, pp. 1–6, 2023.
  * [7] Kutsev Bengisu Ozyoruk, Guliz Irem Gokceler, Taylor L Bobrow, Gulfize Coskun, Kagan Incetan, Yasin Almalioglu, Faisal Mahmood, Eva Curto, Luis Perdigoto, Marina Oliveira, et al.,  “Endoslam dataset and an unsupervised monocular visual odometry and depth estimation approach for endoscopic videos,”  Medical image analysis, vol. 71, pp. 102058, 2021.
  * [8] Anita Rau, PJ Eddie Edwards, Omer F Ahmad, Paul Riordan, Mirek Janatka, Laurence B Lovat, and Danail Stoyanov,  “Implicit domain adaptation with conditional generative adversarial networks for depth prediction in endoscopy,”  International journal of computer assisted radiology and surgery, vol. 14, pp. 1167–1176, 2019.
  * [9] James F Blinn,  “Models of light reflection for computer synthesized pictures,”  in Proceedings of the 4th annual conference on Computer graphics and interactive techniques, 1977, pp. 192–198.
  * [10] Brent Burley and Walt Disney Animation Studios,  “Physically-based shading at disney,”  in Acm Siggraph. vol. 2012, 2012, vol. 2012, pp. 1–7.
  * [11] Pushkar Joshi, Mark Meyer, Tony DeRose, Brian Green, and Tom Sanocki,  “Harmonic coordinates for character articulation,”  ACM Transactions on Graphics (TOG), vol. 26, no. 3, pp. 71–es, 2007.
  * [12] Szymon Rusinkiewicz,  “Estimating curvatures and their derivatives on triangle meshes,”  in Proceedings. 2nd International Symposium on 3D Data Processing, Visualization and Transmission, 2004. 3DPVT 2004. IEEE, 2004, pp. 486–493.
  * [13] Yamid Espinel, Lilian Calvet, Karim Botros, Emmanuel Buc, Christophe Tilmant, and Adrien Bartoli,  “Using multiple images and contours for deformable 3d–2d registration of a preoperative ct in laparoscopic liver surgery,”  International Journal of Computer Assisted Radiology and Surgery, vol. 17, no. 12, pp. 2211–2219, 2022.
  * [14] Bongjin Koo, Erol Özgür, Bertrand Le Roy, Emmanuel Buc, and Adrien Bartoli,  “Deformable registration of a preoperative 3d liver volume to a laparoscopy image using contour and shading cues,”  in International conference on medical image computing and computer-assisted intervention. Springer, 2017, pp. 326–334.
  * [15] Pablo Azagra, Carlos Sostres, Ángel Ferrández, Luis Riazuelo, Clara Tomasini, O León Barbed, Javier Morlana, David Recasens, Víctor M Batlle, Juan J Gómez-Rodríguez, et al.,  “Endomapper dataset of complete calibrated endoscopy procedures,”  Scientific Data, vol. 10, no. 1, pp. 671, 2023.


# Mapping TLE orbital parameters to GNSS ephemeris for LEO PNT mega-
constellation orbit simulations and visibility analysis

Miquel Garcia-Fernandez
Rokubun, Av. Paral·lel 88, 08015, Barcelona, Spain
<EMAIL_ADDRESS>

###### Abstract

The emergence of Low Earth Orbit (LEO) satellite constellations dedicated to
positioning applications holds the promise of improving the capabilities of
existing Global Navigation Satellite Systems (GNSS). However, the absence of
operational systems necessitates a qualitative assessment of potential
improvements through simulation. This paper introduces a methodology to
convert Two Line Element (TLE) orbital parameters—abundantly available for LEO
constellations for communication and Earth Observation—into the widely used
RINEX 4 format employed by GNSS.

The primary goal is to establish a comprehensive database of LEO constellation
orbits directly compatible with the orbit propagation algorithms utilized in
GNSS systems like the Global Positioning System (GPS). This approach enables
seamless integration into simulation tools with minimal adjustments. While TLE
parameters are optimized for the SGP4 propagation model and cautioned against
use in classical Kepler orbit propagation scenarios requiring precision, the
obtained discrepancies, within a few tens of kilometers, suggest that these
representations are realistic for simulation purposes, as demonstrated with
the Spire LEMUR LEO constellation.

As a practical application, the paper conducts a visibility analysis using the
Starlink constellation. Results affirm expectations, showcasing that the
combination of GNSS with a LEO mega-constellation significantly enhances
satellite coverage and reduces Dilution of Precision. This work contributes to
the ongoing discourse on the potential benefits and practicality of
integrating emerging LEO constellations with established GNSS systems,
offering insights into improved navigation and timing capabilities through
simulation-based assessments.

_K_ eywords TLE $\cdot$ orbit mechanics $\cdot$ GNSS $\cdot$ LEO PNT $\cdot$
mega-constellations

##  1 Introduction

The emergence of Low Earth Orbit (LEO) satellite constellations dedicated to
Position, Navigation, and Timing (PNT), exemplified by initiatives such as
Xona or Geely ([1]), signifies a transformative phase in satellite positioning
technologies. While Medium Earth Orbit (MEO) Global Navigation Satellite
Systems (GNSS) like GPS or Galileo have been pivotal in furnishing precise
positioning, the advent of LEO PNT constellations will diversify satellite-
based navigation capabilities. However, these systems are typically under the
control of corporate institutions, limiting access to research data. Despite
ongoing efforts in simulation studies ([2]), the scarcity of available data
hampers comprehensive performance assessments through either real or simulated
datasets.

To address this challenge, this work proposes a novel approach by harnessing
Two Line Element (TLE) sets, widely employed for conveying orbital
information, particularly for LEO satellites. As outlined in [3], TLE
parameters are specifically tailored for orbital propagators like the
Simplified General Perturbations (SGP4) model, offering LEO orbit
representations with errors in the range of kilometers ([4]).

In the domain of GNSS, the Receiver Independent Exchange Format (RINEX) stands
as a universally accepted standard for representing and exchanging satellite
navigation data, encompassing orbits, clocks, and raw measurements ([5]).
Designed primarily for MEO GNSS systems, RINEX ensures compatibility across
diverse receiver manufacturers and analysis tools.

While simulation tools can incorporate SGP4 propagation algorithms to
represent LEO PNT orbits, having this representation in a format widely
adopted by the GNSS community, such as RINEX ([6]), would offer convenience.
This becomes particularly crucial for establishing a consistent satellite
naming convention that seamlessly links satellite orbits with their
corresponding raw measurement sets within the same timescale. Such
standardization expedites development steps, allowing researchers/users to
utilize orbit propagators described in GNSS Interface Control Documents,
compatible with the orbit description in RINEX navigation files.

However, adapting LEO constellation data, especially derived from TLE
parameters, to the RINEX format necessitates a specific mapping strategy, as
RINEX files correspond to classical Keplerian orbital parameters. While
cautionary notes ([7]) advise against this mapping for accurate orbit
representation, this work acknowledges that for simulation purposes, a
realistic (approximated) representation suffices. Thus, a TLE-based RINEX
navigation file is proposed here, with a slight modification to RINEX 4 to
accommodate the large number of LEO PNT satellites (often exceeding 100, as
per [8]). Importantly, this strategy, mapping TLE to RINEX (TLE-RNX),
circumvents the complexities of designing orbital parameters for an entire LEO
mega-constellation, leveraging the abundant TLE sets available on servers like
Celestrak.

This work introduces the methodology for this mapping in the first section,
along with an assessment of the expected orbit accuracy achievable with the
proposed method. Subsequent sections extend the application of this
methodology by conducting a visibility analysis with the Starlink mega-
constellation. A conclusions section wraps this paper.

##  2 Methodology

This section outlines the methodology employed for the conversion of Two Line
Element (TLE) orbital parameters into the RINEX format utilized in Global
Navigation Satellite Systems (GNSS). As previously mentioned, this conversion
aims to provide a broad representation of the orbit. The primary objective is
twofold: firstly, to capitalize on the abundance of publicly available orbit
files for Low Earth Orbit (LEO) constellations accessible through servers like
Celestrak, and secondly, to leverage orbit propagation algorithms utilized by
Medium Earth Orbit (MEO) GNSS systems, exemplified by the Global Positioning
System (GPS) algorithms detailed in Table 20-IV of the GPS Interface Control
Document ([9]).

By doing this, simulation tools could benefit from orbit representation in
currently existing formats (with minimum modifications, as shown in e.g. [8]).

The mapping strategy used is described in the following steps, based on the
necessary parameters of the RINEX file (and its computation using the related
TLE parameter)

  * •

Time of ephemeris ($toe$) is the epoch specified in the TLE field (fields 7
and 8 from line 1) assuming UTC timescale.

  * •

Mean Anomaly at Reference Time ($M_{0}$), is also the mean anomaly of the TLE
(field 7 from line 2) converted to radians

  * •

Eccentricity ($e$), is also the eccentricity of the TLE (field 5 of line 2)

  * •

Square Root of the Semi-Major Axis ($\sqrt{A}$) can be computed from the mean
motion of the TLE ($n$, field 8 of line 2) and the Earth Gravitation constant
($\mu=3.986005e14m^{3}/s^{2}$), using the expression:

$\sqrt{A}=\frac{\sqrt[6]{\mu}}{\sqrt[3]{n}}$ (1)

  * •

Longitude of Ascending Node of Orbit Plane at Weekly Epoch ($\Omega_{0}$) is
the parameter that needs further attention for two reasons: it needs to be
provided in an Earth Centered Earth Fixed reference frame (the TLE parameters
are referred to an Earth Centered Inertial reference frame) and it needs to be
referred to the start of the week (not at the reference epoch $toe$), in order
to be compatible with the computation of the "_Corrected longitude of
ascending node_ " ($\Omega_{k}$) specified in Table 20-IV of the GPS Interface
Control Document ([9]), reproduced here for convenience ($\dot{\Omega}_{e}$ is
the Earth rotation rate):

$\Omega_{k}=\Omega_{0}+(\dot{\Omega}-\dot{\Omega}_{e})\cdot
tk-\dot{\Omega}_{e}\cdot toe$ (2)

Therefore, to compute $\Omega_{0}$, the following steps are necessary:

    1. 1.

Take the TLE ephemeris epoch and compute the weeks and seconds within the week
elapsed since an arbitrary epoch ($t_{0})$. The choice of $t_{0}$ is not
critical as long as it starts on Sunday (due to the specification in e.g. GPS
system). Therefore, a possible date (actually used in this work) is the GPS
reference epoch ($t_{0}$ set to January 6th 1980). The epoch corresponding to
the start of the week ($t_{week}$) is the epoch reconstructed from the number
of weeks since $t_{0}$ (i.e. ignoring the seconds of the week).

    2. 2.

Compute the Right Ascension of the Greenwich meridian at $t_{week}$, i.e.
$\Omega_{G}(t_{week})$. This can be computed with the local sidereal time at
longitude $0^{\circ}$ and $t_{week}$ and then transforming from hours to
degrees/radians. Open source tools such as Poliastro [10] (using the method
poliastro.earth.util.get_local_sidereal_time) can be used for this purpose.

    3. 3.

Subtract the right ascension of the ascending node of the TLE ($\Omega_{TLE}$,
field 4 of line 2) so that:

$\Omega_{0}=\Omega_{TLE}-\Omega_{G}(t_{week})$ (3)

You might also need to compute the modulus operation to ensure that
$\Omega_{0}\in[0,2\pi)$.

  * •

Inclination Angle at Reference Time ($i_{0}$) is the inclination of the TLE
(field 3, line 2), converted to radians.

  * •

Argument of Perigee ($\omega$) is the argument of perigee of the TLE (field 6,
line 2), converted to radians.

  * •

Mean Motion Difference From Computed Value ($\Delta n$), which will be later
multiplied by the time elapsed since the $toe$ corresponds to the first
derivative of the mean motion (i.e. ballistic coefficient, field 9 of line 1).

Other parameters such as rate of Right Ascension ($\dot{\Omega}$), rate of
inclination angle ($idot$) and the amplitude of sine and cosine harmonics
($C*$) are set to 0.

###  2.1 Accuracy of TLE-RNX

As previously articulated, the purpose of generating TLE-RNX files (i.e.,
RINEX navigation files, as delineated in [5] and [8]) for Low Earth Orbit
(LEO) constellations based on Two Line Element (TLE) parameters is not
centered on furnishing highly precise or accurate orbital information.
Instead, the aim is to establish a comprehensive database comprising orbital
parameters that align with the requisites of orbital propagators defined in
the Global Navigation Satellite Systems (GNSS) Interface Control Documents.
Notably, the TLE parameters are specifically tailored for compatibility with
the SGP4 orbit propagation model (see [3]), and are not intended for
interpretation as classical Keplerian orbital parameters. Despite this warning
notice, they have been utilized in this study as a coarse representation of
the satellite orbits (a sample orbit is shown in Figure 1).

Figure 1: Sample orbit for Spire satellite 149 obtained from SP3 and orbit
propagation using TLE-RNX for July 4th 2023.

To assess the disparities achievable through the proposed approach in
constructing a massive database of Low Earth Orbit (LEO) satellite orbits in
RINEX format, the Spire constellation serves as a pertinent case study. The
rationale behind selecting Spire satellites lies in the fact that, on the one
hand, the TLE data for these satellites are included in the Celestrak
database, and on the other hand, precise orbit representation in SP3 format is
also available in the Radio-Occultation COSMIC UCAR product database (i.e.
leoOrb product). The discrepancies between the two orbits (TLE-RNX and SP3)
are shown in Figure 2 for a single day (July 4th 2023) and for the cross-,
along- and radial components.

For the radial and cross components, maximum discrepancies of approximately 10
km can be anticipated, while larger discrepancies, reaching up to 50 km, are
evident in the along component. Notably, these discernible discrepancies
exhibit a well-defined pattern and repeatability, underscoring the inadequacy
of the proposed approach as a precise representation of the orbit.
Nevertheless, these discrepancies, although non-negligible, remain
sufficiently modest to be considered valuable for simulation purposes. This
holds true, at least until operational Low Earth Orbit (LEO) Positioning,
Navigation, and Timing (PNT) constellations are deployed, and their
operational institutions release the accurate orbital parameters.

Figure 2: Single day comparison of Spire satellite positions from SP3 files
and positions computed from orbital parameters derived from TLE files.

##  3 Visibility analysis

A practical application of the proposed strategy for generating a database of
Low Earth Orbit (LEO) orbits in the standard RINEX format involves conducting
a comprehensive satellite visibility analysis for specific locations. This
analysis helps to evaluate the enhancements attainable through the utilization
of LEO constellations in combination with GNSS concerning Dilution of
Precision (DOP) and general coverage. As it is already known, for MEO GNSS,
there exists a satellite deficit in the North azimuths in northern mid-
latitudes, and visibility is notably compromised in high latitudes, where no
satellites are observable in the zenith direction (high elevation).

As shown in Figure 3, augmenting MEO GNSS constellations with LEO ones,
specially if highly populated with a mega-constellation such as Starlink,
improves the satellite visibility and coverage. This is specially true in
regions where MEO GNSS constellations may fall short (notably, addressing gaps
in the north for locations at northern mid-latitudes and enhancing visibility
at high elevations in high-latitude areas). This improvement is also shown in
Figure 4, where the Dilution of Precision (DOP) is shown for two locations and
for the cases with only GNSS as well as the combined GNSS + LEO (Starlink)
case.

|
---|---
Figure 3: Starlink and GNSS combined visibility on 2023-11-02 (between 12h and 16h) for two locations (left panel) Barcelona, Spain (mid-latitude) and (right panel) Tromsø, Norway (high-latitude |
---|---
Figure 4: Dilution of precision for GNSS and GNSS + Starlink on 2023-11-02
(between 12h and 16h) for two locations (left panel) Barcelona, Spain (mid-
latitude) and (right panel) Tromsø, Norway (high-latitude

##  4 Conclusion

This paper introduces a methodology for mapping Two Line Element (TLE)
parameters of Low Earth Orbit (LEO) mega-constellations into orbital
parameters compatible with the RINEX 4 format, enabling seamless integration
with Global Navigation Satellite Systems (GNSS) orbit propagation algorithms.
As it is known, the TLE parameters are finely tuned for the SGP4 propagation
model and not suited for classical Keplerian orbit propagation models.
However, the discrepancies observed with the proposed approach relative to
precise SP3 orbits, for the Spire LEMUR LEO constellation, which amounts to
tens of km, are within reasonable limits to use this strategy as a tool for
simulation.

The visibility analysis using TLE-based RINEX files for the Starlink mega-
constellation, boasting approximately 5000 LEO satellites, affirms the
expected advantages of combining Medium Earth Orbit (MEO) and LEO GNSS. The
results showcase a substantial enhancement in satellite coverage and a
reduction in Dilution of Precision values, reinforcing the assumption that LEO
Positioning, Navigation, and Timing (PNT) systems hold potential for improving
position robustness, particularly in terms of convergence time.

The wealth of TLE data available for LEO constellations in communication and
Earth Observation sets the stage for creating a comprehensive database of
orbits in a slightly modified RINEX 4 format, so that these simulated LEO can
be used for PNT purposes. This lays the foundation for qualitative assessments
through simulation, addressing the current lack of operational systems. In
this context, our work serves as a stepping stone for further advancements in
simulation steps, including the computation of synthetic raw measurements
(pseudorange, carrier phase, Doppler) from LEO PNT constellations.
Incorporating error factors such as ionosphere and troposphere, and advancing
hybridization algorithms between MEO and LEO GNSS, will further enhance key
performance indicators such as accuracy, coverage, convergence time, and
integrity.

## Acknowledgments

The author extend their gratitude to the team at the Earth Intelligence
Department of Spire Global UK for sharing insights into the correlation
between NORAD ID catalog numbers and satellite identification numbers as
specified in the SP3 files hosted at the UCAR Radio Occultation server. Their
contribution has significantly helped to obtain the results presented in this
paper.

The author wish to thank also Celestrack and COSMIC UCAR to make available the
TLE element sets and SP3 orbit files respectively. This data is available at
their servers https://celestrak.org/ and (https://data.cosmic.ucar.edu/gnss-
ro).

## References

  * [1] Fabricio S Prol, R Morales Ferre, Zainab Saleem, Petri Välisuo, Cristina Pinell, Elena-Simona Lohan, Mahmoud Elsanhoury, Mohammed Elmusrati, Saiful Islam, Kaan Çelikbilek, et al.  Position, navigation, and timing (pnt) through low earth orbit (leo) satellites: A survey on current status, challenges, and opportunities.  IEEE Access, 2022.
  * [2] Fabricio S Prol, Sanna Kaasalainen, E Simona Lohan, M Zahidul H Bhuiyan, Jaan Praks, and Heidi Kuusniemi.  Simulations using leo-pnt systems: A brief survey.  In 2023 IEEE/ION Position, Location and Navigation Symposium (PLANS), pages 381–387. IEEE, 2023.
  * [3] Felix R Hoots and Ronald L Roehrich.  SPACETRACK REPORT NO. 3: Models for propagation of NORAD element sets.  Office of Astrodynamics, 1980.  Accessed: December 7th 2023.
  * [4] Jizhang Sang and James C Bennett.  Achievable debris orbit prediction accuracy using laser ranging data from a single station.  Advances in Space Research, 54(1):119–124, 2014.
  * [5] IGS/RTCM RINEX Working Group et al.  The receiver independent exchange format, version 4.00. 1 december 2021\.  https://files.igs.org/pub/data/format/rinex_4.00.pdf, 2021.  Accessed: December 7th 2023.
  * [6] Mikaela Dobbin and Penina Axelrad.  A flexible ephemeris representation for gnss and alternative pnt signal sources using b-splines.  NAVIGATION: Journal of the Institute of Navigation, 70(4), 2023\.
  * [7] David Harding.  Bad science: Converting tles to keplerian elements.  https://blog.hardinglabs.com/tle-to-kep.html, 2019.  Accessed: December 7th 2023.
  * [8] Miquel Garcia-Fernandez.  A possible rinex extension for leo pnt ephemeris.  https://www.rokubun.cat/a-possible-rinex-extension-for-leo-pnt-ephemeris/, 2023\.  Accessed: December 7th 2023.
  * [9] Tony Anthony and Mary Kerns.  Gps interface specification document is-gps-200n.  https://www.gps.gov/technical/icwg/IS-GPS-200N.pdf, 2022.  Accessed: December 7th 2023.
  * [10] Juan Luis Cano Rodríguez, Helge Eichhorn, and Frazer McLean.  Poliastro: an astrodynamics library written in python with fortran performance.  In 6th International Conference on Astrodynamics Tools and Techniques, 2016.