Reconstructing Group Wavelet Transform From Feature Maps With a Reproducing Kernel Iteration

In this article, we consider the problem of reconstructing an image that is downsampled in the space of its SE(2) wavelet transform, which is motivated by classical models of simple cell receptive fields and feature preference maps in the primary visual cortex. We prove that, whenever the problem is solvable, the reconstruction can be obtained by an elementary project and replace iterative scheme based on the reproducing kernel arising from the group structure, and show numerical results on real images.


INTRODUCTION
This article introduces an iterative scheme for solving a problem of image reconstruction, motivated by the classical behavior of the primary visual cortex (V1), in the setting of group wavelet analysis. The mathematical formulation of the problem is that of the reconstruction of a function on the plane which, once represented as a function on the group SE(2) = R 2 ⋊ S 1 of rotations and translations of the Euclidean plane via the group wavelet transform, is known only on a certain twodimensional subset of this three-dimensional group. This problem is equivalent to that of filling in missing information related to a large subset of the SE(2) group and ultimately inquires about the completeness of the image representation provided by feature maps observed in V1.
One of the main motivations for the present study comes indeed from neuroscience and the modeling of classical receptive fields of simple cells in terms of group actions restricted to feature maps such as the orientation preference maps. The attempts to model mathematically the measured behavior of the brain's primary visual cortex (V1) have led to the introduction of the linearnonlinear-Poisson (LNP) model (Carandini et al., 2005), which defines what is sometimes referred to as classical behavior. It describes the activity of a neuron in response to a visual stimulus as a Poisson spiking process driven by a linear operation on the visual stimulus, modeled by the receptive field of the neuron, passed through a sigmoidal nonlinearity. A series of thorough studies of single cell behavior could find a rather accurate description of the receptive fields of a large amount of V1 cells, called simple cells, within the LNP model, in terms of integrals over Gabor functions located in a given position of the visual field (Marcelja, 1980;Ringach, 2002). This description formally reduces the variability in the classical behavior of these cells to few parameters, regulating the position on the visual field, the size, the shape, and the local orientation of a twodimensional modulated Gaussian. A slightly simplified description of a receptive field activity F in response to a visual stimulus defined by a real function f on the plane, representing light intensity, is the following: denoting by x = x 1 x 2 , y = y 1 y 2 ∈ R 2 , f (y)e −2π ip (x 1 −y 1 ) cos θ +(x 2 −y 2 ) sin θ e −2π 2 s 2 |x−y| 2 dy (1) where the parameters s, p ∈ R + define the local (inverse) scale and spatial frequency, the angle θ ∈ [0, 2π) defines the local direction and x ∈ R 2 define the local position of the receptive field in the visual field, while the complex formulation can be formally justified by the prevalence of so-called even and odd cells (Ringach, 2002). We will focus on the sole parameters x, θ . This may be considered restrictive (Hubel and Wiesel, 1962;Sarti et al., 2008;Barbieri, 2015), but it is nevertheless interesting since angles represent a relevant local feature detected by V1 (Hubel and Wiesel, 1962) whose identification has given rise to effective geometrical models of perception (Petitot and Tondut, 1999;Paul C. Bressloff, 2002;John Zweck, 2004;Sarti, 2006, 2015). We will give more details below concerning this restriction. In this case, one can then model the linear part of the classical behavior of simple cells in terms of an object that is well-studied in harmonic analysis: by rephrasing Equation (1) as for fixed values of p, s, we obtain a SE(2) group wavelet transform of f (refer to Section 2). On the other hand, classical experiments with optical imaging have shown that not all parameters θ ∈ S 1 are available for these cortical operations (Blasdel and Salama, 1986;Bonhoeffer and Grinvald, 1991;Weliky et al., 1996;Bosking et al., 1997;White and Fitzpatrick, 2007). Furthermore, with single-cell precision two-photon imaging techniques (Ohki et al., 2006), we can say with a good approximation that in many mammals, a pinwheel-shaped function : R 2 → S 1 determines the available angle at each position. This feature map can be introduced in the Model (1) by saying that the receptive fields in V1 record the data.
Within this setting, a natural question is, thus, whether this data actually contains all the sufficient information to reconstruct the original image f , and how can we reobtain f from these data. This is the main problem we aim to address. Note that one can equivalently consider the graph G = {(x, θ ) ∈ R 2 × S 1 : θ = (x)} of the feature map instead of the feature map itself: this would allow us to address the problem as that of the injectivity of the restriction of the function (Equation 2) to the subset G , which is the point of view taken in Section 3.
Before proceeding, we recall that a severe limitation of the purely spatial Model (Equation 1) is that of disregarding temporal behaviors (DeAngelis et al., 1995;Cocci et al., 2012). Moreover, the classical behavior described by the LNP model is well-known to be effective only to a limited extent: several other mechanisms are present that describe a substantial amount of the neural activity in V1, such as Carandini-Heeger normalization (Carandini and Heeger, 2012), so-called nonclassical behaviors (Fitzpatrick, 2000;Carandini et al., 2005), and neural connectivity (Angelucci et al., 2002). However, spatial receptive fields and the LNP provide relevant insights on the functioning of V1 (Ringach, 2002). Moreover, the ideas behind the LNP model have been a main source of inspiration in other disciplines, notably for the design of relevant mathematical and computational theories, such as wavelets and convolutional neural networks (Marr, 1980;LeCun et al., 2010). We also point out that the use of groups and invariances to describe the variability of the neural activity has proved to be a solid idea to build effective models (Citti and Sarti, 2006;Anselmi and Poggio, 2014;Petitot, 2017), whose influence on the design of artificial learning architectures is still very strong (Anselmi et al., 2019(Anselmi et al., , 2020Montobbio et al., 2019;Lafarge et al., 2021).
We would like to mention that, in addition, some other simplifications already at the level of modeling the receptive fields of V1 cells are assumed when considering (Equation 3) as the data collected by V1, i.e., by restricting the parameters of interest to only (x, θ ) and by considering just one angle θ for each position x. First, in Equation (1), one performs operations x − y, meaning that both x and y are coordinates on the visual plane. On the other hand, it is known that a nontrivial retinotopic map links the cortical coordinates with the visual field (Tootell et al., 1982). Thus, in order for a cortical map to be compared with actual measurements such as those of Blasdel and Salama (1986), Bonhoeffer and Grinvald (1991), Weliky et al. (1996), Bosking et al. (1997), Ohki et al. (2006), White and Fitzpatrick (2007), in Equation (3), we should not consider as computed directly on the visual field coordinate x but rather compose with the retinotopic map. Equation (3) is, thus, to be considered, from the strictly neural point of view, either as a local approximation or as a formulation where also contains the retinotopic map. Moreover, we are assuming yet another simplification of what would be a realistic model by considering the parameters p, s in Equation (1) as fixed. More specifically, it is known that the spatial frequency p of the receptive field is itself a feature that is organized in a columnar fashion, i.e., by a cortical map (Ribot et al., 2013), which is correlated with that of orientations, and recent refined mathematical models are able to reproduce these cortical maps on the basis of the sole structure of the receptive field (Baspinar et al., 2018). Concerning the (inverse) scale parameter s, it is known (Hubel and Wiesel, 1974) that it does not only vary horizontally on the cortex, but it also changes along the transverse direction, hence providing a local analysis at more than one scale. Moreover, considering the variation of the transversal average local scale with respect to the horizontal cortical coordinates, it has been observed that, up to the scatter variations, this scale is also correlated with another cortical map, the retinotopic map (Harvey and Dumoulin, 2011). Finally, the two parameters p and s appear to be statistically correlated in the populations of simple cells whose receptive fields have been measured and fitted with the Gabor model (Ringach, 2002). If these assumptions are considered satisfactory, however, it becomes easy to reconcile a description such as that of Equation (3) with other classical ones based on the concept of tuning curves (refer to e.g., Webster and Valois, 1985. An orientation tuning curve for a given V1 cell is a measured bellshaped curve that shows the response of the cell to different orientations. That is, it describes how a cell that has a given orientation preference θ 0 actually has a nonzero activation when a visual stimulus with similar but not equally oriented stimuli are presented in its receptive field area. In this sense, each cell is not sensitive to just one orientation, in the same way as it is not sensitive to just one point in the visual space. These areas of influence in the parameter space (x, θ ) for a given cell of V1 can actually be related to the structure of receptive fields and the uncertainty principle of the SE(2) group (Barbieri et al., 2014), and orientation tuning curves too can be obtained from Gabor-like receptive fields (Barbieri, 2015).
For the sake of completeness and previously unfamiliar readers, we also recall that simple cells receptive fields are sometimes modeled as scaled directional derivatives of Gaussians. This model, which can actually reproduce quite closely certain Gabor functions, was initially introduced for its theoretical interest and its neural plausibility (Webster and Valois, 1985;Young, 1985;Koenderink and van Doorn, 1987) and later used extensively in the scale space theory in computer vision (Florack et al., 1992;Lindeberg, 1994) as well as for a largely developed axiomatic theory of receptive fields (Lindeberg, 2013(Lindeberg, , 2021. A thorough discussion on the relationship between this model and that of Gabor functions can be found in (Petitot, 2017, Ch. 3).
Another motivation for studying the completeness of the data collected with Equation (3) comes from the relationship of this problem with that of sampling in reproducing kernel Hilbert spaces (RKHS). The RKHS structure, in this case, is that of the range of the group wavelet transform, and will be discussed in detail, together with the basics on the SE(2) wavelet transform, in Section 2. Reducing the number of parameters in a wavelet transform is a common operation, that one typically performs for discretization purposes (refer to e.g., Daubechies, 1992), or e.g., because it is useful to achieve square-integrability (Ali et al., 1991). The main issue is that these operations are typically constrained nontrivially in order to retain sufficient information that allows one to distinguish all interesting signals. For example, when discretizing the short-time Fourier transform to obtain a discrete Gabor Transform, the well-known density theorem (Heil, 2007) imposes a minimum density of points in phase space where the values of the continuous transform must be known in order to retain injectivity for signals of finite energy. Much is known about this class of problems in the context of sampling, from classical Shannon's theorem and uncertaintyrelated signal recovery results (Donoho and Stark, 1989) to general sampling results in RKHS (Fuhr et al., 2017;Grochenig et al., 2018). However, the kind of restriction implemented by feature maps does not seem to fit into these settings, even if some similarities may be found between the pinwheel-shaped orientation preference maps (Bonhoeffer and Grinvald, 1991) and Meyer's quasicrystals, which have been recently used for extending compressed sensing results (Candes et al., 2006;Matei and Meyer, 2010;Agora et al., 2019).
We remark that this article does not focus on dictionary learning techniques for image representation that may give rise to realistic receptive fields and feature maps, such as e.g., those studied in Hyvarinen and Hoyer (2001) and Miikkulainen et al. (2005). We will consider only the problem of reconstruction of the SE(2) wavelet transform from a given feature map.
The plan of the article is as follows. In Section 2, after a brief overview of the SE(2) transform and its RKHS structure, we will formalize in a precise way the problem of functional reconstruction after the restriction (Equation 3) has been performed. In Section 3, we will introduce a technique to tackle this problem, given by an iterative kernel method based on projecting the restricted SE(2) wavelet transform onto the RKHS defined by the group representation. Moreover, we will consider a discretization of the problem and, in the setting of finite dimensional vector spaces, we will give proof that the proposed iteration converges to the desired solution. Finally, in Section 4, we will show and discuss numerical results on real images.

OVERVIEW OF THE SE(2) TRANSFORM
The purpose of this section is to review the fundamental notions of harmonic analysis needed to provide a formal statement of the problem. We will focus on the group wavelet transform defined by the action on L 2 (R 2 ) of the group of rotations and translation of the Euclidean plane, expressed as a convolutional integral transform.
We will denote the Fourier transform of f ∈ L 1 (R 2 ) ∩ L 2 (R 2 ) by and, as customary, we will use the same notation for its extension by density to the whole L 2 (R 2 ). We will also denote by * the convolution on R 2 : Let S 1 be the abelian group of angles of the unit circle, which is isomorphic either to the one dimensional torus T = [0, 2π) or to the group SO(2) of rotations of the plane R 2 . The group SE(2) = R 2 ⋊ S 1 is (refer to e.g., Sugiura, 1990, Ch. IV) the semidirect product group with elements (x, θ ) ∈ R 2 × S 1 and composition law where r θ = cos θ − sin θ sin θ cos θ . Its Haar measure, which is the (Radon) measure on the group which is invariant under group operations, is the Lebesgue measure on R 2 ⋊ S 1 . A standard way to perform a wavelet analysis with respect to the SE(2) group on two-dimensional signals is given by the operator defined as follows. DEFINITION 1. Let us denote by R : S 1 → U(L 2 (R 2 )) the unitary action by rotations of S 1 on L 2 (R 2 ): Let ψ ∈ L 2 (R 2 ), and denote by ψ θ = R(θ )ψ. The SE(2) wavelet transform on L 2 (R 2 ) with the mother wavelet ψ is In terms of this definition, we can then see that if we choose s, p ∈ R + , and let ψ s,p ∈ L 2 (R 2 ) be we can write (Equation 1) as F = W ψ s,p f (x, θ ). Moreover, by making use of the quasiregular representation of the SE(2) group and denoting by ψ † (x) = ψ(−x), we can rewrite (Equation 4) as follows: which is a standard form to write the so-called group wavelet transform (refer to e.g., Führ, 2005;Deitmar and Echterhoff, 2014). Note that, in the interesting case (Equation 5), we have ψ † s,p = ψ s,p .
The SE (2) transform (Equation 4), together with the notion of group wavelet transform, has been studied in multiple contexts (refer to e.g., Weiss and Wilson, 2001;Antoine et al., 2004;Führ, 2005;Duits et al., 2007;Deitmar and Echterhoff, 2014;Dahlke et al., 2015 and references therein), and several of its properties are well-known. In particular, if W ψ f is a bounded operator from L 2 (R 2 ) to L 2 (R 2 × S 1 ), which happens e.g., for ψ ∈ L 1 (R 2 ) ∩ L 2 (R 2 ) by Young's convolution inequality and the compactness of S 1 , it is easy to see that its adjoint reads It is also well-known (Weiss and Wilson, 2001) that W ψ f cannot be injective on the whole L 2 (R 2 ), i.e., we cannot retrieve an arbitrary element of L 2 (R 2 ) by knowing its SE(2) transform. However, as shown in Duits (2005); Duits et al. (2007), and applied successfully in a large subsequent series of works (e.g., Duits and Franken, 2010;Zhang et al., 2016;Lafarge et al., 2021), it is possible to obtain a bounded invertible transform by extending the notion of SE(2) transform to mother wavelets ψ that do not belong to L 2 (R 2 ), or by simplifying the problem and reduce the wavelet analysis to the space of bandlimited functions, that are those functions whose Fourier transform is supported on a bounded set, whenever a Calderón's admissibility condition holds. Since our main point in this article is not the reconstruction of the whole L 2 (R 2 ), we will consider the SE(2) transform with this second, simplified, approach. In this case, the image of the SE(2) transform is a reproducing kernel Hilbert subspace of L 2 (R 2 × S 1 ) whose kernel will play an important role in the next section. For convenience, we formalize these statements with the next two theorems, and provide a sketch of the proof in the Appendix, even if they can be considered standard material.
THEOREM 2. For R > 0, let B R = {ξ ∈ R 2 : |ξ | < R} and let The SE(2) wavelet transform (Equation 4) for a mother wavelet ψ ∈ L 2 (R 2 ) is a bounded injective operator from H R to L 2 (R 2 × S 1 ) if and only if there exist two constants 0 holds for almost every ξ ∈ B R .
Before stating the next result, we repeat the observation of Weiss and Wilson (2001) and show, using this theorem, that the SE(2) transform cannot be a bounded injective operator on the whole L 2 (R 2 ). Indeed, using that we can see that the Calderón's function in condition (Equation 9) is actually a radial function . Hence, the lower bound in condition (Equation 9) cannot be satisfied on the whole R 2 by any ψ ∈ L 2 (R 2 ).
On the other hand, for the mother wavelet (Equation 5), since s 2 cos α dα.
From here, we can see that C ψ (ξ ) > 0 for all ξ ∈ R 2 , so, even if C ψ (ξ ) → 0 as |ξ | → ∞, we have that the lower bound in Equation (9) is larger than zero for any finite R. The next theorem shows how to construct the inverse of the SE(2) wavelet transform on bandlimited functions, and what is the structure of the closed subspace defined by its image. THEOREM 3. Let ψ ∈ L 2 (R 2 ) and R > 0 be such that Equation (9) holds. Let also γ ∈ L 2 (R 2 ) be defined by where χ B R (ξ ) = 1 ξ ∈ B R 0 otherwise , and let H R be as in Equation (8).
is a reproducing kernel Hilbert subspace of continuous functions of L 2 (R 2 × S 1 ), and the orthogonal projection Since W ψ (H R ) is a Hilbert space of continuous functions, it makes sense to consider its values on the zero measure set provided by the graph of a function : R 2 → S 1 , as in Equation (3). We can then provide a formal statement of the problem discussed in Section 1: for f ∈ H R and : R 2 → S 1 , reconstruct f using only the values W ψ f (x, (x)).
For this problem to be solvable, the graph G = {(x, θ ) ∈ R 2 × S 1 : θ = (x)} must be a set of uniqueness for W ψ (H R ). That is, the following condition must hold: if F ∈ W ψ (H R ) and F| G = 0 , then F = 0.
Indeed, if this was not the case, for any nonzero That is, we could not solve the problem (Equation 13).
Condition (Equation 14) is in general hard to be checked, and the formal characterization of the interplay between ψ and that makes it hold true is out of the scope of this article. However, in the next section, we will provide a technique for addressing (Equation 13) in a discrete setting, which will allow us to explore in Section 4, the behavior of this problem for various functions inspired by the feature maps measured in V1.

A RECONSTRUCTION ALGORITHM
In this section, we describe the discretization of the problem (Equation 13) which is used in the next section. Then, we introduce a kernel based iterative algorithm and prove its convergence to the solution whenever the solvability condition (Equation 14) holds.

Discretization of the Problem
The setting described in Section 2 can be discretized in a standard fashion by replacing L 2 (R 2 ) with C N×N , endowed with the usual Euclidean scalar product, circular convolution, and discrete Fourier transform (FFT), which amounts to replacing R with Z N , the integers modulo N. Explicitly, given f , ψ ∈ C N×N , With a uniform discretization of angles, i.e., by replacing S 1 with 2π M Z M = {0, 2π M , 4π M , . . . , 2π M−1 M }, we obtain the following discretization of the SE(2) transform with the mother wavelet (Equation 5): Thus, in particular, W ψ f ∈ C N×N×M . Note that here, for simplicity, we have removed the normalization used in Equation (5).
This allows us to process N × N digital images while retaining the results of Theorems 2 and 3 as statements on finite frames (refer to e.g., Casazza et al., 2013) since Plancherel's theorem and Fourier convolution theorem still holds. In particular, when computing numerically the inverse of W ψ using (i), Theorem 3, one has to choose an R > 0 so that Calderón's condition (Equation 9) for ∈ B R } and, due to the finiteness of the space, now it can be achieved for all R, i.e., without bandlimiting. However, since this is equivalent to the frame inequalities (Equation 24), the quantity B A defines actually the condition number of W ψ . Thus, in order to keep stability for the inversion, the ratio B A cannot be arbitrarily large (refer to e.g., Duits et al., 2007;Casazza et al., 2013). Once the parameters s, p, R are chosen in such a way that this ratio provides an acceptable numerical accuracy, one can then compute the projection P given by (ii), Theorem 3, on F ∈ C N×N×M , by making use of Fourier convolution theorem: We note at this point that this standard discretization, in general (for M different from 2 or 4), retains all of the approaches of Section 2 but the overall semidirect group structure of R 2 ⋊ S 1 . Let us now consider the discretization of the problem (Equation 13) and denote the graph of : If we denote by O the selection operator that sets to zero all the components of an F ∈ C N×N×M that do not belong to G , i.e., We can see that this is now an orthogonal projection of C N×N×M . Hence, problem (Equation 13) can be rewritten in the present discrete setting as follows: given f ∈ H R , find F ∈ C N×N×M that solves the linear problem The meaning of Equation (17) is indeed to recover W ψ f , and hence f , knowing only the values W ψ f (x, (x)). We propose to look for such a solution using the following iteration in C N×N×M : for F 0 = O W ψ f , compute F n = PF n−1 − O PF n−1 + F 0 , n = 1, 2, . . .
The idea behind this iteration is elementary: we start with the values of W ψ f selected by , we project them on the RKHS defined by the image of W ψ , and we replace the values on G of the result with the known values of W ψ f . The convergence of this iteration is discussed in the next section. Before that, we observe that Equation (18) can be seen as a linearized version of the Wilson, Cowan, and Ermentrout equation (Wilson and Cowan, 1972;Ermentrout and Cowan, 1980). Indeed, denoting by K = (1 − O )P, Equation (18) can be seen as a forward Euler scheme (time discretization) for the vector valued ordinary differential equation.
Apart from the absence of a sigmoid, this is indeed a classical model of population dynamics. Here, the "kernel" K is not just the reproducing kernel P, but it also contains the projection on a feature map O . Returning to the model of V1, here, the forcing term F 0 is the data collected by simple cells, while the stationary solution F, if it exists, is the full group representation of the visual stimulus defined as the solution to the Volterra-type equation F = KF + F 0 .

The Project and Replace Iteration
We show here that, whenever the problem (Equation 17) is solvable, the iteration (Equation 18) converges to its solution.
Since the argument is general, we will consider in this subsection the setting of an arbitrary finite dimensional vector space V endowed with a scalar product and the induced norm, and two arbitrary orthogonal projections P, Q. For an orthogonal projection P, we will denote by P ⊥ = 1 − P the complementary orthogonal projection. We will also denote by Ran the range, or image, and by Ker the kernel, of a matrix. We start with a basic observation, which is just a restatement of the solvability condition (Equation 14) as that of a linear system defined by an orthogonal projection, in this case, Q, on a subspace, in this case, characterized as Ran(P). The simple proof is included for convenience, and it can be found in the Appendix.
The problem posed by the system (Equation 19) is a problem of linear algebra: if we know that a vector F belongs to a given subspace Ran(P) ⊂ V, and we know the projection of F on a different subspace Ran(O) ⊂ V, can we recover F? The next theorem shows that, if the system (Equation 19) has a unique solution, such a solution can be obtained by the project and replace iteration (Equation 18). Its proof is given in the Appendix. THEOREM 5. Let V be finite dimensional vector space, and let P, Q be orthogonal projections of V. Given F ∈ Ran(P), set F 0 = Q F, and let {F n } n∈N , {H n } n∈N ⊂ V be the sequences defined by the project and replace the iteration.

NUMERICAL RESULTS
We present in this study the reconstruction results of the project and replace iteration on the restriction to feature maps of the SE(2) transform of real images. We have chosen eight 512 × 512 pixels, 8-bit grayscale digital images {f i } 8 i=1 ⊂ {0, . . . , 255} 512×512 , which are shown in Figure 1 together with their Fourier spectra. Note that, for processing, they have been bandlimited in order to formally maintain the structure described in Section 2. However, this bandlimiting has minimal effects, not visible to the eye, since the spectra have a strong decay: for this reason, the bandlimited images are not shown.
For the discrete SE(2) transform, we have chosen a discretization of S 1 with 12 angles, so that, with respect to the notation of the previous section, we have N = 512 and M = 12. We have also chosen the parameter values s = 51, p = 170, R = 252. The mother wavelet ψ and dual wavelet γ produced by these parameters are shown in Figure 2, top and center, on a crop of the full 512 × 512 space for better visualization. We stress that the stability of the transform and the numerical accuracy of the projection (Equation 16) depend only on the behavior of the Calderón's function, while the accuracy of image representation under bandlimiting with radius R depends on the decay of the spectra of the images. The Calderón's function C ψ , computed as in Equation (15), is shown in Figure 2, bottom left: the chosen parameters define a ratio B/A ≈ 6 · 10 3 , corresponding to a condition number for W ψ of less than 10 2 . In Figure 2, bottom right, we have shown in log 10 scale the multiplier χ B R /C ψ that defines the dual wavelet γ as in Equation (11), and in particular, we can see the bandlimiting radius, to be compared with the spectra of Figure 1. For visualization purposes, in Figure 3, we have shown in spatial coordinates the integral kernel defining the projection (Equation 16), which is the reproducing kernel for the discrete SE(2) transform. Its real and imaginary parts are shown on the same crop used to display the wavelet of Figure 2.
We have implemented iteration (Equation 18) for the restriction of this discrete SE(2) transform to different types of feature maps , shown in Figure 4. We will comment below each one of these cases. We have chosen to illustrate the effect of that iteration as follows. We have computed the sequence {H n } ν n=1 as in Equation (21), for a number ν of iterations, and applied the inverse SE(2) transform W * γ to each H n . This allows us to obtain real images that are directly comparable with the original ones. We then have shown W * γ H 1 , representing the first image that can be directly retrieved from the feature parameters, and W * γ H ν , that is the image obtained when we stopped the iteration. Moreover, as a measure of reconstruction error, we have considered the following rescaling of the Euclidean norm, at each step n ∈ {1, . . . , ν}: This adimensional quantity measures a % error obtained as the average square difference by pixel of an image f i in the dataset from its n-steps reconstruction W * γ H n [f i ], divided by the size of the admissible pixel range for 8 bit images, which is {0, . . . , 255}.

First Feature Map: Purely Random Selection of Angles
The first map that we have considered is a : Z N × Z N → Z M that, for each x ∈ Z N × Z N , simply chooses one value in Z M as a uniformly distributed random variable. This map is shown in the first line of Figure 4, left, and in the first line of Figure 4, right, we have shown an enlargement to the same crop at which the wavelets in Figure 2 and the reproducing kernel in Figure 3 are shown. In Figure 5, we have shown the images resulting from W * γ H 1 and W * γ H 1000 , and the evolution of the error (Equation 22) in log 10 scale, respectively in the left, center, and right column, for the first four images of Figure 1. In this case, we can see that the error n goes beyond 1%, indicated by 0 on the y-axis, in just about 500 iterations. As a remark, feature maps that are similar to such configurations are commonly encountered in rodents (refer to e.g., Ho et al., 2021 and references therein).

Pinwheel-Shaped Feature Maps
Next, we present the results for three selection maps that are pinwheel-shaped, as it is commonly observed for orientation and direction preference maps of V1 in primates and other mammals. These maps can be constructed as follows Petitot (2017): for ρ ∈ R + , let φ ρ : Z N × Z N → C be given by where Ŵ is a purely random process with values in [0, 2π). The maps ρ : Z N × Z N → Z M that we have considered are obtained by where angle(z) is the phase of a complex number z ∈ C, and ⌊t⌋ is the integer part of a real number t. The resulting maps are quasiperiodic, with a characteristic correlation length that corresponds to the fact that the spectrum of φ ρ is concentrated on a ring of radius ρ 2π . The main feature of those maps ρ is that they possess points, called pinwheel points, around which all angles are mapped, and these points are spaced, on average, by a distance of 2π/ρ (refer to e.g., Petitot, 2017 and references therein).
In the second, third, and fourth line of Figure 4, left, we have shown the resulting maps ρ with ρ, respectively, given by ρ = 0.8, ρ = 0.4, and ρ = 0.06. On the right column, we have shown an enlargement to the same crop used before.
The results of the iteration are presented, as described above, in Figures 6-8, whose structure is the same as in Figure 5 with the exception of the number of iterations, which is now larger.
In order to discuss these results, we first recall that the correlation length of orientation preference maps has been often related to the size of receptive fields, as a "coverage" constraint to obtain a faithful representation of the visual stimulus (refer to Swindale, 1991;Swindale et al., 2000;Bosking et al., 2002;Keil and Wolf, 2011;Barbieri et al., 2014 and references therein). However, neither mathematical proof of this principle in terms of image reconstruction has been given so far, nor has the word size received a more quantitative meaning within the models of type (Equation 1), and we have not tried to give any more specific meaning either. However, for the three cases that we present, by comparing the crops of the maps ρ in Figure 4 with the wavelet ψ in Figure 2, we can see that for ρ = 0.8, the correlation length of the map ρ is approximately similar to what we could call effective support of the receptive field, while for ρ = 0.4, we have that the area of influence of the receptive field does not include two different pinwheel points, and for ρ = 0.06, the two scales are very different. Heuristically, one could then be led to think that the reconstruction properties in the three cases may present qualitative differences. For example, that condition (Equation 14), or its discrete counterpart (Equation 20), may hold in the first case and may not hold in the third case.
As can be seen from the numerical results of the proposed algorithm, there is actually a difference in the behavior of the decay. For larger values of the parameter ρ, when the map ρ is more similar to the purely random selection described above, the decay of the error is faster, while for smaller values of ρ the decay is slower, but nevertheless, the error appears to be monotonically decreasing. In the presented cases, for ρ = 0.8, we can see in the right column of Figure 6 that in about 2,000 iterations the error decay appears to enter an exponential regime, which is rectilinear in the log 10 scale. We see in Figure 7 that it takes roughly twice as many iterations for ρ = 0.4 to enter the same regime. On the other hand, for ρ = 0.06, we can see in Figure 8 that after a relatively small number of iterations the decay becomes very small and does not seem to become exponential even after 10,000 iterations. However, visual inspection of the results (which "measures" the error in a different way than the Euclidean norm) in this last case, show that the starting image appears to be qualitatively highly corrupted, while the image obtained after the iteration was stopped is remarkably true to the original one and does not display evident artifacts away from the boundaries.

Selection of a Single Orientation: Deconvolution
The surprisingly good performance in the reconstruction problem for the last map 0.06 has motivated a performance test for an additional feature selection map, given by a function that is independent of x, hence selecting just one angle in the SE(2) transform. In the same fashion as seeing the purely random distribution as a limiting case for large ρ of the pinwheelshaped maps, this can be considered as a limiting case for small ρ. However, keeping only the values of W ψ f for a single angle θ concretely corresponds to performing a convolution with one function ψ θ , and aiming to reconstruct f is actually a deconvolution problem. In this case, the frequency behavior of the convolution filter is that of a Gaussian centered away from the origin, and its shape can be observed in the Calderón's function of Figure 2, which clearly shows the sum of 12 such Gaussians. We have shown the reconstruction results for this problem, using now 15,000 iterations, in Figure 9. The error decay seems to share, qualitatively, much of the same properties of the case ρ = 0.06.

CONCLUSIONS
In this article, we have proposed an elementary iterative technique to address the problem of the reconstruction of images from a fixed reduced set of values of its SE(2) group wavelet transform. We have formally defined these restrictions   in terms of cortical maps, as this problem is inspired by visual perception, since the SE(2) symmetry and the associated integral transform have proved to be relevant for mathematical modeling of V1, refer to e.g., (Petitot and Tondut, 1999;Paul C. Bressloff, 2002;John Zweck, 2004;Sarti, 2006, 2015). Moreover, the presented numerical simulations directly compare with the studies on the relationship between cortical maps and the efficiency of single cell encoding of information in terms of coverage (Swindale, 1991;Swindale et al., 2000;Bosking et al., 2002;Keil and Wolf, 2011;Barbieri et al., 2014). Indeed, as a result of the main theorem, and for cortical maps as formally defined in this article, we can see from the proposed implementation that, when a pinwheel structure is chosen, the relationship between the average distance between pinwheel centers and the size of receptive fields influences the invertibility of the restricted SE(2) transform. A possible interpretation of the proposed iteration with a kernel defined by the SE(2) group as a neural computation in V1 comes from the modeling of the neural connectivity as a kernel operation (Wilson and Cowan, 1972;Ermentrout and Cowan, 1980;Citti and Sarti, 2015;Montobbio et al., 2018), especially if considered in the framework of a neural system that aims to learn group invariant representations of visual stimuli (Anselmi and Poggio, 2014;Anselmi et al., 2020). A direct comparison of the proposed technique with kernel techniques recently introduced with radically different purposes in Montobbio et al. (2018) and Montobbio et al. (2019) shows, however, two main differences at the level of the kernel that is used: here, we need the dual wavelet to build the projection kernel, and the iteration kernel effectively contains the feature maps. On the other hand, a possible application is the inclusion of a solvability condition such as Equation (14) as iterative steps within learning frameworks such as those of Anselmi et al. (2019) and Anselmi et al. (2020). From the point of view of actual neural implementation, even if it was possible to see a formal analogy between the proposed algorithm and a neural field equation, we believe that more has to be done concerning the performances and that the proposed mechanism has to be considered as a much more elementary procedure than the ones that could take place in a real visual cortex. For example, while the proposed iteration is faster than typical convex optimization implementations used for other completion problems, it still requires a large number of iterations even in the exponential regime. Improvements could be searched in several directions, e.g., by including non-classical behaviors in neural modeling.
We would like to observe also that, since the proof of convergence of this technique is general, it could be applied to other problems with a similar structure. The computational cost essentially relies on the availability of efficient methods to implement the two projections that define the problem in the discrete setting, as it happens to be the case for the setting studied in this article. In particular, similar arguments could be applied to other wavelet transforms based on semidirect product groups R d ⋊ G, with G a subgroup of GL d (R) that defines what is sometimes referred to as local features, and to sampling projections obtained for example, but not only, from other types of feature maps : R d → G. From the dimensional point of view, in the discrete setting the selection of local features with a feature map can be seen as downsampling that allows one to maintain in the transformed space the same dimension of the input vector. This is often a desirable property and it is already commonly realized e.g., by the MRA decomposition algorithm of classical wavelets or by the pooling operation in neural networks. Moreover, its apparent stability of convergence seems to suggest that this operation can be performed a priori, without the need for a previous study of the solvability of the problem.
Several questions remain open after this study. Probably the most fundamental one is the characterization of those maps that, for a given mother wavelet ψ, satisfy the solvability condition (Equation 14). In terms of the study of the convergence of the project and replace iteration, it is plausible that one could obtain convergence under weaker conditions than Equation (20), even if maybe to a different solution, as it appears to happen in some of the numerical simulations presented.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and has approved it for publication.

FUNDING
This project has received funding from the European Union's Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie grant agreement No 777822, and from Grant PID2019-105599GB-I00, Ministerio de Ciencia e Innovación, Spain.

APPENDIX: PROOFS OF THE MAIN STATEMENTS
This section contains the proofs of the formal statements introduced in the text.

SKETCH OF THE PROOF OF THEOREM 2. By
Fourier Convolution Theorem, Plancherel's theorem, and the definition of H R , we have Thus, (9) is equivalent to the so-called frame condition for all f ∈ H R , which is equivalent [refer to e.g., Ali et al. (1993)] to W being bounded and invertible on H R .

SKETCH OF THE PROOF OF THEOREM 3. Observe first that
. Thus, recalling (7) and using Fourier Convolution Theorem, we can compute which proves (i). To prove (ii), one needs to show that the elements of W ψ (H R ) are continuous functions, that W ψ (H R ) W ψ W * γ is self adjoint and idempotent, and that (12) holds. The continuity can be obtained as a consequence of the continuity of the unitary representation (6) and, by the same arguments used to prove (i) we can easily see that W γ W * ψ = W ψ W * γ . On the other hand, (i) implies that W ψ W * γ W ψ f = W ψ f for all f ∈ H R , hence W ψ W * γ F = F for all F ∈ W ψ (H R ). Equation (12) can be obtained directly by (7) and the definition of W ψ .
PROOF OF LEMMA 4. The system (19) is of course solved by F = F, so we only need to prove that, for all F ∈ Ran(P), this solution is unique if and only if (20) holds. In order to see this, let S ∈ V be a solution to (19), and denote by E = F − S. Then E ∈ Ran(P) ∩ Ker(Q (PQ ⊥ ) k H 1 F n = Q ⊥ PF n−1 + F 0 = . . . = n k=0 (Q ⊥ P) k F 0 .
If (20) holds, then PQ ⊥ < 1 and Q ⊥ P < 1, because Ran(P) ∩ Ran(O ⊥ ) = Ran(P) ∩ Ker(O ) = {0}. Hence, the existence of the limits H = lim n→∞ H n and F = lim n→∞ F n is due to the convergence of the Neumann series [refer to e.g., Riesz and Sz.-Nagy (1990)]. Thus, getting rid of the notations F 0 = Q F and H 1 = PF 0 , we have that the limits H and F are the solutions to the equations We can now see that F = F, because this equation reads (1 − Q ⊥ P) −1 Q F = F and, using that F ∈ Ran(P), we can rewrite it as which is true by definition of Q ⊥ = 1 − Q. Moreover, by taking the limit in the first equation of (21), and using again that F ∈ Ran(P), we obtain H = PF = P F = F.
Finally, the convergence of (18) is exponential because, using the series expression (A2) for F = F, A similar argument applies to F − H n .