Fast Normalized Cross-Correlation for Template Matching with Rotations

Normalized cross-correlation is the reference approach to carry out template matching on images. When it is computed in Fourier space, it can handle efficiently template translations but it cannot do so with template rotations. Including rotations requires sampling the whole space of rotations, repeating the computation of the correlation each time. This article develops an alternative mathematical theory to handle efficiently, at the same time, rotations and translations. Our proposal has a reduced computational complexity because it does not require to repeatedly sample the space of rotations. To do so, we integrate the information relative to all rotated versions of the template into a unique symmetric tensor template -which is computed only once per template-. Afterward, we demonstrate that the correlation between the image to be processed with the independent tensor components of the tensorial template contains enough information to recover template instance positions and rotations. Our proposed method has the potential to speed up conventional template matching computations by a factor of several magnitude orders for the case of 3D images.


Introduction
A classical problem in image processing and, particularly, in pattern recognition, is to identify if a large image contains copies -and how many, and their locations and orientations-of a small image, named "template".The resulting algorithms are generically known as template matching algorithms [Brunelli, 2009, Forsyth and Ponce, 2002, Gonzalez and Woods, 2017].The most classical solution is based on using cross-correlations, although there are other approaches based, for example, in metaheuristic algorithms [Corona et al., 2023] or on deep learning [Brunelli, 2009, Lamm et al., 2022, Moebel et al., 2021].In this paper, we show the mathematical foundations of the crosscorrelation-based template matching algorithm (TM in all that follows), and we introduce a new fast algorithm that solves the problem using tensors.
The main advantages of TM, when compared to the algorithms based on machine learning, are that TM is a white box model, it is directly applicable when you have just one template and one larger image (not requiring any kind of training, which may be a very difficult task in some applications), and locates rotations with arbitrary precision (Current deep learning based algorithms for template matching in three-dimensional images are not able to estimate rotations accurately [Lamm et al., 2022]).
On the other hand, a major drawback of TM is its computational cost.TM basic idea is to compute the inner product between the (rotated) template and the (translated) image, and normalize the result.These computations are made, for each rotation, in the Fourier domain to efficiently address translations [Lewis, 1995, Böhm et al., 2000, Roseman, 2003].However, this process has to be repeated for every rotation to be investigated, thus the resulting complexity has a dependency with the rotations processed.The computational cost of this process may become restrictive for 3D images since SO(3), the space of rotations of R 3 , is a (compact) manifold of dimension 3.In application domains such as cryo-electron microscopy, there are required more than ten thousand rotations for achieving an angular precision of a few degrees.
We propose an algorithm called tensorial template matching, TTM, which integrates into a unique symmetric tensor the information relative to the template in all rotations.In other words, the tensor template incorporates in a unique object the information about all rotations of the template, thus allowing us to find the position and rotation of instances of the template in any tomogram with just a few correlations with the linearly independent components of the tensor.The tensor template is computed only once per template, and, as soon as it is generated, it enables to process any image.

Classical template matching
Let us introduce some notation.d-dimensional images are just elements of L 2 (R d ), which is a Hilbert space with the inner product f, g = R d f (x)g(x)dx.It is natural to use the inner product to compare two images f, g of the same size.Concretely, we can use that f, g = f 2 g 2 cos θ, where θ is the angle formed by f and g.In particular, f = αg for some positive constant α if f,g f 2 g 2 = 1.
Template matching is typically used to study if instances of a "small" image t (the template) is present in a larger image f , e.g.look for instances of an specific macromolecule in a cryo-electron tomogram (3D volumetric image).The size of the image is connected to the set of points where the image does not vanish, the support of the image.That is, t is meant "small" when the set K = supp(t) = {x : t(x) = 0} R d is small (e.g., is a subset of a small ball D).Let's assume that f and t have quite different sizes, so our interest is to compare t (the template, the small image) with just a part of f .In such case we need to introduce some special operators S : L 2 (R d ) → L 2 (R d ) that fix our attention in just a part of the domain of f .An interesting example of such operators is where U : R → R denotes Heaviside's unit step function and r > 0. If the support of the template t is D 0 (r) = {x : x 2 ≤ r}, the ball of radius r centered at 0 ∈ R d , the normalized inner product S r (f ), t S r (f ) 2 t 2 informs about the similarity between t and the restriction of f to D 0 (r).Moreover, if we introduce the translation operator the result informs about the similarity between t and the restriction of f to D x (r) = {z : x − z 2 ≤ r}.Of course, it may happen that f contains a copy of a rotated version of t, so rotations are also necessary for a complete discussion of the problem.Thus, given R ∈ SO(d), we define the operator The normalized inner product informs about the similarity of t R and the restriction of The operator S r defined by (1) has some special properties.Concretely, it is symmetric, semidefinite positive and commutes with rotations Recall that, given (X, •, • X ) a (real) Hilbert space1 , an operator S : X → X is named: If S : X → X is a symmetric semidefinite positive operator (SSP, in all what follows), then X becomes a pre-Hilbert space with the inner product and the seminorm Observe, for example, that if S is given by (1), then f S = 0 means that f |D 0 (r) = 0 almost everywhere.
) be an SSP operator, and consider the inner product given by (3).Then Moreover, if f, g ∈ L 2 (R d ), g S = 0, the following are equivalent statements: (c) f − f S g S g S = 0. Proof.As S is SSP, we have that, for all α ∈ R, Hence, if g S = 0, the only way that the quadratic polynomial (in α) above is nonnegative everywhere is that On the other hand, if g S = 0, the only way to satisfy (5) is f, g S = 0, in whose case (6) also holds.This proves (a).
Let us now demonstrate (b) ⇔ (c) whenever g S = 0. Indeed, (c) is equivalent to which holds if and only if ) are two images, α > 0, and we take S = S r given by (1), then τ x f −αO R −1 (t) S = 0 means that f has a match with t R in the unit ball centered at x. Indeed, there are many ways to define operators S with the property that f − g S = 0 means that f = g in a neighbourhood of 0, so that τ x f −αO R −1 (t) S = 0 means that f has a match with a rotated version of t in a neighbourhood of 0. Although arbitrary SSP operators may not enjoy this property, they allow the creation of a general way to deal with this kind of operators.
Thus, in all what follows, we assume that S : Rotations and composition of operators will play an important role in this paper.Thus, it is natural to ask how the composition of rotations acts on the images.This is, indeed, a simple computation: 2 1 is not an element of L 2 (R d ), but this can be managed in several ways.In fact, although we use L 2 (R d ) to denote the space of d-dimensional images, in practice we only consider images f with compact support K.Then, when we compute f, 1 , we mean f, 1 = R d f (x)dx = K f (x)dx = f, 1χ K .Moreover, since our interest is on operators S that vanish on functions vanishing outside of a certain neighbourhood D of 0, by 1, 1 S we mean 1χ D , 1χ D S .

Hence
Given an image f , we consider its projection onto the space of images which are S− orthogonal to the constant image 1, Remark 2.2 These projections are important to study invariant properties with respect to constant brightness changes in the images of translations and rotations.Note that there is no "real" difference between an image f and the images of the form f + α1, α ∈ R. When we modify the constant α, what we observe is a uniform change in the density or the brightness, but not the apparition of new structures or forms, in the image f .Thus, f and its projection P S (f ) essentially represent the very same image, since f = P S (f ) + α1 for certain α ∈ R.
Given two images f, t, we have that f = P S (f ) + α1 and t = P S (t) + β1 for certain constants α, β.
Hence f, t S = P S (f ) + α1, P S (t) + β1 S = P S (f ), P S (t) S + αβ 1, 1 S since P S (f ), P S (t) ⊥ S 1. Consequently, if x ∈ R d and R ∈ SO(d), there are two constants ρ = ρ(x) and δ = δ(R) such that Assume that S commutes with rotations, and take x ∈ R d fixed.Then, for each R ∈ SO(d) we have that Hence O R −1 (P S (t)) ⊥ S 1 and this, in conjunction with (10), implies that is an S-orthogonal decomposition of t R , which means that the constant β that multiplies 1 in the S-orthogonal decomposition of t R does not depend on R, and In particular, for each x ∈ R d , the problems: • Maximize τ x (f ), t R S over rotations R.
• Maximize P S (τ x (f )), P S (t) R S over rotations R.
• Maximize P S (τ x (f )), P S (t R ) S over rotations R. are equivalent.
Let us define: Lemma 2.3 If S is an SSP operator that commutes with rotations, the parameter δ that appears in the S-orthogonal decomposition does not depend on R. Consequently, given x ∈ R d , the problems • Maximize P S (f −x,R −1 ), P S (t) S over rotations R. are equivalent.
Proof.We know that δ , so that we only need to prove that f −x,R −1 , 1 S does not depend on R. Indeed, We can now state and demonstrate the following: Theorem 2.4 (Classical template matching) Let S be a SSP operator which commutes with rotations and let x ∈ R d be fixed.Then the following are equivalent problems: Moreover, if t S > 0 and S also has the property that f S = 0 implies f |D = 0 for a certain neighborhood D of 0 ∈ R d which contains the supports of all the rotated templates t Q with Q ∈ SO(d), then a match between f and t R in x is got whenever any one of the following claims hold: The other claims are a direct consequence of Theorem 2.1.✷ In all that follows, we assume that S is an SSP operator that commutes with rotations and t is normalized in the sense that t ⊥ S 1 and t S = 1.Then P S (t R ) = t R and t R S = 1 for all rotation R. Consequently, attains its maximum (= 1) if and only if there is a perfect match between f and t R in x.Moreover, if we define w(x) A perfect match is, in general terms, never attained.This is so because the desired image, represented by the template t, is usually supported on a strict subset Ω of the domain D were the operator S is able to distinguish functions.Thus, the image f may well contain a copy of the image represented by t R but in the neighbourhoods of the support of t R , f will contain some information which is not present in t R .In addition, f is usually corrupted by noise and distortions.This means that the normalized correlations described in items (a * ) − (d * ) of Theorem 2.4, will never equal 1.Consequently, a threshold should be introduced in order to decide if a match has (or has not) been produced.
In order to find the rotation which maximizes c(x, R), the cross-correlation (f ⋆ S(t) R )(x) should be computed for a huge amount of rotations R, which makes classical matching an inefficient approach for template matching.Indeed, for d = 3, the size of the set of rotations R used to sample SO(3) well enough to guarantee a reliable result varies between 10 4 and 5 • 10 5 rotations [Chaillet et al., 2023].
Due to numerical reasons, high frequencies may be altered during rotation transformation.Thus, in practice, we do not apply the operator S to the original images f, t but to a filtered version of them that eliminates these high frequencies.Concretely, we apply an isotropic (i.e.rotation invariant) lowpass filter h to both images and, after that, we apply the template matching algorithm to the resulting images.The idea behind this is that, if there is a match between f and t, there will be a match between f = f * h and t = t * h too.The operator S results from applying a rotationally symmetric mask m(x) = ρ( x ) to the given image.Thus, we substitute f by f = f * h and t by t = t * h.Then we apply the classical (or tensor) matching algorithm to the pair of images f, t using the SSP operator S(f)(x) = m(x)f(x).Usually, the mask m equals 1 within a certain radius around 0 and equals 0 outside a sightly larger radius.In between these radii the mask takes values between 0 and 1.Under these restrictions, it is clear that the operator S is SSP and commutes with rotations.Moreover, if 0 = f S = f, S(f ) ≥ D f 2 (x)dx ≥ 0, we have that f |D = 0 where D is a ball of positive radius centered at 0. Let us compute the inner product (since every filter is translation invariant, and and we use • to denote the standard product of real functions.This means that we would have the same effect just considering the template matching algorithm associated with the operator S applied to the images f, t.Moreover, the following holds: with h defining an isotropic filter and m a rotationally symmetric mask as described above.Then S is SSP. Proof.For the proof, we use the following (well-known) formulae: For functions a, b, c Let us now consider the product f, S(f ) : (since h = h and m ≥ 0).This proves that S is semidefinite positive.Let us show the symmetry: ✷ In all that follows, we assume that the SSP operator S is of the form ( 14) with h, m verifying the hypotheses of Lemma 2.5.Thus, the template matching algorithm is applied with this operator and a fast computation of c(x, R) is needed.
A direct computation leads to: (since h is isotropic, and m is rotationally symmetric).Moreover, Now, using the definition of S (and imposing h * 1 = 1), we can simplify the computation as follows:

S
Note that the FFT algorithm can be used to compute the inner products appearing at the end of the formula above, which helps to fasten the algorithm.
Moreover, the following identities also hold: Thus, Algorithm 1, devoted to classical template matching, has been developed using the formulae above.
An important tool we will use in this paper is the set H of quaternions.In particular, we will use that rotations can be parametrized by unit quaternions (which can be identified with the unit 3-sphere S 3 ), as well as the following formulae (see, e.g.[Ebbinghaus et al., 1991, Pontryagin, 2010]): • Given x ∈ H, x = a + bi + cj + dk, we identify x with a pair (a, v) where a ∈ R and v = (b, c, d) ∈ R 3 , and call a the real part of x, a = Re(x).
Then, if x = (a, v), y = (b, w) ∈ H, we have that We end this section with a well-known result about composition of SSP operators that will be used in the proof of the main theorem of the paper: ) are semidefinite positive symmetric and commute, then T S is also semidefinite positive symmetric.
Proof.If S, T are SSP, then they have (unique) positive square roots √ T , √ S, which are symmetric operators too.Moreover, the composition √ T √ S is also symmetric.Since T, S commute, their positive square roots commute.Hence and the square of any symmetric operator is positive.✷

Tensor template matching
This section introduces a tensorial template matching (TTM) algorithm.
The purpose is to handle translations and rotations efficiently at the same time.First, we introduce some background necessary to understand further mathematical developments.Second, we present the main theorem for TTM, which allows us to determine the optimal rotation of the template, t, on every match in the image f without sampling the SO(3) by computing some tensors.Finally, we explain how to determine match positions (template translations) directly from the computed tensors.
Algorithm 1 Classical template matching with rotations Load and Fourier transform of the image Load mask and pre-process the image Load the mask into m Load and pre-process the template Load template into t

Tensor background
A tensor A ∈ T n (R d ) of order n and dimension d is just an array of the form for every permutation σ ∈ Σ n (the set of permutations of {1, • • • , n}).We denote by S n (R d ) the set of symmetric tensors of order n and dimension d.An important example of symmetric tensor of order n is the so called n-th tensor power of a vector It is well known that T n (R d ) and S n (R d ) are real vector spaces with the natural operations (pointwise sum and multiplication by a scalar), and that dim Moreover, every symmetric tensor is a finite sum of tensor powers, which allows us to introduce the concept of the (symmetric) rank of a symmetric tensor as the minimal number of tensor powers used to represent the tensor with their sum [Comon et al., 2008].
The map •, • : defines an inner product.It is also usual to denote A• B = A, B .Moreover, with this notation, if x, y ∈ R d are d-dimensional vectors, a direct application of the multinomial theorem shows that

we can also consider the inner product
which can be seen as an homogeneous polynomial in d variables, of degree n, which justifies using the notation ) denotes the symmetric tensor whose components are In particular, where ∇ϕ denotes the gradient of the function ϕ : R d → R.
Note that the vector x can be chosen from C d in the definitions above.This justifies the following definition (see [Cui et al., 2014] Using gradients, we can rewrite the equation Au n−1 = λBu m−1 as Hence u is a B-eigenvector of A if and only if it is a critical point of the following optimization problem: Maximize: Ax n under the restriction: Bx m = 1. (23) Two particularly important cases are the H-eigenvectors

Maximize:
Ax n under the restriction: and the Z-eigenvectors

Maximize:
Ax n under the restriction: The optimization problem associated with finding Z-eigenvectors of a given symmetric tensor is particularly important for us since the tensor matching algorithm we propose is reduced to one of these problems in each position, and, fortunately, there are good iterative algorithms to approximate the solutions of (25) (see e.g.[Kofidis andRegalia, 2001, Kolda andMayo, 2010]).These algorithms have a linear rate of convergence.In section 3.4 we show an heuristics that can be used to select the positions where a match is probable, so that solving (25) is necessary.

Defining of the Tensor template
In all that follows in this paper, S : L 2 (R d ) → L 2 (R d ) denotes an SSP operator which commutes with rotations and the template t ∈ L 2 (R d ) is assumed to be normalized by t ⊥ S 1 and t S = 1.In section 2 we proved that, for each x ∈ R d , c(x, R) = w(x)(f ⋆ S(t) R )(x) attains its maximum value on rotation R (and this value equals 1) if and only if there is a match between f and t at (x, R) (i.e., a match between τ x (f ) and t R ).Let us define the symmetric tensor C n (x) ∈ S n (R d ′ ), where d ′ is the number of parameters used to describe the rotations SO(d) (in particular, for d = 3, we get d ′ = 4) by the formula: This means that where is a tensor template (or tensorial needle).
It is of fundamental importance to observe that T (z) is computed only once and contains a reduced number of components, since dim (in particular, dim S 4 (R 4 ) = 7 4 = 35).Indeed, this is one of the main reasons why the tensor template matching algorithm we introduce in this paper is fast.The other reason is that rotations R defining a match between f and t R at x are Z-eigenvectors of the symmetric tensor C n (x), which is really remarkable because the power method used in [Kolda andMayo, 2010, Kofidis andRegalia, 2001] for the solution of the corresponding optimization problem does not require using myriads or even millions of rotations -as is the case with classical matching algorithms-but just a reduced set of them: one by iteration.

Finding the correct rotation
Let us state the main result of this paper: and n ∈ 2N be given.If there is a match between f and t R at x, the function ϕ(Q) = C n (x) • Q ⊙n , defined on rotations of R 3 , when parametrized by unit quaternions Q, attains its global maximum at Q = R.
Proof.We prove the result as a consequence of Theorem 2.4.Thus, our main goal is to represent ϕ(Q) in terms of a scalar product, ϕ(Q) = w(x) τ x (f ), t Q S ′ , for some SSP operator S ′ , which would guarantee that if there is a match between f and t R at x, then ϕ(Q) attains its global maximum (and is equal to 1) at Q = R.
Let us compute ϕ(Q): (by definition of c(x, R)).Hence, dividing by w(x), we get: In other words, Here, must be interpreted as a function defined on R 3 with values on R (indeed, it is an element of L 2 (R 3 )): where S(t)(z) : SO(3) → R is given by Indeed, in general, every element t ∈ L 2 (R 3 ) can be interpreted, for each z ∈ R 3 , as an element of L 2 (SO(3)) just making t(z It follows that S(t) QP (z)K(P −1 )dP (in the last equality, set R = QP and use that |Q| = 1) S(t) P K(P −1 )dP where I d denotes the identity rotation.
Let us now denote by S 2 : L 2 (R 3 ) → L 2 (R 3 ) the operator given by and let S ′ = S 2 • S. Then Thus, the proof ends as soon as we demonstrate that S ′ is an SSP operator, and it is for this that we need to use Lemma 2.6.Indeed, S ′ = S 2 • S is a composition of operators, S is, by hypothesis, symmetric semidefinite positive, and S, S 2 commute because S commutes with rotations and S 2 is defined in terms of convolution in SO(3).Thus, Lemma 2.6 implies that S ′ is SSP whenever S 2 is SSP.
To prove that S 2 is symmetric semidefinite positive, we use the properties of the convolution on SO(3) when interpreted as a hyperspherical convolution on S 3 , the unit sphere of R 4 .Recall that if , which makes of S d−1 a homogeneous space and allows to introduce the convolution of functions defined on S d−1 as follows: where η ∈ S d−1 is the north pole of the sphere and f, g ∈ L 2 (S d−1 ).Thus, if we use that the elements of SO(3) are parametrized by quaternions of norm 1, which that can be identified with the elements of the sphere S 3 = {x ∈ H : xx = |x| = 1} in fourth-dimensional space, then assuming that the north pole of S 3 is given precisely by the identity rotation I d , the convolution of f, g ∈ L 2 (SO(3)) can be interpreted as a hyperspherical convolution on S 3 : Now, as it is well known, L 2 (S 3 ) is a Hilbert space and the so-called hyperspherical harmonics, {Ξ ℓ M }, form an orthonormal basis of this space.Thus, every function f ∈ L 2 (S 3 ) admits a Fourier expansion Moreover, in [Dokmanic and Petrinovic, 2009], it was proven that, if f, g ∈ L 2 (S 3 ) and It follows that, given a template t, for each z ∈ R 3 , the map t(z)(R) = t R (z) belongs to L 2 (S 3 ) (here the rotations R are parametrized as unit quaternions, so that R ∈ S 3 ) and We need the following Lemma, whose proof is included in section 4: This means that convolution with K, which is an operator , is semidefinite positive.Moreover, it is well known that this operator is symmetric (and we will use both things in our computations bellow).
In order to prove that S 2 is SSP, we introduce the operator L : L 2 (R 3 ) → C(SO(3), L 2 (R 3 )) defined by L(t)(R) = t R , as well as the operator L * : where ) and a(w)(R) := a(R)(w).On the other hand, Thus, S 2 (t) is semidefinite positive.Moreover, the same type of computation shows that which proves that S 2 is symmetric.This ends the proof of the theorem.✷ Note that Theorem 3.1 connects the problem of finding, at a given position x, the rotation R which gives a match between f and t R at x with the problem of finding the dominant Z-eigenvalue-eigenvector pair by solving (25) with A = C n (x) ∈ S n (R 4 ) and n even.

Finding the correct position
Although we can find the spatial positions of peaks by running an algorithm to find the dominant Z-eigenvalue-eigenvector pair for each and every voxel, this is fairly expensive using the current decomposition algorithms for higher degree tensors.However, the Frobenius norm of a tensor is related to its spectral norm, and in practice it turns out it can be used as an excellent proxy for finding the spatial locations of peaks.Indeed, we know that C n (x) ∈ S n (R d ′ ) = S n (R 4 ).Now, if T σ denotes the spectral norm of tensor T and T F denotes its Frobenius norm, it is well-known that the largest singular value of T equals its spectral norm, and that (see e.g., [Cao et al., 2023, Kozhasov andTonelli-Cueto, 2022]).
In fact, the connection between T σ and T F is stronger than just this inequality.As is well-known, every tensor is a finite sum of tensors of rank 1 (indeed, if the tensor is symmetric, the tensors of rank one can also be chosen symmetric) [Comon et al., 2008].Moreover, if W 1 is a tensor of rank 1 satisfying T − W F then (see, e.g., [Regalia and Kofidis, 2000]) (T ) is preserved, an increment on the size of T σ ( T F , respectively) is translated into an increment on the size of T F ( T σ , respectively).
Moreover, in 1938 in [Banach, 1938] was demonstrated that, for any symmetric tensor T , Thus, large T F implies large spectral norm of T , and the spectral norm of C n (x) is strongly connected to the optimization problem solved in Theorem 3.1, which justifies using the Frobenius norm of C n (x) as a parameter to select positions x where a match is possible.
For each position x identified as a potential peak, the SS-HOPM algorithm (see [Kofidis and Regalia, 2001, Kolda and Mayo, 2010, Regalia and Kofidis, 2000] for precise definition and implementation of this algorithm) is used to find the exact dominant Z-eigenvalue and its associated Z-eigenvector, which is the rotation R candidate to give a match at x.
We have just explained an heuristics to locate the positions -and, after that, the rotations-where a match is possible.Now, sometimes a false positive may occur.Indeed, in the previous subsection we showed that the tensor-based correlation function C n (x), Q ⊙n can be seen as using a slightly different degenerate inner product, based on S ′ from the proof of Theorem 3.1, rather than S. Concretely, we proved that where w(x) = 1 P S (τx(f )) S , t R S = 1 and t R = P S (t R ).This implies that the relation −1 ≤ C n (x), Q ⊙n ≤ 1 does not necessarily hold because the normalizations were taken in terms of S instead of S ′ .Taking S ′ into account would lead to the equality ≤ 1 (and it is equal to 1 when we have a match).So what is the impact of this?First of all, observe that the operation that is missing from S in the normalization is effectively a kind of convolution, so that its effect on the constant component of an image is to scale it.Consequently, if an image is S-orthogonal to 1, it will also be S ′ -orthogonal to 1.However, the norms are affected.
For the template t, this means that the normalization is off by a certain factor, but this factor is the same everywhere.For the image f , the impact is less benign though, as τ x (f ) S will differ from τ x (f ) S ′ in a nonuniform way.
When will this shortcoming would cause a false positive?For this to happen, the normalization factor used at a non-match position would have to be much higher than the "correct" normalization factor, and/or the normalization factor would have to be too low at a match position.Since the difference between S and S ′ is essentially a smoothing operation, and the normalization factor is the reciprocal of the norm of the projected image, the image would thus have to be (very) smooth at the non-match position, while exhibiting a lot of high frequency energy around the match position.Such a situation would not be impossible, but would at the very least be unusual in the context of a typical application like the analysis of electron microscopy images.

Proof of Lemma 3.2
Let us start recalling the formulae associated to Fourier expansions in hyperspherical harmonics on the sphere S 3 .The parametrization of the sphere we consider is the following one: where (a, b, c, d) ∈ S 3 is identified with the unit quaternion Q = a+bi+cj+dk, which represents a rotation of three dimensional euclidean space R 3 .The volume element (used for integration on S 3 and, henceforth, also in SO( 3)) is then given by dV = sin 2 θ sin φdθdφdϕ.
Then every function f ∈ L 2 (S 3 ) can be decomposed as where {Ξ ℓ (k 1 ,k 2 ) } denotes the orthonormal basis of L 2 (S 3 ) formed by the hyperspherical harmonics and f (ℓ, (k 1 , k 2 )) = f, Ξ ℓ (k 1 ,k 2 ) S 3 are the Fourier coefficients of f in this basis.We want to prove that K(ℓ, (0, 0)) ≥ 0 for all ℓ.Now, where A ℓ (0,0) is a positive constant and C λ ℓ (t) denotes the Gegenbauer polynomial of degree ℓ, which appears as the ℓ-th Taylor coefficient in the expansion: The parity of ℓ and n implies that all factors that appear multiplying the variable θ inside of the cosine functions are even numbers.This makes the corresponding integrals (on [0, π]) equal to 0, except in the case that the factor itself is 0. In such case, cos(0) = 1 implies that only the cosine functions that appear with a minus sign in front of them in the formula can contribute with a negative number to the integral.Now clearly ℓ + 2 > 0 always since ℓ ≥ 0, and ℓ + 2 + n − 2k = 0 implies 2k = n + ℓ + 2 > n, so that k > n/2 which is impossible since the sum's range goes from k = 1 to k = n/2.This means that the term − cos(θ(ℓ + 2 + n−2k)) never contributes with a negative number to the sum.On the other hand, if the cosine function with factor ℓ + 2 + 2k − n contributes, which means that ℓ + 2 + 2k − n = 0, then k = (n − ℓ − 2)/2 < n/2.In particular, taking k * = k + 1, we have that 1 ≤ k * ≤ n/2 so that cos(θ(ℓ + 2k * −n)) = cos(θ(ℓ + 2k + 2 −n)) = cos(0) = 1 and the corresponding term effectively appears in the sum.In particular, adding these two terms of the sum we get since n even and k < n/2.This ends the proof of Lemma 3.2 ✷

Conclusions
We have exposed the maths of classical template matching with rotations.Moreover, an alternative to the classical algorithm, named tensorial template matching (or TTM), has been shown.TTM integrates the information relative to all rotated versions of a template t into a unique symmetric tensor template T , which is computed only once per template.The main theorem of the paper, Theorem 3.1, shows that finding an exact match between an image f and a rotated version t R of the template t at a given position x is equivalent to finding a best rank 1 approximation (in the Frobenius norm) to a certain tensor C n (x).The resulting algorithm has reduced computational complexity when compared to the classical one.TTM finds the position and rotation of instances of the template in any tomogram with just a few correlations with the linearly independent components of T .In particular, Cryo-electron tomography (3D images) for macromolecular detection requires 7112, 45123 and 553680 rotations to achieve an accuracy of 13 • , 7 • and 3 • respectively [Chaillet et al., 2023].Therefore, and considering 4-degree tensors (35 linearly independent components), the potential speed-up of our approach with respect to TM is 203x, 1239x and 184560x in these cases, while the angular accuracy remains constant for TTM and it is limited by the computation of tensorial template.