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) = ℝ² ⋊ 𝕊¹ of rotations and translations of the Euclidean plane via the group wavelet transform, is known only on a certain two-dimensional 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 linear-nonlinear-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 two-dimensional 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=(x1x2),y=(y1y2)∈ℝ2,

where the parameters s, p ∈ ℝ⁺ define the local (inverse) scale and spatial frequency, the angle θ ∈ [0, 2π) defines the local direction and x ∈ ℝ² 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; Citti and 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 θ ∈ 𝕊¹ 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 Θ : ℝ² → 𝕊¹ 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,θ)∈ℝ2×S1:θ=Θ(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 non-classical 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, 2020; Montobbio 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 bell-shaped curve that shows the response of the cell to different orientations. That is, it describes how a cell that has a given orientation preference θ₀ 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, 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 uncertainty-related 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.

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²(ℝ²) 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¹(ℝ²) ∩ L²(ℝ²) by

and, as customary, we will use the same notation for its extension by density to the whole L²(ℝ²). We will also denote by * the convolution on ℝ²:

Let 𝕊¹ be the abelian group of angles of the unit circle, which is isomorphic either to the one dimensional torus 𝕋 = [0, 2π) or to the group SO(2) of rotations of the plane ℝ². The group SE(2) = ℝ² ⋊ 𝕊¹ is (refer to e.g., Sugiura, 1990, Ch. IV) the semidirect product group with elements (x, θ) ∈ ℝ² × 𝕊¹ 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 ℝ² ⋊ 𝕊¹.

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:S1→U(L2(ℝ2)) the unitary action by rotations of 𝕊¹ on L²(ℝ²):

Let ψ ∈ L²(ℝ²), and denote by ψ_θ = R(θ)ψ. The SE(2) wavelet transform on L²(ℝ²) with the mother wavelet ψ is

In terms of this definition, we can then see that if we choose s, p ∈ ℝ⁺, and let ψs,p∈L2(ℝ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²(ℝ²) to L²(ℝ² × 𝕊¹), which happens e.g., for ψ ∈ L¹(ℝ²) ∩ L²(ℝ²) by Young's convolution inequality and the compactness of 𝕊¹, 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²(ℝ²), i.e., we cannot retrieve an arbitrary element of L²(ℝ²) 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²(ℝ²), 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²(ℝ²), 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²(ℝ² × 𝕊¹) 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 BR={ξ∈ℝ2:|ξ|<R} and let

The SE(2) wavelet transform (Equation 4) for a mother wavelet ψ ∈ L²(ℝ²) is a bounded injective operator from HR to L²(ℝ² × 𝕊¹) if and only if there exist two constants 0 < A ≤ B < ∞ such that

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²(ℝ²). Indeed, using that

we can see that the Calderón's function in condition (Equation 9) is actually a radial function

which, by Plancherel's theorem, satisfies ∫0∞Cψ(ρ)ρdρ=||ψ||L2(ℝ2)2. Hence, the lower bound in condition (Equation 9) cannot be satisfied on the whole ℝ² by any ψ ∈ L²(ℝ²).

On the other hand, for the mother wavelet (Equation 5), since ψ^(ξ)=1s2πe-|ξ+(p0)|2/2s2, we have

From here, we can see that C_ψ(ξ) > 0 for all ξ ∈ ℝ², 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²(ℝ²) and R > 0 be such that Equation (9) holds. Let also γ ∈ L²(ℝ²) be defined by

where χBR(ξ)={1ξ∈BR0otherwise, and let HR be as in Equation (8). Then

Since Wψ(HR) 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 Θ : ℝ² → 𝕊¹, as in Equation (3). We can then provide a formal statement of the problem discussed in Section 1:

For this problem to be solvable, the graph GΘ={(x,θ)∈ℝ2×S1:θ=Θ(x)} must be a set of uniqueness for Wψ(HR). That is, the following condition must hold:

Indeed, if this was not the case, for any nonzero F∈Wψ(HR) that vanishes on GΘ, the function fF=Wγ*(Wψf+F)∈HR would be different from f but W_ψf_F would coincide with W_ψf on GΘ. 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.

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.

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 {fi}i=18⊂{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 𝕊¹ 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³, corresponding to a condition number for W_ψ of less than 10². In Figure 2, bottom right, we have shown in log₁₀ scale the multiplier χ_BR/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 {Hn}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γ*H1, 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γ*Hn[fi], divided by the size of the admissible pixel range for 8 bit images, which is {0, …, 255}.

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; Citti and 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 ℝ^d ⋊ G, with G a subgroup of GL_d(ℝ) 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 Θ : ℝ^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.

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

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

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.

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.