Detecting Approximate Reflection Symmetry in a Point Set using Optimization on Manifold

TL;DR

This paper introduces an optimization-based algorithm on manifolds for detecting approximate reflection symmetry in point sets, demonstrating state-of-the-art accuracy and robustness in 2D and 3D cases.

Contribution

It formulates symmetry detection as a combined correspondence and transformation optimization on a Riemannian manifold, independent of descriptors.

Findings

Achieves state-of-the-art performance on benchmark datasets.

Robustly detects symmetry under varying distortions.

Effective in both 2D and 3D symmetry detection.

Abstract

We propose an algorithm to detect approximate reflection symmetry present in a set of volumetrically distributed points belonging to $\mathbb{R}^d$ containing a distorted reflection symmetry pattern. We pose the problem of detecting approximate reflection symmetry as the problem of establishing correspondences between the points which are reflections of each other and we determine the reflection symmetry transformation. We formulate an optimization framework in which the problem of establishing the correspondences amounts to solving a linear assignment problem and the problem of determining the reflection symmetry transformation amounts to solving an optimization problem on a smooth Riemannian product manifold. The proposed approach estimates the symmetry from the geometry of the points and is descriptor independent. We evaluate the performance of the proposed approach on the standard benchmark dataset and achieve the state-of-the-art performance. We further show the robustness of our approach by varying the amount of distortion in a perfect reflection symmetry pattern where we perturb each point by a different amount of perturbation. We demonstrate the effectiveness of the method by applying it to the problem of 2-D and 3-D reflection symmetry detection along with comparisons.

Figures

Tables

Table 1. TABLE I: Precision and recall rates for the methods [ 39 ] , [ 46 ] , and the proposed approach on the dataset [ 60 ] .

	Precision			Recall
	[39]	[46]	Ours	[39]	[46]	Ours
Real Images	0.21	0.30	0.42	0.75	0.95	0.93
Synthetic Images	0.28	0.29	0.73	0.93	1.00	0.96

Equations

$${\mathbf{x}}_{i^{\prime }}$$

$$\mathbf{TE}{\mathbf{T}}^{\top }{\mathbf{x}}_{i^{\prime }}={\mathbf{x}}_i-\mathbf{t}+\mathbf{TE}{\mathbf{T}}^{\top }\mathbf{t}$$

$$\Rightarrow {\mathbf{x}}_i=\mathbf{TE}{\mathbf{T}}^{\top }({\mathbf{x}}_{i^{\prime }}-\mathbf{t})+\mathbf{t}.$$

$${\mathbf{x}}_i-{\mathbf{x}}_{i^{\prime }}=\mathbf{TE}{\mathbf{T}}^{\top }({\mathbf{x}}_{i^{\prime }}-{\mathbf{x}}_i)\Rightarrow ({\mathbf{I}}_d+\mathbf{TE}{\mathbf{T}}^{\top })({\mathbf{x}}_i-{\mathbf{x}}_{i^{\prime }})=0.$$

$${\mathbf{p}}^{\top }({\mathbf{I}}_d+\mathbf{TE}{\mathbf{T}}^{\top })\mathbf{p}=0\Rightarrow {\mathbf{p}}^{\top }{\mathbf{I}}_d\mathbf{p}+({\mathbf{T}}^{\top }\mathbf{p}{)}^{\top }\mathbf{E}({\mathbf{T}}^{\top }\mathbf{p})=0$$

$$\Rightarrow {\mathbf{p}}^{\top }{\mathbf{I}}_d\mathbf{p}+{\mathbf{b}}^{\top }\mathbf{Eb}=0\Rightarrow \sum _{u=1}^d{p}_u^2+\sum _{u=1}^{d-1}{b}_u^2-b_d^2=0$$

$$\Rightarrow \sum _{u=1}^d{p}_u^2+\sum _{u=1}^d{b}_u^2-2b_d^2=0.$$

$$\sum _{u=1}^d{p}_u^2+\sum _{u=1}^d{p}_u^2-2b_d^2=0\Rightarrow \sum _{u=1}^d{p}_u^2=b_d^2\Rightarrow \sum _{u=1}^{d-1}{b}_u^2=0.$$

$$\mathbf{T}\mathbf{E}\mathbf{T}^{\top}\big{(}\mathbf{X}-\mathbf{t}\mathbf{e}^{\top}\big{)}+\mathbf{t}\mathbf{e}^{\top}=\mathbf{X}\mathbf{P}.$$

$$\underset{\mathbf{R},\mathbf{t},\mathbf{P}}{\min }$$

$$\mathbf{Pe}=\mathbf{e},{\mathbf{P}}^{\top }\mathbf{e}=\mathbf{e},\mathbf{P}\in \{0,1{\}}^{n\times n},$$

$$\{(\mathbf{R\Omega },\mathbf{t}):{\mathbf{\Omega }}_u^{\top }=-{\mathbf{\Omega }}_u,{\mathbf{\Omega }}_u\in {\mathbb{R}}^{d\times d},\forall u\in [d-1],\mathbf{t}\in {\mathbb{R}}^d\}.$$

$${⟨(\mathbf{R\Omega },{\mathbf{\eta }}_{\mathbf{t}}),(\mathbf{R}{\mathbf{\Omega }}^{\prime },{\mathbf{\eta }}_{\mathbf{t}}^{\prime })⟩}_{(\mathbf{R},\mathbf{t})}={\mathbf{\eta }}_{\mathbf{t}}^{\top }{\mathbf{\eta }}_{\mathbf{t}}^{\prime }+\sum _{u=1}^{d-1}\text{trace}({\mathbf{\Omega }}_u^{\top }{\mathbf{\Omega }}_u^{\prime }).$$

$$\text{grad }f(\mathbf{R},\mathbf{t})=({\mathbb{P}}_{\mathbf{R}}({\nabla }_{\mathbf{R}}\bar{f}),{\mathbb{P}}_{\mathbf{t}}({\nabla }_{\mathbf{t}}\bar{f}))\in {\mathcal{T}}_{(\mathbf{R},\mathbf{t})}\mathcal{M}.$$

$${\mathbb{P}}_{\mathbf{R}}({\nabla }_{\mathbf{R}}\bar{f})=({\mathbb{P}}_{{\mathbf{R}}_1}({\nabla }_{{\mathbf{R}}_1}\bar{f}),\dots ,{\mathbb{P}}_{{\mathbf{R}}_{d-1}}({\nabla }_{{\mathbf{R}}_{d-1}}\bar{f})).$$

$${\mathbb{P}}_{{\mathbf{R}}_j}({\nabla }_{{\mathbf{R}}_j}\bar{f})={\mathbf{R}}_j\text{skew}({\mathbf{R}}_j^{\top }{\nabla }_{{\mathbf{R}}_j}\bar{f}),$$

$$\text{Hess }f(\mathbf{R},\mathbf{t})[{\mathbf{\eta }}_{\mathbf{R}},{\mathbf{\eta }}_{\mathbf{t}}]=({\mathbb{P}}_{\mathbf{R}}(\text{D}{\mathbf{\xi }}_{\mathbf{R}}(\mathbf{R})[{\mathbf{\eta }}_{\mathbf{R}}]),{\mathbb{P}}_{\mathbf{t}}(\text{D}{\mathbf{\xi }}_{\mathbf{t}}(\mathbf{t})[{\mathbf{\eta }}_{\mathbf{t}}])).$$

$$({\mathbb{P}}_{{\mathbf{R}}_1}(\text{D}{\mathbf{\xi }}_{{\mathbf{R}}_1}({\mathbf{R}}_1)[{\mathbf{\eta }}_{{\mathbf{R}}_1}]),\dots ,{\mathbb{P}}_{{\mathbf{R}}_{d-1}}(\text{D}{\mathbf{\xi }}_{{\mathbf{R}}_{d-1}}({\mathbf{R}}_{d-1})[{\mathbf{\eta }}_{{\mathbf{R}}_{d-1}}])),$$

$$\text{D}{\mathbf{\xi }}_{\mathbf{x}}(\mathbf{x})[{\mathbf{\eta }}_{\mathbf{x}}]=\underset{t\to 0}{\lim }\frac{\mathbf{\xi }(\mathbf{x}+t{\mathbf{\eta }}_{\mathbf{x}})-\mathbf{\xi }(\mathbf{x})}{t}$$

$$\nabla_{\mathbf{t}}\bar{f}=2\big{(}\mathbf{I}_{d}-\mathbf{T}\mathbf{E}\mathbf{T}^{\top}\big{)}(2\mathbf{e}^{\top}\mathbf{e}\mathbf{t}-\mathbf{Xe}-\mathbf{XPe}).$$

$$\nabla_{\mathbf{R}_{j}}\bar{f}=-2\big{(}\prod_{u=1}^{j-1}\mathbf{R}_{u})^{\top}\mathbf{A}\big{(}\prod_{u=1}^{d-1}\mathbf{R}_{u}\big{)}\mathbf{E}\big{(}\prod_{u=j+1}^{d-1}\mathbf{R}_{u}\big{)}^{\top}.$$

$$\mathbf{A}=(\mathbf{XP}-\mathbf{t}{\mathbf{e}}^{\top })(\mathbf{X}-\mathbf{t}{\mathbf{e}}^{\top }{)}^{\top }+(\mathbf{X}-\mathbf{t}{\mathbf{e}}^{\top })(\mathbf{XP}-\mathbf{t}{\mathbf{e}}^{\top }{)}^{\top }$$

$$\boldsymbol{\xi}_{\mathbf{t}}(\mathbf{t})=2\big{(}\mathbf{I}_{d}-\mathbf{T}\mathbf{E}\mathbf{T}^{\top}\big{)}(2\mathbf{e}^{\top}\mathbf{e}\mathbf{t}-\mathbf{Xe}-\mathbf{XPe}),$$

$$\boldsymbol{\xi}_{\mathbf{R}_{j}}(\mathbf{R}_{j})=-\mathbf{R}_{j}\big{(}\prod\limits_{u=1}^{j}\mathbf{R}_{u}\big{)}^{\top}\mathbf{A}\big{(}\prod\limits_{u=1}^{d-1}\mathbf{R}_{u}\big{)}\mathbf{E}\big{(}\prod\limits_{u=j+1}^{d-1}\mathbf{R}_{u}\big{)}^{\top}$$

$$+\mathbf{R}_{j}\big{(}\prod\limits_{u=j+1}^{d-1}\mathbf{R}_{u}\big{)}\mathbf{E}\big{(}\prod\limits_{u=1}^{d-1}\mathbf{R}_{u}\big{)}^{\top}\mathbf{A}^{\top}\big{(}\prod\limits_{u=1}^{j}\mathbf{R}_{u}\big{)}.$$

$$\frac{1}{2}{\mathbf{R}}_j([{\mathbf{B}}_1,[{\mathbf{R}}_j^{\top }{\mathbf{B}}_2{\mathbf{R}}_j,{\mathbf{\Omega }}_j]]+[[{\mathbf{\Omega }}_j,{\mathbf{B}}_1],{\mathbf{R}}_j^{\top }{\mathbf{B}}_2{\mathbf{R}}_j]).$$

$$\mathbf{B}_{1}=\big{(}\prod_{u=j+1}^{d-1}\mathbf{R}_{u})\mathbf{E}\big{(}\prod_{u=j+1}^{d-1}\mathbf{R}_{u})^{\top},$$

$$\mathbf{B}_{2}=\big{(}\prod_{u=1}^{j-1}\mathbf{R}_{u}\big{)}^{\top}\mathbf{A}\big{(}\prod_{u=1}^{j-1}\mathbf{R}_{u}\big{)}.$$

$$\text{Hess}_{\mathbf{t}}(f(\mathbf{R,t}))[\boldsymbol{\eta}_{\mathbf{t}}]=4n\big{(}\mathbf{I}_{d}-\mathbf{T}\mathbf{E}\mathbf{T}^{\top}\big{)}\boldsymbol{\eta}_{\mathbf{t}}.$$

$$\Vert {\mathbf{X}}_{\text{m}}-\mathbf{XP}{\Vert }_{\text{F}}^2=\text{trace}(({\mathbf{X}}_{\text{m}}-\mathbf{XP}{)}^{\top }({\mathbf{X}}_{\text{m}}-\mathbf{XP}))$$

$$=\text{trace}({\mathbf{X}}_{\text{m}}^{\top }{\mathbf{X}}_{\text{m}}-2{\mathbf{X}}_{\text{m}}^{\top }\mathbf{XP}+{\mathbf{X}}^{\top }\mathbf{XP}{\mathbf{P}}^{\top }).$$

Full text

I Introduction

Symmetry present in natural and man-made objects enriches the objects to be physically balanced, beautiful, easy to recognize, and easy to understand. Characterizing and finding the symmetry has been an active topic of research in computer vision and computer graphics as physical objects form the basis for these research areas. The digitized objects are mainly represented in the form of meshes, volumes, sets of points, and images. The primary objective has been to detect symmetry in objects depicted through these different representations. We particularly aim to detect reflection symmetry present in objects represented by a set of finite number of points belonging to $\mathbb{R}^{d}$. In Fig. 1, we present an example result of the proposed approach for illustration.

The motivation behind detecting symmetry in higher dimensional spaces ($d>3$) is inspired by the fact that many physical data points reside in the space of dimensions greater than three. For example, an RGB-D image captured using a Kinect sensor, which has become a major tool for interaction of human with machine, has four dimensions at each pixel location. Another example is the embedding of feature points or shapes into a higher-dimensional space. In the scale invariant feature transform (SIFT) algorithm, each keypoint is represented in a 128-dimensional space [1]. We not only target data residing in 2-D (image) and 3-D (point cloud), but also develop a generic framework to detect symmetry in higher dimensional data.

The problem of establishing correspondences between reflection symmetry points and determining the hyperplane of reflection symmetry has been extensively studied due to its astounding applications such as compression of objects, symmetrization, shape matching, and symmetry aware segmentation of shapes [2]. Most of the existing algorithms attempt this problem by using surface signatures such as Gaussian curvature, eigenbases of the Laplace-Beltrami operator, and heat kernels for the points sampled on a given surface ([3, 4, 2]). The challenge we face is that, we only have a set of discrete points in $\mathbb{R}^{d}$. We can not take benefits from local surface signatures by fitting a surface over these points. For the case $d=2$, an explanation could be that the prominent surface signatures, such as Gaussian curvatures, are meaningful only if the surface is non-linear. For the case $d\geq 3$, an explanation could be that if the point set represents a volumetric shape, fitting a surface could be hard and eigenbases of Laplace-Beltrami operator are not defined for a set of finite points since it is not a compact manifold without the boundary [5]. Prominent methods such as [6] and [7] are independent of surface features and employ randomized algorithms to establish correspondences between the reflection symmetry points. However, they require fine tuning of a hyper-parameter to handle the reflection symmetry patterns perturbed by an unknown source of noise and an improper choice of this parameter could lead to higher time complexity.

Both these categories of algorithms are sequential in the sense that they first establish the correspondences between the reflection symmetry points and then determine the reflection symmetry hyperplane. Therefore, many outlier correspondences could be detected along with the correct correspondences. In summary, detecting symmetry in a set containing a finite number of points is a non-trivial problem. In this work, we propose an optimization framework where we jointly establish correspondences between reflection symmetry points and determine the reflection symmetry hyperplane in a set of points containing a distorted reflection symmetry pattern. In order to design the cost function, we introduce an affine transformation to obtain the reflection point of a point in $\mathbb{R}^{d}$. The main intuition behind forming this cost function is that the reflection point of a point obtained through the optimal reflection hyperplane should be present closest to its ground-truth reflection point.

The primary contributions of this work are listed below.

We propose an optimization based algorithm to establish correspondences between the reflection symmetry points and determine the reflection symmetry transformation in a set of discrete points residing in $\mathbb{R}^{d}$ containing a distorted reflection symmetry pattern. 2. 2.

We show that the proposed optimization framework is convex in translation and correspondences matrix, and locally convex in each of the rotation matrices. 3. 3.

The proposed approach is shown to not use any shape descriptors and can be applied to point sets obtained by sampling any shape residing in $\mathbb{R}^{d}$. 4. 4.

We demonstrate the effectiveness of the proposed approach by detecting symmetry in 2-D images and 3-D point clouds.

We organize the remainder of the paper as follows. In §II, we present the related works to our approach. In §III-A, we formulate the energy minimization problem. In §III-B, we find the optimal rotations and translation. In §III-C, we find the optimal mirror symmetric correspondences. In §III-D, we prove the convergence properties. In §IV, we report the computational complexity of our algorithm. In §V, we report the results and the evaluation of the proposed approach. In §VI, we conclude the work with future directions.

II Related Works

The problem of characterizing and detecting the reflection symmetry in digitized objects has been extensively studied. The works [8] and [2] provide a survey of the symmetry detection algorithms. The symmetry detection algorithms can be categorized based on either the form of the input data or whether the algorithm is dependent or independent of the surface features. General forms of the input data are: set of points, mesh, volume, and image. Most of the methods for symmetry detection in meshes first extract salient keypoints on the surface and then describe each point using local surface features. The prominent surface features are: Gaussian curvatures, slippage features, moments, geodesic distances, and extended Gaussian images ([8, 2]).

Symmetry detection in a set of points without features. These algorithms detect reflection symmetry in a set of points without using surface features. Our work also falls in this category. In the work by Zabrodsky et al., the authors find the closest shape to a given shape represented by a set of points in $\mathbb{R}^{2}$ and it requires point correspondences [9]. However, our goal is different in the sense that we find reflection correspondences within the given set of points in $\mathbb{R}^{d}$. In the work by Lipman et al., the authors propose the concept of symmetry factored embedding where they represent pairs of points which are in the same orbit in a new space and propose the concept of symmetry factored distance to find the distance between such pairs [6]. In the work by Xu et al., the authors detect multi-scale symmetry [7]. The authors use a randomized algorithm to detect the correspondences efficiently. However, performance degrades as the perfect pattern gets perturbed due to noisy measurements. We compare the correspondences established by our method to that of this method and show that our method performs better than this method when the patterns are perturbed. It is fair to compare with this method on the perturbed patterns because most of the real world patterns are not perfectly symmetric, e.g., human face and butterfly wings. In the works by Combès et al. [10], Speciale et al. [11], Ecins et al. [12], Cicconet et al. [13], Li et al. [14], and Sipiran et al. [15], the authors automatically detect the symmetry plane in a point cloud. But, the methods in [10], [11], [13], [14], and [15] do not establish correspondences. However, correspondences are an important aspect as shown in ([7, 6]). Ecins et al. proposed an ICP based approach [12] where they used the normals at each point to determine the symmetry. However, this method is applicable only to non-volumetric point clouds, i.e., points sampled from a surface.

Symmetry detection in meshes using surface features. These algorithms either directly use surface patches described using local features or first detect the salient keypoints on the surface and describe them using the local surface features. Here, we review only the salient works to give an idea of these algorithms. Mitra et al. detect partial and approximate symmetries in 3D models [4]. They start with sampling salient keypoints on the surface and pair them up using their local principal curvatures. Then using the Hough transformation, they find the pairs of reflection symmetry points. Then in the Hough transformation space, they perform the clustering of the pairs to determine all the partial symmetries. Martinet et al. detect symmetries by generalized moment functions where the shape symmetry gets inherited as symmetry in these functions [16]. Raviv et al. detect symmetry in non-rigid shapes by observing that the intrinsic geometry of a shape is invariant under non-rigid shape transformations [17]. Berner et al. start with constructing a graph based on the similarity of slippage features detected on the surface [3]. Then, they detect the structural regularities by matching the sub-graphs. Cohen et al. detect symmetry using geometric and image cues [18]. They use it to reconstruct accurate 3D models. We refer the reader to some of the pioneering works for more details on this category ([19, 20, 21, 22, 23, 24, 25, 26, 27]). There exist algorithms which find symmetry in meshes and volume without sampling keypoints. The works described in ([28, 29, 5, 30, 31, 32, 33, 34, 35]) belong to this category.

Symmetry detection algorithm for real images. These algorithms primarily rely on the local image features such as edge orientations, curvatures, and gradients. The recent works such as ([36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47]) present excellent algorithms for reflection symmetry detection in images. Given the accurate detection of keypoints, the algorithm developed in this work can be used to detect reflection symmetry in images without using local features.

Our algorithm is similar to Iterative Closest Point (ICP) algorithm ([48, 49]) only in the sense that we also follow the alternation between the optimization of reflection transformation (rotation and translation in ICP) and the correspondences between the mirror symmetric points (correspondences between the points of two different shapes in ICP). Our algorithm differs from the ICP algorithm since ICP has a different error function in the transformation parameters than the error function of our problem. Furthermore, our matching is bijective since we impose the bijectivity constraints in our optimization framework. These constraints ensure that each point has exactly one mirror image point.

III Proposed Approach

Consider a set $\mathcal{S}=\{\mathbf{x}_{i}\}_{i=1}^{n}$ of points, where $\mathbf{x}_{i}\in\mathbb{R}^{d}$, containing a distorted reflection symmetry pattern. Our goal is to determine the reflection symmetry transformation and establish the correspondences between the reflection symmetry points. In Fig. 2, we show the graphical representation of our problem. We formulate an optimization framework in which both the correspondences between reflection symmetry points and the reflection symmetry transformation are variables as described below. We use the notation $[k]$ for the set $\{1,2,\ldots,k\}$, where $k$ is a natural number.

III-A Problem Formulation

We introduce reflection transformation in $\mathbb{R}^{d}$ in order to obtain the reflection of a point through a hyperplane $\boldsymbol{\pi}$, not necessarily passing through the origin. The intuition is based on the fact that any hyperplane in $\mathbb{R}^{d}$ is a $d-1$ dimensional subspace. Therefore, it can be made parallel to the subspace spanned by any $d-1$ coordinate axes by translating the origin of the coordinate system on the hyperplane $\boldsymbol{\pi}$ and then rotating these $d-1$ axes sequentially (by the angle between the hyperplane $\boldsymbol{\pi}$ and the axis). In this new coordinate system, the reflection of a point through the hyperplane $\boldsymbol{\pi}$ can be obtained by multiplying the coordinate corresponding to the remaining axis of the point by $-1$. Then the reflection in the original coordinate system is obtained by applying the inverse procedure.

Definition 1**.**

The reflection point $\mathbf{x}_{i^{\prime}}\in\mathbb{R}^{d}$ of a point $\mathbf{x}_{i}\in\mathbb{R}^{d}$ through the reflection symmetry hyperplane $\boldsymbol{\pi}$ is determined by an affine transformation as shown in Equation 1.

[TABLE]

Here, $i,i^{\prime}\in[n]$, $\mathbf{t}\in\mathbb{R}^{d}$ is the translation vector which translates the origin of the coordinate system on the hyperplane $\boldsymbol{\pi}$, $\mathbf{R}_{u}$ is a rotation matrix of size $d\times d$ that rotates the $u^{\text{th}}$ axis about the origin such that it becomes perpendicular to the normal of the hyperplane $\boldsymbol{\pi}$, and the matrix $\mathbf{E}$ is defined as $\mathbf{E}=\begin{bmatrix}\mathbf{I}_{d-1}&\mathbf{0}_{d-1}\\ \mathbf{0}_{d-1}^{\top}&-1\end{bmatrix}$ and satisfies $\mathbf{E}^{\top}=\mathbf{E}$, $\mathbf{E}^{\top}\mathbf{E}=\mathbf{I}_{d}$. The matrix $\mathbf{R}_{u}$ is an orthogonal matrix ($\mathbf{R}_{u}^{\top}\mathbf{R}_{u}=\mathbf{R}_{u}\mathbf{R}_{u}^{\top}=\mathbf{I}_{d}$) with determinant equal to +1, $\forall u\in\{1,2,\ldots,d-1\}$. Here, $\mathbf{0}_{d-1}$ is a vector of size $(d-1)\times 1$ with all the coordinates equal to zero and $\mathbf{I}_{d-1}$ is an identity matrix of size $(d-1)\times(d-1)$.

Now, we introduce the essential properties of this transformation in order to formulate the problem. We show that the rotation matrices $(\mathbf{R}_{1},\ldots,\mathbf{R}_{d-1})$ and the translation vector $\mathbf{t}$ uniquely determine the reflection hyper-plane $\boldsymbol{\pi}$. We let $\mathbf{T}=\prod_{u=1}^{d-1}\mathbf{R}_{u}$ throughout this paper and note that it is again an orthogonal matrix with determinant equal to +1.

Theorem 1: *The point $\mathbf{x}_{i^{\prime}}$ is the reflection of the point $\mathbf{x}_{i}$ through the hyperplane $\boldsymbol{\pi}$ if and only if the point $\mathbf{x}_{i}$ is the reflection of the point $\mathbf{x}_{i^{\prime}}$ through the hyperplane $\boldsymbol{\pi}$.

*Proof: We prove the forward direction of the Theorem 1, since the backward direction can be proved in a similar way. Let us assume that the point $\mathbf{x}_{i^{\prime}}$ is the reflection of the point $\mathbf{x}_{i}$. Therefore, Equation (1) holds true. Now, we multiply Equation (1) by $\mathbf{T}\mathbf{E}\mathbf{T}^{\top}$ from left and use the identities $\mathbf{E}^{\top}=\mathbf{E},\mathbf{E}\mathbf{E}=\mathbf{I}_{d},\mathbf{T}^{\top}\mathbf{T}=\mathbf{T}\mathbf{T}^{\top}=\mathbf{I}_{d}$ to achieve,

[TABLE]

[TABLE]

∎

Theorem 2: The normal vector of the reflection hyper-plane $\boldsymbol{\pi}$ lies in the null space of the matrix $\mathbf{I}_{d}+\mathbf{T}\mathbf{E}\mathbf{T}^{\top}$, the hyper-plane $\boldsymbol{\pi}$ passes through $\mathbf{t}$, and the null space of the matrix $\mathbf{I}_{d}+\mathbf{T}\mathbf{E}\mathbf{T}^{\top}$ is an one-dimensional subspace of $\mathbb{R}^{d}$.

Proof: We subtract Equation (1) from Equation (2) to achieve

[TABLE]

Therefore, the normal to the reflection hyperplane $\boldsymbol{\pi}$, which is in the direction of the vector $(\mathbf{x}_{i}-\mathbf{x}_{i^{\prime}})$, lies in the null space of the matrix $(\mathbf{I}_{d}+\mathbf{T}\mathbf{E}\mathbf{T}^{\top})$. It is easy to show that the reflection hyperplane $\boldsymbol{\pi}$ passes through the translation $\mathbf{t}$ by noting that the reflection point of the point $\mathbf{t}$ is $\mathbf{t}$. This is possible only if the point $\mathbf{t}$ lies on the reflection hyperplane $\boldsymbol{\pi}$. We prove that the null space of the matrix $\mathbf{I}_{d}+\mathbf{T}\mathbf{E}\mathbf{T}^{\top}$ is an one-dimensional subspace of $\mathbb{R}^{d}$ in order to show that there exists an unique hyperplane $\boldsymbol{\pi}$. The nullspace of a matrix is the space spanned by the eigenvectors corresponding to the zero eigenvalue. Let $\mathbf{p}=\begin{bmatrix}p_{1}&p_{2}&\ldots&p_{d}\end{bmatrix}^{\top}\in\mathbb{R}^{d}$ be any vector. If $\mathbf{p}$ is an eigenvector corresponding to the zero eigenvalue of the matrix $\mathbf{I}_{d}+\mathbf{T}\mathbf{E}\mathbf{T}^{\top}$, then we must have

[TABLE]

[TABLE]

[TABLE]

Here, $\mathbf{b}=\mathbf{T}^{\top}\mathbf{p}$. We note that $\|\mathbf{b}\|_{2}^{2}=(\mathbf{T}^{\top}\mathbf{p})^{\top}(\mathbf{T}^{\top}\mathbf{p})=\mathbf{p}^{\top}\mathbf{p}=\|\mathbf{p}\|_{2}^{2}$. Therefore, from Equation (3) we have

[TABLE]

Therefore, $b_{1}=b_{2}=\ldots=b_{d-1}=0$ and $b_{d}\in\mathbb{R}$. Hence, the vector $\mathbf{b}$ lies in the one dimensional space $\{\mathbf{q}:q_{1}=q_{2}=\ldots=q_{d-1}=0,q_{d}\in\mathbb{R}\}$. Since $\mathbf{b}=\mathbf{T}^{\top}\mathbf{p}\Rightarrow\mathbf{p}=\mathbf{T}\mathbf{b}$. Since the rotation does not change the dimension of a linear space, the vector $\mathbf{p}$ also lies in one dimensional space. ∎

Given the set $\mathcal{S}$, our goal is to find all the correct reflection correspondences $(i,i^{\prime})\in[n]\times[n]$ and the matrices $(\mathbf{R}_{1},\mathbf{R}_{2},\dots,\mathbf{R}_{d-1},\mathbf{t})$ which define the reflection symmetry hyperplane $\boldsymbol{\pi}$. We represent all the correspondences by a permutation matrix $\mathbf{P}\in\{0,1\}^{n\times n}$, such that $\mathbf{P}_{ii^{\prime}}=1$ if the point $\mathbf{x}_{i^{\prime}}$ is the reflection point of the point $\mathbf{x}_{i}$ and $\mathbf{P}_{ii^{\prime}}=0$, otherwise. Here, we note from Theorem 1 that $\mathbf{P}_{ii^{\prime}}=1\Leftrightarrow\mathbf{P}_{i^{\prime}i}=1$.

Now, we let $\mathbf{R}=(\mathbf{R}_{1},\mathbf{R}_{2},\dots,\mathbf{R}_{d-1})\in\mathbb{V}$. Here, $\mathbb{V}=\mathbb{R}^{d\times d}\times\mathbb{R}^{d\times d}\times\ldots\times\mathbb{R}^{d\times d}$. Let $\mathbf{X}=\begin{bmatrix}\mathbf{x}_{1}&\mathbf{x}_{2}&\ldots&\mathbf{x}_{n}\end{bmatrix}\in\mathbb{R}^{d\times n}$ be the matrix containing all the points of the set $\mathcal{S}$ as its columns. Since the $i^{\text{th}}$ column of the matrix $\mathbf{X}\mathbf{P}$ is the reflection point of the point $\mathbf{x}_{i}$, the reflection transformation $(\mathbf{R},\mathbf{t})$ maps the matrix $\mathbf{X}$ to the reflected points matrix $\mathbf{X}\mathbf{P}$. Using Equation 1, we write the reflected points in the form of the matrix $\mathbf{T}\mathbf{E}\mathbf{T}^{\top}(\mathbf{X}-\mathbf{t}\mathbf{e}^{\top})+\mathbf{t}\mathbf{e}^{\top}$, where $\mathbf{e}=\begin{bmatrix}1&1&\ldots&1\end{bmatrix}^{\top}$ is a vector of size $n\times 1$. Therefore, Equation (4) holds true when the input set contains a perfect reflection symmetry pattern.

[TABLE]

In practice, a reflection symmetry pattern might have been distorted. Therefore, we would be able to find only the approximate reflection symmetry. We find the reflection transformation $(\mathbf{R},\mathbf{t})$ and the correspondences matrix $\mathbf{P}$ in such a way that the symmetry error, which we define as $\|\mathbf{T}\mathbf{E}\mathbf{T}^{\top}(\mathbf{X}-\mathbf{t}\mathbf{e}^{\top})+\mathbf{t}\mathbf{e}^{\top}-\mathbf{X}\mathbf{P}\|_{\text{F}}^{2}$, is minimized. Here $\|.\|_{\text{F}}$ is the Frobenius norm operator. We frame this problem in an optimization framework as shown in Equation (5).

[TABLE]

By imposing the constraints $\mathbf{P}\mathbf{e}=\mathbf{e}$ and $\mathbf{P}^{\top}\mathbf{e}=\mathbf{e}$, we ensure that each point has only one reflection point. We adopt an alternating optimization approach to solve the problem defined in Equation (5). We start with initializing the reflection transformation $(\mathbf{R,t})$ and solve for the optimal correspondences $\mathbf{P}$ and then for this optimal $\mathbf{P}$, we solve for optimal the $\mathbf{(R,t)}$. We continue to alternate till convergence.

Once $\mathbf{P}$ is fixed, if we minimize the cost over the set $\mathbb{V}\times\mathbb{R}^{d}$, then we have to make sure that the orthogonality and the unit determinant constraints hold true for the matrices $\mathbf{R}_{u},\forall u\in[d-1]$. One approach could be the Lagrange augmentation which requires us to handle $3d-3$ additional Lagrange multipliers. However, we observe that the set $\mathcal{M}=\{(\mathbf{R}_{1},\ldots,\mathbf{R}_{d-1},\mathbf{t}):\mathbf{R}_{u}^{\top}\mathbf{R}_{u}=\mathbf{R}_{u}\mathbf{R}_{u}^{\top}=\mathbf{I}_{d},\det(\mathbf{R}_{u})=1,\mathbf{R}_{u}\in\mathbb{R}^{d\times d},\forall u\in[d-1],\mathbf{t}\in\mathbb{R}^{d}\}$ of constraints is a smooth Riemannian product manifold over which the optimization algorithms are well studied [50].

We solve the sub-optimization problem for optimal $\mathbf{(R,t)}$ on a manifold which we discuss in Section III-B. We observe that the optimization of Equation (5) for $\mathbf{P}$ is a standard linear assignment problem for which we formulate an integer linear program which we discuss in Section III-C.

III-B Optimizing reflection transformation ($\mathbf{R,t}$)

In this step, we fix the correspondences matrix $\mathbf{P}$ and find the optimal reflection transformation $(\mathbf{R},\mathbf{t})$ by taking advantages from the differential structure of the set $\mathcal{M}$. We shall now briefly introduce the differential geometry of the set $\mathcal{M}$.

Differential geometry of the set $\mathcal{M}$ of constraints. In order to introduce the essential differential geometry of the set $\mathcal{M}$, we follow [50]. The elements of the set $\mathcal{M}$ are of the form $(\mathbf{R},\mathbf{t})\simeq(\mathbf{R_{1}},\ldots,\mathbf{R_{d-1}},\mathbf{t})$. All the orthogonal matrices (each for rotation along a single axis) with determinant +1 form a Lie group, also known as special orthogonal group, which is a smooth Riemannian manifold. The Euclidean space $\mathbb{R}^{d}$ is also a smooth Riemannian manifold. Therefore the set $\mathcal{M}$ is a product manifold, $\mathcal{SO}(2,d)\times\ldots\times\mathcal{SO}(2,d)\times\mathbb{R}^{d}$, the product of $d-1$ special orthogonal groups $\mathcal{SO}(2,d)$ and an Euclidean space $\mathbb{R}^{d}$. Each rotation matrix performs rotation about a single axis. Therefore, all the possible rotation matrices about a particular axis form a $\mathcal{SO}(2)$ embedded in the Euclidean space $\mathbb{R}^{d\times d}$. We denote this group as $\mathcal{SO}(2,d)$.

The tangent space $\mathcal{T}_{(\mathbf{R},\mathbf{t})}\mathcal{M}$ at the point $(\mathbf{R,t})\in\mathcal{M}$ is

[TABLE]

Here, $\mathbf{R\Omega}=(\mathbf{R}_{1}\mathbf{\Omega}_{1},\ldots,\mathbf{R}_{d-1}\mathbf{\Omega}_{d-1})$. The Riemannian metric $\left\langle.,.\right\rangle_{\mathbf{(R,t)}}$ on the product manifold $\mathcal{M}$, which gives the intrinsic distance between two elements $(\mathbf{R\Omega},\boldsymbol{\eta}_{\mathbf{t}})$ and $(\mathbf{R}\mathbf{\Omega}^{\prime},\boldsymbol{\eta}^{\prime}_{\mathbf{t}})$ of the tangent space at the point $\mathbf{(R,t)}$ of the manifold $\mathcal{M}$, is defined in Equation (7).

[TABLE]

Let $\bar{f}:\mathbb{V}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a scalar function. Let the function $f=\bar{f}\mid_{\mathcal{M}}$ be the restriction of the function $\bar{f}$ on the product manifold $\mathcal{M}$. Since the product manifold $\mathcal{M}$ is a submanifold of the Riemannian manifold $\mathbb{V}\times\mathbb{R}^{d}$, the Riemannian gradient of the function $f$ at the point $(\mathbf{R,t})$ is obtained by projecting the Riemannian gradient of the function $\bar{f}$ at the point $(\mathbf{R,t})\in\mathbb{V}\times\mathbb{R}^{d}$ on the tangent space at the point $(\mathbf{R,t})\in\mathcal{M}$. Therefore, the Riemannian gradient of the function $f$ at the point $(\mathbf{R,t})$ is defined in Equation (8).

[TABLE]

Since the tangent space at a point in an Euclidean space is again an Euclidean space, the second component is given by $\mathbb{P}_{\mathbf{t}}(\nabla_{\mathbf{t}}\bar{f})=\nabla_{\mathbf{t}}\bar{f}$. The first component is defined as

[TABLE]

Here,

[TABLE]

where $\text{skew}(\mathbf{A})=0.5(\mathbf{A}-\mathbf{A}^{\top})$. We define $\boldsymbol{\xi}_{\mathbf{R}_{j}}(\mathbf{R}_{j})=\mathbb{P}_{\mathbf{R}_{j}}(\nabla_{\mathbf{R}_{j}}\bar{f})$. The Riemannian Hessian of the function $f$ at a point $(\mathbf{R},\mathbf{t})$ is a linear map, $\text{Hess }f:\mathcal{T}_{(\mathbf{R},\mathbf{t})}\mathcal{M}\rightarrow\mathcal{T}_{(\mathbf{R},\mathbf{t})}\mathcal{M}$ and is defined as shown in Equation (9).

[TABLE]

Here, the first component $\mathbb{P}_{\mathbf{R}}(\text{D}\boldsymbol{\xi}_{\mathbf{R}}(\mathbf{R})[\boldsymbol{\eta}_{\mathbf{R}}])$ is equal to

[TABLE]

where $\boldsymbol{\eta}_{\mathbf{R}_{j}}=\mathbf{R}_{j}\mathbf{\Omega}_{j}$. The term

[TABLE]

is the classical derivative of the vector field $\boldsymbol{\xi}(\mathbf{x})$ in the direction $\boldsymbol{\eta}_{\mathbf{x}}$.

The Riemannian trust region method. Our goal is to minimize the function $f(\mathbf{R,t})$ over the product manifold $\mathcal{M}$. There exists a generalization of the popular optimization methods on the Riemannian manifolds. Since our problem is locally convex in each variable $\mathbf{R}_{j}$, which we prove in Theorem 6, we employ the Riemannian trust region approach [51]. It requires the Riemannian gradient and the Riemannian Hessian operator for the function $f$, which we find as follows. Let $\bar{f}$ be a function from the set $\mathbb{V}\times\mathbb{R}^{d}$ to $\mathbb{R}$ and defined as $\bar{f}(\mathbf{R},\mathbf{t})=\|\mathbf{T}\mathbf{E}\mathbf{T}^{\top}(\mathbf{X}-\mathbf{t}\mathbf{e}^{\top})+\mathbf{t}\mathbf{e}^{\top}-\mathbf{X}\mathbf{P}\|_{\text{F}}^{2}.$ Its classical gradients with respect to both the variables are given in the Equations (10) and (11). The detailed derivation is given in the Appendices §A1 and §A2.

[TABLE]

[TABLE]

Here,

[TABLE]

which satisfies $\mathbf{A}^{\top}=\mathbf{A}$. Now let the function $f=\bar{f}\mid_{\mathcal{M}}$ be the restriction of the function $\bar{f}$ on the set $\mathcal{M}$. We obtain the Riemannian gradient of the function $f$ at a point $(\mathbf{R,t})$ by projecting the Riemannian gradient of the function $\bar{f}$ over the tangent space $\mathcal{T}_{(\mathbf{R,t})}$ at the point $(\mathbf{R,t})$. Since the manifold $\mathbb{V}\times\mathbb{R}^{d}$ is an Euclidean space, the Riemannian gradient of the function $\bar{f}$ is equal to its classical gradient. Therefore, we apply the definition given in Equation (8) in order to find the Riemannian gradient $\text{grad}f(\mathbf{R,t})$ of the function $f$ which we denote as $(\boldsymbol{\xi}_{\mathbf{R}_{1}}(\mathbf{R}_{1}),\ldots,\boldsymbol{\xi}_{\mathbf{R}_{d-1}}(\mathbf{R}_{d-1}),\boldsymbol{\xi}_{\mathbf{t}})$ and define in Equations (12) and (13). The detailed derivation is given in the Appendices §A3 and §A4.

[TABLE]

[TABLE]

[TABLE]

We determine the Riemannian Hessian of the function $f$ using the definition given in Equation (9). In order to determine the $j^{\text{th}}$ component $\text{Hess}_{\mathbf{R}_{j}}(f(\mathbf{R,t}))[\mathbf{R}_{j}\Omega_{j}]$ of the Riemannian Hessian, which is equal to $\mathbb{P}_{\mathbf{R}_{j}}(\text{D}\boldsymbol{\xi}_{\mathbf{R}_{j}}(\mathbf{R}_{j})[\mathbf{R}_{j}\mathbf{\Omega}_{j}])$, we first find the classical derivative $\text{D}\boldsymbol{\xi}_{\mathbf{R}_{j}}(\mathbf{R}_{j})[\mathbf{R}_{j}\mathbf{\Omega}_{j}]$ of the Riemannian gradient $\boldsymbol{\xi}_{\mathbf{R}_{j}}(\mathbf{R}_{j})$ in the direction $\mathbf{R}_{j}\mathbf{\Omega}_{j}$ and then apply the projection operator $\mathbb{P}_{\mathbf{R}_{j}}$. Therefore, the $j^{\text{th}}$ component $\text{Hess}_{\mathbf{R}_{j}}(f(\mathbf{R,t}))[\mathbf{R}_{j}\Omega_{j}]$ of the Riemannian Hessian is equal to

[TABLE]

The detailed derivation is given in the Appendix §A5. Here $[.,.]$ is the Lie bracket and defined as $[\mathbf{U},\mathbf{V}]=\mathbf{UV}-\mathbf{VU}$ for any two matrices $\mathbf{U}$ and $\mathbf{V}$,

[TABLE]

and

[TABLE]

In a similar way, we determine the component, $\mathbb{P}_{\mathbf{t}}(\text{D}\boldsymbol{\xi}_{\mathbf{t}}(\mathbf{t})[\boldsymbol{\eta}_{\mathbf{t}}])$, of the Riemannian Hessian which is shown in Equation (15).

[TABLE]

The detailed derivation is given in the Appendix §A6. Now, we apply the Riemannian-trust-region method using the Riemannian gradient and Hessian defined in Equations (12), (13), (14), and (15) in order to obtain the optimal solution. We use the manopt toolbox in order to implement the optimization problem given in Equation (5) for a fixed $\mathbf{P}$ [52].

Determining the reflection symmetry hyperplane $\boldsymbol{\pi}$. In order to determine the reflection hyperplane $\boldsymbol{\pi}$, we use Theorem 2 which states that the normal vector of $\boldsymbol{\pi}$ lies in the null space of the matrix \mathbf{I}_{d}+\big{(}\prod_{u=1}^{d-1}\mathbf{R}_{u}\big{)}\mathbf{E}\big{(}\prod_{u=1}^{d-1}\mathbf{R}_{u}\big{)}^{\top} and the optimal translation $\mathbf{t}$ lies on the hyperplane.

III-C Optimizing Correspondences $\mathbf{P}$

After obtaining the current estimate of the reflection transformation $(\mathbf{R,t})$, we improve the correspondences matrix $\mathbf{P}$ by solving the problem given in Equation (5) while fixing $(\mathbf{R,t})$. We show that this sub-problem is equivalent to a linear assignment problem, where an assignment is a pair $(i,i^{\prime})$ of reflection symmetry points.

Claim 1: The optimization problem given in Equation (5) is a linear assignment problem in $\mathbf{P}$, for a fixed $(\mathbf{R,t})$.

Proof: Let us consider the cost function in Equation (5) and let $\mathbf{X}_{\text{m}}=\mathbf{T}\mathbf{E}\mathbf{T}^{\top}(\mathbf{X}-\mathbf{t}\mathbf{e}^{\top})+\mathbf{t}\mathbf{e}^{\top}$. We have

[TABLE]

[TABLE]

Since, the first and the third terms (using the fact that the permutation matrices are orthogonal) are not the functions of $\mathbf{P}$, the problem of finding the point of minimum of the function $\|\mathbf{X}_{\text{m}}-\mathbf{X}\mathbf{P}\|_{\text{F}}^{2}$ is identical to the problem of finding the point of maximum of the function $\text{trace}(\mathbf{X}_{\text{m}}^{\top}\mathbf{X}\mathbf{P})$. Using the identity $\text{trace}(\mathbf{A}^{\top}\mathbf{B})=\mathtt{vec}(\mathbf{A})^{\top}\mathtt{vec}(\mathbf{B})$, we have that $\text{trace}(\mathbf{X}_{\text{m}}^{\top}\mathbf{X}\mathbf{P})=\mathtt{vec}(\mathbf{X}^{\top}\mathbf{X}_{\text{m}})^{\top}\mathtt{vec}(\mathbf{P})$, where the operator $\mathtt{vec}$ vectorizes a matrix by stacking all the columns successively in a column vector.

Therefore, for a fixed reflection transformation, the problem defined in Equation (5) is equivalent to the problem defined in Equation (16).

[TABLE]

which is a standard linear assignment problem.$\square$

**Claim 2: **The problem defined in Equation (16) is an integer linear program.

Proof: Let $\mathbf{v}_{1}$ be a vector of size $n^{2}\times 1$ with the first $n$ coordinates equal to one and the last $n(n-1)$ coordinates equal to zero. Let $\mathbf{e}_{1}$ be a vector of size $n\times 1$ with all the coordinates equal to zero except the first coordinate which is equal to one. Let $\mathbf{v}_{2}=\begin{bmatrix}\mathbf{e}_{1}^{\top}&\mathbf{e}_{1}^{\top}&\ldots&\mathbf{e}_{1}^{\top}\end{bmatrix}^{\top}$ be a vector of size $n^{2}\times 1$. Now let us construct the matrices $\mathbf{A}_{1}$ and $\mathbf{A}_{2}$, each of size $n\times n^{2}$, such that the $i^{\text{th}}$ row of the matrix $\mathbf{A}_{1}$ is equal to the row vector $\mathtt{cs}(\mathbf{v}_{1}^{\top},n(i-1))$ and the $i^{\text{th}}$ row of the matrix $\mathbf{A}_{2}$ is equal to the row vector $\mathtt{cs}(\mathbf{v}_{2}^{\top},i-1)$. Here $\mathtt{cs}(\mathbf{v}^{\top},i)$ is a row vector obtained by circularly shifting any row vector $\mathbf{v}^{\top}$ right by $i$ coordinates.

Now, it is trivial to verify that the constraint $\mathbf{P}^{\top}\mathbf{e}\leq\mathbf{e}$ is equivalent to $\mathbf{A}_{1}\mathtt{vec}(\mathbf{P})\leq\mathbf{e}$ and the constraint $\mathbf{P}\mathbf{e}\leq\mathbf{e}$ is equivalent to $\mathbf{A}_{2}\mathtt{vec}(\mathbf{P})\leq\mathbf{e}$. Therefore, the problem defined in Equation (16) is equivalent to the problem defined in Equation (17).

[TABLE]

which is an integer linear program with $\mathbf{a}=\mathtt{vec}(\mathbf{P}).\square$

Solving the ILP. Since ILP is an NP-complete problem, there may not exist a polynomial time algorithm to find the optimal solution. We relax this ILP to a linear program by converting the constraint $\mathbf{a}\in\{0,1\}^{n^{2}\times 1}$ into $\mathbf{a}\in[0,1]^{n^{2}\times 1}$. Now, the above ILP becomes a linear program. We first solve this LP using the Karmarkar’s algorithm in [53] which takes $O(n^{3.5})$ time. The solution $\mathbf{a}^{\star}=\begin{bmatrix}a^{\star}_{1}&a^{\star}_{2}&\ldots&a^{\star}_{n^{2}}\end{bmatrix}^{\top}$ of this LP belongs to $[0,1]^{n^{2}\times 1}$ which is a continuous solution. However, our final solution $\mathbf{a}^{f}=\begin{bmatrix}a^{f}_{1}&a^{f}_{2}&\ldots&a^{f}_{n^{2}}\end{bmatrix}^{\top}$ of the proposed ILP should be a discrete solution. We follow the rounding approach, as explained in ([54], ch. 11). The $i$-th element $a^{f}_{i}$ of the final solution is equal to $1$, if $a^{\star}_{i}\geq 0.5$ and equal to 0, if $a^{\star}_{i}<0.5$. This solution $\mathbf{a}^{f}$ may not be the optimal solution because according to [54], $\mathtt{vec}(\mathbf{X}_{\text{m}}^{\top}\mathbf{X})^{\top}\mathbf{a}^{f}\geq\frac{1}{2}\times\mathtt{vec}(\mathbf{X}_{\text{m}}^{\top}\mathbf{X})^{\top}\mathbf{a}^{OPT}$. Here, $\mathbf{a}^{OPT}$ is the optimal solution of the above ILP.

III-D Convergence Analysis

We derive the essential results in order to prove that the alternating optimization framework converges.

**Theorem 3: ** *The cost function $f(\mathbf{R,t,P})$ is convex in the variable $\mathbf{t}$. *

Proof: In order to prove this, we prove that the Riemannian Hessian of the function $f$ with respect to the variable $\mathbf{t}$ is a positive semi-definite (PSD) matrix. Let $\boldsymbol{\eta}_{\mathbf{t}}=\begin{bmatrix}\eta_{1}&\eta_{2}&\ldots&\eta_{d}\end{bmatrix}^{\top}\in\mathbb{R}^{d}$. Then using the definition of Riemannian metric, we have

[TABLE]

Now, using the Riemannian Hessian $\text{Hess}_{\mathbf{t}}(f)[\boldsymbol{\eta}_{\mathbf{t}}]$ defined in Equation (15), we have that

[TABLE]

Now let $\mathbf{q}=\mathbf{T}^{\top}\mathbf{\boldsymbol{\eta}_{\mathbf{t}}}$. Then, we obtain

[TABLE]

[TABLE]

Now, we know that $\mathbf{T}\mathbf{T}^{\top}=\mathbf{I}$. Hence, we have

[TABLE]

Therefore,

[TABLE]

**Theorem 4: ** At the critical point, the matrix $\mathbf{T}^{\star}=\prod_{u=1}^{d}\mathbf{R}_{u}^{\star}$ contains the eigenvectors of the matrix $\mathbf{A}$ as columns.

Proof: At the critical point, the Riemannian gradient given in Equation (13) vanishes. Therefore, $\boldsymbol{\xi}_{\mathbf{R}_{j}}(\mathbf{R}_{j})=\mathbf{0}_{d\times d}$. Now pre-multiplying it with \big{(}\prod_{u=1}^{j}\mathbf{R}_{u}\big{)}\mathbf{R}_{j}^{\top} and then post-multiplying with \big{(}\prod_{u=j+1}^{d-1}\mathbf{R}_{u}\big{)}, we achieve

[TABLE]

Now, let $\mathbf{Q}=\begin{bmatrix}\mathbf{Q}_{1}&\mathbf{q}_{2}\\ \mathbf{q}^{\top}_{3}&q_{4}\end{bmatrix}=(\mathbf{T}^{\star})^{\top}\mathbf{A}\mathbf{T}^{\star}$ be a matrix. Then, we have $\mathbf{QE}=\mathbf{EQ}$. Therefore,

[TABLE]

[TABLE]

Since, $\mathbf{I}_{d-1}$ is a diagonal matrix and the equality $\mathbf{Q}_{1}\mathbf{I}_{d-1}=\mathbf{I}_{d-1}\mathbf{Q}_{1}$ holds true, it is easy to prove that $\mathbf{Q}_{1}$ is a diagonal matrix. Therefore, the matrix $\mathbf{Q}$ is also diagonal. The spectral theorem states that every real symmetric matrix has eigenvalue decomposition with real eigenvalues and orthogonal eigenvectors. Here, we have observed that the matrix $\mathbf{A}$ is a real symmetric matrix and satisfies $\mathbf{Q}=(\mathbf{T}^{\star})^{\top}\mathbf{A}\mathbf{T}^{\star}$, where the matrix $\mathbf{Q}$ is a diagonal matrix and the matrix $\mathbf{T}^{\star}$ is an orthogonal matrix. Therefore, the matrix $\mathbf{T}^{\star}$ is the matrix containing the eigenvectors of the matrix $\mathbf{A}$. In Theorem 5, we prove that the order of stacking eigenvectors of $\mathbf{A}$ as columns of $\mathbf{T}^{\star}$ affects the convexity of the problem. $\square$

Theorem 5: *The cost function $f(\mathbf{R,t,P})$ is locally convex in each rotation matrix $\mathbf{R}_{j}$. *

Proof: In order to show the local convexity in $\mathbf{R}_{j}$, we have to show that the value $\left\langle\mathbf{R}_{j}\mathbf{\Omega}_{j},\mathbf{H}[\mathbf{R}_{j}\mathbf{\Omega}_{j}]\right\rangle_{\mathbf{R}_{j}}\geq 0$ in the neighborhood of the optimal angle $\theta_{j}^{\star}$. Here, $\mathbf{H}[\mathbf{R}_{j}\mathbf{\Omega}_{j}]=\text{Hess}_{\mathbf{R}_{j}}(f(\mathbf{R,t}))[\mathbf{R}_{j}\mathbf{\Omega}_{j}]$. By using the Riemannian metric defined in Equation (7), we have

[TABLE]

By using Equation (14), the matrix $\mathbf{R}_{j}^{\top}\mathbf{H}[\mathbf{R}_{j}\mathbf{\Omega}_{j}]$ is equal to

[TABLE]

In the Appendix §A7, we show that the $\text{trace}(\mathbf{\Omega}_{j}^{\top}\mathbf{R}_{j}^{\top}\mathbf{H}[\mathbf{R}_{j}\mathbf{\Omega}_{j}])$ is equal to

[TABLE]

We visualize this term for $d=2$. For $d=2$, the matrix $\mathbf{\Omega}=\begin{bmatrix}0&-\theta\\ \theta&0\end{bmatrix}$, $\mathbf{E}=\begin{bmatrix}1&0\\ 0&-1\end{bmatrix}$, and let $\mathbf{A}=\begin{bmatrix}a_{1}&a_{2}\\ a_{2}&a_{3}\end{bmatrix}$ and $\mathbf{R}=\begin{bmatrix}\cos\theta&-\sin\theta\\ \sin\theta&\cos\theta\end{bmatrix}$. We have that

[TABLE]

In Fig. 3, we plot the value $\frac{\left\langle\mathbf{R\Omega},\mathbf{H}[\mathbf{R\Omega}]\right\rangle_{\mathbf{R}}}{\theta^{2}}$ against the initialization angle $\theta$ for six reflection symmetry patterns having different orientations for symmetry axis. We observe that the PSD values are positive in the proximity of the optimal angles. Therefore, it is locally convex. We further observe that this quantity is maximum at the optimal angle. We also observe that, if $\theta$ is the symmetry axis orientation, then the PSD value becomes positive in the proximity of $\theta$ and $\theta+180^{\circ}$. The reason for the second range is that, if $\theta$ is the slope of a line, then $\theta+180^{\circ}$ is also the slope of the same line.

In Theorem 4, we claimed that the order in which the eigenvectors are stacked as columns of the matrix $\mathbf{R}$ affects the local convexity. We prove it as follows. At the critical point, we have that $\mathbf{R}^{\top}\mathbf{A}\mathbf{R}=\text{diag}(d_{1},d_{2})$. We note that $\mathbf{\Omega}_{j}\mathbf{B}_{1}\mathbf{\Omega}_{j}-\mathbf{\Omega}_{j}\mathbf{\Omega}_{j}\mathbf{B}_{1}=\mathbf{E}$ for $d=2$. Now from Equation (18), we achieve

[TABLE]

Therefore, the first column of the matrix $\mathbf{R}^{\star}$ should be the eigenvector corresponding to the maximum eigenvalue and the second column of the matrix $\mathbf{R}^{\star}$ should be the eigenvector corresponding to the minimum eigenvalue of the matrix $\mathbf{A}$.$\square$

**Theorem 6: ** The proposed alternating framework converges to the global minimum if the initialization of the rotation matrices $\mathbf{R}_{1}$, $\ldots$, $\mathbf{R}_{d-1}$ are within the proximity of the optimal rotation matrices and initialization of the translation $\mathbf{t}$ is any random vector.

Proof: We observe that the proposed alternation framework is basically the block coordinate descent (BCD) method, where $(\mathbf{R}_{1},\ldots,\mathbf{R}_{d-1},\mathbf{t})$ and $\mathbf{P}$ are two blocks of coordinates. According to [55], the BCD method converges if the cost function is convex in each block of coordinates. Here, we have seen that the cost function is convex in the coordinates $\mathbf{t}$ (Theorem 3), convex in the coordinates $\mathbf{P}$ on the relaxed domain $[0,1]^{n\times n}$, and locally convex in the coordinates $(\mathbf{R}_{1},\ldots,\mathbf{R}_{d-1})$ (Theorem 5). This implies that if the initialization of $(\mathbf{R}_{1},\ldots,\mathbf{R}_{d-1})$ is within the proximity of the optimal solution, then the alternating framework converges to the global minimum. We experimentally show this theorem for the case $d=2$. We use the dataset for $d=2$ with $\sigma=0$ as mentioned in §V-B. In Fig. 4, we plot the error (averaged over all optimal angles) at the convergence point against the initialization angles for the case $d=2$ (we shift the error vectors for different optimal angles so that the optimal angle is always $90^{\circ}$). We observe that the variance becomes zero for initialization angle $\theta_{0}\in(90^{\circ}-12^{\circ},90^{\circ}+9^{\circ})$ and $\theta_{0}\in(270^{\circ}-12^{\circ},270^{\circ}+9^{\circ})$. The reason for the second range is that, if $\theta$ is the slope of a line, then $\theta+180^{\circ}$ is also the slope of the same line.$\square$

In summary, in order to obtain the optimal $\mathbf{(R^{\star},t^{\star},\mathbf{P^{\star}})}$, we follow Algorithm 1.

Initialization Strategy: In the Theorem 5, we have shown that $f(\mathbf{R,t,P})$ is locally convex in rotation matrix $\mathbf{R}$. Therefore, Algorithm 1 converges to the global minimum if we initialize the rotation matrix in the proximity of the global solution. Hence, we approximate the initial $\mathbf{R}$ by finding a small set of candidate pairs of mirror symmetric points. We discuss the proposed approach for finding a small set of candidate pairs of mirror symmetric points as follows.

Let us consider the input set $\mathcal{S}=\{\mathbf{x}_{i}\}_{i=1}^{n}$. We propose a randomized approach to find a small set of candidate pairs of mirror symmetric points. We select two points, $\mathbf{x}_{p}$ and $\mathbf{x}_{q}$, uniformly at random from the set $\mathcal{S}$. Let $\mathbf{x}_{p^{\prime}}$ and $\mathbf{x}_{q^{\prime}}$ be their actual mirror images, respectively. We then construct two sets, $\mathcal{P}=\{(\mathbf{x}_{p},\mathbf{x}_{i})\}_{i=1,i\neq p,q}^{n}$ and $\mathcal{Q}=\{(\mathbf{x}_{q},\mathbf{x}_{i})\}_{i=1,i\neq q,p}^{n}$ of pairs of points. Given the sets $\mathcal{P}$ and $\mathcal{Q}$, our goal is to find the pairs $(\mathbf{x}_{p},\mathbf{x}_{p^{\prime}})$ and $(\mathbf{x}_{q},\mathbf{x}_{q^{\prime}})$. It is trivial to observe that $(\mathbf{x}_{p},\mathbf{x}_{p^{\prime}})\in\mathcal{P}$ and $(\mathbf{x}_{q},\mathbf{x}_{q^{\prime}})\in\mathcal{Q}$. We note that each pair of points define its own symmetry plane, the one which is perpendicular to the line segment joining the two points and passing through the mid-point of this line segment. Now, if the pairs $(\mathbf{x}_{p},\mathbf{x}_{p^{\prime}})$ and $(\mathbf{x}_{q},\mathbf{x}_{q^{\prime}})$ are true pairs then both the reflection planes, defined by these two pairs, should be the same. For each pair $(\mathbf{x}_{p},\mathbf{x}_{i})\in\mathcal{P}$, we keep sampling a pair $(\mathbf{x}_{q},\mathbf{x}_{j})\in\mathcal{Q}$ uniformly at random without replacement until the reflection planes defined by these two pairs are the same. We determine whether the two reflection planes, defined by these two pairs, $\boldsymbol{\pi}_{pi}:\boldsymbol{\eta}_{pi}^{\top}\mathbf{x}-c_{pi}=0$ and $\boldsymbol{\pi}_{qj}:\boldsymbol{\eta}_{qj}^{\top}\mathbf{x}-c_{qj}=0$ are the same if the conditions, $\cos^{-1}(\boldsymbol{\eta}_{pi}^{\top}\boldsymbol{\eta}_{qj})\leq\epsilon_{\theta}$ and $\frac{\min\{d_{q},d_{j}\}}{\max\{d_{q},d_{j}\}}\geq 1-\epsilon_{d}$ are true. Here, $\boldsymbol{\eta}_{pi}=\frac{\mathbf{x}_{p}-\mathbf{x}_{i}}{\|\mathbf{x}_{p}-\mathbf{x}_{i}\|_{2}}$ is the normal vector to the plane $\boldsymbol{\pi}_{pi}$, $c_{pi}=\boldsymbol{\eta}_{pi}^{\top}(\frac{\mathbf{x}_{p}+\mathbf{x}_{i}}{2})$ is the distance of the origin from the plane $\boldsymbol{\pi}_{pi}$, $\boldsymbol{\eta}_{qj}=\frac{\mathbf{x}_{q}-\mathbf{x}_{j}}{\|\mathbf{x}_{q}-\mathbf{x}_{j}\|_{2}}$ is the normal vector to the plane $\boldsymbol{\pi}_{qj}$, $c_{qj}=\boldsymbol{\eta}_{qj}^{\top}(\frac{\mathbf{x}_{q}+\mathbf{x}_{j}}{2})$ is the distance of the origin from the plane $\boldsymbol{\pi}_{qj}$, $d_{q}=|\boldsymbol{\eta}_{pi}^{\top}\mathbf{x}_{q}-c_{pi}|$, and $d_{j}=|\boldsymbol{\eta}_{pi}^{\top}\mathbf{x}_{j}-c_{pi}|$.

We repeat the above procedure ten times. With this, we get a set of 20 (2 for each run) candidate pairs of mirror symmetric points. Since we consider the case where only a single symmetric object is present in the input set, we consider the median plane of the 20 planes defined by the above computed 20 candidate pairs. Now, we use the normal $\boldsymbol{\eta}$ to this median plane for initialization. We also initialize the initial translation vector $\mathbf{t}$ as the median of the mid-points of the line segment joining the points of the candidate pairs of the mirror symmetric points.

First, we subtract each data point of the point cloud from the estimated $\mathbf{t}$ of the point cloud. This ensures that the reflection symmetry plane passes through the origin. Now, we know the unit normal to the reflection symmetry plane. Therefore, we use the Householder transform to reflect each point which is $\mathbf{x}_{i^{\prime}}=(\mathbf{I}-2\boldsymbol{\eta}\boldsymbol{\eta}^{\top})\mathbf{x}_{i}$. Therefore, we have the matrix $\mathbf{X}$ containing the original point cloud and the matrix $\mathbf{X}_{\text{m}}$ containing the reflected point cloud about the estimated reflection symmetry plane. Now, using $\mathbf{X}$ and $\mathbf{X}_{\text{m}}$, we solve the linear assignment problem, defined in Equation (17) to find the matrix $\mathbf{P}$. Now, we use these approximate correspondences to estimate the reflection symmetry plane as step 4 of Algorithm 1.

IV Time Complexity

There are two main steps involved in our algorithm. The first one is to solve for reflection symmetry transformation matrices $\mathbf{R}_{1},\mathbf{R}_{2},\ldots,\mathbf{R}_{d-1},\mathbf{t}$ using the Riemannian trust region [51]. The second step is to find the pairs of reflective symmetric points using an integer linear program. The time complexity of Riemannian trust region method is $O(nd^{2})$. Since solving integer linear program is an $NP$-complete problem, we first relax it to a linear program (as discussed at the end of §III-C). The time complexity of solving a linear program is polynomial in the number of points in the point cloud. We use the Karmarkar’s algorithm in [53] which has the time complexity of $O(n^{3.5})$. Therefore, the overall complexity of our approach is polynomial in the number of points in the point cloud which is equal to $O(nd^{2})+O(n^{3.5})\approx O(n^{3.5})$, since $d<<n$. It takes around 38.4 seconds ($d=3$) to find the symmetry plane and all the pairs of mirror symmetric points in a point cloud with 500 points using MATLAB on a Linux machine with i7-7500U CPU @ 2.70GHz, and 16GB RAM.

V Evaluation and Results

V-A Evaluation of Reflection Symmetry Plane

In order to evaluate the performance of reflection symmetry plane detection, we compare the performance of our approach with the performance of the methods in [12], [13], and [11]. We compare the detected plane of reflection symmetry to that of these methods on the dataset in [56] with F-score as the evaluation metric proposed in [56]. The dataset given in [56] contains models of 1354 3D real world objects in which the ground-truth plane of reflection symmetry is provided for all the objects.

Speciale et al. proposed a Hough transform voting based approach [11]. Ecins et al. proposed an ICP based approach [12]. First, they initialize the reflection symmetry plane and then iteratively update the reflection symmetry plane using the Levenberg-Marquardt solver till convergence. They have further used the normals at each point to reject outliers points. Therefore, they need oriented point clouds, i.e., normal at each point be given. Cicconet et al. first reflected the original point cloud about an arbitrary reflection plane and then used the ICP algorithm to align the original point cloud and the reflected point cloud [13]. Then, they determine the reflection symmetry plane.

In order to evaluate the accuracy of detecting reflection symmetry plane for each method, we find the precision and recall rates and the $F$-score. According to [56], the precision and the recall rates are defined as $P=\frac{TP}{TP+FP}$, $R=\frac{TP}{TP+FN}$, respectively. The $F$-Score is defined as $F=\frac{2RP}{R+P}$. According to [56], $TP$ is equal to the number of correctly estimated reflection symmetry planes, $FP$ is equal to the number of incorrectly estimated reflection symmetry planes, and $FN$ is equal to the number of ground-truth reflection symmetry planes which are not detected. According to [56], a detected plane of reflective symmetry is declared to be correct or incorrect as follows. Let $\mathbf{x}_{1}^{\text{e}}$, $\mathbf{x}_{2}^{\text{e}}$, and $\mathbf{x}_{3}^{\text{e}}$ be three points on the detected plane of reflection symmetry. Let $\mathbf{x}_{1}^{\text{g}}$, $\mathbf{x}_{2}^{\text{g}}$, and $\mathbf{x}_{3}^{\text{g}}$ be three points on the ground truth plane of reflection symmetry of the underlying symmetric object. These three points on the plane of reflection symmetry planes are any three points from the four points of intersection of the plane of reflection symmetry with the bounding box of the underlying reflective symmetric object. Now, according to [56], the detected plane of reflection symmetry is declared correct if the angle between the normal of the detected plane of reflection symmetry, which is defined as $\boldsymbol{\eta}_{\text{e}}=(\mathbf{x}_{1}^{\text{e}}-\mathbf{x}_{2}^{\text{e}})\times(\mathbf{x}_{1}^{\text{e}}-\mathbf{x}_{3}^{\text{e}})$, and the normal of the ground truth plane of reflection symmetry, which is defined as $\boldsymbol{\eta}_{\text{g}}=(\mathbf{x}_{1}^{\text{g}}-\mathbf{x}_{2}^{\text{g}})\times(\mathbf{x}_{1}^{\text{g}}-\mathbf{x}_{3}^{\text{g}})$, is less than a predefined threshold, i.e., $\cos^{-1}\left(\frac{|\boldsymbol{\eta}_{\text{e}}^{\top}\boldsymbol{\eta}_{\text{g}}|}{\|\boldsymbol{\eta}_{\text{e}}\|_{2}\|\boldsymbol{\eta}_{\text{g}}\|_{2}}\right)<t_{\theta}$. Furthermore, according to [56], the distance between the center of the detected plane of reflection symmetry, which defined as $\mathbf{c}_{\text{e}}=\frac{\mathbf{x}_{1}^{\text{e}}+\mathbf{x}_{2}^{\text{e}}}{2}$, from the ground truth plane of reflection symmetry is less than a predefined threshold, i.e., $\frac{|\mathbf{c}_{\text{e}}^{\top}\boldsymbol{\eta}_{\text{g}}-\boldsymbol{\eta}_{\text{g}}^{\top}\mathbf{x}_{1}^{\text{g}}|}{\|\boldsymbol{\eta}_{\text{g}}\|_{2}}<t_{d}$. In order to find the precision vs. recall curve, we change the threshold for angle as $t_{\theta}\in[0,45^{\circ}]$ and the threshold for distance as $t_{d}\in[0,2s]$. Here, $s=\min\{\|\mathbf{x}_{1}^{\text{e}}-\mathbf{x}_{2}^{\text{e}}\|_{2},\|\mathbf{x}_{1}^{\text{e}}-\mathbf{x}_{3}^{\text{e}}\|_{2},\|\mathbf{x}_{1}^{\text{g}}-\mathbf{x}_{2}^{\text{g}}\|_{2},\|\mathbf{x}_{1}^{\text{g}}-\mathbf{x}_{3}^{\text{g}}\|_{2}\}$. In Fig. 5, we plot the recall vs. precision curves for the methods in [13], [12], [11], and the proposed approach on the dataset given in [56]. We show the maximum $F$-score for each method in the legends. The maximum $F$-score for [12] is equal to 0.83, for [13] is equal to 0.67, for [11] is equal to 0.73, and for the proposed approach is equal to 0.86.

V-B Robustness to Perturbations

In order to measure the qualitative performance of the proposed approach, we investigate the following two errors which are functions of the perturbation radius $\sigma^{2}$:

[TABLE]

The error $e_{\text{d}}$ represents how well the vectors, along the line segments joining the estimated reflection symmetry points, align with the normal to the hyperplane $\boldsymbol{\pi}$ at convergence. The error $e_{\text{m}}$ represents how well the mid-points of line segments joining reflection symmetry points lie on the estimated hyperplane $\boldsymbol{\pi}$. Here, $\hat{\mathbf{z}}_{i}=\frac{\mathbf{x}_{i}-\mathbf{x}_{i^{\prime}}}{\|\mathbf{x}_{i}-\mathbf{x}_{i^{\prime}}\|_{2}}$, $\hat{\mathbf{v}}$ is the unit normal to the hyperplane $\boldsymbol{\pi}$, $\mathbf{x}_{i}^{\text{m}}=\frac{\mathbf{x}_{i}+\mathbf{x}_{i^{\prime}}}{2}$, and $w_{0}$ is the distance of the hyperplane $\boldsymbol{\pi}$ from the origin. In Fig. 6, we show the errors $e_{\text{d}}$ and $e_{\text{m}}$ against the perturbation radius $\sigma^{2}$. We observe that the values $e_{\text{m}}$ and $e_{\text{d}}$ for the proposed approach are close to that of the ground-truth reflection symmetry even as the value of $\sigma^{2}$ increases. We construct the following dataset to perform the above experiment. Let $\{\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{\frac{n}{2}}\}$ be the randomly chosen $\frac{n}{2}$ points. Given the reflection transformations $\{\mathbf{R}_{1},\ldots,\mathbf{R}_{d-1},\mathbf{t}\}$, we reflect these points using the Definition 1 in order to get the final symmetric set $\mathcal{S}=\{\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{\frac{n}{2}},\mathbf{x}_{1}^{\prime},\mathbf{x}_{2}^{\prime},\ldots,\mathbf{x}_{\frac{n}{2}}^{\prime}\}$. Then, we perturb each point with random noise as $\mathbf{x}\leftarrow\mathbf{x}+\mathcal{N}(\mathbf{0}_{d},\text{diag}(\sigma^{2},\sigma^{2},\ldots,\sigma^{2})),\forall\mathbf{x}\in\mathcal{S}$. Here, $\sigma^{2}$ is the perturbation radius and the perturbation is different for each point. For the case $d=2$, we create sets containing reflection symmetry patterns with $n\in\{50,100,150,200,250,300\}$ with $0\leq x,y\leq 1$. For each $n$, we take 19 different symmetry axis orientations in the range from $-90^{\circ}$ to $90^{\circ}$ with step size of $10^{\circ}$. We choose $\sigma^{2}\in\{0,0.01,0.02,\ldots,0.1\}$ to get 11 different perturbations. In total, we have 1254 sets for the evaluation. In Fig. 7, we show an example point set from this dataset. For the case $d=3$, we create reflective symmetric sets with $n\in\{50,100,150,200,250,300\}$ with $0\leq x,y\leq 1$. For each $n$, we take 16 different symmetry plane orientations by considering $\theta_{1}\in\{-30^{\circ},0^{\circ},35^{\circ},80^{\circ}\}$ and $\theta_{2}\in\{-30^{\circ},0^{\circ},35^{\circ},80^{\circ}\}$. We choose $\sigma^{2}\in\{0,0.01,\ldots,0.1\}$. In total, we obtain 1056 point sets.

V-C Evaluation in Higher Dimensional Data

Datasets. Since datasets for higher dimensions ($d>3$) are not available with ground-truth reflection symmetry, we synthetically create datasets as follows. For the case $d=6$ and $d=8$, we create mirror symmetric point clouds using Definition 1, with $n\in\{50,100,150,200,250,300\}$ and $0\leq x,y\leq 1$. For each $n$, we take 20 random symmetry plane normals. We choose $\sigma^{2}\in\{0,0.02,0.04,\ldots,0.1\}$ to get 6 different perturbations. In total, we have 720 sets for evaluation. For all these point clouds, we have the ground-truth correspondences between the symmetric points and the normals to the ground-truth symmetry planes.

Evaluation of correspondences. In order to evaluate the performance, we measure the correspondence rate which is the number of correct correspondences out of the estimated correspondences. Let $(i,i_{\text{e}}^{\prime})$ be the estimated correspondence and let $(i,i_{\text{g}}^{\prime})$ be the ground-truth correspondence. Then, we decide if the estimated correspondence $(i,i_{\text{e}}^{\prime})$ is correct based on a distance threshold $\tau_{d}$. If the distance $\|\mathbf{x}_{i_{\text{e}}^{\prime}}-\mathbf{x}_{i_{\text{g}}^{\prime}}\|_{2}$ between the points $\mathbf{x}_{i_{\text{e}}^{\prime}}$ and $\mathbf{x}_{i_{\text{g}}^{\prime}}$ is less than the distance threshold $\tau_{d}$, then the correspondence $(i,i_{\text{e}}^{\prime})$ is correct and otherwise, incorrect. For a given threshold $\tau_{d}$, we count the correspondences $(i,i_{\text{e}}^{\prime})$ for which the condition $\|\mathbf{x}_{i_{\text{e}}^{\prime}}-\mathbf{x}_{i_{\text{g}}^{\prime}}\|_{2}<\tau_{d}$ holds true. In Fig. 8, we show the correspondence rate vs the distance threshold curves for the different perturbation radius $\sigma^{2}\in\{0,0.02,0.04,\ldots,0.1\}$ and for $d=6$ and $d=8$. We vary the distance threshold as $\tau_{d}\in\{0,0.01,0.02,\ldots,0.34\}$. We observe that the correspondence rate increases as the distance threshold increases and the correspondence rate decreases as the perturbation radius increases for both $d=6$ and $d=8$.

Evaluation of symmetry plane. To evaluate the performance of the reflection plane detection in higher dimensional point clouds ($d>3$), instead of finding $d-1$ points on the estimated hyperplane (since finding $d-1$ points could be difficult), we measure the distance between their normals. Without loss of generality, we prepare the dataset such that the reflection symmetry plane passes through the origin. Now, let $\boldsymbol{\eta}^{\text{g}}$ and $\boldsymbol{\eta}^{\text{e}}$ be the unit normals to the ground-truth and the estimated reflection symmetry planes, respectively. Then, we declare the estimated reflection symmetry plane to be correct, if $\cos^{-1}(|(\boldsymbol{\eta}^{\text{g}})^{\top}\boldsymbol{\eta}^{\text{e}}|)<\tau_{\theta}$. We vary the angle threshold $\tau_{\theta}$ in the range $[0^{\circ},5^{\circ}]$ with a step size of $0.01^{\circ}$. In Fig. 9, we show the precision rate vs the angle threshold $\tau_{\theta}$ curves for different perturbation radius $\sigma^{2}\in\{0,0.02,0.04,\ldots,0.1\}$ and for $d=6$ and $d=8$. We observe that the precision rate increases as the angle threshold increases and decreases as the perturbation radius increases for both $d=6$ and $d=8$.

V-D Results

In Fig. 10, we show the detected reflection symmetry for two real 3D scans of buildings from the dataset [56]. In Fig. 11, we present the results for the case $d=3$. The point cloud in Fig. 10(a) contains 912045 points and the point cloud in Fig. 10(b) contains 767474 points. Since the computational complexity is $O(n^{3.5})+O(nd^{2})$, the computation time and space requirement (storing the matrices $\mathbf{A}_{1}$ and $\mathbf{A}_{2}$) are very high. Therefore, in order to compute the reflection symmetry in these scans, we randomly sample around 600 points. In both cases, we show the reflection symmetry plane by the blue color and estimated pairs of reflective symmetric points by the red colored line segment joining them. In order to make our algorithm robust to the part removal, we simply put the extra constraint $\mathbf{e}^{\top}\mathbf{Pe}\leq 2k$ in ILP defined in Equation (17) which limits the number of pairs to at most $k$. For $d=2$, we detect reflection symmetry in the set of corner points in a real image. In order to determine the symmetry axis, we use Theorem 2. For $d=3$, we use existing standard 3D models dataset [57]. In order to calculate the symmetry axis in an image using the proposed approach, we first find the set of corner points [58]. This set may contain the corners not lying on the symmetric object. Therefore, we apply the proposed approach with RANSAC [59]. We compare the proposed results with the results of two descriptor based methods [39] and [46]. We evaluate on real and synthetic images containing single symmetric object from the dataset [60]. In TABLE I, we present the precision and the recall rates. We observe that for synthetic images, the precision rate is very high for the proposed approach because most of the corner points lie on the symmetric object. Whereas, in real images, the set of corner points contains many outlier corners which leads to the degraded performance. Precision rates for the proposed approach are higher than that for the methods [39] and [46]. The recall rates are better than that of the method [39] and comparable to that of the method [46]. This leads to the conclusion that symmetry detection can be performed even when the feature descriptors are not available. In Fig. 12, we show the results on the datasets [61, 62], and [60]. The last two images show the failure cases from the datasets [60]. The reason could be that the pixels which are responsible for symmetry detection such as pixels on eyes and ear tips in the second image are not detected in the corner point detection step.

Influence of Different Initializations. We first create the following dataset of 3D point clouds. We create 5000 point clouds $\{S_{i}\}_{i=1}^{5000}$ with known ground-truth symmetries as discussed in §5.2. We keep 500 points in each point cloud. Without loss of generality, we choose the reflection symmetry plane such that it makes $90^{\circ}$ angle with $x$-axis and $y$-axis, i.e., the $x$-$y$ plane. For each point cloud, we initialize the variable $\mathbf{t}_{i}^{0}=\text{mean}(S_{i})$ and $(\theta^{0}_{x},\theta^{0}_{y})$ on every point of the grid domain $\{-90^{\circ},-80^{\circ},\ldots,+80^{\circ},+90^{\circ}\}\times\{-90^{\circ},-80^{\circ},\ldots,+80^{\circ},+90^{\circ}\}$. We then run our approach and measure the error at the convergence $e_{i}(\theta_{x}^{0},\theta_{y}^{0})=\left\|\mathbf{R}^{\star}_{x}\mathbf{R}_{y}^{\star}\mathbf{E}(\mathbf{R}^{\star}_{x}\mathbf{R}_{y}^{\star})^{\top}(\mathbf{X}_{i}-\mathbf{t}^{\star}\mathbf{e}^{\top})+\mathbf{t}^{\star}\mathbf{e}^{\top}-\mathbf{X}_{i}\mathbf{P}^{\star}\right\|_{\text{F}}^{2}$ for each initialization $(\theta_{x}^{0},\theta_{y}^{0})$. Then, we find the average error $e(\theta_{x}^{0},\theta_{y}^{0})=\frac{1}{5000}\sum_{i=1}^{5000}e_{i}(\theta_{x}^{0},\theta_{y}^{0})$. Here, $\mathbf{R}_{x}$ and $\mathbf{R}_{y}$ are defined as follows.

[TABLE]

In Fig. 13, we show the average error $e(\theta_{x}^{0},\theta_{y}^{0})$. We observe that if the initialization $(\theta_{x}^{0},\theta_{y}^{0})$ is far away from the global optimum $(0^{\circ},0^{\circ})$, then the error is very high. As the distance between the initialization angles $(\theta_{x}^{0},\theta_{y}^{0})$ and the global optimum angles $(0^{\circ},0^{\circ})$ decrease, the error $e(\theta_{x}^{0},\theta_{y}^{0})$ remains approximately constant and suddenly drops to near zero after a particular distance. This indicates that, if the initialization angles are within a particular distance from the global optimum, then our approach always find the global optimum solution. This empirical result concurs with the result we already proved in Theorems 5 and 6.

VI Conclusion

In this work, we have developed a theory for establishing the correspondences between the mirror symmetric points in $\mathbb{R}^{d}$. We, further, determine the reflection symmetry transformation in a volumetric set of points in $\mathbb{R}^{d}$ containing a perturbed reflection symmetry pattern using optimization on Riemannian manifold. We have shown that our method is robust to a significant amount of perturbation and achieves 100% accuracy for no perturbation. We have further shown the significance of this work by detecting reflection symmetry in real images and comparing with state-of-the-art methods. The proposed approach is particularly suitable for detecting reflection symmetry of objects in applications where obtaining a robust local descriptor is highly challenging. The linear assignment problem is a time consuming step which restricts us to apply it on the large point sets. However, a proper sampling method can be employed to reduce the size of the point set without losing the symmetry present in the point set. We believe that the fundamental theory and algorithm developed in this work will pave the way for researchers to exploit them for scenarios where estimating feature descriptors is a challenging task.

Our approach detects single reflection symmetry plane of an object. Consider the third row of Fig. 11 in which there are multiple reflection symmetry planes present. In such cases, the detected reflection symmetry plane will be the one to which the initialized plane is the closest. For example, in the third row of Fig. 11, we have shown both the reflection symmetry planes detected depending on different initializations. This may not be a proper way of detecting multiple symmetries, though this is an interesting direction. We would like to extend our approach for the detection of multiple reflection symmetry planes of a symmetric object exhibiting multiple symmetries or a point cloud containing more than one symmetric objects.

VII Appendix

A1. Euclidean gradient of the function $\bar{f}$ with respect to the variable $\mathbf{t}$ (Equation (10))

We write the cost function as follows.

[TABLE]

We note that

[TABLE]

Therefore, we have (the terms which are not functions of $\mathbf{t}$ are not shown)

[TABLE]

[TABLE]

Now taking the derivative with respect to $\mathbf{t}$ we have,

[TABLE]

[TABLE]

[TABLE]

Here we have that $(\mathbf{I}_{d}-\mathbf{T}\mathbf{E}\mathbf{T}^{\top})\mathbf{T}\mathbf{E}\mathbf{T}^{\top}=-(\mathbf{I}_{d}-\mathbf{T}\mathbf{E}\mathbf{T}^{\top})$. Therefore,

[TABLE]

A2. Euclidean gradient of the function $\bar{f}$ with respect to the variable $\mathbf{R}_{j}$ (Equation (11))

Let us consider the cost function as defined in equation (5) (in main manuscript):

[TABLE]

Now, let us define $\mathbf{T}=\prod\limits_{u=1}^{d-1}\mathbf{R}_{u}$, $\mathbf{U}=\mathbf{X}-\mathbf{t}\mathbf{e}^{\top}$ and $\mathbf{V}=\mathbf{X}\mathbf{P}-\mathbf{t}\mathbf{e}^{\top}$.

Then the cost function becomes.

[TABLE]

Here we note that $\mathbf{T}\mathbf{E}\mathbf{T}^{\top}\mathbf{T}\mathbf{E}\mathbf{T}^{\top}=\mathbf{I}_{d}$, therefore

[TABLE]

Now taking the classical gradient of $\bar{f}$ with respect to $\mathbf{R}_{j}$ we have. (We follow [63] for the necessary properties.)

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Where

[TABLE]

A3. The Riemannian gradient of the function $f$ with respect to the variable $\mathbf{t}$ (Equation (12))

Using the definition, as defined in main paper, of Riemannian gradient $\boldsymbol{\xi}_{\mathbf{t}}(\mathbf{t})$ of the function $f$ with respect to the variable $\mathbf{t}$ we have

[TABLE]

[TABLE]

A4. The Riemannian gradient of the function $f$ with respect to the variable $\mathbf{R}_{j}$ (Equation (13))

Using the definition, as defined in main paper, of Riemannian gradient $\boldsymbol{\xi}_{\mathbf{R}_{j}}(\mathbf{R}_{j})$ of the function $f$ with respect to the variable $\mathbf{R}_{j}$ we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Therefore,

[TABLE]

[TABLE]

[TABLE]

A5. The Riemannian Hessian of the function $f$ with respect to $\mathbf{R}_{j}$ (Equation (14))

Next, we determine the Riemannian Hessian of the function $f$. In order to determine the $j^{\text{th}}$ component $\text{Hess}_{\mathbf{R}_{j}}(f(\mathbf{R,t}))[\mathbf{R}_{j}\Omega_{j}]=\mathbb{P}_{\mathbf{R}_{j}}(\text{D}\boldsymbol{\xi}_{\mathbf{R}_{j}}(\mathbf{R}_{j})[\mathbf{R}_{j}\mathbf{\Omega}_{j}])$, of the Riemannian Hessian, we first find the classical derivative $\text{D}\boldsymbol{\xi}_{\mathbf{R}_{j}}(\mathbf{R}_{j})[\mathbf{R}_{j}\mathbf{\Omega}_{j}]$ of the Riemannian gradient $\boldsymbol{\xi}_{\mathbf{R}_{j}}(\mathbf{R}_{j})$ in the direction $\mathbf{R}_{j}\mathbf{\Omega}_{j}$ and then we apply the projection operator $\mathbb{P}_{\mathbf{R}_{j}}$. Now using Equation 19 we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here,

[TABLE]

[TABLE]

Now

[TABLE]

[TABLE]

[TABLE]

The first term is equal to

[TABLE]

The second term is equal to

[TABLE]

Therefore,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here $[.,.]$ is the Lie bracket and defined as $[\mathbf{U},\mathbf{V}]=\mathbf{UV}-\mathbf{VU}$ for any two matrices $\mathbf{U}$ and $\mathbf{V}$.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A6. The Riemannian Hessian of the function $f$ with respect to $\mathbf{t}$ (Equation (15))

In a similar way, we determine the second component, $\mathbb{P}_{\mathbf{t}}(\text{D}\boldsymbol{\xi}_{\mathbf{t}}(\mathbf{t})[\boldsymbol{\eta}_{\mathbf{t}}])$, of the Riemannian Hessian. Since $\mathbb{R}^{d}$ is a vector space we have $\mathbb{P}_{\mathbf{t}}(\text{D}\boldsymbol{\xi}_{\mathbf{t}}(\mathbf{t})[\boldsymbol{\eta}_{\mathbf{t}}])=\text{D}\boldsymbol{\xi}_{\mathbf{t}}(\mathbf{t})[\boldsymbol{\eta}_{\mathbf{t}}]$

[TABLE]

[TABLE]

Therefore

[TABLE]

A7. Steps of Theorem 5

Showing the fact

[TABLE]

Now

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Since,

[TABLE]

[TABLE]

[TABLE]

we have

[TABLE]

[TABLE]

[TABLE]