Registering Image Volumes using 3D SIFT
and Discrete SP-Symmetry

Laurent Chauvin, William Wells III and Matthew Toews
L. Chauvin and M. Toews are with the Ecole de Technologie Superieure, Montreal, QC H3C 1K3, Canada (email<EMAIL_ADDRESS>matt.toews@gmail.com).
W. Wells is with the Harvard Medical School, Boston, MA 02115, USA and the Massachusetts Institute of Technology, Boston, MA 02139, USA (email: sw@bwh.harvard.edu).
Manuscript submitted May 15, 2022.

Journal of  Class Files, Vol. 14, No. 8, May 2022
Chauvin et al.: Image Registration based on 3D Feature Properties

This paper proposes to extend local image features in 3D to include invariance to discrete symmetry including inversion of spatial axes and image contrast. A binary feature sign $s \in \{-1,+1\}$ is defined as the sign of the Laplacian operator $\nabla^2$, and used to obtain a descriptor that is invariant to image sign inversion $s \rightarrow -s$ and 3D parity transforms $(x,y,z)\rightarrow(-x,-y,-z)$, i.e. SP-invariant or SP-symmetric. SP-symmetry applies to arbitrary scalar image fields $I: R^3 \rightarrow R^1$ mapping 3D coordinates $(x,y,z) \in R^3$  to scalar intensity $I(x,y,z) \in R^1$, generalizing the well-known charge conjugation and parity symmetry (CP-symmetry) applying to elementary charged particles. Feature orientation is modeled as a set of discrete states corresponding to potential axis reflections, independently of image contrast inversion. Two primary axis vectors are derived from image observations and potentially subject to reflection, and a third axis is an axial vector defined by the right-hand rule. Augmenting local feature properties with sign in addition to standard (location, scale, orientation) geometry leads to descriptors that are invariant to coordinate reflections and intensity contrast inversion. Feature properties are factored in to probabilistic point-based registration as symmetric kernels, based on a model of binary feature correspondence. Experiments using the well-known coherent point drift (CPD) algorithm demonstrate that SIFT-CPD kernels achieve the most accurate and rapid registration of the human brain and CT chest, including multiple MRI modalities of differing intensity contrast, and abnormal local variations such as tumors or occlusions. SIFT-CPD image registration is invariant to global scaling, rotation and translation and image intensity inversions of the input data.

Local Features, 3D SIFT, Discrete Symmetry, CP-symmetry, Charge conjugation, Sign inversion, Parity transform, Contrast Inversion, Image Registration, Multiple Modalities, right-hand rule.

§ INTRODUCTION

3D structure is often described in terms of localized phenomena or features. In classical physics, fine-scale matter in the atomic model is characterized via the Schrodinger equation <cit.>, where the symmetry of the Laplacian operator ensures that atomic properties remain invariant to changes in the viewpoint of the observer. Similarly, elementary 3D image features may be identified in generic image volumes via fundamental mathematical operators, e.g. the 3D Laplacian-of-Gaussian operator <cit.>, following from scale-space theory and the 2D scale-invariant feature transform (SIFT) <cit.>. These are used in highly robust and efficient memory-based applications, e.g. detection <cit.>, registration <cit.>, segmentation <cit.>, indexing <cit.> in a manner invariant to transformations including translation, rotation and scaling, and are applicable in generic contexts with no training or calibration procedures.

Our work here seeks to advance the understanding of 3D keypoints observed in image volumes from the notion of symmetry, specifically discrete symmetry in feature orientation. A geometrical object may be said to exhibit a symmetry if its observed properties remain unaffected by (or are invariant to) a group of transforms, i.e., a symmetry group. In particle physics, discrete symmetries include axis reflections and also simultaneous charge conjugation (inversion), 3D parity transformation and time reversal (CPT-symmetry). The CPT transform models the transition of a charge-associated particle to its anti-particle, e.g. a proton to an anti-proton, and is currently a topic of high interest <cit.> as symmetry violations have been observed in special cases including charge conjugation and parity (CP) transforms <cit.>. Nevertheless in image analysis, discrete shape symmetries including reflections <cit.> and image contrast inversion <cit.> are typically represented as separate, unrelated phenomena, typically in 2D coordinate space where 3D parity does not apply. Indeed, image processing and graphics research focuses largely on representations and challenges stemming from 2D data, including photographs, surface meshes or projections, and achieving invariance to viewpoint changes. Despite the compelling analogy between local features and the atomic model, the notion of a local feature sign, analogous to a particle charge, has not been proposed nor linked to the 3D parity transform in the context of image processing.

Illustrating the geometrical properties of 3D features (colored spheres) extracted from an MRI image of the human brain. Properties include location $\vect{x}$, orientation $\Theta$, scale $\sigma$, we introduce a binary feature sign $s = \{-,+\}$ based on the sign of the Laplacian-of-Gaussian.

Our primary contribution is the first model of discrete symmetry including intensity inversion in the context of local feature processing as shown in Figure <ref>. We note that a CP-like transform occurs in the context of image processing, where a local feature observed in different imaging modalities may exhibit simultaneous contrast inversion and gradient reversal. We refer to this as a sign inversion and parity transform (SP-transform), which applies generally to arbitrary scalar image fields $I: R^3 \rightarrow R^1$ mapping 3D coordinates $(x,y,z) \in R^3$ to a scalar image intensity $I(x,y,z) \in R^1$, and we seek SP-invariance. Local features are endowed with a binary sign $s \in \{-1,+1\}$ and a set of discrete orientation states, allowing feature descriptors and image registration to be computed in a manner invariant to reflections and intensity contrast inversions, thus achieving SP-invariance. Volumetric data allows physical matter and local image features to be observed in situ within isotropic 3D space and free from projective distortion, thus facilitating the analogy between local image features and atoms in terms of charge. Local feature orientation is defined by a 3D rotation matrix or coordinate reference frame, where two primary axes are estimated from dominant image gradient directions and the third vector is defined by the cross product, ensuring feature reference frames are restricted to $SO(3)$ despite potential parity transform due to contrast inversion. Four discrete orientation states are defined, corresponding to inversions of the two primary axis vectors, due to a combination of symmetric image gradient patterns and intensity contrast inversion. Gradient-based descriptors may be inverted via sign change, independently of orientation state.

Our secondary contribution is to incorporate invariant feature properties into a probabilistic point-based registration framework. The variability of feature scale, orientation and location between images is modeled via an exponential kernel function, approximating a zero-mean noise model, invariant to feature reflection and intensity contrast inversion. Experiments adapt our model to the well-known Coherent Point Drift (CPD) <cit.> algorithm, and demonstrate that our SIFT-CPD model of enhanced feature geometry leads to faster and more accurate inter-subject registration, in a diverse variety of image data including multiple T1 and T2 weighted MRI modalities of the human brain, tumors and CT images of the human torso. We also demonstrate that our model may be used to achieve SP-invariance independently of the method used to estimate feature orientation, e.g. via principal component analysis or maximum gradient directions.

The remainder of this paper describes related work, method and experiments. Throughout the paper, we aim to provide an accessible yet complete presentation of related concepts in 3D geometry, image processing and particle physics.

§ RELATED WORK

Our work focuses on characterizing 3D features extracted from volumetric images, and extending invariance to account for discrete symmetry due to both axis reflection and intensity contrast inversion.

§.§ Local Image Features

In computer vision and image processing, informative image structures are typically identified via mathematically-defined filters or convolution operators. First order gradient $\nabla I$ operators <cit.> may be used to identify corner-like structures, in a manner invariant to translation and rotation. Second-order derivative operators such as the Laplacian $\nabla^2 I$ may be used to identify blob-like structures, and to extend feature invariance to scale or size based on Gaussian scale-space theory<cit.>, as in the seminal scale-invariant feature transform (SIFT) <cit.>. The gradient and Laplacian operators may be regarded as minimal filter representations required to localize singular patterns in a scalar image $I: \mathcal{R}^D\rightarrow\mathcal{R}^1$ defined over a $D$-dimensional coordinate space. Extended representations may be used to achieve affine invariance <cit.>, or make use of alternative representations, eg. non-linear scale spaces <cit.>, information theory <cit.>, integral images for efficiency<cit.>. Over-complete filters banks, ie. convolutional neural networks (CNNs), may be trained for the purpose of identifying point correspondences <cit.>, possibly from training labels generated from fundamental operators <cit.>.

Local features incorporate descriptors of the image content in a region or receptive field. The Gaussian scale-space <cit.> may be viewed as a CNN with a single Gaussian channel per layer <cit.>, where the scale or receptive field size is defined explicitly by the Gaussian standard deviation $\sigma$, and is sampled in 3 logarithmic increments between $[\sigma,2\sigma]$ prior to subsampling. The receptive fields of deep network activations are a function of the filter sizes and depth <cit.>, and may be relatively large and dependent on the image contents due to non-linear operations including max pooling and ReLu <cit.>. Traditional descriptors have made use of gradients $\nabla I$, including gradient orientation histograms SIFT <cit.>, binary comparisons <cit.> and variants such as rank ordering <cit.>, RootSIFT <cit.>. Descriptors may be learned via triplet loss and gradient analysis <cit.>. While GPU-based deep learning is the defacto state-of-the-art for image classification <cit.>, variants of the traditional gradient representations remain competitive with learned variants <cit.>, specifically for matching images of non-planar objects <cit.> and image retrieval <cit.>.

The SIFT algorithm has been generalized to 3D volumetric images. A body of work investigates SIFT-like detectors in 3D video coordinates $(x,y,t)$ defined by 2D space $x,y$ and 1D time $t$ <cit.>. In the context of 3D space $x,y,z$, applications include object detection <cit.>, medical image analysis <cit.>, segmentation <cit.>, alignment <cit.>, image stitching <cit.>, large-scale indexing <cit.>, population studies <cit.>. Most recently, the Jaccard distance between feature sets was introduced to automatically flag errors in large public training MRI datasets <cit.>, and to identify family members from brain MRI <cit.>, capabilities facilitated by highly efficient feature indexing. We adopt the 3D SIFT-Rank method <cit.> using a GPU-optimized implementation <cit.> first used in the context of brain MRI analysis <cit.>.

§.§ Image Registration

Registration is a fundamental image processing task, and seeks to identify a coordinate transform $T: \Omega_1 \rightarrow \Omega_2$ between the coordinate systems $\Omega_1 \in R^3$, $\Omega_2 \in R^3$ of the same object or scene, often in three spatial dimensions, based on a pair of observed images $(I_1$,$I_2)$. Registration generally makes use of intensity and geometry information <cit.>, our work focuses on properties of local feature geometry, specifically location, orientation, scale and sign.

Point-based registration approaches are most generally applicable to local features, where the loss function minimizes the distance between pairs of points or points and a model. Examples of point cloud registration algorithms <cit.> include iterative algorithms such as the Iterative Closest Point (ICP) algorithm <cit.> where points contribute uniformly to a solution, or the Coherent Point Drift (CPD) algorithm <cit.> where points contribute according to a probabilistic weighting. 3D SIFT keypoint correspondences have been used to achieve point-based registration the context of image-guided neurosurgery <cit.>, including non-rigid image registration via thin plate splines <cit.>, finite element methods<cit.> and the CPD algorithm <cit.>, however these have made use of keypoint locations and not orientation and scale properties as we propose.

Feature-based registration methods generally consider points with properties beyond simple location. Points on a surface model may be endowed with properties such as mass <cit.> or charge <cit.>, however these properties tend to be assigned algorithmically, eg. uniformly assigned across points, and not derived from the image content itself. Local properties including orientation, scale and affine deformation may be used to enhance point-based registration <cit.>. Smooth deformations may be computed in a manner consistent with local rigid or affine reference frames <cit.>. The variability of geometrical misalignment is related to feature scale $\sigma$ <cit.>, consistent with uncertainty due to Gaussian blur. The localization accuracy of SIFT correspondences has been compared to that of manual human labeling in 2D and 3D  <cit.>, showing similar accuracy, where experts preferred automatic SIFT correspondences in 80% cases <cit.> in the case of brain imaging. Our 3D SIFT-CPD algorithm proposed in the following section accounts for feature properties of orientation and scale, in addition to point locations.

Multi-modal image registration is a major challenge, where intensities may vary locally in a non-linear manner between modalities. Training may be used to approximate an intensity mapping for a specific domains, e.g. MRI and CT <cit.>, however this mapping may not be functional or stationary throughout the image. In the general case with no specific training domain or data, intensity-based registration must adopt statistical similarity measures such as mutual information <cit.>. In the case of 2D image keypoints, invariance to local contrast inversion may be achieved by transforming the image intensity into a contrast-invariant format, including phase congruency <cit.> or Laplacian <cit.> images. As contrast inversion leads to image gradient reversal, attempts have been made to reverse aspects linked to the gradient, including descriptor orientation
<cit.>. Combinations of multiple keypoints <cit.> including self-similarity <cit.> may be adopted. In preliminary work, we proposed to account for 3D gradient reversal <cit.>, here we extend this to include reflections and SP-symmetry via the notion of a binary sign $s \in \{-1,+1\}$.

§.§ Discrete Symmetry

In a physical system, a symmetry refers to a property that remains unchanged under set of transforms, i.e. an invariant property. The symmetry group is the Lie group of transforms under which a geometrical object such as a particle is invariant <cit.>. A symmetry may be described as continuous or discrete. Continuous symmetry pertains to continuous parameters of pose and scale, eg. translation and orientation relative to the reference frame of an image acquisition device, and are analogous to a consequence of Noether's theorem for Lagrangian mechanical systems <cit.>. Continuous symmetry has been investigated in the computer vision literature in terms of differential invariance <cit.>. Classical neural networks such as CNNs are generally invariant to translations but not to reflections or scale changes <cit.>. Robustness to certain transforms may be improved via regularization <cit.>. Recently models of invariance or equivariance in trainable neural networks have emerged, including SO(3) <cit.> rotations, 3D rigid transforms<cit.> non-linear transforms <cit.>, probabilistic symmetries <cit.> and gauge theory <cit.>. Much of this work has been limited to 2D image modalities, i.e. surfaces or projective images, and has not considered the 3D notion of SP-symmetry including a binary sign $s \in \{-1,+1\}$.

A discrete symmetry generally involves discrete displacements in location and/or orientation, e.g. translations, rotations, reflections, glide reflections or repetitive crystal structure. Charge conjugation and parity symmetry (CP-symmetry) <cit.>, and generally space-time reversal in CPT-symmetry including time relating to the transition of a particle to its anti-particle <cit.> are important aspects of discrete symmetry relating to 3D reflection of elementary physical particles. CP-symmetry is closely related to reflection symmetry, which exists in a function exhibiting one or more planes of symmetry, i.e. a  reflection about a single axis $x=0$ such that $f(x)=f(-x)$, or in a 3D system, the parity coordinate transform $f(x,y,z)=f(-x,-y,-z)$, both characterized in 3D by a rotation matrix $R$ with determinant $det(R)=-1$. The effect of the parity transform may be seen by first representing $f(x)$ in terms of $f(x) = f_s(x) + f_a(x)$, the sum purely symmetric and anti-symmetric image components $f_s(x)$ and $f_a(x)$ defined by $f_s(-x)=f_s(x)$ and $f_a(-x)=-f_a(x)$. Symmetric or even functions $f_s(x)$ include the Gaussian and Laplacian-of-Gaussian operators, wave functions associated with bosonic particles such as the photon, the cosine function, etc. Anti-symmetric or odd functions $f_a(x)$ include the gradient operator $\nabla I$, wave functions associated with fermionic particles such as electrons, and the sine function, etc. As the gradient is an anti-symmetric operator, an image contrast inversion $-\nabla I(x)=\nabla I(-x)$ is equivalent to an axis reflection. This phenomenon is analogous to charge conjugation and parity (CP) symmetry, and describes the transition of a charge-associated particle to its antiparticle, e.g. electrons and positrons. Notably, CP-symmetry was thought to be an inviolable in nature until the discovery of special cases beginning with the composite meson particle <cit.>. In our work here, charge is equivalent to image sign or contrast, which may be inverted, and we propose SP-invariance to register images of differing modality.

Image processing methods specifically relating to our work have analyzed discrete shape symmetry (ie. reflections) <cit.> and image contrast inversion <cit.> as separate, unrelated phenomena. Our work is the first to consider these within the unified framework of discrete SP-symmetry, including a binary sign and discrete states of local reference frame orientation. Our method generalizes to various contexts where orientation is estimated from dominant local image gradients $\nabla I$. For example, Toews et al. identify maxima in a 3D spherical gradient histogram <cit.> such that $\vect{\hat{\theta}}_1 =  \underset{\vect{\hat{\theta}}}{\mathrm{argmax}}~|\nabla I\cdot \hat{\theta}|$,  $\vect{\hat{\theta}}_2 =  \underset{\vect{\hat{\theta}}}{\mathrm{argmax}}~|\vect{\hat{\theta}} \times (\nabla I \times \vect{\hat{\theta}}_1)|$, with a third axis defined by the cross product $\vect{\hat{\theta}}_3 = \vect{\hat{\theta}}_1 \times \vect{\hat{\theta}}_2$. Rister et al. compute the eigenvectors of the local gradient structure tensor, ie. the $3\times3$ gradient correlation matrix <cit.>. We demonstrate that both such representations may be modeled in terms of discrete SP-symmetry.

§ METHOD

Our method provides a means of characterizing 3D image features in a manner invariant to discrete symmetry transforms, including symmetric image patterns and sign inversion and parity (SP-symmetry) transforms due to image contrast inversion. We begin by describing the properties of a single local image feature, which we augment to include a binary sign and a set of discrete axis reflections. This allows us to compute a descriptor that is invariant to reflections and changes in sign due inversion of intensity contrast. We then define a kernel function $\mathcal{K}(f_n,f_m)$ to quantify the similarity of a pair of features $(f_n,f_m)$ potentially arising from the same underlying structure in two different images. Finally, this kernel is  incorporated into probabilistic image registration in order to estimate a transform $T: f_m \rightarrow f_n$ aligning the coordinate systems of two images.

§.§ Single-Feature Properties

A feature is a distinctive spherical region localized in 3D image space, as illustrated in Figure <ref>. An individual feature $f=\{g,\vect{a},s\}$ is characterized within the image by its geometry $g$, a descriptor of the image appearance $\vect{a}$ surrounding the feature and a binary sign $s$. The feature geometry $g=\{\Theta,\sigma,\vect{x}\} \in SO(3)\times R^+ \times R^3$ is a scaled coordinate reference frame in 3D defined by 7 parameters. These include a scale $\sigma \in \mathbb{R}^+$, a coordinate location $\vect{x} \in \mathbb{R}^3$ and a reference frame $\Theta \in SO(3)$. The reference frame may be defined as a set $\Theta = \{\vect{\hat{\theta}}_1,\vect{\hat{\theta}}_2,\vect{\hat{\theta}}_3\}$ of orthonormal unit vectors $\vect{\hat{\theta}}_i, i \in [1,2,3]$ where $\vect{\hat{\theta}}_i \cdot \vect{\hat{\theta}}_j = 0, i \ne j$. Note equivalently 3D rotation could make use of alternative representations, eg. $SU(2)$ Pauli matrices or unit quaternions. Feature appearance $\vect{a}\in \mathbb{R}^M$ is an M-dimensional descriptor of the image, resampled according to geometry $g$. The sign $s \in \{-1,1\}$ is used to achieve invariance to image contrast inversion.

Let $I: \mathbb{R}^3 \rightarrow \mathbb{R}^1$ be a scalar image sampled over 3D coordinates where $I(\vect{x})$ represents a voxel measurement at location $\vect{x}$. Let $I(\vect{x},\sigma)=I(\vect{x})*G(\vect{x},\sigma)$ be a scale-space defined by convolution of the image with a Gaussian filter $G(\vect{x},\sigma)$ of parameter $\sigma$ <cit.>. Features are detected as points in scale-space $\{\vect{x_m},\sigma_m,s_m\}$ at which the magnitude of the scale-normalized Laplacian operator is maximized, similarly to the original SIFT method <cit.>, here $m$ is a feature index. Furthermore, they are endowed with a binary sign as follows.

\begin{align}
\{\vect{x_m},\sigma_m\} &=  \underset{\vect{x},\sigma}{\mathrm{ local~argmax}}~|\nabla^2 I(\vect{x},\sigma)|, \notag \\
s_m &= {\mathrm{sign}} \left( \nabla^2 I(\vect{x}_m,\sigma_m) \right).
\end{align}

Scale-space loci $\{\vect{x_m},\sigma_m,s_m\}$ represent generic spherical regions in which the divergence of the local gradient field $\nabla I(\vect{x},\sigma)$ is maximized. Feature sign $s_m$ is novel to our characterization, and is defined as the sign of the Laplacian-of-Gaussian.

Feature orientation $\Theta$ is derived from the image gradient $\nabla I$ surrounding the feature origin $\vect{x}$. Orientation axis vectors $\vect{\hat{\theta}}_1,\vect{\hat{\theta}}_2,\vect{\hat{\theta}}_3$ may be determined according dominant gradient orientations within a window $w(\vect{x},\sigma) \in [0,1]$ centered upon $\vect{x}$ with spatial extent proportional to scale $\sigma$. As $\Theta \in SO(3)$, the determinant $det(\Theta)=1$ is positive, ie. no reflections. Axes are thus defined by two primary axis vectors $\vect{\hat{\theta}}_1$ and $\vect{\hat{\theta}}_2$, with the third constrained according to a handedness convention, ie. the right-hand rule via the cross product $\vect{\hat{\theta}}_3 = \vect{\hat{\theta}}_1 \times \vect{\hat{\theta}}_2$. Furthermore, we assume that axis vectors may be ordered uniquely according to the gradient magnitude along their directions, i.e. $[| \nabla I\cdot\vect{\hat{\theta}}_1|] > [| \nabla I\cdot\vect{\hat{\theta}}_2|]$, where $[ | \nabla I\cdot\vect{\hat{\theta}}_i|] = \int w(\vect{x},\sigma)| \nabla I(\vect{u})\cdot\vect{\hat{\theta}}_i| du$ is the expected value of the gradient magnitude along axis $\vect{\hat{\theta}}_i$.

Our primary contribution is a descriptor that is invariant to discrete SP-transforms, i.e. is SP-symmetric, based on binary sign $s$. We first note that the sign of the orientation axes $\vect{\hat{\theta}}_1,\vect{\hat{\theta}}_2$ may generally be ambiguous due to a variety of factors including bilateral symmetry of the image pattern or its local gradient distribution (e.g. an ellipsoid or rectangular box), an image contrast inversion (e.g. different imaging modalities), or the algorithm used to estimate orientation (e.g. principal component analysis has an inherent eigenvector sign ambiguity). Sign ambiguity may be represented by maintaining a discrete set of possible orientation reference frames $\Theta$. We consider four potential reference frame states $\{\pm \vect{\hat{\theta}_1},\pm \vect{\hat{\theta}_2}\}$ (with $\vect{\hat{\theta}_3}=\vect{\pm\hat{\theta}_1} \times \pm\vect{\hat{\theta}_2}$) resulting from axis inversion, i.e. due to either intensity contrast inversion or pattern symmetry, which are equivalent to the identity in addition to rotations of $\pi$ about each of three axes. By defining the third axis vector $\vect{\hat{\theta}_3}$ using the cross product, we ensure feature orientation follows the right-hand rule. Figure <ref> illustrates the four possible orientation states, including an initial reference frame and rotations of $\pi$ about each of the 3 axes.

Feature sign $s$ is used to generate a novel contrast inversion-invariant descriptor $\vect{a}_m$ for each orientation state, based on the image content surrounding location $\vect{x_m}$ at scale $\sigma_m$ and orientation $\Theta_m$. Note that the geometry $g$ defines a 7-parameter similarity transform from the local feature reference frame to a canonical reference frame, including rotation, translation and scaling. This this is used to resample the image according to a characteristic reference frame $\hat{I}(\vect{u}) = I(\sigma_m\Theta_m\vect{u}+\vect{x_m})$ that is invariant to global similarity transforms of the image. A variety of descriptors may be computed, eg. spherical harmonics <cit.>, we adopt a computationally efficient variant of the gradient orientation histogram <cit.>. As in <cit.>, 3D space surrounding the feature origin $\vect{x}_m$ is quantized uniformly into eight discrete octants $\vect{r}=[\pm 1,\pm 1,\pm 1]$, and the image gradient at each point $\vect{u}$ into eight discrete symmetric directions $\vect{\phi}=[\pm 1,\pm 1,\pm 1]$. These are combined into an $8\times8=64$-element appearance descriptor, by accumulating the gradient magnitude into a histogram over spatial location and orientation

\begin{align}
a(\vect{r},\vect{\phi}) = \sum_{\vect{u} \in \vect{r}} \| \nabla \hat{I}(\vect{u})\cdot \vect{\phi}  \|[\vect{\phi}=s\vect{\phi}_{\vect{u}}],
\label{eq:descriptor}
\end{align}

where $\vect{\phi}_{\vect{u}}= \underset{\vect{\phi}}{\mathrm{argmax}}
\{ \nabla \hat{I}(\vect{u}) \cdot s\vect{\phi} \}$ is the direction of maximum gradient at image location $\vect{u}$, which is inverted in the case of $s=-1$ by the Iverson bracket $[\vect{\phi}=s\vect{\phi}_{\vect{u}}]$  evaluating to 1 upon equality and 0 otherwise. Thus descriptor orientation elements $\vect{\phi}$ are invariant to reversal due to sign change $a(\vect{r},\vect{\phi}) = a(\vect{r},s\vect{\phi})$, and parity transforms of orthants $\vect{r} \rightarrow -\vect{r}$ are accounted for by feature orientation states $\{\pm \vect{\hat{\theta}_1},\pm \vect{\hat{\theta}_2}\}$, as illustrated in Figure <ref>.

Illustrating an appearance descriptor $a(r,\phi)$ and four discrete orientation states given principal axes $\hat{\theta}_1$, $\hat{\theta}_2$ (left to right) and feature sign $s$ (upper and lower). The descriptor content is shown for an example octant ($r_{+++}$) and gradient orientation ($\phi_{++-}$). Axis orientation is defined independently of sign (left to right), and negative sign $s=-1$ inverts the descriptor gradient (lower row) achieving invariance to contrast inversion.

§.§ Two-Feature Observations

Image registration requires identifying pairs of features $(f_n,f_m)$ arising from the same underlying structure in different images. To do this, we propose a kernel function $K(f_n,f_m)\in [0,1]$ to quantify the similarity of two features $f_n$ and $f_m$ potentially sampled from the same distribution. We construct $K(f_n,f_m)$ as a product of squared exponential kernels, ensuring that it is positive, symmetric $K(f_m,f_n)=K(f_n,f_m)$, and proportional to a Gaussian density. It may be expressed as the product of three factors

\begin{align}
\mathcal{K}(g_{n},g_{m}) = &\mathcal{K}_{\sigma}(\sigma_{n},\sigma_{m})\mathcal{K}_{\Theta}(\Theta_{n},\Theta_{m})\mathcal{K}_{x}(\vect{x_{n}},\vect{x_{m}}), \label{eq:geo_kernel}
\end{align}

quantifying the variability in scale, orientation and location, respectively. As in the notions of linear and angular momentum from classical physics, deviations of location and orientation are modeled as orthogonal and independent.

The first factor in Equation (<ref>),  $\mathcal{K}_{\sigma}(\sigma_{n},\sigma_{m})$, is a kernel penalizing the log scale difference:

\begin{equation}
\mathcal{K}_{\sigma}(\sigma_{n},\sigma_{m}) = \exp \left( - \| \log \sigma_n - \log \sigma_m \|^2 \right),
\label{eq:scale_kernel}
\end{equation}

and is proportional to a log normal density about mean $\log 1=0$ with unit variance $1$. The log maps multiplicative variations in feature scale on the range $\sigma \in [0,\infty]$ to additive variations on the range $\log \sigma \in [-\infty,\infty]$, which may be modelled as symmetric and Gaussian.

The second factor $\mathcal{K}_{\Theta}(\Theta_{n},\Theta_{m})$ quantifies the variability of feature orientations based on the angular displacement of the axes:

\begin{equation}
\mathcal{K}_{\Theta}(\Theta_{n},\Theta_{m}) = \exp \left(-3 + \sum_{i=1,2,3} \vect{\hat{\theta}_{in}}\cdot \vect{\hat{\theta}_{im}} \right).
\label{eq:orient_kernel}
\end{equation}

In Equation (<ref>), the scalar product
$\vect{\hat{\theta}_{in}}\cdot \vect{\hat{\theta}_{im}}=\cos(\vect{\hat{\theta}_{in}}- \vect{\hat{\theta}_{im}})$ quantifies the angular separation between axis unit vectors $\vect{\hat{\theta}_{in}}$ and  $\vect{\hat{\theta}_{im}}$ on the range $[-1,1]$. This is equivalent to modeling the angular difference between features independently in each axis $i$ as proportional to a von Mises density over the unit circle.

The third factor $\mathcal{K}_{x}(\vect{x_{n}},\vect{x_{m}})$ is a kernel that quantifies the similarity of feature locations based on their squared displacement:

\begin{equation}
\mathcal{K}_{x}(\vect{x_{n}},\vect{x_{m}}) = \exp \left( -\frac{||\vect{x_n} - \vect{x_m}||^2}{k~\sigma_{n}\sigma_{m}+\sigma_{T}^2} \right).
\label{eq:loc_kernel}
\end{equation}

where the denominator normalizes the spatial displacement $||\vect{x_n} - \vect{x_m}||^2$ via a linear function of the product of feature scales $\sigma_{n}\sigma_{m}$.
The first term $k~\sigma_{n}\sigma_{m}$ embodies intrinsic variance due to the scale of the observed feature, where $\sigma_{n}\sigma_{m}$ represents the square of the geometric mean of feature scales and $k$ is a positive proportionality constant. This allows quantifying the deviation in location in a manner invariant to keypoint scale. The second term $\sigma_{T}^2$ represents the minimum achievable variance in mapping $T$ given the particular imaging context, which dominates in the case of small-scaled features, i.e. $k~\sigma_{n}\sigma_{m} < \sigma_{T}^2$. Note that in previous work <cit.> $k=1$ and $\sigma_{T}=0$, here we use $k$ and $\sigma_{T}$ to estimate a more accurate linear relationship between scale and observed variability. These are empirically set to provide a generally useful weighting, e.g. here we use $k=12$, $\sigma_T^2=200$.

§.§ Multi-Feature Registration

Registration seeks a transform $\mathcal{T}$ mapping the coordinate system of a fixed image $I_f$ to that of a moving image $I_m$, here from sets of features extracted in the fixed $\mathcal{F}=\{f_{n}\}$ and moving $\mathcal{M}=\{f_{m}\}$ images. The transform maps the properties of features from one coordinate space to the next $\mathcal{T}:\vect{x_m} \rightarrow \vect{x_n$}, and here is taken to be a global similarity transform, followed by independent zero-mean deviations for individual features. Note that the similarity transform is a 7-parameter representation  $\mathcal{T}=\{d\Theta,d\sigma,d\vect{x}\} \in SO(3)\times R^+ \times R^3$ equivalent that of an individual feature.

Let $g_m'$ represent the geometry of feature $f_m$ as transformed by $\mathcal{T}$, ie. $g_m'=T \circ g_m$, note that $T$ may act on location, scale and orientation parameters in a consistent manner. We propose using our kernel function $\mathcal{K}(g_n,g_m')$ to account for feature properties including scale and orientation in a standard point-based registration framework. Our kernel lends itself naturally to a probabilistic algorithm such as Coherent Point Drift (CPD) in Algorithm~\ref{alg:cpd}, which seeks a maximum likelihood solution based on the Expectation-Maximization (EM) algorithm~\cite{dempster1977maximum}. The expectation step (E) intend estimates the log-likelihood based on the current parameters, while the maximization step (M) maximizes the log-likelihood over the parameters using the result of the E step. CPD registration is based on a probability map $p_{mn}$ between each pair of points $\{f_n,f_m\}$, and takes the form of a softmax function. We propose to bias this probability map based on a kernel function $\mathcal{K}(g_n,g_m')$ as follows:
\begin{equation}
\label{eq:pij}
\end{equation}
where $\eta = (2\pi\lambda^2)^{D/2}\frac{w}{1-w}\frac{M}{N}$ is a constant ensuring the denominator $p_{mn}$ is non-zero, where $w$ is parameter accounting for the relative probability of a uniform background density not associated with any specific feature. In Equation~\eqref{eq:pij}, the exponential expression with parameter $\lambda$ is a variance parameter that is reduced iteratively to converge to a solution, and a constant kernel function $\mathcal{K}(g_n,g_m')=1$ is identical to the original CPD algorithm~\cite{myronenko2010point}. Thus in our proposed SIFT-CPD, two exponential factors operate on point coordinates $\vect{x}$ including Equation~\ref{eq:loc_kernel}.

Iterative algorithms such as EM must be initialized within a neighborhood of the solution in order to converge correctly. We initialize registration via a global 3D Hough transform between feature sets $\mathcal{F}$ and $\mathcal{M}$ in 3D~\cite{Toews2013a} analogously to the classic 2D SIFT algorithm~\cite{Lowe2004DistinctiveKeypoints}, where all features vote independently as to the predicted transform, after which a most likely transform is identified\footnote{The Hough transform method was originally proposed to track particles in bubble chamber photographs~\cite{hough1959machine}}. Here, each feature $f_{n} \in \mathcal{F}$ votes for a transform $\mathcal{T}_{nm}: g_n \rightarrow g_m$ mapping the geometry of feature $f_{n}$ to that of feature $f_{m}$, where $f_{m}=\underset{f}{\mathrm{argmin}}~\| \vect{a}_n-\vect{a} \|$ is the nearest neighbor of $f_{n}$ in terms of the Euclidean distance between appearance descriptors. Matches are identified between all features and discrete orientation states $\{\pm\hat{\theta}_1,\pm\hat{\theta}_2\}$ including parity transforms, and in a manner invariant to intensity contrast inversions due to the binary sign $s \in \{-1,+1\}$ in Equation~\eqref{eq:descriptor}. A dominant transform $\mathcal{T}^* \in \{\mathcal{T}_{nm}\}$ is then identified such that it is consistent with the largest number of matches, i.e. inlier correspondences. A pair of features $(f_{n},f_{m})$ represents an inlier of the transform $\mathcal{T}^*$ if their transform $\mathcal{T}_{nm}: g_n \rightarrow g_m$ differs by less than a threshold $\|\mathcal{T}_{nm}-\mathcal{T}^*\| < Thres$ as follows:
\begin{align}
\mathcal{T}^* =
\underset{T}{\mathrm{argmax}}~
|\{T: \| \mathcal{T}_{nm}-\mathcal{T} \| < Thres\} |.
\label{eq:thres}
\end{align}
The thresholding operation in Equation~\eqref{eq:thres} may be defined as the logical conjunction of thresholds independently applied in rotation, scaling and squared displacement
\begin{align}
\| \mathcal{T}_{nm}-\mathcal{T} \| < Thres
\rightarrow &  \\
&\cos(\vect{\hat{\theta}}_i-\vect{\hat{\theta}}_{inm}) < ~\epsilon_{cos\theta} \notag \\
& \wedge~|\log \sigma-\log\sigma_{nm}| < ~\epsilon_{log\sigma}  \notag \\
&\wedge~\|\vect{x}-\vect{x_{nm}}\|^2/\sigma\sigma_{nm} < \epsilon_{x/\sigma}, \notag
\end{align}
using thresholds $(\epsilon_{cos\theta},\epsilon_{log\sigma},\epsilon_{x/\sigma})$. These may be set generously to apply in a wide variety contexts, here we use
$(\epsilon_{cos\theta},\epsilon_{log\sigma},\epsilon_{x/\sigma}) = (0.7,log 1.5, 0.25)$. The Hough transform process may be implemented efficiently using hash tables and the mean-shift clustering algorithm~\cite{comaniciu2002mean}. Our modified CPD algorithm is provided in Algorithm~\ref{alg:cpd}, and the $solve$ functions for similarity transforms is found in Algorithm~\ref{alg:solve_rigid} in Appendix~\ref{app:solve}.

\begin{algorithm}
\SetKwInOut{InHere}{Inputs}
\SetKwInOut{OutHere}{Outputs}
\SetKwInOut{Init}{Initialization}
\SetKwInOut{Opti}{EM Optimization}
\DontPrintSemicolon

\InHere{ $\mathcal{F}=\{f_{n}\}$~\qquad\text{Fixed Features} \\  $\mathcal{M}=\{f_{m}\}$\qquad\text{Moving Features}
\OutHere{ $\mathcal{T}:\mathcal{M} \rightarrow \mathcal{F}$~\quad\text{Transform}

\Init{$0\leq w \leq 1$\;
$\lambda^2 = \frac{1}{DNM}\sum_{n=1}^{N}\sum_{m=1}^{M}||\vect{x_n}-\vect{x_m}||^2$\;
$\mathcal{T} \leftarrow \mathcal{T}^*$}

\EMOpti{
\EStep{
\MStep{
$\{\mathcal{T},\lambda^2\} = solve(\mathcal{F},\mathcal{M},P)$
The aligned point set is $\mathcal{T}(\mathcal{M})$
\\ The correspondence probability is $P$.

\caption{Probabilistic SIFT-CPD registration algorithm, adapted from~\cite{myronenko2010point} to include kernel $\mathcal{K}(g_n,g'_m)$.}
\label{alg:cpd}
\end{algorithm}
\vspace{\baselineskip}
\noindent
%The $solve$ function in Algorithm~\ref{alg:cpd} can estimate a rigid (Algorithm~\ref{alg:solve_rigid}) or affine (Algorithm~\ref{alg:solve_affine}) transform, described in Appendix~\ref{app:solve}.

%However, given the density of keypoints in each dataset, most keypoints won't have a counterpart in the second dataset, increasing the noise in the CPD algorithm, hence reducing accuracy and computational speed. By preserving only matching keypoints between two datasets, we can greatly increase the signal-to-noise ratio, speeding up the converge to a solution. To do so, the search of an initial Hough transform can be leveraged to create a sparse subset of inliers through a simple threshold on keypoints spatial location, orientation and scale, hence reducing the computational complexity of the subsequent SIFT-CPD algorithm, as seen in Figure~\ref{fig:synth_data_brain}.

\section{Experiments}

We hypothesize that our kernel function $\mathcal{K}(g_n,g_m)$ incorporating feature scale and orientation will improve the accuracy of registration beyond point information alone, and that feature orientation state and binary sign will provide invariance to intensity contrast inversion and reflections, independently of the method used to identify feature orientation $\Theta$. Experiments validate these hypotheses in the context of 3D medical imaging data, which provide the opportunity to investigate features arising from realistic, diverse anatomical structure observed across subjects and 3D imaging modalities.

Three sets of experiments are performed. The first is based on synthetic similarity transforms, and establishes baseline accuracy of methods in the case of known ground truth. The second consists of inter-subject image registration trials, first between brain MRI and then between chest CT volumes of different people, with confounds including brain tumors and contrast variations due to multiple T1 and T2-weighted MRI modalities, demonstrating superior algorithm performance in a diverse set of contexts. Variability is estimated by repeating trials following synthetic transforms of images and multiple subjects. The third investigates orientation state changes in registration trials for different feature orientation estimation methods, i.e., principal components~\cite{rister2017volumetric} or maximum orientation directions~\cite{Toews2013a}.

Note that the true mapping $T$ between the anatomies of different subjects is generally non-linear and may not exist throughout the image, due to aspects of anatomy specific to individuals or in the case of occlusion or missing structure. As mentioned, features are assumed to follow a transform $T$ that is globally linear followed by random feature-specific deviations as specified by our kernel $\mathcal{K}(g_n,g_m)$. For inter-subject registration, ground truth is established from inlier correspondences identified via the Hough transform, which are visually validated for correctness and rejected or manually adjusted if necessary. This follows the work of~\cite{toews2007statistical,Machado2018Non-rigidMatching}, where automatic and manually labeled correspondences exhibit the same error range. The average Point Registration Error (PRE) measure is used to quantify registration performance, based on the sum of 3D point differences between the registration solution and ground truth.

In all experiments, feature sets are extracted using a GPU implementation of the SIFT-Rank algorithm for computational efficiency~\cite{carluer2021gpu}, and registered using five point-cloud registration methods: ICP with 20 and 100 iterations (respectively ICP20 and ICP100), the original CPD algorithm using feature centers $\vect{x}$, SIFT-CPD using full feature geometry and our kernel $\mathcal{K}(g_n,g_m)$, and SIFT-CPD* using only inlier features from the Hough transform. All algorithms are initialized to the Hough transform solution prior to iterative registration in order to ensure a fair comparison.

\subsection{Synthetic Image Registration}

A preliminary experiment was first performed establish the baseline performances of algorithms against known ground truth, here synthetic similarity transformations of a single image. A T1w brain image from the Human Connectome Project (HCP)~\cite{VanEssen2012} dataset was selected (see Table~\ref{tab:experiments} for more details), and 100 different synthetic transforms were applied, where each was generated by a random rotation about each axis on an angle range of $\rho \in \pm [10^{\circ},30^{\circ}]$ followed by a translation $\delta \in \pm[0~mm,10~mm]$. The registration error in rotation and translation was evaluated relative to the known transform applied, and the results are presented in Table~\ref{tab:samesubject_registration}. As expected, the lowest error is achieved for SIFT-CPD* (inliers alone), validating this method as a baseline error. SIFT-CPD and CPD perform similarly, with slightly less error for SIFT-CPD. ICP results are generally poor, failing to converge in 8 and 21 cases for 20 and 100 iterations.

%Note that CPD, SIFT-CPD and SIFT-CPD Optimized performed well in term of accuracy compared to ICP, where SIFT-CPD Optimized lead to the lowest error while running respectively 5x and 22x faster than SIFT-CPD and CPD.\\~\\

\begin{table*}
\large
\begin{center}
\resizebox{2\columnwidth}{!}{
\bgroup
\def\arraystretch{1.8}
\begin{tabular}{|cc|ccc|ccc|cc|cc}
\hline
\multicolumn{2}{|P{3.5cm}|}{\multirow{2}{*}{\backslashbox[3.9cm]{Algorithm}{Error}}}&\multicolumn{3}{|c|}{\textit{Rotation Error (deg)}}&\multicolumn{3}{|c|}{\textit{Translation Error (mm)}}&\multicolumn{2}{|c|}{\textit{PRE (mm)}}&\multicolumn{2}{|c|}{\textit{SSD ($e^{11}$)}}\\
\cline{3-12}
\multicolumn{2}{|P{3.5cm}|}{}&\multicolumn{1}{|c|}{X-Axis}&\multicolumn{1}{|c|}{Y-Axis}&\multicolumn{1}{c|}{Z-Axis}&\multicolumn{1}{|c|}{X-Axis}&\multicolumn{1}{|c|}{Y-Axis}&\multicolumn{1}{|c|}{Z-Axis}&\multicolumn{2}{|c|}{Mean $\pm$ Std.Dev.}&\multicolumn{2}{|c|}{Mean $\pm$ Std.Dev.}\\
\thickhline
\multicolumn{2}{|P{3.5cm}|}{ICP20} & \multicolumn{1}{|c|}{$25.7 \pm 22.0$} & \multicolumn{1}{|c|}{$16.8 \pm 14.8$} & \multicolumn{1}{|c|}{$24.8 \pm 17.2$} & \multicolumn{1}{|c|}{$14.2 \pm 19.6$} & \multicolumn{1}{|c|}{$12.8 \pm 29.1$} & \multicolumn{1}{|c|}{$17.4 \pm 29.0$} & \multicolumn{2}{|c|}{$16.37 \pm 8.41$} &
\multicolumn{2}{|c|}{$7.62 \pm 1.03$}\\
\hline
\multicolumn{2}{|P{3.5cm}|}{ICP100} & \multicolumn{1}{|c|}{$8.1 \pm 18.9$} & \multicolumn{1}{|c|}{$6.3 \pm 18.8$} & \multicolumn{1}{|c|}{$8.5 \pm 18.5$} & \multicolumn{1}{|c|}{$9.0 \pm 24.0$} & \multicolumn{1}{|c|}{$9.2 \pm 32.8$} & \multicolumn{1}{|c|}{$8.5 \pm 23.2$} & \multicolumn{2}{|c|}{$8.16 \pm 17.15$} &
\multicolumn{2}{|c|}{$4.89 \pm 2.31$}\\
\hline
\multicolumn{2}{|P{3.5cm}|}{CPD} & \multicolumn{1}{|c|}{$0.06 \pm 0.03$} & \multicolumn{1}{|c|}{$0.06 \pm 0.03$} & \multicolumn{1}{|c|}{$0.06 \pm 0.03$} & \multicolumn{1}{|c|}{$0.8 \pm 0.6$} & \multicolumn{1}{|c|}{$1.3 \pm 1.1$} & \multicolumn{1}{|c|}{$0.7 \pm 0.5$} & \multicolumn{2}{|c|}{$1.81 \pm 0.88$} &
\multicolumn{2}{|c|}{$1.38 \pm 0.39$}\\
\hline
\multicolumn{2}{|P{3.5cm}|}{\textbf{SIFT-CPD}} & \multicolumn{1}{|c|}{$\bf{0.02 \pm 0.01}$} & \multicolumn{1}{|c|}{$\bf{0.01 \pm 0.01}$} & \multicolumn{1}{|c|}{$\bf{0.02 \pm 0.01}$} & \multicolumn{1}{|c|}{$\bf{0.9 \pm 0.5}$} & \multicolumn{1}{|c|}{$\bf{0.8 \pm 0.6}$} & \multicolumn{1}{|c|}{$\bf{0.5 \pm 0.4}$} & \multicolumn{2}{|c|}{$\bf{1.05 \pm 0.32}$} &
\multicolumn{2}{|c|}{$\bf{1.43 \pm 0.33}$}\\
\hline
\multicolumn{2}{|P{3.5cm}|}{\textbf{SIFT-CPD*}} & \multicolumn{1}{|c|}{$\bf{0.02 \pm 0.01}$} & \multicolumn{1}{|c|}{$\bf{0.01 \pm 0.01}$} & \multicolumn{1}{|c|}{$\bf{0.02 \pm 0.01}$} & \multicolumn{1}{|c|}{$\bf{0.9 \pm 0.4}$} & \multicolumn{1}{|c|}{$\bf{0.9 \pm 0.4}$} & \multicolumn{1}{|c|}{$\bf{0.4 \pm 0.3}$} & \multicolumn{2}{|c|}{$\bf{0.68 \pm 0.15}$} &
\multicolumn{2}{|c|}{$\bf{0.96 \pm 0.39}$}\\
\hline
\end{tabular}
\egroup
\end{center}
\caption{Comparison of five registration algorithms based on 100 synthetic transforms of a single T1w MRI brain volume. Error is listed terms of individual rotation and translation axes, overall point registration error (PRE), and the sum of squared intensity differences (SSD) following registration.
%CPD, SIFT-CPD and SIFT-CPD Optimized performed well, but SIFT-CPD Optimized performed the registration 5x faster than SIFT-CPD, and 22x faster than CPD.
%\\ * 21 values have been removed as they failed to converge.\\ ** 8 values have been removed as they failed to converge.
\label{tab:samesubject_registration}
\end{table*}

\subsection{Inter-subject Registration}

Four sets of experiments quantified the accuracy of inter-subject registration, i.e. registration of images of different individuals. The goal was to evaluate the accuracy and speed of registration in a diverse set of contexts, including different aspects of anatomy (brain, chest), abnormalities (tumor) and modalities (T1w and T2w MRI, CT). Information regarding datasets is provided in Table~\ref{tab:experiments}. Note that minimum registration error is generally non-zero due to differences in anatomy, and the goal is to identify a robust transform minimizing the PRE for sets of inlier correspondences.

\begin{table*}
\begin{center}
\resizebox{1.5\columnwidth}{!}{
\bgroup
\def\arraystretch{1.8}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\textbf{Dataset} & \textbf{Anatomy} & \textbf{Modality} & \textbf{Resolution} & \textbf{Voxel Size (mm)}& \textbf{Features} & \textbf{Inliers}\\
\thickhline
\cite{VanEssen2012} & Brain & T1w MRI &$260\times260\times311$&$0.7\times0.7\times0.7$&$6000$&$350$\\
\hline
\cite{regan2011genetic} & Chest & CT &$280\times280\times235$&$1.2\times1.2\times1.2$&$3000$&$70$\\
\hline
\cite{Menze2015TheBRATS} & Brain Tumor & T1w MRI &$240\times240\times155$&$1.0\times1.0\times1.0$&$1200$&$70$\\
\hline
\cite{VanEssen2012} & Brain & T1w-T2w MRI &$260\times260\times311$&$0.7\times0.7\times0.7$&$5000$&$115$\\
\hline
\end{tabular}
\egroup
\end{center}
\caption{Inter-subject registration experimental data.}
\label{tab:experiments}
\end{table*}

\subsubsection{T1w Brain Images}

The following experiment quantified inter-subject registration accuracy between brain images of a pair of different subjects. As in the previous experiment, one image was fixed, while the other was subjected to a random transform from the previously mentioned ranges $\rho$ and $\delta$. A pair of monozygotic twin subjects was selected, which generally share more inlier correspondences in comparison to unrelated subjects due to higher neuroanatomical similarity~\cite{Chauvin2020NeuroimageRelatives}. The results for 100 random transforms are shown in Figure~\ref{fig:synth_data_brain}, these are generally consistent across all registration experiments. ICP variants which weight all points uniformly do not correctly converge (red and orange). CPD and SIFT-CPD (green and light blue) both converge consistently to acceptable results in term of PRE, both lower than that the initial robust Hough transform (purple). SIFT-CPD is consistently lower in terms of error and computation time in comparison to CPD, demonstrating that geometrical properties improve estimation, in the case where the majority of features may be outliers with no valid correspondence between images. Overall, the SIFT-CPD* (dark blue) leads to the lowest error and computation time, as it refines the original set of Hough transform inliers.

\begin{figure}[!h]
\centering
\includegraphics[width=0.5\textwidth]{figures/fig_Registration_Merged_Twin_iter100_Hough_ICP_CPD_AllPoints_SIFTCPD_InliersPoints2.pdf}
\caption{Inter-subject registration performance for healthy brain T1w MRI.
%The estimated original error (black cross) represents the estimation of the registration error from original HCP images (registered via a rigid transform through FSL) and the transformed error (red cross) the error after applying a randomly generated transform.
%ICP has been evaluated with 20 and 100 iterations, respectively ICP20 and ICP100. The registration error has been evaluated using the Keypoint-based Registration Error (KRE) between two keypoint sets, representing the two images. The estimated original error (black cross) represents the estimation of the registration error from original HCP images (registered via a rigid transform through FSL), the transformed error (red cross) the error after applying a randomly generated transform, and color circles the KRE after registration with different methods. We notice that ICP does not perform well (with ICP20 relatively unstable in term of performances), while SIFT-CPD performs very well with a low KRE and a much greater stability.
\label{fig:synth_data_brain}
\end{figure}

% TODO: Add label x-axis

% \begin{figure*}[!h]
%     \centering
%     \includegraphics[width=0.8\textwidth]{figures/fig_Boxplot_Registration_Time_Twin_iter100_ICP20_ICP100_CPD_SIFTCPD_SIFTCPDOpti.pdf}
%     \caption{Computational time in milliseconds (log scale) necessary for different algorithms to register two point sets from monozygotic twins from HCP for 100 random transforms. }
%     \label{fig:computation_time_brain}
% \end{figure*}

% \begin{table}[h]
%     \setlength\extrarowheight{3pt}
%     \begin{center}
%      \resizebox{\columnwidth}{!}   {
%         \begin{tabular}{ccccccc}
%             \hline
%             & \textbf{No Registration} & \textbf{Hough} & \textbf{ICP20} & \textbf{ICP100} & \textbf{CPD} & \textbf{SIFT-CPD} \\
%             \hline
%             \textbf{MZ} & $\bm{0.9999}$ & $0.9998$ & $0.9991$ & $0.9998$ & $\bm{0.9999}$ & $\bm{0.9999}$ \\
%             \textbf{DZ} & $0.9391$ & $0.9314$ & $0.9338$ & $\bm{0.9434}$ & $0.9371$ & $0.936$ \\
%             \textbf{FS} & $0.8989$ & $0.9021$ & $0.8883$ & $\bm{0.9048}$ & $\bm{0.9048}$ & $0.9043$ \\
%             \hline
%         \end{tabular}
%         }
%     \end{center}
%     \caption{AUC with Appearance+Geometry kernel for each family relationship, k=200, after different registration method}
%     \label{tab:my_label}
% \end{table}

% \begin{table}[h]
%     \setlength\extrarowheight{3pt}
%     \begin{center}
%      \resizebox{\columnwidth}{!}   {
%         \begin{tabular}{cccccc}
%             \hline
%             & \textbf{Hough} & \textbf{ICP20} & \textbf{ICP100} & \textbf{CPD} & \textbf{SIFT-CPD} \\
%             \hline
%             \textbf{X-Axis (deg)} & $1.5741$ & $15.0531$ & $0.6906$ & $\bm{0.0264}$ & $0.0315$ \\
%             \textbf{Y-Axis (deg)} & $1.4151$ & $11.0593$ & $0.2870$ & $\bm{0.0355}$ & $0.0476$ \\
%             \textbf{Z-Axis (deg)} & $1.6011$ & $13.0147$ & $0.7310$ & $\bm{0.0227}$ & $0.0294$ \\
%             \textbf{Average (deg)} & $1.5301$ & $13.0423$ & $0.5695$ & $\bm{0.0282}$ & $0.0362$ \\
%             \textbf{Time (ms)} & $1249.33$ & $\bm{70.27}$ & $348.81$ & $26670.83$ & $200.09$ \\
%             \hline
%         \end{tabular}
%         }
%     \end{center}
%     \caption{Rotation error (deg) and computation time (ms). CPD leads to better accuracy but at the cost of computation time, while SIFT-CPD leads to similar accuracy but approximately 130 times faster. ICP with 100 iteration is less accurate, but require 1.74 the computation time.}
%     \label{tab:my_label}
% \end{table}

\subsubsection{CT Chest Images}
In order to show that our method is not tied to specific contexts or imaging modalities, we performed inter-subject registration of CT chest images from the Chronic Obstructive Pulmonary Disorder Genetics (COPDGene) dataset~\cite{regan2011genetic}. Two different subject images in expiration breathing state were randomly selected (see Table~\ref{tab:experiments} for more details).
As in the previous experiment, we generated 100 random transform with the same $\rho$ and $\delta$ parameters, applied the transform to one image, and registered it using different registration approaches. The PRE for each registration is shown on Figure~\ref{fig:synth_data_lungs}. Although the computational complexity of SIFT-CPD is comparable to the complexity of CPD in Figure~\ref{fig:synth_data_lungs}, SIFT-CPD* is significantly faster, closer to the ICP20 performances.

% \begin{figure}
%     \centering
%     \includegraphics[width=0.75\textwidth]{figures/fig_Lung_Illustration.pdf}
%     \caption{Particles extracted from lung images with Chronic Obstructive Pulmonary Disorder (COPD). }
%     \label{fig:lungs_illustration}
% \end{figure}

\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{figures/fig_Registration_Merged_Lungs_iter100_CPD_ICP_SIFT_CPD_AllPoints_SIFTCPD_InliersPoints_HoughTimeAdded2.pdf}
\caption{Inter-subject registration performance for chest CT images. %The estimated original error (black cross) represents the estimation of the registration error from original COPD images and the transformed error (red cross) the error after applying a randomly generated transform.
%ICP has been evaluated with 20 and 100 iterations, respectively ICP20 and ICP100. The registration error has been evaluated using the Keypoint-based Registration Error (KRE) between two keypoint sets, representing the two images. The estimated original error (black cross) represents the estimation of the registration error from original COPD images, the transformed error (red cross) the error after applying a randomly generated transform, and color circles the KRE after registration with different methods. We notice that ICP does not perform well (with ICP20 relatively unstable in term of performances), while SIFT-CPD performs very well with a low KRE and a much greater stability. Circles around 175mm are registration that failed to converge and gave irrelevant registration transforms.
\label{fig:synth_data_lungs}
\end{figure}

% TODO: Add label x-axis
% \begin{figure*}
%     \centering
%     \includegraphics[width=1\textwidth]{figures/fig_Boxplot_Registration_Time_Lungs_iter100_ICP20_ICP100_CPD_SIFTCPD_SIFTCPDOpti.pdf}
%     \caption{Computational time in milliseconds (log scale) necessary for different algorithms to register two point sets from two different subjects from COPD after random transforms. }
%     \label{fig:computation_time_lungs}
% \end{figure*}

% \begin{figure}
%     \centering
%     \includegraphics[width=1\textwidth]{figures/fig_Boxplot_Family_Registrations_Geo_NoReg_Hough_ICP_CPD_CPDScale.pdf}
%     \caption{Geometrical weight before and after registration}
%     \label{fig:HCP_Registration}
% \end{figure}

\subsubsection{Abnormal Variations and Occlusion}
Practical image analysis requires robustly coping with occlusions, including a lack of one-to-one correspondence due to inter-subject variations including pathological tumor structure. Our method is based on local feature properties, and is thus particularly robust to occlusion and missing structure. We used 10 T1w MRI images of different subjects with tumors from the BraTS dataset~\cite{Menze2015TheBRATS} in order to evaluate registration in the presence of occlusion (see Table~\ref{tab:experiments} for more details). Nine subject images are randomly selected, and transformed 5 times each, and registered to the remaining subject image. Results are presented in Figure~\ref{fig:occlusion_results}.

%generated 100 random transforms (same parameters as in previous experiments) in addition of an occlusion at a random location with a radius $r \in [3mm, 45mm]$. A registration with and without the occlusion has been performed to evaluate the degradation in performances relatively to the size of the occlusion. An example of generated images with the largest possible occlusion (almost the size of an hemisphere) can be seen in Figure~\ref{fig:example_occlusion}. Results are presented in Figure~\ref{fig:occlusion_results}, and we observe virtually no degradation in the registration performances, up to 45mm radius occlusion.

\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{figures/fig_Registration_Merged_Brats_Hough_ICP_CPD_SIFTCPD_SIFTCPDOpti2.pdf}
\caption{Inter-subject registration performance for brains in T1w MRI, with occlusions due to tumors~\cite{Menze2015TheBRATS}.}
\label{fig:occlusion_results}
\end{figure}

\subsubsection{Multi-Modal Images}
To evaluate the registration performance of our approach on multi-modal images, a random pair of brain images was selected from HCP dataset, one with a T1w image and the other one with a T2w (see Table~\ref{tab:experiments} for more details). T1w and T2w images are acquired with different MRI parameters and highlight different tissue properties, where T2w imaging involves longer repetition time (TR) and time-to-echo (TE) parameters. The T2w image was randomly transformed and registered as in previous experiments, with 100 trials. Results are presented in Figure~\ref{fig:synth_data_multimodal_brain}. An example of corresponding features identified in different modalities is shown in Figure~\ref{fig:multimodal_matches} and are investigated in greater detail in the following experiment section.

\begin{figure}[!h]
\centering
\includegraphics[width=0.5\textwidth]{figures/fig_Registration_Merged_Twin_iter100_MultiModal2.pdf}
\caption{Inter-subject registration performance for healthy brain images between T1w and T2w MRI modalities~\cite{VanEssen2012}, exhibiting intensity contrast inversion.
%The estimated original error (black cross) represents the estimation of the registration error from original COPD images, the transformed error (red cross) the error after applying a randomly generated transform, and color circles the KRE after registration with different methods. We notice that ICP does not perform well (with ICP20 relatively unstable in term of performances), while SIFT-CPD performs very well with a low KRE and a much greater stability.
\label{fig:synth_data_multimodal_brain}
\end{figure}

\begin{figure}[!h]
\centering
\includegraphics[width=0.5\textwidth]{figures/fig_Matches_MultiModal.pdf}
\caption{Examples of features matched between T1w and T2w images, showing feature signs, note the sign changes do not occur everywhere.}
\label{fig:multimodal_matches}
\end{figure}

% \begin{figure*}
%     \centering
%     \includegraphics[width=1\textwidth]{figures/fig_Boxplot_Registration_Time_MultiModal_ICP20_ICP100_CPD_SIFTCPD_SIFTCPDOpti.pdf}
%     \caption{Computational time in milliseconds (log scale) necessary for different algorithms to register two point sets from T1w and T2w from two HCP random monozygotic twins.}
%     \label{fig:computation_time_multimodal}
% \end{figure*}

\subsection{Feature Sign and Orientation State}

Here we investigated the binary sign $s$ and orientation state $(\theta_1,\theta_2)$ of local features in registration experiments, using two different methods for estimating feature orientation $\Theta$: eigenvectors of the local gradient structure tensor matrix~\cite{rister2017volumetric} or maximum orthogonal gradient directions~\cite{Toews2013a}. Both methods have been demonstrated to achieve effective local feature correspondences between images of the same modality, but not different modalities where discrete orientation state changes may occur.

Figures~\ref{fig:states} show distributions of orientation state transitions of inlier correspondences between images. For same modality registration (a, b, c) there are relatively few state changes (main diagonal), however in the case of multi-modal registration (d) noticeable differences exist. In the case of orientation determined via principal component analysis (PCA) of gradients~\cite{rister2017volumetric}, approximately 40\% of all correspondences exhibit single-axis reflections (Figure~\ref{fig:states} d) rather than SP transforms, this is due to eigenvector sign ambiguity, which is accounted for in our algorithm. In the case of orientation determined by maximum gradient directions, correspondences generally exhibit full SP-transforms, i.e. sign inversion and parity (Figure~\ref{fig:states} d). Figure~\ref{fig:symmetries} shows examples of corresponding features and state transitions in the case of sign inversion, including identity (0-0), single axis reflections of major (0-1) and minor axes (0-2), and parity transforms (0-3).

\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{figures/fig_MatchingStates_All.pdf}
\caption{Orientation state transition histograms for principal gradient orientation estimation (see~\cite{rister2017volumetric}) and max gradient orientation estimation (see~\cite{Toews2013a}) for: a) T1w brain, b) T1w brain tumor~\cite{Menze2015TheBRATS}, c) CT chest~\cite{regan2011genetic}, d) Multi-modal T1w/T2w brain. States
described in Figure~\ref{fig:workflow}, are combinations of discrete axis reflections $\{\pm\hat{\theta}_1,\pm\hat{\theta}_2\}$ as 0:$\{++\}$,
3:$\{--\}$. Bins along the diagonal indicate no state change, off-diagonals indicate single-axis reflections and inverse diagonal indicates parity transform.
\label{fig:states}
\end{figure}

% \begin{figure}
%     \centering
%     \includegraphics[width=0.33\textwidth]{figures/fig_MatchingStates_T1_Brats_COPD_MultiModal_Normalized_Canonical.pdf}
%     \caption{Orientation state change histograms for  for: (a) T1 brain, (b) T1 brain tumor ~\cite{Menze2015TheBRATS}, (c) CT chest~\cite{regan2011genetic}, (d) T1 and T2 brain. The four states are the binary combinations of axis reflections $\{\pm\hat{\theta}_1,\pm\hat{\theta}_2\}$ as 0:$\{++\}$,
%     1:$\{+-\}$,
%     2:$\{-+\}$,
%     3:$\{--\}$, bins along the ascending diagonal indicate no state change, multi-modal contrast inversion generally leads to an axis inversion (d).}
%     \label{fig:states_max}
% \end{figure}

\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{figures/fig_MatchingStates_Symmetries2.pdf}
\caption{Examples of feature orientation state changes observed between T1w and T2w modalities. Red and blue arrows indicate the primary $\hat{\theta}_1$ and secondary $\hat{\theta}_2$ orientation axes projected onto an image slice plane. Circles indicate the feature location and scale in a slice plane, sign $s \in \{-1,+1\}$ is indicated by the white-black or black-white transitions. }
\label{fig:symmetries}
\end{figure}
% \begin{figure*}
%     \centering
%     \includegraphics[width=0.8\textwidth]{figures/fig_Occlusion.pdf}
%     \caption{Occlusion}
%     \label{fig:example_occlusion}
% \end{figure*}

\section{Discussion}

This paper proposed to extend local feature methods to account for SP-symmetry and to achieved 3D multi-modal image registration. We introduced the notion of discrete SP-symmetry in the of case features extracted from a scalar intensity image $I(\vect{x}) \in R^1$ over $\vect{x} \in R^3$ space. We proposed to model the geometrical properties of 3D image features in a manner invariant to SP-transforms observed in multi-modal image registration. We proposed to include a feature sign $s \in \{-,+\} = {\text sign}~\nabla^2~I(\vect{x},\sigma)$ analogous to electric charge, and maintain a discrete set of four orientation states including reflections of primary axes and parity transforms. Sign is used to invert gradient descriptor elements in the case of intensity contrast inversion, i.e. charge conjugation. Feature geometry and appearance descriptors are thus invariant to SP transforms, in addition to continuous rotations, scalings and translations as in the 3D SIFT method~\cite{Toews2013a,rister2017volumetric}.

%Features are extracted in scale and space extrema of Laplacian-of-Gaussian operator $\nabla^2~I(\vect{x},\sigma)$, and defined geometrically by a scale, 3D location and rotation matrix, such as $g=\{\sigma,\vect{x},\Theta\} \in \times R^+ \times R^3 \times SO(3)$. Rotation matrix axes $\Theta$ and appearance descriptors are derived from the gradient field $\nabla~I(\vect{x},\sigma)$ surrounding the feature.

We integrated feature properties into a well-known probabilistic point-based image registration framework, the CPD algorithm~\cite{myronenko2010point}, via a kernel function $\mathcal{K}(f_n,f_m)$, leading to a highly robust and stable SIFT-CPD registration algorithm. Image registration experiments are performed using a range of volumetric image modalities of the human brain and chest, showing that additional properties consistently improve the accuracy of registration in comparison to point locations alone. Registration error and computation time is consistently lower for our proposed SIFT-CPD vs. standard point-based CPD frameworks, indicating that feature scale and orientation improve the estimation of image-to-image transforms. An optimized version of SIFT-CPD functioning solely on sparse inlier feature correspondences exhibits the lowest error, refining the solution achieved via a robust but coarse initial Hough transform. Experiments also demonstrated how SP-symmetry may be applied despite the method used to estimate feature orientation estimation, here principal gradient orientations~\cite{rister2017volumetric} and maximum gradient orientations~\cite{Toews2013a}.  Invariant feature appearance descriptors are used to identify inliers and to initialize registration, there were not explicitly incorporated into the SIFT-CPD kernel function based solely on geometrical properties here, however this possibility is left for future work.

%While other 3D point registration methods have focused primarily on generic points constrained to lie on surface models of 3D objects~\cite{Golyanik2016GravitationalRegistration,Jauer2019EfficientClouds}, we focus on features located in volumetric image structure, strengthening the analogy between invariant image features and the atomic model.

%Image registration trials are performed on a variety of volumetric medical imaging contexts including MRI and CT modalities, human brain and chest anatomies, tumors. These demonstrate the general robust nature of the algorithm in the context of inter-subject registration, where the mapping between the anatomies of different people may be non-linear or non-existent. Our results here are based on estimating a linear similarity transform, this could be extended affine or deformable transforms.

%Our approach models multi-scale image information via an isotropic Gaussian diffusion process $I(\vect{x},\sigma)=I(\vect{x})*G(\sigma)$, where elementary Laplacian $\nabla^2 I(\vect{x},\sigma)$ and gradient $\nabla I(\vect{x},\sigma)$ operators may be evaluated over a range of spatial scales $\sigma$ related to diffusion time $t=\sigma^2$. As the Laplacian-of-Gaussian (LoG) operator shares the symmetric center-surround electric potential of atoms of the periodic table, it could thus be used to identify and characterize atoms in molecular arrangements, i.e. in volumes produced via density functional theory (DFT)~\cite{grimme2007density,hohenberg1964inhomogeneous}, and generally to model 3D structure at a variety of spatial scales.

We note that the theory of SP-symmetry we present applies generally to individual channels of convolutional neural networks, and could be used within the general neural network learning framework to achieve SP-invariance, in addition to continuous invariance~\cite{bardes2021vicreg,esteves2018learning,spezialetti2019learning,cohen2019gauge}. Here, 3D invariant features have been identified via unbiased, symmetric operators with no explicit training procedure, which could be used as is in general imaging contexts, or to train or validate domain-specific detectors via neural networks~\cite{Yi2016LIFT:Transform,Detone2018SuperPoint:Description}.\\~\\
All code required to reproduce our results may be obtained at \href{https://github.com/3dsift-rank/SIFT-CPD}{https://github.com/3dsift-rank/SIFT-CPD}.

%where the ReLu~\cite{nair2010rectified} activation is typically used to suppress negativity.

%Overall, all results presented here in the different experiments are consistent, with ICP20 achieving best results in terms of computational speed but with inaccurate precision in term of registration. While ICP100 achieve better slightly better accuracy, it's computational complexity is also higher. Precision-wise, the Hough transform usually performs equally or better than ICP100, but at a fraction of the computational speed. While CPD converge globally to an accurate solution, its complexity makes it unsuitable for large scale image registration. On the other hand, while SIFT-CPD performs better than CPD computationally and accuracy-wise, the gain achieved are marginal. However, the optimization performed in SIFT-CPD allows it to significantly outperform CPD in terms of precision, at a speed equivalent to the one of ICP20, making this algorithm particularly suitable for large scale image registration for it's accuracy, speed and robustness to local variations as well as contract inversion.

%Experiments show that these modifications achieved faster and more accurate registration that the original implementation of the CPD algorithm. Applied in the context of medical imaging, it allows to register a large number of volumetric images with a reduced complexity, while maintaining a high robustness to occlusions. Our method also has the advantages of being independent of the type of images (e.g. brain, lungs) or modalities (e.g. T1w, T2w, CT).

%Although this work focused on rigid transform, it can easily be extended to affine transform, as well as other point-cloud registration algorithm such as the Iterative Closest Point (ICP), the Robust Point Matching (RPM). With a bit more modification, it should also be possible to extend the Bayesian CPD algorithms to features to achieve deformable registration, or using a thin-plate deformation as in~\cite{frisken2020comparison} . The sign of features could also be used to normalize the gradient orientation to achieve multi-modality registration as in~\cite{toews2013feature}.

% use section* for acknowledgment
\ifCLASSOPTIONcompsoc
% The Computer Society usually uses the plural form
\section*{Acknowledgments}
\else
% regular IEEE prefers the singular form
\section*{Acknowledgment}
\fi
This work was supported by NIH grant P41EB015902 (NAC), the Quebec (FRQNT) New Researchers Startup Program and the Canadian National Sciences and Research Council (NSERC) Discovery Grant. Data were provided by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.

% Can use something like this to put references on a page
% by themselves when using endfloat and the captionsoff option.
\ifCLASSOPTIONcaptionsoff
\newpage
\fi

% trigger a \newpage just before the given reference
% number - used to balance the columns on the last page
% adjust value as needed - may need to be readjusted if
% the document is modified later
% The "triggered" command can be changed if desired:

% references section

% can use a bibliography generated by BibTeX as a .bbl file
% BibTeX documentation can be easily obtained at:
% http://mirror.ctan.org/biblio/bibtex/contrib/doc/
% The IEEEtran BibTeX style support page is at:
% http://www.michaelshell.org/tex/ieeetran/bibtex/
\bibliographystyle{IEEEtran}
% argument is your BibTeX string definitions and bibliography database(s)
\bibliography{mendeley,cpd}
% <OR> manually copy in the resultant .bbl file
% set second argument of \begin to the number of references
% (used to reserve space for the reference number labels box)

% biography section
% If you have an EPS/PDF photo (graphicx package needed) extra braces are
% needed around the contents of the optional argument to biography to prevent
% the LaTeX parser from getting confused when it sees the complicated
% \includegraphics command within an optional argument. (You could create
% your own custom macro containing the \includegraphics command to make things
% simpler here.)
%\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{mshell}}]{Michael Shell}
% or if you just want to reserve a space for a photo:

% if you will not have a photo at all:

% insert where needed to balance the two columns on the last page with
% biographies

% You can push biographies down or up by placing
% a \vfill before or after them. The appropriate
% use of \vfill depends on what kind of text is
% on the last page and whether or not the columns
% are being equalized.

% Can be used to pull up biographies so that the bottom of the last one
% is flush with the other column.

\newpage
\appendices

\section{Probabilistic Registration Details}
\label{app:solve}

The probabilistic CPD algorithm is based on a $solve$ function defined for rigid and affine transforms. The $solve\_rigid$ algorithm used in this work is presented in Algorithm~\ref{alg:solve_rigid}, while the $solve\_affine$ can be found in Algorithm~\ref{alg:solve_affine}, from~\cite{myronenko2010point}. The $solve\_rigid()$ function takes as inputs a fixed and a moving point sets $\mathcal{F}$ and $\mathcal{M}$ along with their corresponding probability map $P=\{p_{ij}\}$, and outputs the covariance $\lambda^2$, and a scaled rigid transform with parameters rotation matrix $R$, translation vector $t$, and magnification b. The $solve\_affine()$ function returns an affine transform including invertible affine matrix $B$ and translation vector $t$.

\begin{algorithm}
\DontPrintSemicolon
\SetKwFunction{FMain}{solve\_rigid}
\SetKwProg{Fn}{}{:}{}
\Fn{\FMain{$\mathcal{F}$, $\mathcal{M}$, $P$}}{
$\mathit{N}_{\bf{P}}=\bf{1}^{\mathit{T}}\bf{P1},~ \mu_\mathcal{F}=\frac{1}{\mathit{N}_{\bf{P}}}\bm{\mathcal{F}}^{\mathit{T}}\bf{P}^{\mathit{T}}\bf{1},~\mu_\mathcal{M}=\frac{1}{\mathit{N}_{\bf{P}}}\bm{\mathcal{M}}^{\mathit{T}}\bf{P1},$\;
$\bf{A}=\bm{\hat{\mathcal{F}}}^{\mathit{T}}\bf{P}^{\mathit{T}}\bm{\hat{\mathcal{M}}}$, compute SVD of $\bf{A}$ such as $\bf{A}=\bf{U}\textit{SS}\bf{V}^{\mathit{T}},$\;
$\bf{R}=\bf{UCV}^{\mathit{T}}$, where $\bf{C}=\mathrm{d(1,..,1,det(\bf{UV}^{\mathit{T}}))},$\;
$\lambda^2 = \frac{1}{\mathit{N}_{\bf{P}}D}(tr(\bm{\hat{\mathcal{F}}}^{\mathit{T}}\mathrm{d(\bf{P}^{\mathit{T}}\bf{1})}\bm{\hat{\mathcal{F}}})-\mathnormal{b~tr(\bf{A}^{\mathit{T}}\bf{R}))}.$\;
\KwRet $\{b,\bf{R},\bf{t}\},\lambda^2$
\label{alg:solve_rigid}
\caption{Solve for rigid transform.}
\end{algorithm}

% that's all folks
\end{document}

%In the algorithm, $t$ is a translation vector, $R$ is a rotation matrix, $s$ is an isotropic scaling parameter.

\begin{algorithm}
\DontPrintSemicolon
\SetKwFunction{FMain}{solve\_affine}
\SetKwProg{Fn}{}{:}{}
\Fn{\FMain{$\mathcal{F}$, $\mathcal{M}$, $P$}}{
$\mathit{N}_{\bf{P}}=\bf{1}^{\mathit{T}}\bf{P1},~ \mu_\mathcal{F}=\frac{1}{\mathit{N}_{\bf{P}}}\bm{\mathcal{F}}^{\mathit{T}}\bf{P}^{\mathit{T}}\bf{1},~\mu_\mathcal{M}=\frac{1}{\mathit{N}_{\bf{P}}}\bm{\mathcal{M}}^{\mathit{T}}\bf{P1},$\;
$\lambda^2 = \frac{1}{\mathit{N}_{\bf{P}}D}(tr(\bm{\hat{\mathcal{F}}}^{\mathit{T}}\mathrm{d(\bf{P}^{\mathit{T}}\bf{1})}\bm{\hat{\mathcal{F}}})-\mathnormal{tr(\bm{\hat{\mathcal{F}}}^{\mathit{T}}\bf{P}^{\mathit{T}}\bm{\hat{\mathcal{M}}}\bf{B}^{\mathit{T}}))}.$\;
\KwRet $\{\bf{B},\bf{t}\},\lambda^2$
\label{alg:solve_affine}
\caption{Solve for affine transform.}
\end{algorithm}


22mm22mm25mm25mm

# Harrison center and products of sums of powers

Hua-Lin Huang, Lili Liao, Huajun Lu and Chi Zhang  School of Mathematical
Sciences, Huaqiao University, Quanzhou 362021, China<EMAIL_ADDRESS>School of Mathematical Sciences, Huaqiao University, Quanzhou 362021, China
<EMAIL_ADDRESS>School of Mathematical Sciences, Huaqiao University,
Quanzhou 362021, China<EMAIL_ADDRESS>Department of Mathematics,
Northeastern University, Shenyang 110819, China<EMAIL_ADDRESS>
###### Abstract.

This paper is mainly concerned with identities like

$(x_{1}^{d}+x_{2}^{d}+\cdots+x_{r}^{d})(y_{1}^{d}+y_{2}^{d}+\cdots
y_{n}^{d})=z_{1}^{d}+z_{2}^{d}+\cdots+z_{n}^{d}$

where $d>2,$ $x=(x_{1},x_{2},\dots,x_{r})$ and $y=(y_{1},y_{2},\dots,y_{n})$
are systems of indeterminates and each $z_{k}$ is a linear form in $y$ with
coefficients in the rational function field $\mathbbm{k}(x)$ over any field
$\mathbbm{k}$ of characteristic $0$ or greater than $d.$ These identities are
higher degree analogue of the well-known composition formulas of sums of
squares of Hurwitz, Radon and Pfister. We show that such composition
identities of sums of powers of degree at least $3$ are trivial, i.e., if
$d>2,$ then $r=1.$ Our proof is simple and elementary, in which the crux is
Harrison’s center theory of homogeneous polynomials.

###### Key words and phrases:

sum of powers, composition

###### 2010 Mathematics Subject Classification:

11E76, 15A69

Supported by NSFC 11971181 and 12131015.

##  1\. Introduction

Throughout the paper, let $d\geq 3$ be an integer and $\mathbbm{k}$ be a field
of characteristic $0$ or $>d.$

In 1898, Hurwitz [10] proved the well-known “$1$, $2$, $4$, $8$ Theorem” which
states that the only values of $n$ for which there is an identity of the type

(1.1)
$(x_{1}^{2}+x_{2}^{2}+\cdots+x_{n}^{2})(y_{1}^{2}+y_{2}^{2}+\cdots+y_{n}^{2})=z_{1}^{2}+z_{2}^{2}+\cdots+z_{n}^{2}$

where the $z_{k}$ are bilinear forms in $x=(x_{1},x_{2},\dots,x_{n})$ and
$y=(y_{1},y_{2},\dots,y_{n})$ are $n=1,2,4,8.$ Twenty years later, Hurwitz
[11] and Radon[16] proved independently that there is an identity

(1.2)
$(x_{1}^{2}+x_{2}^{2}+\cdots+x_{r}^{2})(y_{1}^{2}+y_{2}^{2}+\cdots+y_{n}^{2})=z_{1}^{2}+z_{2}^{2}+\cdots+z_{n}^{2}$

satisfying the previous conditions if and only if $r\leq\rho(n)$ where
$\rho(n)$ is the Hurwitz-Radon function defined as $\rho(n)=8a+2^{b}$ if
$n=2^{4a+b}n_{0}$ where $n_{0}$ is odd and $0\leq b\leq 3.$ In the 1960s,
Pfister [14] showed that if one allows $z_{k}$ to be linear in $y$ with
coefficients in $\mathbbm{k}(x)=\mathbbm{k}(x_{1},x_{2},\dots,x_{r}),$ then
there is an identity of form (1.1) if and only if $n=2^{m}$ is a power of $2.$
These beautiful square identities are nice motivation for the algebraic theory
of quadratic forms, and have close connection to various topics in algebra,
arithmetics, combinatorics, geometry and topology.

For a long time, people are wondering if there are similar identities of
products of sums of powers in higher degrees. For example, in 1975 Nathanson
[13] considered the following identity

(1.3)
$(x_{1}^{d}+x_{2}^{d}+\cdots+x_{n}^{d})(y_{1}^{d}+y_{2}^{d}+\cdots+y_{n}^{d})=z_{1}^{d}+z_{2}^{d}+\cdots+z_{n}^{d}$

requiring that the $z_{k}$ are polynomials in $x$ and $y$ with integer
coefficients. He applied simple arithmetic technique to show that such
identities cannot hold for $(n,d)=(3,3)$ and $n=2,d>2.$ He also conjectured
that there exist no identities like (1.3) for any $n>1$ and $d>2.$ In the
1980s, from his work on higher pythagoras numbers, Becker [1] proposed a much
stronger conjecture than the one of Nathanson’s: if the following identity

(1.4)
$(x_{1}^{d}+x_{2}^{d}+\cdots+x_{r}^{d})(y_{1}^{d}+y_{2}^{d}+\cdots+y_{s}^{d})=z_{1}^{d}+z_{2}^{d}+\cdots+z_{n}^{d}$

holds, where each $z_{k}\in\mathbbm{k}(x,y),$ then $rs\leq n.$ To the best of
our knowledge, so far very little seems to be known about this conjecture, but
Gude made some progress in his thesis [4]. See also Shapiro’s book [18,
Chapter 16] for a brief summary of Gude’s results.

In this paper, we are mainly concerned about identities of the following type

(1.5)
$(x_{1}^{d}+x_{2}^{d}+\cdots+x_{r}^{d})(y_{1}^{d}+y_{2}^{d}+\cdots+y_{n}^{d})=z_{1}^{d}+z_{2}^{d}+\cdots+z_{n}^{d}$

where each $z_{k}$ is a linear form in $y$ with coefficients in
$\mathbbm{k}(x).$ These identities are higher degree analogue of the well-
known composition formulas of sums of squares of Hurwitz, Radon and Pfister.
We show that such composition identities of sums of higher degree powers are
trivial, namely, if (1.5) holds for $d>2,$ then $r=1.$ In particular, it
follows that there are no nontrivial higher degree $n$-power formulas like
(1.3). As a matter of fact, we confirm a stronger version of the
aforementioned Nathanson’s conjecture. Our proof is simple and elementary, in
which the crux is Harrison’s center theory of higher degree forms [5]. We
remark that nonexistence of (1.5) for $r>1,d>2$ is one of Gude’s main results
and his proof involves some sophisticated methods such as $p$-adic analysis,
Jordan’s theorem, Hurwitz curves, Kodaira dimension, etc. By the same idea, we
also prove that if $r>1$ and $(n,d)\neq(2,3),$ then there is no identity of
type

(1.6)
$(x_{1}^{d}+x_{2}^{d}+\cdots+x_{r}^{d})(y_{1}^{d}+y_{2}^{d}+\cdots+y_{n}^{d})=z_{1}^{d}+z_{2}^{d}+\cdots+z_{n+1}^{d}$

where each $z_{k}$ is a linear form in $y$ with coefficients in
$\mathbbm{k}(x).$ This confirms some special case, namely $r=n=2$ and $d>3,$
of Becker’s conjecture.

The notion of centers of homogeneous polynomials was introduced in the 1970s
by Harrison in his seminal work [5] of the algebraic theory of higher degree
forms. The centers of homogeneous polynomials were rediscovered independently
by several other authors [2, 12, 15, 17] and found various applications in
commutative algebra, algebraic geometry, multilinear algebra and invariant
theory, see e.g. [3, 7, 8, 9]. One major application of centers is to the
problem of direct sum decomposition, or in plain words, separation of
variables of polynomials. There is one-to-one correspondence between the
direct sum decompositions of a homogeneous polynomial of degree $>2$ and
complete sets of orthogonal idempotents of the center algebra of the
polynomial. Moreover, the direct sum decomposition into indecomposable
summands for a higher degree homogeneous polynomial is essentially unique.
This is the crucial ingredient used in our argument.

The paper is organized as follows. In Section 2 we recall Harrison’s theory of
centers and its application to direct sum decompositions of homogeneous
polynomials. Then we apply the center theory to prove our main results in
Section 3.

##  2\. Harrison centers of homogeneous polynomials

Let $f(x_{1},x_{2},\dots,x_{n})\in\mathbbm{k}[x_{1},\dots,x_{n}]$ be a
homogeneous polynomial, also called form, of degree $d.$ Harrison [5] defined
the center of $f=f(x_{1},x_{2},\dots,x_{n})$ as

$Z(f):=\\{X\in\mathbbm{k}^{n\times n}\mid HX=X^{T}H\\},$

where $H=(\frac{\partial^{2}f}{\partial x_{i}\partial x_{j}})_{1\leq i,\ j\leq
n}$ is the hessian matrix.

A homogeneous polynomial $f$ is said to be nondegenerate, if no variable can
be removed by an invertible linear change of variables. It is not hard to
observe that $f$ is nondegenerate if and only if the first-order differentials
$\frac{\partial f}{\partial x_{i}}$ are linearly independent. To avoid
triviality, usually it is enough to consider nondegenerate homogeneous
polynomials. A homogeneous polynomial $f(x_{1},x_{2},\dots,x_{n})$ is called a
direct sum if, after an invertible linear change of variables $y=Px,$ it can
be written as a sum of $t\geq 2$ forms in disjoint sets of variables as
follows

$f=f_{1}(y_{1},\dots,y_{a_{1}})+f_{2}(y_{a_{1}+1},\dots,y_{a_{2}})+\cdots+f_{t}(y_{a_{t-1}+1},\dots,y_{n}).$

If this is not the case, then $f$ is said to be indecomposable. On the other
extreme, if the $f_{i}$’s are forms in only one variable, then $f$ is said to
be diagonalizable.

Corresponding to the homogeneous polynomial $f(x_{1},x_{2},\dots,x_{n}),$
there is an associated symmetric $d$-linear space as follows. First write $f$
in the symmetric way

$f(x_{1},\dots,x_{n})=\sum_{1\leq i_{1},\cdots,i_{d}\leq n}a_{i_{1}\cdots
i_{d}}x_{i_{1}}\cdots x_{i_{d}}$

where the $a_{i_{1}\cdots i_{d}}$’s are symmetric, i.e., they remain unchanged
under any permutation of their subscripts. Then take a vector space $V$ over
$\mathbbm{k}$ of dimension $n$ with a basis $e_{1},\dots,e_{n}$ and define a
symmetric $d$-linear function $\Theta\colon V\times\cdots\times
V\longrightarrow\mathbbm{k}$ by
$\Theta(e_{i_{1}},\dots,e_{i_{d}})=a_{i_{1}\cdots i_{d}}.$ The pair
$(V,\Theta)$ is called the associated symmetric $d$-linear space of $f$ under
the basis $e_{1},\dots,e_{n}.$ One can recover the polynomial $f$ from
$(V,\Theta)$ as

$f(x_{1},\dots,x_{n})=\Theta\left(\sum_{1\leq i\leq n}x_{i}e_{i},\ \dots,\
\sum_{1\leq i\leq n}x_{i}e_{i}\right).$

Clearly, a change of basis for $V$ results in a change of variables for $f.$
The center can also be equivalently defined in terms of symmetric $d$-linear
space $(V,\Theta)$ as

$Z(V,\Theta):=\\{\phi\in\operatorname{End}(V)\mid\Theta(\phi(v_{1}),v_{2},\dots,v_{d})=\Theta(v_{1},\phi(v_{2}),\dots,v_{d}),\
\forall v_{1},v_{2},\dots,v_{d}\in V\\}.$

Nonzero subspaces $V_{1},\dots,V_{t}$ of $(V,\Theta)$ of a symmetric
$d$-linear space $(V,\Theta)$ are said to be orthogonal, if
$\Theta(v_{1},\dots,v_{d})=0$ for all $v_{1},\dots,v_{d}\in
V_{1}\cup\cdots\cup V_{t}$ unless all the $v_{i}$’s are in the same $V_{s}$
for some $1\leq s\leq t.$ If $V=V_{1}\oplus\cdots\oplus V_{t}$ for $t\geq 2$
nonzero orthogonal subspaces, then call $(V,\Theta)$ decomposable. Otherwise,
call $(V,\Theta)$ indecomposable. Clearly, the orthogonal decompositions of
$(V,\Theta)$ are in bijection with the direct sum decompositions of its
associated homogeneous polynomial $f.$ Suppose $V=V_{1}\oplus\cdots\oplus
V_{t}$ is an orthogonal decomposition. Choose a basis
$\epsilon_{1},\dots,\epsilon_{a_{1}}$ for $V_{1},$ a basis
$\epsilon_{a_{1}+1},\dots,\epsilon_{a_{2}}$ for $V_{2},$ and so on and so
forth. Assume $(e_{1},\dots,e_{n})=(\epsilon_{1},\dots,\epsilon_{n})P$ with
$P\in\operatorname{GL}_{n}(\mathbbm{k}).$ Take the change of variables $x=Py,$
then we have the direct sum decomposition for the polynomial
$f=f_{1}(y_{1},\dots,y_{a_{1}})+f_{2}(y_{a_{1}+1},\dots,y_{a_{2}})+\cdots+f_{t}(y_{a_{t-1}+1},\dots,y_{n}).$
In fact, $f_{i}$ is the associated polynomials of $(V_{i},\Theta|_{V_{i}}).$
More precisely,

(2.1) $f_{i}=\Theta\left(\sum_{a_{i-1}+1\leq i\leq a_{i}}y_{i}\epsilon_{i},\
\dots,\ \sum_{a_{i-1}+1\leq i\leq a_{i}}y_{i}\epsilon_{i}\right).$

Conversely, given a direct sum decomposition
$f=f_{1}(y_{1},\dots,y_{a_{1}})+f_{2}(y_{a_{1}+1},\dots,y_{a_{2}})+\cdots+f_{t}(y_{a_{t-1}+1},\dots,y_{n})$
with change of variables $x=Py,$ let
$(\epsilon_{1},\dots,\epsilon_{n})=(e_{1},\dots,e_{n})P.$ Then let $V_{1}$ be
the subspace spanned by $\epsilon_{1},\dots,\epsilon_{a_{1}},$ $V_{2}$ be the
subspace spanned by $\epsilon_{a_{1}+1},\dots,\epsilon_{a_{2}},$ etc. It is
clear that $V=V_{1}\oplus\cdots\oplus V_{t}$ is an orthogonal decomposition.

The following are some useful facts about centers and direct sum
decompositions of homogeneous polynomials obtained in [5], see also [6, 7, 8,
9].

######  Proposition 2.1.

Suppose $f\in\mathbbm{k}[x_{1},x_{2},\dots,x_{n}]$ is nondegenerate and is of
degree $d.$ Then

  * (1)

The center $Z(f)$ is a commutative subalgebra of the full matrix algebra
$\mathbbm{k}^{n\times n}.$

  * (2)

If $f=f_{1}+f_{2}+\cdots+f_{t}$ is a direct sum decomposition, then $Z(f)\cong
Z(f_{1})\times Z(f_{2})\times\cdots\times Z(f_{t}).$

  * (3)

There is a one-to-one correspondence between direct sum decompositions of $f$
and complete sets of orthogonal idempotents of $Z(f).$

  * (4)

If $f=g_{1}+g_{2}+\cdots+g_{s}$ and $f=h_{1}+h_{2}+\cdots+h_{t}$ are direct
sum decompositions, where the $g_{i}$ and $h_{j}$ are indecomposable, then
$s=t$ and after a reordering $g_{i}=h_{i}$ for all $1\leq i\leq s.$

  * (5)

$f$ is diagonalizable if and only if
$Z(f)\cong\mathbbm{k}\times\mathbbm{k}\times\cdots\times\mathbbm{k}\ ($n$\
copies).$

###### Proof.

For the convenience of the reader, we include a sketchy proof for item (4)
which is the key ingredient used to prove our main results. More details can
be found in the aforementioned references. Let $(V,\Theta)$ be the associated
symmetric $d$-linear space of $f.$ Suppose
$f=g_{1}+g_{2}+\cdots+g_{s}=h_{1}+h_{2}+\cdots+h_{t}$ as in the proposition.
Then there are two associated complete sets of primitive orthogonal
idempotents in $Z(V,\Theta)$ corresponding to the given two direct sum
decompositions of $f.$ However, the center algebra $Z(V,\Theta)$ is finite-
dimensional and commutative, so it has only one unique complete set of
primitive orthogonal idempotents [9]. Thus for $V,$ there is a unique
orthogonal decomposition $V=V_{1}\oplus\cdots\oplus V_{r}$ where each $V_{i}$
is indecomposable. As (2.1), each $g_{i}$ or $h_{j}$ is the associated
polynomial of some $(V_{k},\Theta|_{V_{k}}).$ It follows directly that $s=t=r$
and after a reordering $g_{i}=h_{i}$ for all $1\leq i\leq s.$ ∎

##  3\. Nonexistence of some products of sums of powers

In this section, we apply Harrison’s center theory, in particular the
uniqueness of the direct sum decomposition of indecomposable summands, to
prove the nonexistence of some products of sums of powers.

######  Theorem 3.1.

Suppose $d>2.$ There is an identity of product of sums of powers

$(x_{1}^{d}+x_{2}^{d}+\cdots+x_{r}^{d})(y_{1}^{d}+y_{2}^{d}+\cdots+y_{n}^{d})=z_{1}^{d}+z_{2}^{d}+\cdots+z_{n}^{d}$

where each $z_{k}$ is a linear form in $y$ with coefficients in
$\mathbbm{k}(x)$ if and only if $r=1.$

###### Proof.

View both sides of the previous identity as polynomials in $y$ with
coefficients in $\mathbbm{k}(x).$ First note that if the $z_{k}$ are linearly
dependent as linear forms in $y,$ then the polynomial
$z_{1}^{d}+z_{2}^{d}+\cdots+z_{n}^{d}$ is degenerate in $\mathbbm{k}(x)[y]$ as
it can be presented by fewer variables, namely a maximally linearly
independent subset of the $z_{k}.$ Now assume the $z_{k}$ are linearly
independent. Then $z_{1}^{d}+z_{2}^{d}+\cdots+z_{n}^{d}$ is a direct sum, and
obviously each $z_{k}^{d}$ is indecomposable. Thus by item (4) of Proposition
2.1, after a reordering we may assume
$(x_{1}^{d}+x_{2}^{d}+\cdots+x_{r}^{d})y_{j}^{d}=z_{j}^{d}$ for all $1\leq
j\leq n.$ By assumption $z_{k}\in\mathbbm{k}(x)[y],$ hence
$z_{k}=\frac{p(x)}{q(x)}y_{k}$ where $p(x),q(x)\in\mathbbm{k}[x].$ It follows
that $x_{1}^{d}+x_{2}^{d}+\cdots+x_{r}^{d}=\frac{p(x)^{d}}{q(x)^{d}}.$ Since
$\mathbbm{k}[x]$ is a unique factorization domain, this forces $q(x)|p(x).$
Therefore we can assume in the beginning $q(x)=1.$ Now we have
$x_{1}^{d}+x_{2}^{d}+\cdots+x_{r}^{d}=p(x)^{d}.$ This means $p(x)$ is a
$\mathbbm{k}$-linear form in $x.$ Obviously, this happens if and only if
$r=1.$ ∎

######  Remark 3.2.

_It follows easily by the previous theorem that there are no nontrivial
identities like ( 1.3) even if the $z_{k}$ are allowed to be chosen from
$\mathbbm{k}(x)[y].$ Therefore, we confirm a stronger version of Nathason’s
conjecture._

By the same idea, Theorem 3.1 can be slightly generalized to the following

######  Theorem 3.3.

Suppose $r>1,\ n>1$ and $(n,d)\neq(2,3).$ Then there is no identity of the
following type

$(x_{1}^{d}+x_{2}^{d}+\cdots+x_{r}^{d})(y_{1}^{d}+y_{2}^{d}+\cdots+y_{n}^{d})=z_{1}^{d}+z_{2}^{d}+\cdots+z_{n+1}^{d}$

where each $z_{k}$ is a linear form in $y$ with coefficients in
$\mathbbm{k}(x).$

###### Proof.

Assume the contrary and we will seek for a contradiction. Just as in the proof
of the previous theorem, we may assume without loss of generality that
$z_{1},z_{2},\dots,z_{n}$ are linearly independent and $z_{n+1}=\sum_{1\leq
k\leq n}a_{k}z_{k}$ with $a_{k}\in\mathbbm{k}(x).$ If some $a_{k}$ are zero,
say $a_{l}=\cdots=a_{n}=0,$ then it is clear that $z_{n}$ is an indecomposable
summand and so $z_{n}^{d}=(x_{1}^{d}+x_{2}^{d}+\cdots+x_{r}^{d})y_{j}^{d}$ for
some $j.$ For the same reason as in Theorem 3.1, this is impossible. If
$a_{1}a_{2}\cdots a_{n}\neq 0,$ then we compute the center of
$g=z_{1}^{d}+z_{2}^{d}+\cdots+z_{n+1}^{d}$ and come up with a contradiction.
Indeed, let $H$ be the hessian matrix of $g,$ then
$H_{ii}=d(d-1)(z_{i}^{d-2}+a_{i}^{2}z_{n+1}^{d-2})$ for all $1\leq i\leq n$
and $H_{ij}=d(d-1)a_{i}a_{j}z_{n+1}^{d-2}$ for all $i\neq j$. Suppose
$X=(x_{ij})\in Z(g).$ According to the definition of centers, compare the
$ij$-entry and $ji$-entry of $HX$ for all $i\neq j,$ then we have

$z_{i}^{d-2}x_{ij}+z_{n+1}^{d-2}\sum_{1\leq k\leq
n}a_{i}a_{k}x_{kj}=z_{j}^{d-2}x_{ji}+z_{n+1}^{d-2}\sum_{1\leq k\leq
n}a_{j}a_{k}x_{ki}.$

It is not hard to observe that if $(n,d)\neq(2,3),$ then
$z_{i}^{d-2},z_{j}^{d-2},z_{n+1}^{d-2}$ are linearly independent. It follows
easily that $x_{ij}=0$ and $x_{ii}=x_{jj}$ for all $i\neq j,$ that is,
$Z(g)\cong\mathbbm{k}(x).$ This is absurd, since the center of
$(x_{1}^{d}+x_{2}^{d}+\cdots+x_{r}^{d})(y_{1}^{d}+y_{2}^{d}+\cdots+y_{n}^{d})$
is isomorphic to $\mathbbm{k}(x)^{n}$ by item (5) of Proposition 2.1. The
proof is completed. ∎

######  Remark 3.4.

_If $(n,d)=(2,3),$ then $z_{1},z_{2},z_{3}$ are obviously linearly dependent.
Thus the previous argument does not work in this situation. For the moment, we
do not know if one can deal with this special case by the theory of centers._

######  Corollary 3.5.

Suppose $d>3.$ There holds an identity

(3.1)
$(x_{1}^{d}+x_{2}^{d})(y_{1}^{d}+y_{2}^{d})=z_{1}^{d}+z_{2}^{d}+\cdots+z_{n}^{d}$

where each $z_{k}$ is a linear form in $y$ with coefficients in
$\mathbbm{k}(x)$ if and only if $n>3.$

###### Proof.

Consider the special case $r=n=2$ of Theorem 3.3. It follows that the identity
(3.1) holds only if $n>3.$ Conversely, if $n>3,$ let $z_{1}=x_{1}y_{1},\
z_{2}=x_{1}y_{2},\ z_{3}=x_{2}y_{1},\ z_{4}=x_{2}y_{2}$ and the other
$z_{k}=0,$ then obviously 3.1 holds. ∎

######  Remark 3.6.

_This corollary confirms a special case of Becker’s conjecture, namely the
case of $r=s=2$ and $d>3$ with $z_{k}$ assumed to be polynomials in $y.$_

## References

  * [1] E. Becker, The real holomorphy ring and sums of $2^{n}$-th powers, In: Géometrie Algébrique Réelle et Formes Quadratiques, (J.-L. Colliot-Thélène, M. Coste, L. Mahé and M.-F. Roy, eds.), Lecture Notes in Math. 959 (1981), Springer, Berlin, 139–181.
  * [2] J. Carlson, P. Griffiths, Infinitesimal variations of Hodge structure and the global Torelli problem, In: Beauville, A. (Ed.) Journées de Géometrie Algébrique d’Angers, pp. 51–76. Sijthoff and Noordhoff (1980).
  * [3] M. Fedorchuk, Direct sum decomposability of polynomials and factorization of associated forms, Proc. London Math. Soc. (3) 120 (2020) 305–327.
  * [4] U.-G. Gude, Über die Nichtexistenz rationaler Kompositionsformeln bei Formen höheren Grades, doctoral dissertation, Univ. Dortmund, 1988.
  * [5] D. Harrison, A Grothendieck ring of higher degree forms, J. Algebra 35 (1975) 123–138.
  * [6] D. Harrison, B. Pareigis, Witt rings of higher degree forms, Comm. Algebra 16 (6) (1988) 1275–1313.
  * [7] H.-L. Huang, H. Lu, Y. Ye, C. Zhang, Diagonalizable higher degree forms and symmetric tensors, Linear Alg. Appl. 613 (2021) 151–169.
  * [8] H.-L. Huang, H. Lu, Y. Ye, C. Zhang, On centers and direct sum decompositions of higher degree forms, Linear Multilinear Algebra, DOI: 10.1080/03081087.2021.1985057.
  * [9] H.-L. Huang, H. Lu, Y. Ye, C. Zhang, Centers of multilinear forms and applications, arXiv:2110.07915.
  * [10] A. Hurwitz, Über die Komposition der quadratischen Formen von beliebig vielen Variablen, Nahr. Ges. Wiss. Göttingen (1898) 309–316.
  * [11] A. Hurwitz, Über die Komposition der quadratischen Formen, Math. Ann. 88 (1923) 1–25.
  * [12] J. Kleppe, Additive splittings of homogeneous polynomials, Ph.D. thesis, U. Oslo, 2005, arXiv:1307.3532.
  * [13] M. Nathanson, Products of sums of powers, Math. Mag. 48 (1975) 112–113.
  * [14] A. Pfister, Quadratische Formen in beliebigen Körpern, Invent. Math. 1 (1966) 116–132.
  * [15] A. Prószyński, On orthogonal decomposition of homogeneous polynomials, Fund. Math. 98 (1978) 201–217.
  * [16] J. Radon, Lineare scharen orthogonale Matrizen, Abh. Math. Sem. Univ. Hamburg 1 (1922) 1–14.
  * [17] B. Reichstein, On symmetric operators of higher degree and their applications, Linear Alg. Appl. 75 (1986) 155–172.
  * [18] D. Shapiro, Compositions of quadratic forms, Walter de Gruyter & Co., Berlin, 2000.