# 6D Object Pose Regression via Supervised Learning on Point Clouds

Ge Gao1, Mikko Lauri1, Yulong Wang2, Xiaolin Hu2, Jianwei Zhang1 and Simone
Frintrop1 *This work is partially funded by the German Science Foundation
(DFG) in project Crossmodal Learning, TRR 169, and the National Natural
Science Foundation of China (NSFC) under Grant No. 61621136008.1Department of
Informatics, University of Hamburg, Germany<EMAIL_ADDRESS><EMAIL_ADDRESS>Key Laboratory of
Intelligent Technology and Systems, Department of Computer Science and
Technology, Tsinghua University, China<EMAIL_ADDRESS>
###### Abstract

This paper addresses the task of estimating the 6 degrees of freedom pose of a
known 3D object from depth information represented by a point cloud. Deep
features learned by convolutional neural networks from color information have
been the dominant features to be used for inferring object poses, while depth
information receives much less attention. However, depth information contains
rich geometric information of the object shape, which is important for
inferring the object pose. We use depth information represented by point
clouds as the input to both deep networks and geometry-based pose refinement
and use separate networks for rotation and translation regression. We argue
that the axis-angle representation is a suitable rotation representation for
deep learning, and use a geodesic loss function for rotation regression.
Ablation studies show that these design choices outperform alternatives such
as the quaternion representation and L2 loss, or regressing translation and
rotation with the same network. Our simple yet effective approach clearly
outperforms state-of-the-art methods on the YCB-video dataset.

##  I Introduction

The problem of 6 degrees of freedom (6D) object pose estimation is to
determine the transformation from a local object coordinate system to a
reference coordinate system (e.g., camera or robot coordinate) [1]. The
transformation is composed of 3D location and 3D orientation. Robust and
accurate 6D object pose estimation is of primary importance for many robotic
applications such as grasping and dexterous manipulation. The recent success
of convolutional neural networks (CNNs) in visual recognition has inspired
methods that use deep networks for learning features from RGB images [2, 3,
4]. These learned features are used for inferring 6D object poses. Similarly,
CNNs can also be applied to RGB-D images and treat depth information as an
additional channel for feature learning [5, 6, 7, 8]. However, in some
scenarios, color information may not be available, and depth information is
not in the 2-dimensional matrix format (e.g., laser range finder data), which
can be easily processed with CNN-based systems. Depth information can also be
represented by a point cloud, which is an unordered set of points in a metric
3D space. In existing methods, point clouds are mainly used for pose
refinement [9, 10, 11, 12] or template matching with hand-crafted features
extracted from point clouds [13, 14]. Using point clouds in the registration
stage confines its usage scope and hand-crafted features are usually less
effective compared to deep-learned features.

In this work, we investigate how to accurately estimate the 6D pose of a known
object represented by a point cloud segment containing only geometric
information using deep networks. Our approach is inspired by PointNet [15], a
deep network for object classification and segmentation operating on point
clouds. We adapt the system to the problem of pose estimation. PointNet
provides a method to apply deep learning to unordered point sets, and it is a
suitable architecture for our purpose. Our 6D pose regression method can be
applied to any type of range data, e.g., data from laser range finders.

For developing our system, we investigate three open questions. The first is
how to efficiently use depth information in a deep learning-based system.
Although it has been shown by many applications that CNNs can extract powerful
features from RGB-D information for specific tasks, due to the inherent
difference between color and depth information, it is unclear whether this is
an efficient way to treat depth information. We argue that a point cloud is a
more suitable structure and should be used in the scope of both deep networks
and geometry-based optimization.

The second question is whether translation and orientation should be estimated
with separate networks or a single network in a supervised learning system.
During a supervised learning process, a network learns the mapping from its
input to the desired output guided by a loss function. Since the metric units
for translation and orientation are different (i.e., meters and radians), we
argue that regressing them using separate networks and loss functions is a
more suitable choice. Our experiments show that an architecture with separate
networks outperforms those with shared layers.

Another question is the choice of rotation representation and the loss
function for measuring the distance between two rotations. Quaternions have
been a popular choice for many learning-based systems [11, 16]. However,
quaternions have the unit-norm constraint which imposes a limit on the network
output range. We argue that axis-angle is a more suitable choice because it is
a constraint-free representation. Concerning loss functions, L2 loss is a
popular choice for measuring the distance between two rotations [16, 17]. We
argue that the geodesic distance is a more suitable choice since it is well
justified mathematically, and provides a clearer learning goal compared to the
L2 loss. We show experimentally that our arguments are valid.

Figure 1: System overview. A point cloud is created using the depth data and
the output from a semantic segmentation method. This segment is processed with
farthest point sampling to obtain a down-sampled segment with consistent
surface structure. The segment with object class information is fed into two
networks for rotation and translation prediction. The geometry-based iterative
closest point algorithm is used for pose refinement.

Our contributions are thus as follows:

  * •

We present a simple yet effective system that infers the 6D pose of known
objects represented by point cloud segments containing geometric information.
This system is based on PointNet and exploits the point cloud structure by
utilizing it as the input to both deep networks and an iterative closest point
(ICP) refinement process. To the best of our knowledge, this is the first deep
learning architecture based system that regresses 6D poses from only unordered
point sets.

  * •

We demonstrate that the proposed method outperforms state-of-the-art methods
on a public dataset. Experimental results show that our system outperforms
methods that use both color and depth information during pose inference stage.
This evaluation result indicates that the proposed system is an efficient way
to use depth information for pose estimation.

  * •

Ablation studies provide an evaluation of each system component. We show
experimentally that each design choice has an impact on system performance.

##  II Related Work

Pose estimation has been well studied in the past using both color [2, 3, 4]
and depth data [18, 19]. Since in this work we focus on investigating how to
use geometric information during the pose inference stage, we mainly review
works that use depth information.

With the common usage of depth cameras and laser range finders (e.g. LiDAR),
methods using depth information have been proposed [18, 20, 21, 1, 22]. LINE-
MOD [18, 23] is one of the first works that use hand-crafted features from
depth information for pose estimation. It uses surface normals as part of its
local patch features. This patch representation is adapted and used with
random forest in [20, 19]. Surface normal is also used as an additional
modality in [6]. Another way of using depth information is to treat it as an
extra image depth channel (RGBD) and feed it into a CNN [5, 24, 6, 7, 8], or
random forest [21, 1, 22] or a fully connected sparse autoencoder [25] for
feature extraction. Depth is also used to create point clouds, which are used
for generating pose hypotheses with 3D-3D correspondences and ICP refinement
in [12]. Point clouds are used to facilitate Point-to-Point matching in [13,
14]. Geometry embeddings are extracted from point cloud segments with a deep
network in [16]. Approaches such as [11] use color information to provide an
initial pose estimate, then refine it with ICP using depth information. Point
cloud segments are also used for pose estimation in [26]. However, they only
report experimental data for pose estimation for a single angle as opposed to
all three as we do. Furthermore, they formulate rotation estimation as
classification and discretize the rotation angles to bins. It is not clear if
this approach would scale up to three rotation angles. In our previous work
[27], we only considered rotation regression for single objects. Here, we
regress the full 6D pose in a multi-class setting.

For learning-based systems, [11] propose to predict translation and rotation
with separate networks sequentially. [16] propose to regress translation and
rotation with the same network. Quaternions are used as the rotation
representation for regression [11, 16, 8]. Bui et al. [8] propose to use L2
loss function for rotation learning. Axis-angle representation is also used in
[28]. However, they only address estimating 3D rotation from RGB information,
while we address both 3D rotation and 3D translation from point clouds. Works
that use deep networks for extracting features and do nearest neighbor search
pose retrieval or classification into discretized bins to obtain object poses
[5, 6, 7, 24, 26] are not in the scope of this work and are not discussed in
detail here.

Our work is most similar to DenseFusion [16] as we both use PointNet [15] to
extract features from point cloud segments. However, there are two significant
differences: first, during the pose estimation stage, we only use coordinate
information from point clouds while DenseFusion also uses color information
extracted by a CNN; second, we design regression targets and loss functions
for rotation and translation individually while DenseFusion uses one
regression target for the 6D pose.

##  III System Architecture

Figure 1 shows an overview of our system. The proposed system for object 6D
pose estimation is a multi-class system, i.e., we use the same system to
predict poses for objects from different classes. Hence, an object segment, as
well as the corresponding class information, is required as input to the pose
estimation networks. As semantic segmentation is a well-studied problem, we
assume the object segment and class information is provided by an off-the-
shelf method111Here, we use the semantic segmentation from [11]. and focus on
the object pose estimation from a point cloud segment in this work. A point
cloud segment is created using depth and the target object segment. This
segment is processed with Farthest Point Sampling (FPS) [29] to obtain a down-
sampled segment with a consistent surface structure representation. This down-
sampled segment and class information are combined as the input for two
separate networks for rotation estimation in axis-angle representation and
translation prediction through translation residual regression. The 6D pose is
refined with a geometry-based optimization process to produce the final pose
estimate.

Figure 2a illustrates BaseNet, which is the basic building block of our
system. BaseNet is an adapted version of PointNet [15]. Given a point cloud
with $n$ points as input, PointNet is invariant to all $n!$ possible
permutations. Each point is processed independently using multi-layer
perceptrons (MLPs) with shared weights. Compared to PointNet, we remove the
spatial transformer blocks and adapt the dimension for the output layer to be
$3$. For each input point with class information, a feature vector is learned
with shared weights. These feature vectors are max-pooled to create a global
representation of the input point cloud. Finally, we use a three-layer
regression MLP on top of this global feature to predict the poses.

Figure 2b shows a more detailed diagram of our system. We use two separate
networks to handle translation and rotation estimation. Input to the rotation
network is a point cloud with $n$ points concatenated with the one-hot encoded
class information. In total, the input is a $n$ by $(3+k)$ array where $k$ is
the total number of classes. The output of the rotation network is the
estimated rotation in axis-angle representation. The translation network takes
normalized point coordinates concatenated with class information and estimates
the translation residual. The full translation is obtained by adding back the
coordinate mean.

(a)

(b)

Figure 2: (a) The architecture of BaseNet. The numbers in parentheses are
numbers of MLP layers. Numbers not in parentheses indicate the dimensions of
intermediate feature vectors. A feature vector for each point is learned with
shared weights. A max pooling layer aggregates the individual features into a
global feature vector. A regression network with 3 fully-connected layers
outputs the pose prediction. (b) Diagram for input and output of our pose
networks. For the rotation network, the input is point coordinate information
concatenated with class information per point, the output is rotation in axis-
angle representation. For the translation network, the input coordinates are
normalized by removing the mean. It outputs translation residual. The full
translation is obtained by adding back the coordinate mean. The number of
input points is $n$ and $k$ is the total number of classes.

##  IV Supervised Learning for 6D Pose Regression

This section describes how the loss functions for 6D pose regression are
formulated in our supervised learning framework. Given a set of points
$\mathcal{X}=\left\\{\mathbf{x}_{i}\in\mathbb{R}^{3}\mid
i=1,\ldots,n\right\\}$ on the surface of a known object in the camera
coordinate, the aim of pose estimation is to find a transformation that
transforms $x_{i}$ from the object coordinates to the camera coordinates. This
transformation is composed of a translation and a rotation. A translation
consists of the displacements along the three coordinate axes. A rotation
specifies the rotation around the three coordinate axes, and it has different
representations such as axis-angle, Euler angles, quaternion, and rotation
matrix. For supervised learning, suitable loss functions are required to
measure the differences between predicted poses and ground truth poses. For
rotation learning, we argue that the axis-angle representation is the best
suited for the learning task. Geodesic distance is used as the loss function
for rotation regression. For translation learning, we predict the residual of
translation.

###  IV-A Rotation estimation with axis-angle and geodesic distance loss

In the axis-angle representation, a vector $\mathbf{r}\in\mathbb{R}^{3}$
represents a rotation of $\theta=\left\lVert\mathbf{r}\right\rVert_{2}$
radians around the unit vector
$\frac{\mathbf{r}}{\left\lVert\mathbf{r}\right\rVert_{2}}$ [30]. Given an
axis-angle representation
$\mathbf{r}=\begin{bmatrix}r_{1}&r_{2}&r_{3}\end{bmatrix}^{T}$, the
corresponding rotation matrix $R$ is obtained via the exponential map $\exp:$

$R=\exp(\mathbf{r}_{\times})=\mathrm{I}_{3\times
3}+\frac{\sin\theta}{\theta}\mathbf{r}_{\times}+\frac{1-\cos\theta}{\theta^{2}}\mathbf{r}^{2}_{\times},$
(1)

where $\mathrm{I}_{3\times 3}$ is the identity matrix and
$\mathbf{r}_{\times}$ is the skew-symmetric matrix

$\mathbf{r}_{\times}=\begin{bmatrix}0&-r_{3}&r_{2}\\\ r_{3}&0&-r_{1}\\\
-r_{2}&r_{1}&0\end{bmatrix}.$ (2)

For rotation learning, we regress to a predicted rotation $\hat{\mathbf{r}}$.
Prediction $\hat{\mathbf{r}}$ is compared with ground truth rotation
$\mathbf{r}$ via a rotation loss function $l_{r}$, which is the geodesic
distance between $\hat{R}$ and $R$ [31, 30]:

$l_{r}(\hat{\mathbf{r}},\mathbf{r})=\arccos\left(\frac{\mathrm{trace}(\hat{R}R^{T})-1}{2}\right),$
(3)

where $\hat{R}$ and $R$ are the two rotation matrices corresponding to
$\hat{\mathbf{r}}$ and $\mathbf{r}$, respectively.

This loss function directly measures the magnitude of rotation difference
between $\hat{R}$ and $R$, so it is convenient to interpret. Furthermore, the
network can make constraint-free predictions with axis-angle representations,
in contrast to e.g. quaternion representations which require normalization.

###  IV-B Translation residual estimation

To simplify the learning task by reducing the variance in regression space for
translation prediction, the learning target is chosen to be the residual of
translation. Given a translation residual $\widehat{\Delta\mathbf{t}}$, full
translation prediction $\hat{\mathbf{t}}$ is obtained via

$\hat{\mathbf{t}}=\widehat{\Delta\mathbf{t}}+\mathbf{{\mu}_{t}},$ (4)

where $\mathbf{{\mu}_{t}}$ is the mean of $\mathcal{X}$. $L2$-norm is used to
measure the distance between prediction $\hat{\mathbf{t}}$ and ground truth
$\mathbf{t}$, resulting in the translation loss function
$l_{t}(\hat{\mathbf{t}},\mathbf{t})$:

$l_{t}(\hat{\mathbf{t}},\mathbf{t})=\left\lVert\mathbf{t}-\hat{\mathbf{t}}\right\rVert_{2}.$
(5)

###  IV-C Total loss function for 6D pose regression

The total loss is defined as the combination of the translation and the
rotation loss:

$l(\mathbf{t},\hat{\mathbf{t}},\mathbf{r},\hat{\mathbf{r}})=\alpha
l_{t}(\hat{\mathbf{t}},\mathbf{t})+l_{r}(\hat{\mathbf{r}},\mathbf{r}),$ (6)

where $\alpha$ is a scaling factor. The total loss is used for training the
pose estimation networks.

##  V Experiments

We evaluate the proposed system on the YCB-Video dataset [11] and compare the
performance with the state-of-the-art methods PoseCNN [11] and DenseFusion
[16]. We also compare the performance on a subset of the object classes with a
state-of-the-art RGB-based method DOPE [4]. Besides prediction accuracy and
performance under occlusions, we also investigate the impact of using
different network structures, as well as the influence of different rotation
representations. The implementation of our system is available on-
line222https://github.com/GeeeG/CloudPose.

###  V-A Experiment setup

The YCB video dataset [11] contains 92 video sequences with total 133,827
frames of 21 objects selected from the YCB object set [32] with 6D pose
annotations. We follow the official train/test split and use 80 video
sequences for training. Testing is performed on the $2,949$ key frames chosen
from the remaining 12 sequences [11]. 80,000 frames of synthetic data are also
provided by YCB-Video dataset as an extension to the training set. During
training, Adam optimizer is used with a learning rate of $0.0008$. The batch
size is $128$. For the total loss, we use $\alpha=10$, which is given by the
ratio between the expected error of translation and rotation at the end of the
training [33]. The number of points of the input point cloud segment is
$n=256$. Batch normalization is applied to all layers. No dropout is used. All
of our networks are trained for $90$ epochs. For refinement, we use the Point-
to-Point ICP registration provided by Open3D [34] and refine for 10
iterations. The initial search radius is $0.01$m and is reduced by $10\%$
after each iteration. For a fair comparison, all methods use object
segmentation provided by PoseCNN during testing.

###  V-B Evaluation metrics

We use the average distance (AD) of model points and the average distance for
a rotationally symmetric object (ADS) proposed in [23] as evaluation metrics.
Given a 3D model represented as a set $\mathcal{M}$ with $m$ points, ground
truth rotation $R$ and translation $\mathbf{t}$, as well as estimated rotation
$\hat{R}$ and translation $\mathbf{\hat{t}}$, the AD is defined as:

$\displaystyle\mathrm{AD}=\frac{1}{m}\displaystyle\sum_{\mathbf{x}\in\mathcal{M}}\left\lVert(R\mathbf{x}+\mathbf{t})-(\hat{R}\mathbf{x}+\hat{\mathbf{t}})\right\rVert_{2}.$
(7)

ADS is computed using closest point distance. It provides a distance measure
that considers possible pose ambiguities caused by rotational symmetry:

$\displaystyle\mathrm{ADS}=\frac{1}{m}\displaystyle\sum_{\mathbf{x}_{1}\in\mathcal{M}}\min_{\mathbf{x}_{2}\in\mathcal{M}}\left\lVert(R\mathbf{x}_{1}+\mathbf{t})-(\hat{R}\mathbf{x}_{2}+\hat{\mathbf{t}})\right\rVert_{2}.$
(8)

A 6D pose estimate is considered to be correct if AD and ADS are smaller than
a given threshold. We report the area under error threshold-accuracy curve
(AUC) for AD and ADS. The maximum thresholds for both curves are set to
$0.1$m. Furthermore, we also provide ADS accuracy with threshold $0.01$m
($<$1cm) to illustrate the performance accuracy under a smaller error
tolerance.

###  V-C Prediction accuracy

Evaluation results averaged for all 21 objects in the YCB-Video dataset are
shown in Table I. PoseCNN [11] uses RGB information to provide an initial pose
estimate (PC w/o ICP), then uses depth information with a highly customized
ICP for pose refinement (PC). DenseFusion [16] (DF) uses both color and point
cloud features extracted by deep networks to give per-pixel pose estimate for
final pose voting, and iterative pose refinement is performed with an extra
network module. Ours w/o ICP is the estimated pose from the proposed system
architecture (Section III), and Ours is the result after ICP refinement. We
also perform the ICP refinement on DF results (DF+ICP). For the overall
performance in Table I, we highlight the best performance in bold font.
Details regarding the data type used by pose regression networks and the post
process are also presented.

Our method achieves state-of-the-art performance using only depth information.
In terms of AD, we outperform both PC an DF. We observe that DF+ICP shows
small improvement compared to DF. One possible reason is the sensitivity of
ICP to the initial pose guess, if the method already performs well without
refinement, ICP is able to provide further gains. If the initial guess is
poor, ICP can even make the results worse. This result indicates that features
learned from depth information represented by unordered point clouds are
sufficient for accurately regressing 6D pose. Furthermore, this also shows
that the proposed approach is an efficient way to use depth information in a
deep learning framework for pose regression.

Performance for individual objects is shown in Table II. We use the trained
network for six objects provided by the authors of DOPE [4] and report the
results. The AD results are not available because the object coordinate frames
used in the YCB object dataset [32], YCB video dataset for PoseCNN [11] and
DOPE are different. As our method uses the frames from [11], and the
transformation between [32] and [11] is not publicly available, we can not
find the correspondence between model points required for AD. We also applied
ICP to DOPE pose estimates, but the performance was not improved. A possible
reason is the sensitivity of ICP to the initial pose estimate.

Some qualitative results are shown in Figure 3. Pose estimates from PC, DF and
our method are used for projecting object models onto 2D images. More
qualitative results are available in the supplementary video.

TABLE I: Quantitative evaluation of 6D pose on the YCB-Video Dataset [11].
Best performance is in bold font.

| RGB | Depth | ICP | AD | ADS | <1 cm
---|---|---|---|---|---|---
PC w/o ICP [11] | ✓ |  |  | 51.5 | 75.6 | 26.1
PC [11] | ✓ | ✓ | ✓ | 77.8 | 93.6 | 88.4
DF [16] | ✓ | ✓ |  | 74.7 | 93.9 | 87.6
DF [16] \+ ICP333we apply ICP to DF results after iterative refinement. | ✓ | ✓ | ✓ | 76.3 | 94.7 | 89.0
Ours w/o ICP |  | ✓ |  | 76.0 | 91.3 | 80.9
Ours |  | ✓ | ✓ | 82.7 | 94.7 | 90.3

TABLE II: Pose estimation accuracy per object class on the YCB-Video Dataset
[11]. Best per class performance for AD(S) is in bold font. Ours achieves the
best performance on a majority of object classes.

| DOPE [4] | PC [11] | DF [16] | Ours
---|---|---|---|---
|  ADS | <1cm | AD | ADS | <1cm | AD | ADS | <1cm | AD | ADS | <1cm
02_master_chef |  — | — | 68.1 | 95.8 | 99.5 | 73.2 | 96.4 | 100 | 46.9 | 95.4 | 95.4
03_cracker_box |  62.7 | 29.6 | 83.4 | 92.7 | 84.8 | 94.2 | 95.8 | 97 | 76.7 | 93 | 80.4
04_sugar_box |  85.0 | 33.4 | 97.1 | 98.2 | 100 | 96.5 | 97.6 | 100 | 97.5 | 98.5 | 99.7
05_tomato_soup |  88.5 | 74.5 | 83.6 | 96.6 | 99 | 87.4 | 96.6 | 99.1 | 72.7 | 96.5 | 96.8
06_mustard_bottle |  90.7 | 65.3 | 98 | 98.6 | 98.9 | 94.8 | 97.3 | 97.8 | 79.2 | 97.7 | 94.1
07_tuna_fish_can |  — | — | 83.9 | 97.1 | 97.6 | 81.8 | 97.1 | 99.5 | 72 | 97.7 | 100
08_pudding_box |  — | — | 96.6 | 97.9 | 100 | 93.2 | 95.9 | 98.6 | 94.4 | 97.3 | 91.1
09_gelatin_box |  84.6 | 36.9 | 98.1 | 98.8 | 100 | 96.7 | 98 | 100 | 98.6 | 99 | 100
10_potted_meat |  32.0 | 3.7 | 86 | 94.3 | 87.5 | 87.8 | 95 | 92 | 90.6 | 95.7 | 93.7
11_banana |  — | — | 91.9 | 97.1 | 95 | 83.6 | 96.2 | 98.2 | 95.1 | 97.7 | 95.5
19_pitcher_base |  — | — | 96.9 | 97.8 | 99.6 | 96.6 | 97.5 | 99.5 | 96.1 | 97.9 | 100
21_bleach_clean |  — | — | 92.5 | 96.9 | 95.1 | 89.7 | 95.8 | 99.4 | 95.4 | 97.4 | 98.4
24_bowl |  — | — | 14.4 | 81 | 42.9 | 5.9 | 89.5 | 55.7 | 83.9 | 97.7 | 99.3
25_mug |  — | — | 81.1 | 94.9 | 97.6 | 88.8 | 96.7 | 98 | 93.9 | 97.8 | 99.7
35_power_drill |  — | — | 97.7 | 98.2 | 99.3 | 93 | 96.1 | 97.8 | 94.9 | 97.7 | 96.7
36_wood_block |  — | — | 70.9 | 87.6 | 74.4 | 30.9 | 92.8 | 88.8 | 90 | 94.9 | 97.5
37_scissors |  — | — | 78.4 | 91.7 | 68 | 77.4 | 91.9 | 71.3 | 75.8 | 91.3 | 63
40_large_marker |  — | — | 85.3 | 97.2 | 97.1 | 93 | 97.6 | 100 | 92.2 | 98 | 100
51_large_clamp |  — | — | 52.2 | 75.3 | 67.4 | 26.4 | 72.6 | 33.3 | 68.5 | 77.4 | 69.6
52_e_large_clamp |  — | — | 25.9 | 74.9 | 48.2 | 16.6 | 77.4 | 10.9 | 25.3 | 66.4 | 22
61_foam_brick |  — | — | 48.1 | 97.2 | 99.7 | 59 | 92 | 100 | 92.9 | 98 | 99.3

PC [11] | DF [16] | Ours | PC [11] | DF [16] | Ours
---|---|---|---|---|---
|  |  |  |  |
|  |  |  |  |

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 3: Qualitative results for 6D pose estimation. From left to right:
PoseCNN (PC) [11], DenseFusion (DF) [16], and ours. The colored overlay
indicates the predicted pose of the target object. Our method gives more
accurate translation estimates, and also is able to give accurate rotation
estimation for texture-less object (e.g. the red bowl). More qualitative
results are available in the supplementary video.

###  V-D Occlusion

For a given target object in a frame, the occlusion factor $O$ of the object
is defined as [27]

$\displaystyle O=1-\frac{\lambda}{\nu},$ (9)

where $\lambda$ is the number of pixels in the 2D ground truth segmentation,
and $\nu$ is the number of pixels in the projection of the 3D object model
onto the image plane using the camera intrinsic parameters and the ground
truth 6D pose, when we assume the object would be fully visible. The occlusion
factor of the YCB-Video dataset ranges from $0.8\%$ to $81\%$. We divide this
range into $8$ bins with a bin width of $10\%$ and report the prediction
accuracy (ADS) with a threshold of $1$ cm. Figure 4 illustrates the results.
It can be observed that our method (Ours) has competitive performance when the
occlusion is lower than $40\%$, then both ours and PC start to suffer as the
amount of occlusion increases. One possible reason is that DF outputs per-
pixel prediction with confidence scores, while ours and PC provide only one
pose prediction. This per-pixel prediction may have helped to provide better
performance when the amount of occlusion is higher than $40\%$.

Figure 4: Effect of occlusion compared to PC [11] and DF [16]. The horizontal
axis denotes the upper limit (occlusion in $\%$) of each bin. Width for each
bin is $10\%$. Numbers in parentheses denote the number of samples in
corresponding bin. Ours is competitive with state-of-the-art methods when
occlusion is lower than $40\%$.

###  V-E Network architecture: ablation study

To investigate whether translation and rotation should be regressed with the
same or separate networks, we compare the performance of different
architectures. We alter the network architecture by incrementally sharing the
layers between translation and rotation networks. Table III shows the result
in terms of AD, ADS, and accuracies for translation and rotation under certain
thresholds. None denotes the proposed architecture which regresses translation
and rotation with two separate networks. The numbers in the first column
denote the number of shared layers between translation and rotation BaseNet
(Figure 2). We compare performance without ICP refinement. When sharing
layers, the performance is worse than using two separate networks.

We also tested an architecture that shares all the layers while having the
same amount of parameters as the proposed structure with a doubled layer
width. The performance is similar to the architecture with the single width,
and this verifies that the performance deterioration is not caused by
insufficient network capacity. This result verifies that using separated
networks for translation and rotation is a more suitable design choice.

TABLE III: Accuracy with different network structures. Best performance is in
bold font. Using two separate networks performs best.

shared layers | AD | ADS | rot_err<10° | tran_err<1 cm
---|---|---|---|---
none | 76.0 | 91.3 | 41.3 | 73.0
1 | 75.2 | 91.1 | 39.5 | 72.0
2 | 75.5 | 91.1 | 38.6 | 69.8
3 | 75.1 | 91.1 | 41.2 | 69.6
4 | 75.5 | 91.2 | 39.4 | 70.3
5 | 75.2 | 91.1 | 39.1 | 69.7
all | 63.0 | 87.8 | 23.3 | 69.3

TABLE IV: Prediction accuracy with different rotation representations and loss
functions. Best performance is in bold font.

Rotation representation | <10° | <15° | <20°
---|---|---|---
Axis-angle | 41.3 | 52.7 | 59.6
Quaternion | 37.1 | 48.5 | 56.7

(a) Different rotation representations. Geodesic distance is used as the loss
function for both cases.

Loss function | <10° | <15° | <20°
---|---|---|---
Geodesic | 41.3 | 52.7 | 59.6
L2 | 40.8 | 52.5 | 58.4

(b) Different loss functions for rotation regression. Axis-angle
representation is used for both cases.

###  V-F Rotation representation and loss function: ablation study

We investigate the impact of different rotation representations and loss
functions. For comparing quaternion to axis-angle, we adapted our rotation
network to have $4$-dimensional output instead of $3$. The output is
normalized and then converted to the axis-angle representation. We use the
same loss function as described in Section IV. For comparing L2 loss with
Geodesic distance, we keep the rotation representation in axis-angle format
and apply different loss functions. Table IV shows the accuracy of rotation
prediction with different thresholds. With the same loss function, using axis-
angle yields a better result than quaternion. This indicates that axis-angle
is a better choice for rotation learning. With the same rotation
representation, L2 loss slightly underperforms geodesic loss. Since geodesic
distance also has a better mathematical justification, this makes it a better
choice.

###  V-G Time performance

We measure the time performance on a Nvidia Titan X GPU. The system is
implemented with Tensorflow. Pose estimation by a forward pass through our
network takes 0.11 seconds for a single object. The 10 iterations of ICP
refinement require an additional 0.3 seconds.

##  VI Conclusion

We propose a system for fast and accurate 6D pose estimation of known objects.
We formulate the problem as a supervised learning problem and use two separate
networks for rotation and translation regression, and use point clouds as
input for the regression. We use axis-angle as rotation representation and
geodesic distance as the loss function for rotation regression. Ablation
studies show that these design choices outperform the commonly used quaternion
representation and L2 loss. Experimental results show that the proposed system
outperforms two state-of-the-art methods on a public benchmark.

To the best of our knowledge, this is the first deep learning system that
regresses 6D object poses from only depth information represented by unordered
point clouds. Features extracted from point clouds with deep networks can be
used for accurately regressing object pose. Our pose regression system can be
applied to range data from other sensors such as laser range finders. In the
future work, we will investigate aspects such as pose estimation for
rotational symmetry objects using only geometric information.

## References

  * [1] A. Krull, E. Brachmann, F. Michel, M. Y. Yang, S. Gumhold, and C. Rother, “Learning analysis-by-synthesis for 6D pose estimation in RGB-D images,” in _ICCV_ , 2015.
  * [2] W. Kehl, F. Manhardt, F. Tombari, S. Ilic, and N. Navab, “SSD-6D: making RGB-based 3D detection and 6D pose estimation great again,” in _ICCV_ , 2017.
  * [3] M. Oberweger, M. Rad, and V. Lepetit, “Making deep heatmaps robust to partial occlusions for 3D object pose estimation,” in _ECCV_ , 2018.
  * [4] J. Tremblay, T. To, B. Sundaralingam, Y. Xiang, D. Fox, and S. Birchfield, “Deep object pose estimation for semantic robotic grasping of household objects,” in _CoRL_ , 2018.
  * [5] W. Kehl, F. Milletari, F. Tombari, S. Ilic, and N. Navab, “Deep learning of local RGB-D pacthes for 3D object detection and 6D pose estimation,” in _ECCV_ , 2016.
  * [6] S. Zakharov, W. Kehl, B. Planche, A. Hutter, and S. Ilic, “3D object instance recognition and pose estimation using triplet loss with dynamic margin,” in _IROS_ , 2017.
  * [7] C. Li, J. Bai, and G. D. Hager, “A unified framework for multi-view multi-class object pose estimation,” in _ECCV_ , 2018.
  * [8] M. Bui, S. Zakharov, S. Albarqouni, S. Ilic, and N. Navab, “When regression meets manifold learning for object recognition and pose estimation,” in _ICRA_ , 2018.
  * [9] R. B. Rusu, N. Blodow, and M. Beetz, “Fast point feature histograms (FPFH) for 3D registration,” in _ICRA_ , 2009.
  * [10] C. Sahin, R. Kouskouridas, and T. Kim, “A learning-based variable size part extraction architecture for 6D object pose recovery in depth images,” _Image and Vision Computing_ , vol. 63, no. C, pp. 38–50, 2017.
  * [11] Y. Xiang, T. Schmidt, V. Narayanan, and D. Fox, “PoseCNN: a convolutional neural network for 6D object pose estimation in cluttered scenes,” in _RSS_ , 2018.
  * [12] O. H. Jafari, S. K. Mustikovela, K. Pertsch, E. Brachmann, and C. Rother, “iPose: Instance-aware 6D pose estimation of partly occluded objects,” in _ACCV_ , 2018.
  * [13] B. Drost, M. Ulrich, N. Navab, and S. Ilic, “Model globally, match locally: Efficient and robust 3D object recognition,” in _CVPR_ , 2000.
  * [14] S. Hinterstoisser, V. Lepetit, N. Rajkumar, and K. Konolige, “Going further with point pair features,” in _ECCV_ , 2016.
  * [15] C. R. Qi, H. Su, K. Mo, and L. J. Guibas, “PointNet: Deep learning on point sets for 3D classification and segmentation,” in _CVPR_ , 2017.
  * [16] C. Wang, D. Xu, Y. Zhu, R. Martín-Martín, C. Lu, L. Fei-Fei, and S. Savarese, “DenseFusion: 6D object pose estimation by iterative dense fusion,” in _CVPR_ , 2019.
  * [17] A. Kendall and R. Cipolla, “Geometric loss functions for camera pose regression with deep learning,” in _CVPR_ , 2017.
  * [18] S. Hinterstoisser, S. Holzer, C. Cagniart, S. Ilic, K. Konolige, N. Navab, and V. Lepetit, “Multimodal templates for real-time detection of texture-less objects in heavily cluttered scenes,” in _ICCV_ , 2011.
  * [19] A. Tejani, R. Kouskouridas, A. Doumanoglou, D. Tang, and T. K. Kim, “Latent-class hough forests for 6 DoF object pose estimation,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 40, no. 1, pp. 119–132, 2018.
  * [20] A. Tejani, D. Tang, R. Kouskouridas, and T. K. Kim, “Latent-class hough forests for 3D object detection and pose estimation,” in _ECCV_ , 2014\.
  * [21] E. Brachmann, A. Krull, F. Michel, J. S. S. Gumhold, and C. Rother, “Learning 6D object pose estimation using 3D object coordinates,” in _ECCV_ , 2014\.
  * [22] F. Michel, A. Kirillov, E. Brachmann, A. Krull, S. Gumhold, B. Savchynskyy, and C. Rother, “Global hypothesis generation for 6D object pose estimation,” in _CVPR_ , 2017.
  * [23] S. Hinterstoisser, V. Lepetit, S. Ilic, S. Holzer, G. Bradski, K. Konolige, and N. Navab, “Model based training, detection and pose estimation of texture-less 3D objects in heavily cluttered scenes,” in _ACCV_ , 2012\.
  * [24] V. Balntas, A. Doumanoglou, C. Sahin, J. Sock, R. Kouskouridas, and T. K. Kim, “Pose guided RGBD feature learning for 3D object pose estimation,” in _ICCV_ , 2017.
  * [25] A. Doumanoglou, R. Kouskouridas, S. Malassiotis, and T. K. Kim, “Recovering 6D object pose and predicting next-best-view in the crowd,” in _CVPR_ , 2016.
  * [26] C. R. Qi, W. Liu, C. Wu, H. Su, and L. J. Guibas, “Frustum PointNets for 3D object detection from RGB-D data,” in _CVPR_ , 2018.
  * [27] G. Gao, M. Lauri, J. Zhang, and S. Frintrop, “Occlusion resistant object rotation regression from point cloud segments,” in _ECCV 4th International Workshop on Recovering 6D Object Pose_ , 2018.
  * [28] S. Mahendran, H. Ali, and R. Vidal, “3D pose regression using convolutional neural networks,” in _CVPR workshops_ , 2017.
  * [29] Y. Eldar, M. Lindenbaum, M. Porat, and Y. Y. Zeevi, “The farthest point strategy for progressive image sampling,” _IEEE Transactions on Image Processing_ , vol. 6, no. 9, pp. 1305–1315, 1997.
  * [30] R. Hartley, J. Trumpf, Y. Dai, and H. Li, “Rotation averaging,” _International Journal of Computer Vision_ , vol. 103, no. 3, pp. 267–305, 2013.
  * [31] D. Q. Huynh, “Metrics for 3D rotations: Comparison and analysis,” _Journal of Mathematical Imaging and Vision_ , vol. 35, no. 2, pp. 155–164, 2009.
  * [32] B. Calli, A. Walsman, A. Singh, S. Srinivasa, and P. Abbeel, “Benchmarking in manipulation research using the Yale-CMU-Berkeley object and model set,” _Robotics & Automation Magazine, IEEE_, vol. 22, no. 3, pp. 36–52, 2015\.
  * [33] A. Kendall, M. Grimes, and R. Cipolla, “PoseNet: a convolutional network for real-time 6-DOF camera relocalization,” in _ICCV_ , 2015.
  * [34] Q. Zhou, J. Park, and V. Koltun, “Open3D: a modern library for 3D data processing,” _arXiv:1801.09847_ , 2018.


In 1926–1927, Tarski designed a set of axioms for Euclidean geometry which reached its final form in a manuscript by Schwabhäuser, Szmielew and Tarski in 1983.
The differences amount to simplifications obtained by Tarski and Gupta.
Gupta presented an independent version of Tarski's system of geometry, thus establishing that his version could not be further simplified without modifying the axioms.
To obtain the independence of one of his axioms, namely Pasch's axiom, he proved the independence of one of its consequences: the previously eliminated symmetry of betweenness.
However, an independence model for the non-degenerate part of Pasch's axiom was provided by Szczerba for another version of Tarski's system of geometry in which the symmetry of betweenness holds.
This independence proof cannot be directly used for Gupta's version as the statements of the parallel postulate differ.

In this paper, we present our progress towards obtaining an independent version of a variant of Gupta's system.
Compared to Gupta's version, we split Pasch's axiom into this previously eliminated axiom and its non-degenerate part and change the statement of the parallel postulate.
We verified the independence properties by mechanizing counter-models using the Coq proof-assistant.

§ INTRODUCTION

The independence[We recall that an axiom is said to be independent from a set of axiom if it is not derivable from the axioms in this set.] of axioms for geometry has often been an important topic in the field of geometry.
For centuries, many mathematicians believed that Euclid’s fifth postulate was rather a theorem which could be derived from the first four of Euclid’s postulates.
History is rich with incorrect proofs of Euclid’s fifth postulate.
In 1763, Klügel provided, in his dissertation, a survey of about 30 attempts to prove Euclid’s parallel postulate [13].
The question was finally settled in 1832 and 1840, when Bolyai [5] and Lobachevsky [14] exhibited models of hyperbolic geometry, thus establishing that this postulate was independent.
Later, Hilbert dedicated the second section of his famous Grundlagen der Geometrie [12] to independence properties.
Then, when working out the final version [22] of the axioms for Metamathematische Methoden in der Geometrie [19], commonly referred to as SST, independence results proved very helpful.

So it is a surprise that, now that we have access to tools like proof assistants which we believe to be perfectly suited for the task, the only independence to be mechanized was the one for Euclid’s fifth postulate [18].
To the best of our knowledge, the most recent work on the topic is the formalization of the Poincaré disk model in Isabelle/HOL [20].

In this paper, we study independence properties linked to SST [19].
SST has the advantage of being expressed in the first-order language rather than natural language which leaves room for interpretation leading to possible problems [8].[A possible interpretation of Hilbert's axiom could lead to a degenerate model for first two groups of Hilbert’s axioms.]
There are several ways to prove independence results [4].
Here we focus on independence through counter-model, i.e. constructing a model where the axiom to be proven independent will not hold while all the others will.

In 1965, Gupta presented an independent version of Tarski's system of geometry [11].
To obtain the independence of one of his axioms, namely Pasch's axiom, he proved the independence of one of its consequences: the previously eliminated symmetry of betweenness.
However, an independence model for the non-degenerate part of Pasch's axiom was provided by Szczerba for another version of Tarski's system of geometry in which the symmetry of betweenness holds [21].
This independence proof cannot be directly used for Gupta's version as the statements of the parallel postulate differ.
This can be remedied by carefully choosing the statement of the parallel postulate amongst the ones known to be equivalent [7].
We aim to verify that splitting Pasch's axiom into its non-degenerate part and the symmetry of betweenness in addition to changing the statement of the parallel postulate allows to obtain a system that is still independent.
So we go in the opposite direction of what Makarios did by removing the need of the reflexivity properties for congruence thanks to a modification to the five-segment property [17].
Indeed, our view is that an axiom should capture a limited and well-defined property, instead of trying to minimize the needed number of axioms at any cost.

We remark that a small change in the statement of an axiom can change whether or not it holds in a specific model.
This makes a computer very well suited to the verification that an axiom holds in a model.
So we chose to mechanize the various counter-models needed for this task in the Coq proof-assistant [23].

The rest of the paper is structured as follows.
First, in Sec. <ref>, we present the system we will be working on throughout the rest of this paper.
Then, in Sec. <ref>, we show how to build a model of Tarski's axiom.
Finally, before concluding on the achieved results, we present an example of independence proof in Sec. <ref>.

§ A VARIANT OF TARSKI'S SYSTEM OF GEOMETRY

In this section, we start by recalling the axioms of Tarski's system of geometry.
Then, we present the modification Gupta made to it to obtain a fully independent system [11].
Finally, we describe how to modify his system to combine his results and the ones from Szczerba [21].

§.§ Tarski's System of Geometry

Tarski's axiom system is based on a single primitive type depicting points and two predicates, namely congruence and betweenness.
$\congT{AB}{CD}$ states that the segments $\overline{AB}$ and $\overline{CD}$ have the same length.
$\bT{A}{B}{C}$ means that $A$, $B$ and $C$ are collinear and $B$ is between $A$ and $C$ (and $B$ may be equal to $A$ or $C$).
For an explanation of the axioms and their history see [22].
Table <ref> lists the axioms for planar Euclidean geometry.

A1     Symmetry     $\congT{AB}{BA}$

A2     Pseudo-Transitivity     $\congT{AB}{CD} \land \congT{AB}{EF} \imp \congT{CD}{EF}$

A3     Cong Identity     $\congT{AB}{CC} \imp A=B$

A4     Segment construction    $\exists E, \bT{A}{B}{E} \land \congT{BE}{CD}$

A5     Five-segment     $\congT{AB}{A'B'} \land \congT{BC}{B'C'} \land$

$\congT{AD}{A'D'} \land \congT{BD}{B'D'} \land$

$\bT{A}{B}{C} \land \bT{A'}{B'}{C'} \land A \neq B \imp \congT{CD}{C'D'}$

A6     Between Identity     $\bT{A}{B}{A} \imp A=B$

A7     Inner Pasch     $\bT{A}{P}{C} \land \bT{B}{Q}{C} \imp \exists X,\bT{P}{X}{B} \land \bT{Q}{X}{A}$

A8     Lower Dimension    $\exists A B C, \lnot \bT{A}{B}{C} \land \lnot \bT{B}{C}{A} \land \lnot \bT{C}{A}{B}$

A9     Upper Dimension    $\congT{AP}{AQ} \land \congT{BP}{BQ} \land \congT{CP}{CQ} \land P \neq Q \imp$

$\bT{A}{B}{C} \lor \bT{B}{C}{A} \lor \bT{C}{A}{B}$

A10     Euclid     $\bT{A}{D}{T} \land \bT{B}{D}{C} \land A \neq D \imp$

$\exists X Y, \bT{A}{B}{X} \land \bT{A}{C}{Y} \land \bT{X}{T}{Y}$

A11     Continuity     $(\exists A, (\forall X Y, \Xi (X) \land \Upsilon (Y) \imp \bT{A}{X}{Y})) \imp$

$~\exists B, (\forall X Y, \Xi (X) \land \Upsilon (Y) \imp \bT{X}{B}{Y})$

Tarski's axiom system for planar Euclidean geometry.

§.§ Gupta's Contribution

The problem of the independence of Tarski's axiom system, as defined in Table <ref>, remains open.
Let us introduce the modifications Gupta made to it to obtain an independent system.
He reintroduced the inner transitivity of betweenness A15[We number them as in [22].] in Table <ref>.
Having added this axiom, the identity axiom for betweenness A6 became a theorem and could then be removed from the system.
Finally A2, A9 and A11 are replaced by A2', A9' and A11'.
We omit the details of how to mechanize in Coq that this system, consisting of A1, A2', A3-A5, A7, A8, A9', A10, A11' and A15,[Actually the statement for A7 differs in [22] but the change is not important here.] and Tarski's system are equivalent.
Gupta proves that this system is independent.
To prove that A7 is independent in this system, he shows that A14, a consequence of A7 in this system, does not hold.
However, Szczerba [21] found that A7 does not hold when A14 and all the other axioms, with the exception of A10, in Gupta's system do.
So this would suggest that A7 can be split into A14 and a variant of A7 while still having an independent system.

A0     Point equality decidability     $X = Y \lor X \not = Y$

A2'     Pseudo-Transitivity     $\congT{AB}{EF} \land \congT{CD}{EF} \imp \congT{AB}{CD}$

A7'     Inner Pasch     $\bT{A}{P}{C} \land \bT{B}{Q}{C} \land$

$A \neq P \land  P \neq C \land B \neq Q \land Q \neq C \land$

$\lnot \left( \bT{A}{B}{C} \lor \bT{B}{C}{A} \lor \bT{C}{A}{B} \right) \imp$

$\exists X,\bT{P}{X}{B} \land \bT{Q}{X}{A}$

A9'     Upper Dimension    $\congT{AP}{AQ} \land \congT{BP}{BQ} \land \congT{CP}{CQ} \land$

$P \neq Q \land A \neq B \land A \neq C \land B \neq C \imp$

$\bT{A}{B}{C} \lor \bT{B}{C}{A} \lor \bT{C}{A}{B}$

A10'     Proclus     $\para{AB}{CD} \land \col{A}{B}{P} \land \lnot \col{A}{B}{Q} \imp$

$\exists Y, \col{C}{D}{Y} \land \col{P}{Q}{Y}$

A11'     Continuity     $(\exists A, (\forall X Y, \Xi (X) \land \Upsilon (Y) \imp \bT{A}{X}{Y})) \imp$

$~\exists B, (\forall X Y, \Xi (X) \land \Upsilon (Y) \imp$

$~~~~~~~~X = B \lor B = Y \lor \bT{X}{B}{Y})$

A14     Between Symmetry     $\bT{A}{B}{C} \imp \bT{C}{B}{A}$

A15     Between Inner Transitivity     $\bT{A}{B}{D} \land \bT{B}{C}{D} \imp \bT{A}{B}{C}$

Added axioms to Tarski's system of geometry.

§.§ An Independent System for Planar Geometry

The system that we want to prove independent is very close to the one Gupta studied in his thesis [11].
We split Pasch's axiom A7 into its non-degenerate part A7' and A14, change the version of the parallel postulate A10 and add one axiom (for reasons explained later).
A7' excludes from A7 the degenerate cases where the triangle $ABC$ is flat or when $P$ or $Q$ are respectively not strictly between $A$ and $C$ or $B$ and $C$.
We cannot use A10 as it does not hold in the counter-model found by Szczerba [21].
We chose Proclus postulate,[$\col{A}{B}{C}$ and $\para{AB}{CD}$ denotes that $A$, $B$ and $C$ are collinear and that lines $AB$ and $CD$ are parallel according to the definitions given in SST [19].] denoted as A10' in Table <ref>, verified to be equivalent to it when assuming A0-A9, using Coq [7], as it holds in all the counter-models provided by Gupta as well as in the one found by Szczerba, thanks to Theorem 1 in [21].
Finally, the formal development found in SST [19] is essentially classical due to the many case distinctions found in the proofs of its lemmas.
However, the decidability of point equality is sufficient to obtain the arithmetization of geometry in an intuitionistic setting [6].
So we add the decidability of point equality A0 so that we can work in an intuitionistic setting.
The reader not familiar with the difference between classical and intuitionistic logic may refer to [1].
This system, consisting of A0, A1, A2', A3-A5, A7', A8, A9', A10', A11', A14 and A15, and Tarski's system are equivalent.
Again, we do not detail how to mechanize this fact in Coq.

§ A MODEL OF TARSKI'S SYSTEM OF GEOMETRY

In this section, we present our proof that Cartesian planes over a Pythagorean ordered field form a model of the variant of Tarski's system of geometry that we have introduced in the previous section.
First, we present the structure that we used to define this model.
Then we define the model that we used, that is, the way we instantiated the signature of this system.
Finally, we detail the proofs of some of the more interesting axioms.

§.§ The Real Field Structure

The structure that was used to define this model was built by Cohen [9].
The real field structure results of the addition of operators to a discrete[Discrete fields are fields with a decidable equality.] field: two boolean comparison functions (for strict and non-strict order) and a norm operator.
Elements of this real field structure verify the axioms listed in Table <ref>.
Finally, the elements of a real field structure are all comparable to zero.
We should remark that this field is not necessarily Pythagorean.
In fact, there is no defined structure in the Mathematical Components library [15] for Pythagorean fields.
This can however be added much more easily than before thanks to the recent modification of the Mathematical Components library to make use of the Hierarchy Builder [10].
However, the Pythagorean property is only required for the proof of the segment construction axiom A4.
So we chose to prove that this axiom holds in our model by admitting an extra axiom which was defined in this library: the real closed field axiom.
It states that intermediate value property holds for polynomial with coefficients in the field.
While it is much stronger than Pythagoras' axiom, we only used it to be able to define the square root of a number which is a sum of squares and would therefore have a square root in a Pythagorean field.
Finally, we did not yet prove that $A11$ holds in our model since it would require a much more involved effort.
Indeed, this is similar to verifying that Tarski's system of geometry admits a quantifier elimination procedure.

Subadditivity of the norm operator     $| x + y | \le |x| + |y|$

Compatibility of the addition with the strict comparison     $0 < x \land 0 < y \imp 0 < x + y$

Definiteness of the norm operator     $| x | = 0 \imp x = 0$

Comparability of positive numbers     $ 0 \le x \land 0 \le y \imp (x \le y) || (y \le x)$

The norm operator is a morphism for the multiplication     $| x * y | = |x| * |y|$

Large comparison in terms of the norm     $(x \le y) = (| y - x | == y - x)$

Strict comparison in terms of the large comparison     $(x < y) = (y ~ != x) \&\& (x \le y)$

Axioms of the real field structure.

$==$ denotes the boolean equality test for the elements of the field.

§.§ The Model

Let us now define our model.
Being based on a single primitive type and two predicates, the signature of Tarski's system of geometry is rather simple.
However, this system has the advantage of having a $n$-dimensional variant.
To obtain this variant, one only needs to change the dimension axioms.
So far, we have restricted ourselves to the planar version of this system.
With a view to extend the GeoCoq library to its $n$-dimensional variant, we wanted to define a model in which we could prove all but the dimension axioms in an arbitrary dimension to be able to construct a model of the $n$-dimensional variant by only proving the new dimension axioms.
Hence we chose to define  as a vector of dimension $n+1$ with coefficient in the real field structure $\mathbb{F}$ (we used the real field structure for all the development with the exception of the proof of the segment construction axiom) for a fixed integer $n$, that is .
We adopted Gupta's definition [11] for the congruence , namely that $\congT{AB}{CD}$ if the squares of the Euclidean norms of $B-A$ and $D-C$ are equal.
Actually Gupta also proved that any model of the $n$-dimensional variant of Tarski's system of geometry is isomorphic to his model.
He defined that $\bT{A}{B}{C}$ holds if and only if there exists a $k \in \mathbb{F}$ such that $0 \le k \le 1$ and $B-A = k (C-A)$.
In fact, if such a $k$ exists, it can be computed.
By letting $A = \left( a_i \right)_{1 \le i \le n+1}$, $B = \left( b_i \right)_{1 \le i \le n+1}$ and $C = \left( c_i \right)_{1 \le i \le n+1}$, if $A \neq C$ then there exists a $i \in \mathbb{N}$ such that $1 \le i \le n+1$ and $a_i \neq c_i$ and in this case we set $k$ to $\frac{b_i - a_i}{c_i - a_i}$ and if $A = C$ we set $k$ to zero.
Therefore we defined a function  that computes the possible value for $k$, thus allowing us to define the betweenness by the boolean equality test.
This was actually important as it permitted to directly manipulate the definition for betweenness by rewriting since we defined it as a boolean test.
Finally, as it was often necessary to distinguish whether or not $\bT{A}{B}{C}$ holds due to a degeneracy, we split the definition  of the betweenness into two predicates: the first one, , capturing the general case of $k$ being strictly between $0$ and $1$ and the second one, , capturing the three possible degenerate cases, namely either $A = B$, $B = C$ or $A = B$ and $B = C$.

Formally, we consider the following model:

Variable R : realFieldType.
Variable n : nat.

Implicit Types (a b c d : 'rV[R]_(n.+1)).

Definition cong a b c d := (b - a) *m (b - a)^T == (d - c) *m (d - c)^T.

Definition betE a b c := [ || [        a == b     b == c ], a == b | b == c ].

Definition ratio v1 v2 :=
if [pick k : 'I_(n.+1) | v2 0 k != 0] is Some k
then v1 0 k / v2 0 k else 0.

Definition betR a b c := ratio (b - a) (c - a).

Definition betS a b c (r := betR a b c) :=
[        b - a == r *: (c - a), 0 < r     r < 1].

Definition bet a b c := betE a b c || betS a b c.

§.§ Proof that the Axioms hold in the Model

Now that we have defined the model, we focus on the proof that the axioms of the system from Sec. <ref> hold in this model.
However, we omit the details of the proofs for axioms A1, A2', A3 and A14 since they are rather straightforward.
For the same reason, we do not cover the decidability of point equality A0.

Let us start by focusing on axioms A7' and A15 as the proofs that they hold in our model are quite similar.
In the case of axiom A15 we know that $\bT{A}{B}{D}$ and $\bT{B}{C}{D}$ so let $k_1 \in \mathbb{F}$ be such that $0 < k_1 < 1$ and $B - A = k_1 (D - A)$ (the degenerate case of this axiom is trivial so we only consider the general case) and $k_2 \in \mathbb{F}$ be such that $0 < k_2 < 1$ and $C - B = k_2 (D - B)$.
In order to prove that $\bT{A}{B}{C}$ we need to find a $k \in \mathbb{F}$ such that $0 < k < 1$ and $B - A = k (C - A)$.
By calculation we find that $k = \frac{k_1}{k_1 + k_2 - k_1 k_2}$ and we can verify that $0 < k < 1$.
In a similar way, for axiom A7', we know that $\bT{A}{P}{C}$ and $\bT{B}{Q}{C}$ so let $k_1 \in \mathbb{F}$ be such that $0 < k_1 < 1$ and $P - A = k_1 (C - A)$ (the hypotheses imply that $0 < k_1 < 1$ because $A \neq P$ and $P \neq C$) and $k_2 \in \mathbb{F}$ be such that $0 < k_2 < 1$ and $Q - B = k_2 (C - B)$.
In order to prove that there exists a point $X$ such that $\bT{P}{X}{B}$ and $\bT{Q}{X}{A}$ we need to find a $k_3 \in \mathbb{F}$ and a $k_4 \in \mathbb{F}$ such that $0 < k_3 < 1$, $0 < k_4 < 1$ and $k_3 (B - P) + P = k_4 (A - Q) + Q$.
By calculation we find that $k_3 = \frac{k_1 (1 - k_2)}{k_1 + k_2 - k_1 k_2}$ and $k_4 = \frac{k_2 (1 - k_1)}{k_1 + k_2 - k_1 k_2}$ and we can verify that $0 < k_3 < 1$ and $0 < k_4 < 1$.
In both of these proof, the ratios are almost identical to the point that it suffices to prove the following lemma:

Lemma ratio_bet a b c k1 k2 k3 :
0 < k1 -> 0 < k2 -> k1 < 1 -> 0 < k3 -> k3 < k1+k2-k1*k2 ->
b - a == ((k1+k2-k1*k2)/k3)^-1 *: (c - a) -> bet a b c.

It allows to prove quite easily both of these axioms.
For axiom A4, we proceeded in a analogous way: it suffices to set the point $E$ that can be constructed using this axiom to $\frac{\| D - C \|}{\| B - A \|} (B - A) + A$ and to verify this point satisfies the desired properties by calculation.

We now turn to axiom A5.
We followed Makarios' approach for the proof that this axiom holds in our model [16].
In his proof he used the cosine rule: in a triangle whose vertices are the vectors $A$, $B$ and $C$ we have

$$\| C - B \|^2 = \| C - A \|^2 + \| B - A \|^2 - 2 (B - A) \cdot (C - A).$$

As noted by Makarios, using the cosine rule allows to avoid defining angles and properties about them.
Applying the cosine rule for the triangles $BCD$ and $B'C'D'$ allows to prove that $\| D - C \|^2 = \| D' - C' \|^2$ by showing that

$$(C - B) \cdot (D - B) = (C' - B') \cdot (D' - B')$$

which can be justified, by applying the cosine rule again, this time in the triangles $ABD$ and $A'B'D'$, if

$$\| D - A \| - \| D - B \| - \| A - B \| = \| D' - A' \| - \| D' - B' \| - \| A' - B' \|$$

which we know from the hypotheses and if the ratios corresponding to the betweenness $\bT{A}{B}{C}$ and $\bT{A'}{B'}{C'}$ are equal which can be obtained by calculation.

Next, let us consider axiom A10.[
Similarly to A7, when we were proving Euclid's axiom, we realized that the same kind of distinctions was also needed.
The degenerate cases are implied by the other betweenness axioms so it suffices to show that A10 holds when the angle $\angle BAC$ is non-flat and when $D$ is different from $T$.]
From the hypotheses we have two ratios $k_1 \in \mathbb{F}$ and $k_2 \in \mathbb{F}$ such that $0 < k_1 < 1$, $0 < k_2 < 1$, $D - A = k_1 (T - A)$ and $D - B = k_2 (C - B)$.
Using these ratios, it suffices to define $X$ such that $B - A = k_1 (X - A)$ and $Y$ such that $C - A = k_1 (Y - A)$.
So we know by construction that $\bT{A}{B}{X}$ and $\bT{A}{C}{Y}$ and we easily get that $T - X = k_2 (Y - X)$ by calculation, thus proving that $\bT{X}{T}{Y}$.
Since A10 and A10' are equivalent when A0, A1, A2', A3-A5, A7', A8, A9', A11', A14 and A15 hold, this allows to prove that A10' holds in our model.

Finally the remaining two axioms are treated in a slightly different setting since they are the dimension axioms.
Formally we fix the value of $n$ to $1$.
In order to simplify the many rewriting steps needed for these proofs we started by establishing the following two lemmas:

Definition sqr_L2_norm_2D a b :=
(b 0 0 - a 0 0) ^+ 2 + (b 0 1 - a 0 1) ^+ 2.

Lemma congP a b c d :
reflect (sqr_L2_norm_2D a b = sqr_L2_norm_2D c d) (cong a b c d).

Lemma betSP' a b c (r := betR a b c) :
reflect ([ / b 0 0 - a 0 0 = r * (c 0 0 - a 0 0),
b 0 1 - a 0 1 = r * (c 0 1 - a 0 1), 0 < r     r < 1])
(betS a b c).

The reader familiar with SSReflect will have recognized the  predicate, described in [9] for example.
In practice, these lemmas allowed to spare many steps that would have been repeated in almost every proof concerning the dimension axioms.
It was much more straightforward to prove that axiom A8 holds in our model than for axiom A9'.
In fact, it is enough to find three non-collinear points.
We simply took the points $(0, 0)$, $(0, 1)$ and $(1, 0)$:

Definition row2 R : ringType (a b : R) : 'rV[R]_2 :=
_p [eta  with 0 |-> a, 1 |-> b] p.

Definition a : 'rV[R]_(2) := row2 0 0.
Definition b : 'rV[R]_(2) := row2 0 1.
Definition c : 'rV[R]_(2) := row2 1 0.

It was then an easy matter to verify that axiom A8 holds in our model.
For axiom A9, the idea of the proof that we formalized was to first show that, by letting $M$ be the midpoint of $P$ and $Q$, the equation $(x_P - x_M) (x_M - x_X) + (y_P - y_M) (y_M - y_X) = 0$, capturing the property that the points $P$, $M$, and $X$ form a right angle with the right angle at vertex $M$, was verified when $X$ would be equal to $A$, $B$ or $C$:

Lemma cong_perp (a p q : 'rV[R]_(2)) (m := (1 / (1 + 1)) *: (p + q)) :
cong a p a q ->
(p 0 0 - m 0 0) * (m 0 0 - a 0 0) +
(p 0 1 - m 0 1) * (m 0 1 - a 0 1) = 0.

Next, we demonstrated that for three points $A$, $B$ and $C$ verifying $(x_A - x_B) (y_B - y_C) - (y_A - y_B) (x_B - x_C) = 0$ are collinear in the sense that $\bT{A}{B}{C} \lor \bT{B}{C}{A} \lor \bT{C}{A}{B}$:

Lemma col_2D a b c :
(a 0 0 - b 0 0) * (b 0 1 - c 0 1) ==
(a 0 1 - b 0 1) * (b 0 0 - c 0 0) ->
(bet a b c  bet b c a  bet c a b).

Using the equations implied by  we could derive that

$$(x_P - x_M) (y_M - y_P) \left( (x_A - x_B) (y_B - y_C) - (y_A - y_B) (x_B - x_C) \right) = 0.$$

We were then left with three cases: either the abscissas of $P$ and $M$ are equal in which case the ordinate of $A$, $B$ and $C$ were equal thus sufficing to complete the proof, or the ordinates of $P$ and $M$ are equal in which case the abscissas of $A$, $B$ and $C$ were equal thus completing the proof, or $(x_A - x_B) (y_B - y_C) - (y_A - y_B) (x_B - x_C) = 0$ corresponding to the lemma that we had proved and again allowing to conclude.

Putting everything together, we could prove that Cartesian planes over a Pythagorean ordered field form a model of the variant of Tarski's system of geometry, thus proving the satisfiability of the theory.[ is the type class that captures the theory consisting of axioms A0-A10.]

Global Instance Rcf_to_T2D : Tarski_2D Rcf_to_T_PED.

Global Instance Rcf_to_T_euclidean : Tarski_euclidean Rcf_to_T_PED.

§ AN EXAMPLE OF INDEPENDENCE PROOF

To illustrate how we obtain formal proofs of independence we present an example.
We start by defining the counter-model we will use to prove the independence of axiom A10'.
We then provide the sketch of the formal proof.

§.§ Klein's Model

To prove Euclid’s Parallel Postulate independent from the other axiom we work in Klein's model as defined in SST [19]:

Variable R : realFieldType.
Variable n : nat.

Definition Vector := 'rV[R]_(n.+1).
Definition Point : Type := p : Vector | (p *m p^T) 0 0 < 1.
Notation "#" := proj1_sig.

Implicit Types (a b c d : Point).
Implicit Types (v w x y : Vector).

Definition bet' a b c := bet (#a) (#b) (#c).

Definition omd_v v w := (1 - (v *m (w)^T) 0 0).
Definition cong_v v w x y :=
(omd_v v w)^+2/(omd_v v v * omd_v w w) ==
(omd_v x y)^+2/(omd_v x x * omd_v y y).
Definition cong' a b c d := cong_v (#a) (#b) (#c) (#d).

Here,  is the type of , vectors of dimension $n+1$ with coefficient in the real field structure, lying inside the unit disk and  the projection allowing to recover the coordinate part of this dependent type.
In Klein's model, $b$ is said to be between $a$ and $c$ iff their coordinate parts can be said to be  in the model from Sec. <ref> and line-segments $\overline{ab}$ and $\overline{cd}$ are said to be congruent iff

$$\frac{\left( 1 - \# a \cdot \# b \right)^2}{\left( 1 - \# a \cdot \# a \right) \left( 1 - \# b \cdot \# b \right)} = \frac{\left( 1 - \# c \cdot \# d \right)^2}{\left( 1 - \# c \cdot \# c \right) \left( 1 - \# d \cdot \# d \right)}$$

where $\cdot$ denotes the dot product of two vectors.

§.§ Independence of Euclid’s Parallel Postulate via Klein's Model

Here we only detail the proof that A10' does not hold in this model.
Mechanizing the following proof sketch allows to derive.

Lemma euclid :   euclidP (@Point R 1) (@bet' R 1).

To make sure that we did not introduce any change in the axioms between the various models we relied on predicates such as , which depend on possibly the type for points and the predicate(s) for betweenness and/or congruence.

Axiom A10' does not hold in Klein's model.

Since Klein's model forms a model of neutral geometry,[Neutral geometry is defined by the set of axioms of Euclidean geometry from which the parallel postulate has been removed.] it suffices to prove that any version of the parallel postulate, proven equivalent to A10' in Coq when assuming A0-A9, does not hold.
We choose to work with A10.
Picking $a$, $b$, $c$, $d$ and $t$ to be of coordinates $(0, 0)$, $(0, \frac{1}{2})$, $(\frac{1}{2}, 0)$, $(\frac{1}{4}, \frac{1}{4})$ and $(\frac{1}{2}, \frac{1}{2})$, some computations allow to verify that $\bT{\# a}{\# d}{\# t}$, $\bT{\# b}{\# d}{\# c}$ $b \neq d$, $d \neq c$ and $\lnot \col{\# a}{\# b}{\# c}$.
So, to prove that this version does not hold, it is enough to show that for any $x$ and $y$ such that $x$ lies inside the unit disk, $\bT{\# a}{\# b}{\# x}$, $\bT{\# a}{\# c}{\# y}$ and $\bT{\# x}{\# t}{\# y}$, it holds that $y$ is not a , meaning that it lies outside the unit disk.
Let us first eliminate the case where $b = x$ as it would lead to a contradiction.
Here, we use the algebraic characterization of collinearity[Here we use the converse of .] to obtain that, if $b = x$, the ordinate of x would need to be equal to both $0$ and $\frac{1}{2}$ which is impossible.
Now let us pose $b'$ to be the vector $x + a - b$.
It is an easy matter to check that $\bT{\# a}{\# b'}{\# x}$ so let us pose $k_1$ to be the ratio associated to this betweenness.
We can verify that $k_1 \le \frac{1}{2}$ since $x$ is supposed to belong to the unit disk.
We can then take $d'$ at ratio $k_1$ from $a$ to $t$.
Applying what was proven to show that A10' holds in Cartesian planes over a Pythagorean ordered field, we can show that $y'$ at ratio $\frac{1}{k_1}$ from $a$ to $c$ is such that $\bT{\# a}{\# c}{\# y'}$ and $\bT{\# x}{\# t}{\# y'}$.
If we can prove that $y = y'$ we will be done as $y'$ lies outside of the unit disk because $k_1 \le \frac{1}{2}$ so $2 \le \frac{1}{k_1}$.
Finally, to prove that $y = y'$ we can reason by uniqueness of the intersection of lines which is valid in neutral geometry.

§ CONCLUSION

We defined ten out of the eleven counter-models present in Gupta's thesis [11], thus obtaining the Coq formal proof of the independence of ten out of the thirteen axioms of the system presented in Sec. <ref>.
This seems to indicate that Pasch's axiom could indeed be split into two meaningfully different parts as done in this paper while still having an independent system.
However, we will only be sure of this once we will have formalized the missing three counter-models.
These can be found in Gupta's thesis [11], Szczerba's paper [21], and Beeson's section The recursive model in [2].

Five of the formalized models are finite and the other five are modifications of the model presented in Sec. <ref>.
We highlight that, for the latter five, A11' is not verified for the same reason as for the model from Sec. <ref>..
All these models are available in the GeoCoq library[<http://geocoq.github.io/GeoCoq/>] and represent about 4k lines of formal proof.

We are currently extending this work by proving the independence of a more constructive version[We replace point equality decidability by point equality stability, namely $\forall X Y, \lnot \lnot X = Y \imp X = Y$, which allows to prove equality of points by contradiction but does not allow case distinctions. We do not go as far as in [3] where not even stability is assumed. We also apply the same modifications made to obtain what is called continuous Tarski geometry in [1].] of the axioms which would also allow to capture $n$-dimensional geometry.
For this extension we could not rely on A9$^{(n)}$ from [22].
Indeed, we found that it can only be assumed as an upper $n$-dimensional axiom when $n =2$ or $3$.
A9$^{(n)}$ is stated as follows.

$$ \bigwedge\limits_{1 \le i \le j \le n} P_i \neq P_j \land \bigwedge\limits_{i = 2}^n \congT{A P_1}{A P_i} \land \bigwedge\limits_{i = 2}^n \congT{B P_1}{B P_i} \land \bigwedge\limits_{i = 2}^n \congT{C P_1}{C P_i} \imp \col{A}{B}{C}$$

By taking $P_i = (\cos \frac{2i \pi}{n}, \sin \frac{2i \pi}{n}, 0, \ldots, 0)$ for $1 \le i \le n$ then $(0, 0, x_3, x_4, \ldots, x_n)$ satisfies the premises for any $x_3, x_4, \ldots, x_n$ in the standard n-dimensional model while triplets of points of this form are not necessarily collinear.
The various modifications did not allow to reuse some of the counter-models already mechanized, so new ones are necessary.

We are convinced that using a proof-assistant is crucial when proving the independence of a system, where small changes in a statement are critical.
Actually, there was a typo in Gupta's counter-model for A2 and we just exhibited a problem with axiom A9$^{(n)}$ from [22].
The GeoCoq library also proved very useful as it allowed us to combine the algebraic and geometric[Such as the use of the uniqueness of the intersection in our proof of Theorem <ref>.] reasoning.

Acknowledgments: We would like to thank Marius Hinge for his contribution to the early stage of this work.

[1]
Beeson, M.: A Constructive Version of Tarski's Geometry. Annals of Pure and
Applied Logic  166(11),  1199–1273 (2015).

[2]
Beeson, M.: Constructive Geometry and the Parallel Postulate. Bulletin of
Symbolic Logic  22(1),  1–104 (2016). 10.1017/bsl.2015.41

[3]
Beeson, M.: Brouwer and Euclid. Indagationes Mathematicae  29(1),
483–533 (2018). 10.1016/j.indag.2017.06.002

[4]
Beeson, M., Boutry, P., Narboux, J.: Herbrand’s theorem and non-Euclidean
geometry. The Bulletin of Symbolic Logic  21(2),  111–122 (2015).

[5]
Bolyai, J.: Appendix, Scientiam Spatii absolute veram exhibens: a veritate aut
falsitate Axiomatis XI. Euclidei (a priori haud unquam decidenda)
independentem; adjecta ad casum falsitatis, quadratura circuli geometrica.
Auctore Johanne Bolyai de eadem, Geometrarum in Exercitu Caesareo Regio
Austriaco Castrensium Capitaneo. Coll. Ref. (1832)

[6]
Boutry, P., Braun, G., Narboux, J.: Formalization of the Arithmetization of
Euclidean Plane Geometry and Applications. Journal of Symbolic Computation
(2018). 10.1016/j.jsc.2018.04.007

[7]
Boutry, P., Gries, C., Narboux, J., Schreck, P.: Parallel Postulates and
Continuity Axioms: A Mechanized Study in Intuitionistic Logic Using Coq.
Journal of Automated Reasoning  (2017),

[8]
Braun, G., Boutry, P., Narboux, J.: From Hilbert to Tarski. In: Narboux, J.,
Schreck, P., Streinu, I. (eds.) Proceedings of the Eleventh International
Workshop on Automated Deduction in Geometry. pp. 78–96 (2016),

[9]
Cohen, C.: Formalized algebraic numbers: construction and first-order theory.
Theses, Ecole Polytechnique X (Nov 2012),

[10]
Cohen, C., Sakaguchi, K., Tassi, E.: Hierarchy Builder: Algebraic hierarchies
Made Easy in Coq with Elpi (System Description). In: 5th International
Conference on Formal Structures for Computation and Deduction (FSCD 2020).
vol. 167, pp. 34:1–34:21. Paris, France (Jun 2020).

[11]
Gupta, H.N.: Contributions to the Axiomatic Foundations of Geometry. Ph.D.
thesis, University of California, Berkeley (1965)

[12]
Hilbert, D.: Les fondements de la géométrie. Dunod, Jacques Gabay
edn. (1971), edition critique avec introduction et compléments
préparée par Paul Rossier

[13]
Klugel, G.S.: Conatuum praecipuorum theoriam parallelarum demonstrandi
recensio. Ph.D. thesis, Schultz, Göttingen (1763), German translation
available: <http://www2.math.uni-wuppertal.de/ volkert/versuch.html>

[14]
Lobatschewsky, N.: Geometrische Untersuchungen zur Theorie der
Parallellinien, pp. 159–223. Springer Vienna, Vienna (1985).

[15]
Mahboubi, A., Tassi, E.: Mathematical Components. Zenodo (Sep 2022),

[16]
Makarios, T.J.M.: A Mechanical Verification of the Independence of Tarski's
Euclidean Axiom. Master's thesis, Victoria University of Wellington (2012)

[17]
Makarios, T.J.M.: A further simplification of Tarski's axioms of geometry.
arXiv: Logic  (2013). 10.1285/i15900932v33n2p123.

[18]
Narboux, J., Janicic, P., Fleuriot, J.: Computer-assisted Theorem Proving in
Synthetic Geometry. In: Sitharam, M., John, A.S., Sidman, J. (eds.)
Handbook of Geometric Constraint Systems Principles. Discrete Mathematics
and Its Applications, Chapman and Hall/CRC  (2018).
10.1201/9781315121116-2, <https://inria.hal.science/hal-01779452>

[19]
Schwabhäuser, W., Szmielew, W., Tarski, A.: Metamathematische Methoden in
der Geometrie. Springer-Verlag (1983)

[20]
Simić, D., Marić, F., Boutry, P.: Formalization of the Poincaré
Disc Model of Hyperbolic Geometry. Journal of Automated Reasoning
65(1),  31–73 (Apr 2020). 10.1007/s10817-020-09551-2

[21]
Szczerba, L.W.: Independence of Pasch's axiom. Bull. Acad. Polon. Sci. Sér.
Sci. Math. Astronom. Phys.  18,  491–498 (1970)

[22]
Tarski, A., Givant, S.: Tarski's System of Geometry. The Bulletin of Symbolic
Logic  5(2),  175–214 (1999). 10.2307/421089

[23]
Team, T.C.D.: The Coq Proof Assistant, version 8.15 (2022),


# The N-Gene Model for Evolutionary Interaction

Tyler Clark11footnotemark: 1<EMAIL_ADDRESS>
###### Abstract

The evolution of altruistic behavior, which confers benefits to others at a
cost to the individual, remains a fundamental question in evolutionary
biology. While previous models have investigated the conditions favoring the
emergence and stability of altruism, they have often employed simplifying
assumptions, such as single-gene inheritance and discrete strategies. In this
study, we develop a novel evolutionary model that incorporates Mendelian
genetics, continuous strategies, and the potential for multiple genes to
contribute to a single phenotypic trait. We employ a modified dictator game as
the framework for evolutionary interactions and explore the stability of
altruistic behavior under various conditions. Our primary result demonstrates
that when considering heterozygous genes, altruism can be evolutionarily
stable at cost-to-benefit ratios exceeding unity, even with initially low
frequencies of altruists in the population. This finding contrasts with the
case of homozygous genes, where altruism is only stable at cost-to-benefit
ratios greater than 2. We further illustrate that assortative matching, noise
in the cost-to-benefit ratio, and the presence of ”taking” behavior can
influence the evolutionary dynamics of altruism within our model. The
generality of our approach allows for its application to a diverse range of
evolutionary games and interactions, providing a powerful tool for
investigating the emergence and maintenance of social behaviors and
personality traits. Our results contribute to the understanding of the
evolutionary mechanisms underlying altruism and underscore the importance of
incorporating genetic complexity in evolutionary models. This work has
implications for the study of social evolution and the genetic architecture of
complex behavioral phenotypes.

###### keywords:

Evolutionary game theory, Altruism, Dictator Game

††journal: Journal of Theoretical Biology

##  1 Introduction

Up until the early twentieth century, Darwin’s notion of natural selection was
synonymous with ’survival of the fittest’: to eat, or be eaten. In the popular
imagination at least, this association, or conflation, still persists today.
We tend to think of evolution in competitive, selfish terms.

In the later half of the twentieth century, however, evolutionary theories
also began to explain non-selfish behaviours by introducing the terms
”inclusive fitness” [1] and ”reciprocal altruism” [2]. In this way, these
studies engendered a new way of looking at evolution, one in which survival
depended on cooperation instead of precluding it. Work in the late twentieth
century also brought into question how much one’s genes even impact traits. To
what extent are we a tabula rasa waiting to be written on?

With the wide spread use of twin studies, the nature vs. nurture debate has
been studied scientifically, and the findings are that, on average, genetics
and upbringing are equally significant. This is also roughly the case for
altruism and self-interest [3][4][5], although the exact proportions are not
necessary for my purposes here; it is simply notable that a significant
percentage of the variance is explained via genetics.

One approach to determining the stability of evolutionary strategies is the
field of evolutionary game theory. This field relies on the concept of an
evolutionarily stable strategy (ESS) [6]. An ESS is a modified version of a
strict nash equilibrium [7]. In brief, this a strategy that cannot be improved
after considering the possible strategical space of the other players.

ESSs are motivated entirely differently than those of traditional game theory.
It is presumed that these strategies are biologically encoded and heritable.
Individuals have no conscious intention in their strategy and are not aware of
the game. They reproduce and are subject to the forces of natural selection,
with the payoffs of the game representing reproductive biological fitness.
Alternative game strategies occasionally occur via a process like mutation. To
be an ESS, a strategy must be resistant to these alternatives.

One topic rarely discussed is that there are possible disparities between an
ESS within an individual’s lifetime and the ESS of a gene over many
generations. Often, an analytic solution is found for the ESS of one lifetime,
and assumed to be interchangeable with an intergenerational ESS.

This disparity has been noted in parent-offspring conflict [8]. In a 1971
paper, for example, Trivers discusses how the optimal parental investment (PI)
is different for parents and for offspring. Offspring will want more PI then
their siblings, but parents will want to allocate evenly.

The godfather of evolutionary altruism research, W.D. Hamilton, has set the
benchmark for research on evolutionary cooperation over the past half a
century. He found that within a single lifetime the optimal strategy for
altruistic giving is based on the inequality $rb>c$ [1], where r is wright’s
coefficient of relatedness, b is the net fitness benefit to the recipient, and
c is the cost to the donor.

More recently, Alger and Weibull (AW) have investigated this topic at the
level of traits and within a modified ESS framework [9]. However, they again
look for an ESS within a single lifetime. Both Hamilton’s and AW’s work only
looks at genes that have one allele per loci, also known as homozygeous genes.

In this study I consider a model with a class of pairwise interactions that
can have gradations of strategy, Mendelian genetics, 222Mendelian genetics
refers to meiotic division obeying Mendel’s laws and the possibility of many
genes being correlated to a single trait. This allows me to model
distributions of traits within a population. I am able to reproduce the
results of Hamilton and AW when using homozygeous genes, but when looking at
heterozygeous genes I find that the optimal strategy entails a cost-to-benifit
ratio of $>1$ for full siblings, counter to their findings.

To explain this intuitively, consider a situation where an individual has a
chance to increase their parent’s expected number of offspring rather then
maximize their own (Helping their brother at a net loss to themselves, for
example). If the cost-to-benefit ratio is $>2$; then both AW and Hamilton
would agree that this is not a stable strategy. However, from the parent’s
point of view, it is optimal for their offspring to help one another at a
cost-to-benefit $\alpha>1$. When an individual has offspring, the initially
negative altruism cost pays off through the generations as they enjoy a higher
expected number of children then their short-term selfish counterparts.

###  1.1 Model Introduction

Consider a population of sexually reproducing individuals who are matched into
pairs. Each pair is engaged in some interaction. Each individual carries some
heritable trait $\theta\in T$, where T is a set of potential traits. Here we
use the word ’trait’ in a general sense, allowing it to mean anything that may
influence the outcome of the interaction, such as intelligence, size, decision
rule, etc.

This is the game, and after it is played the individuals survive until
replication with a certain probability which is impacted during the game.

####  1.1.1 Fitness

Many models deal with the fitness associated with a given trait in an abstract
sense, as being the expected number of offspring an individual creates. If we
assume that an individual’s life can be broken into two parts, immaturity and
maturity, once the individual lives until maturity, they will give birth
without possibility of death.

When looked at in a probabilistic frame, we see that changes in genotype
fitness without evolutionary forces/interactions can either act on the number
of offspring ($n$) an individual has or the chances the individual has to live
until reproduction/maturity ($p$).

######  Definition 1.1.

Mature Fitness is the number of offspring (n) an individual about to reproduce
has that live until reproduction without evolutionary forces, where each
$i^{th}$ offspring has a probability ($p_{i}$) of living until this point.

$\omega=\sum_{i=1}^{n}p_{i}$ (1)

######  Definition 1.2.

Immature Fitness is the expected number of offspring an individual that lives
until reproductive age has, taking into account that the individual has not
yet reached reproductive age, and will do so with a probability (p).

$\omega^{\prime}=p\omega$ (2)

The above formula for mature fitness is the expectation of the binomial
poisson distribution, which in the limiting case that all offspring have an
equal probability of living until replication reduces to $\omega=np$.

######  Definition 1.3.

Evolutionary Interaction is the immature fitness resulting from the
interaction of an individual with trait ($\theta$) with another individual
with trait ($\theta^{\prime}$) .

$I(\theta,\theta^{\prime})=(p+\delta p)(\omega+\delta\omega)$ (3)

Because trait inheritance is a distribution, not deterministic, in this model
it is easiest to specify a constant mature fitness $\omega_{c}$ which all
individuals in the gene pool possess.

We can define the change in probability of living until replication after the
interaction as

$f(\theta,\theta^{\prime})=p+\delta p$ (4)

Figure 1: This figure shows how two first generation parents undergo meiotic
reproduction in order to produce two offspring. The offspring then play an
evolutionary game which changes their chance of living until reproduction. The
second generation then undergoes meiotic reproduction to create the third
generation.

####  1.1.2 Traits as Macrostates

Let each individual have a set of n genes that are correlated to the trait
$\theta$. Let the trait property be a function of the individual’s total
number of correlated genes, which we will call the genomic macrostate. This
can be contrasted with the individual’s specific set of genes, or the genomic
microstate.

######  Definition 1.4.

Trait $\theta(k)$ is the phenotypic expression of the genomic macrostate $k$
of the individual; which effects the strategy of an evolutionary interaction.

####  1.1.3 Modified Dictator Game

Let us consider an asymmetric two player game that is a modified version of
the dictator game. Let Player 1 be the dictator who has the ability to give
Player 2 a part of their total fitness with a benefit-to-cost ratio of
$\alpha$. The strategy set for the game would be

$(p,p)\mapsto(f(\theta^{\prime},\theta),f(\theta,\theta^{\prime}))=(p-\delta
p,p+\delta p\alpha)$

Where $\delta p\alpha\leq(1-p)$ and $\delta p<p$

We can define the mapping of the genomic macrostate to the trait as

$\theta(\Omega)=a\Omega+b$ (5)

Where $a$ specifies the magnitude of impact the interaction has on fitness and
$b$ changes the minimum amount of giving required.

###  1.2 Assumptions of the Model

  1. 1.

Let there be $n$ genes

  2. 2.

Let the gene-pool for each gene have two competing alleles. One allele has a
positive influence on the trait, and the other is normalized to have no
effect.

  3. 3.

Let each gene have two alleles, giving rise to the possibility of
heterozygosity.

  4. 4.

Let the frequency of the correlated allele in each of the $n$ genes be equal.

  5. 5.

Mendelian inheritance applies.

  6. 6.

Let the correlated gene be dominant with a probability $p$.

  7. 7.

Let each gene be equally correlated with the trait.

In this application let the genomic microstate of an individual
$s=[g_{1},g_{2}...g_{n}]$ be the specific set of alleles that an individual
possesses. Due to the fact that all genes have an identical effect on the
trait, the relevant factor is not which genes have the correlated allele, but
the total number of correlated alleles.

Let $\theta(\Omega)$ be the trait as a function of the individual’s macrostate
$\Omega=\sum^{n}_{i=1}g_{i}$ where $g_{i}=(1,0)$ is a set indicating whether
the gene is correlated to the trait or not. If $g_{i}=(1,1)$ the trait is
correlated, if $g_{i}=(0,0)$ uncorrelated, and if $g_{i}=(1,0)$ then the trait
is correlated with a probability $p$, which is the chance of being dominant.

Instead of calculating the macrostate directly, one can also do this with
knowledge of the frequencies of the homozygeous correlated alleles
$\rho_{BB}$, the heterzygeous set of alleles $\rho_{BW}$, and the homozygeous
uncorrelated set of alleles $\rho_{WW}$. Because the genes are identical, the
previous frequencies are equal for each gene.

###  1.3 Mating

Figure 1 shows the intuition behind the process of mating in this model. Given
parents of macrostate $i,j$, and assuming Mendel’s law of segregation holds,
we can calculate the macrostate of the gamete $\Omega_{G}$ that contains only
a single allele per gene.

Using the equations shown below, one can derive probability or frequency
estimate of how many individuals will have a certain macrostate, and therefor
how many parents of macrostates $i,j$ will mate together and the probability
of their child having macrostate $k_{s}$

####  1.3.1 Heterozygeous populations

We can create a tensor

$\Gamma=\sum_{i,j,k...}P(\Omega_{G}=i|\rho_{BB},\rho_{BW})P(\Omega_{G}=j|\rho_{BB},\rho_{BW})\prod_{s=1}^{N}P(\Omega=k_{s}|i,j)$
(6)

Which specifies the frequency of N offspring with ordered macrostates
$k_{1,2..N}$, given parents with gamete macrostates i, j.

First consider the last term, $P(\Omega_{G}=k_{s}|i,j)$. This is the
probability of having an offspring with macrostate $k_{s}$ given parents of
gamete macrostates i, j.

We find that when the correlated gene is completely recessive:

$P(\Omega=k_{s}|i,j)_{recess}={n-j\choose i-k_{s}}{j\choose
k_{s}}\frac{k_{s}!(n-i)!}{n!}\ \ \forall i>j$ (7)

Where n is the number of correlated genes. We can see that if $j>i$ the
equation is the same with j and i exchanged.

####  1.3.2 Dominance

When the correlated gene is dominant with a certain probability $p$ we arrive
at the following equation:

$P(\Omega=k_{s}|i,j,p)=\sum_{v=0}^{k_{s}}P(\Omega=k_{s}-v|i,j)_{recess}{j+i-2k_{s}+2v\choose
v}p^{v}(1-p)^{i+j-2k_{s}+2v}$ (8)

$\forall$
a) $k_{s}-v\leq j$
b) $i+j-2k_{s}+2v\geq v$
c) $n\geq i+j-2k_{s}+2v$
d) $k_{s}-v\geq i+j-n$

If each individual has n correlated genes; and each gene loci has the
possibility of being either heterozygeous or homozygous333In the case that the
individual’s gene is homozygeous, it can either be correlated or uncorrelated,
the expected gamete macrostate from choosing one allele from each gene in
which to reproduce can be expressed using the binomial formula below.

$Pr(\Omega_{G}=k|\rho_{BB},\rho_{BW})={n\choose
k}(\rho_{BB}+\frac{\rho_{BW}}{2})^{k}(1-\rho_{BB}-\frac{\rho_{BW}}{2})^{n-k}$
(9)

As you can see, the exact macrostate of each parent does not have to be known:
just the macrostate probability of the single-allele gamete that is exchanged
during mating.

####  1.3.3 Homozygeous and Heterozygeous Frequencies

After inspecting $P(\Omega=k_{s}|i,j)_{recess}$ we can see this implies that
there are $H_{cor}=k_{s}$ Homozygeous correlated genes, $H_{t}=i+j-2k$
Heterozygeous genes, and $H_{un}=n-i-j+k$ Homozygeous uncorrelated genes.

Therefor, the expected frequency of correlated or uncorrelated genes going
into the next generation for two offspring without evolutionary pressures is:

$E[H]=\sum_{i,j,k_{1},k_{2}}P(\Omega_{G}=i|\rho_{BB},\rho_{BW})P(\Omega_{G}=j|\rho_{BB},\rho_{BW})P(\Omega=k_{1}|i,j,p)P^{\prime}(\Omega=k_{2}|i,j,p)$
(10)

Where

$P^{\prime}(\Omega=k_{s}|i,j,p)=\sum_{v=0}^{k_{s}}P(\Omega=k_{s}-v|i,j)_{recess}{j+i-2k_{s}+2v\choose
v}p^{v}(1-p)^{i+j-2k_{s}+2v}H(n,i,j,k_{s})$ (11)

The above notation can be simplified to

$E[H]=\sum_{i,j,k_{1},k_{2}}\Psi(i,j,k_{1},k_{2},p)$ (12)

We would expect this result to replicate the Hardy-Weinburg equilibrium which,
using simulations, we find it does.

###  1.4 The Pairwise Interaction

Due to the symmetry between both offspring in equation 12 we can define

$E[H|F]=\sum_{i,j,k_{1},k_{2}}\Psi(i,j,k_{1},k_{2},p)F(\theta(k_{1}),\theta(k_{2}))$
(13)

Where

$F(\theta(k_{1}),\theta(k_{2}))=f(\theta(k_{1}),\theta(k_{1}))+f(\theta(k_{1}),\theta(k_{2}))$
(14)

###  1.5 Modeling

We can use equation 13 to determine the frequencies of the correlated and
uncorrelated genes of the $t+1$ season. Then we can re-derive the distribution
using equation 9.

##  2 Results

Using simulations, the following important results were found:

###  2.1 Main Results

  1. 1.

Homozygeous genes (single allele per loci) found altruism to be evolutionarily
stable for cost-to-benefit ratios $>2$.

  2. 2.

Heterozygeous genes (two alleles per loci) found altruism to be evolutionarily
stable for cost-to-benefit ratios of $>1$ with very low initial populations of
altruists.

  3. 3.

When considering portions of random assortative matching, the cost-to-benefit
ratio was equal to the inverse of the probability of being matched with the
random individual in the population compared to a sibling.

###  2.2 Secondary Results

  1. 1.

Hardy-Weinberg Equilibrium was found for no evolutionary pressures with
heterozygeous genes, but not for homozygeous genes. Initial conditions only
effected the mean of the distribution.

  2. 2.

Altruism was also stable when taking was allowed, although the process was
slower and the dynamics changed.

  3. 3.

When mean centered stochastic noise is added to a cost-to-benefit ratio equal
to 1, then altruism is preferred.

##  3 Discussion

The present study extends previous work on the evolution of altruistic
behavior by introducing a novel model that incorporates Mendelian genetics,
gradations of strategy, and the potential for multiple genes to contribute to
a single phenotypic trait. The most significant finding is that when
considering heterozygous genes, altruistic behavior can be evolutionarily
stable at cost-to-benefit ratios exceeding unity, even with initially low
frequencies of altruists in the population. This result stands in contrast to
the findings for homozygous genes, which only support the stability of
altruism at cost-to-benefit ratios greater than 2.

The approach employed in this study, which models traits as macrostates
determined by the frequencies of correlated alleles across multiple genes,
provides a useful framework for capturing the distribution of traits within a
population. The linear mapping from genomic macrostate to expressed trait
serves as a reasonable first approximation.

The use of a modified dictator game as the evolutionary interaction
effectively models altruistic versus selfish behavior with adjustable costs
and benefits. However, the model’s flexibility allows for its application to a
wide range of evolutionary games and interactions. For example, the model
could be adapted to study the evolution of cooperation in the prisoner’s
dilemma, the emergence of fairness norms in the ultimatum game, or the
dynamics of reciprocity in the repeated prisoner’s dilemma. By adjusting the
payoff structure and the mapping from macrostates to traits, researchers can
investigate the evolutionary stability of various strategies in different
contexts.

Moreover, the model’s ability to incorporate multiple genes and gradations of
traits makes it well-suited for studying the evolution of complex social
behaviors and personality traits. For instance, the model could be used to
explore the evolutionary basis of empathy, trust, or aggression by considering
the interplay of multiple genetic factors and environmental influences. Such
investigations could shed light on the origins and maintenance of individual
differences in social behavior and contribute to our understanding of the
genetic architecture underlying human psychology.

The derivations of offspring macrostates arising from the macrostates of
parental gametes, considering both recessive and dominant inheritance
patterns, appear rigorous. The model’s ability to replicate the Hardy-Weinberg
equilibrium in the absence of evolutionary pressures serves as a valuable
validation.

Additional notable results include the finding that assortative matching
enables the stability of higher cost-to-benefit ratios in proportion to the
ratio of matching with siblings versus the general population, the
confirmation that altruism can still evolve even when ”taking” behavior is
allowed, albeit at a slower pace, and the demonstration that mean-centered
noise added to the cost-to-benefit ratio favors altruism.

Several avenues for further discussion and extension of this work present
themselves. First, investigating the sensitivity of the results to the assumed
linear mapping between macrostate and trait, and exploring the consequences of
nonlinear mappings, could yield valuable insights. Second, relaxing the
assumption that all genes are equally correlated with the trait may provide a
more nuanced understanding of the system. Third, considering more complex
interactions or multi-round games could reveal additional factors influencing
the evolution of altruism. Finally, conducting larger simulations with more
genes, finer gradations of traits, and larger population sizes could help
validate the scalability of the results.

In conclusion, this study represents a valuable contribution to our
understanding of the evolutionary mechanisms underlying the emergence of
altruistic behaviors. The modeling of traits as genomic macrostates with
Mendelian inheritance provides a useful framework, and the key result that
heterozygosity allows altruism to evolve under a wider range of conditions is
both intriguing and warrants further investigation. The model’s flexibility
and potential for application to various games and interactions make it a
promising tool for future research on the evolution of social behavior and
personality traits. With additional expansion and more detailed analysis, this
work has the potential to make a significant impact in the field.

Figure 2: This figure shows how the distribution of traits in a population
changes depending on pressures from the evolutionary game. The x-axis
represents the genomic macrostate, the y-axis represents the season of the
population. The 0 macrostate represents no giving in the dictator game,
whereas 3 represents giving to the highest possible amount. Specifically, the
cost-to-benefit ratio $\alpha$ was 4, and the initial frequency of genes was
$97.5\%$ selfish.  Figure 3: This figure shows how the distribution of traits
in a population changes during pressures from the evolutionary game. The
x-axis represents the genomic macrostate frequency, the y-axis represents the
season of the population. The 0 macrostate represents no giving in the
dictator game; whereas 3 represents giving to the highest possible amount.
Specifically the cost to benefit ratio $\alpha$ was 4, and the initial
frequency of genes was $97.5\%$ selfish

## References

  * [1] W. D. Hamilton, The genetical evolution of social behaviour, Journal of Theoretical Biology 7 (1964a) 1–16.
  * [2] R. Trivers, The evolution of reciprocal altruism, The Quarterly Review of Biology 46 (1) (1971) 35–57.
  * [3] J. Rushton, Genetic and environmental contributions to pro-social attitudes: a twin study of social responsibility., Proceedings of the Royal Society London 271 (2004) 2583–5.
  * [4] F. D. N. M. N. D. E. H. Rushton, J.P., Altruism and aggression: the heritability of individual differences. journal of personality and social psychology, Journal of Personality and Social Psychology 50 (1986) 1192–8.
  * [5] F. C. W. N. T. M. S. Reuter, M., C. Montag, Investigating the genetic basis of altruism: the role of the comt val158met polymorphism, Soc Cogn Affect Neurosci. (2011) 662–8doi:doi:10.1093/scan/nsq083.
  * [6] J. Maynard Smith, G. Price, The logic of animal conflict, Nature 246 (5427) (1973) 15–18.  doi:doi:10.1038/246015a0.
  * [7] J. F. Nash, Equilibrium points in n-person games, Proceedings of the National Academy of Sciences 36 (1950) 48–49.  doi:DOI:10.1073/pnas.36.1.48.
  * [8] R. L. Trivers, Parental investment and sexual selection, Proceedings of the National Academy of Sciences 36 (1972) 48–49.  doi:DOI:10.1073/pnas.36.1.48.
  * [9] I. Alger, J. W. Weibull, A generalization of hamilton’s rule—love others how much?, Journal of Theoretical Biology 299 (2012) 42–54.