Topology-preserving smoothing of retinotopic maps

Retinotopic mapping, i.e., the mapping between visual inputs on the retina and neuronal activations in cortical visual areas, is one of the central topics in visual neuroscience. For human observers, the mapping is obtained by analyzing functional magnetic resonance imaging (fMRI) signals of cortical responses to slowly moving visual stimuli on the retina. Although it is well known from neurophysiology that the mapping is topological (i.e., the topology of neighborhood connectivity is preserved) within each visual area, retinotopic maps derived from the state-of-the-art methods are often not topological because of the low signal-to-noise ratio and spatial resolution of fMRI. The violation of topological condition is most severe in cortical regions corresponding to the neighborhood of the fovea (e.g., < 1 degree eccentricity in the Human Connectome Project (HCP) dataset), significantly impeding accurate analysis of retinotopic maps. This study aims to directly model the topological condition and generate topology-preserving and smooth retinotopic maps. Specifically, we adopted the Beltrami coefficient, a metric of quasiconformal mapping, to define the topological condition, developed a mathematical model to quantify topological smoothing as a constrained optimization problem, and elaborated an efficient numerical method to solve the problem. The method was then applied to V1, V2, and V3 simultaneously in the HCP dataset. Experiments with both simulated and real retinotopy data demonstrated that the proposed method could generate topological and smooth retinotopic maps.


Introduction
Retinotopic mapping, i.e., the mapping between visual input on the retina to neuronal activations in visual cortical areas, is one of the central topics in visual neuroscience [1]. Although initially discovered in neurophysiology [2], Blood Oxygenation Level Dependent (BOLD) functional magnetic resonance imaging (fMRI) provides an excellent noninvasive tool to perform in vivo retinotopic mapping on human observers [3,4]. The procedure typically consists of (1) recording BOLD fMRI cortical activations generated by carefully designed slowly moving visual stimuli on the retina to uniquely encode the visual space [5,6], (2) co-registering the fMRI activations with structural MRI, (3) flattening the 3D structural MRI data to obtain visual cortical surface, and (4) decoding the retinal visual coordinates underlying the observed fMRI time series at each location on the cortical surface [3,7]. In the last two decades, the development of a variety of retinotopic mapping techniques has greatly extended our understanding of visual cortical organization in normal and abnormal populations [3,[8][9][10][11][12][13].
Although the exact shape of visual inputs on the retina is not preserved on the cortical surface, neurophysiology studies [14][15][16][17] have shown that the mapping preserves local neighborhood geometric relationships. More specifically, neighboring points on the cortical surface should have neighboring retinal visual coordinates [18]. We adopted the term topology preserving (or topological) to refer to the preservation of neighborhood relationships of the observations in the output space [19]. In this work, the topological condition or topology-preserving are not related to the surface genus, handles or holes. The concept of topological condition is illustrated in Fig 1 for retinotopic mapping. Fig 1A and 1B show an example of topological mapping that preserves the structure of the source polygon (on the cortical surface) in the target (retina) domain ( Fig 1B). Fig 1A and 1C show an example of non-topological mapping that breaks the source polygon structure after mapping: f i is out of the polygon edge � f j f k , resulting in what is called a flipped triangle Δf i f j f k (the orientations of the triangle are different in the source and target domains).
The topological condition is a natural consequence of the hierarchical organization of the visual system: there cannot be duplicated neuronal representations of the same retinal location within each visual area (otherwise the visual area should be further divided into more subareas). Although neurophysiological studies on animals have found that the retinotopic maps in lower visual areas are topological [2], the "raw" retinotopic maps derived from BOLD fMRI recordings are usually not because of the low signal-to-noise ratio and spatial resolution of the technology [3]. Although the quality of fMRI data has improved during the last three decades, the signal-to-noise ratio is still relatively low even with ultra-high field MRI systems [20]. A typical raw retinotopic map of V1 from the state of the art dataset and method [20] still contains about 20% flipped triangles ( Fig 1D; Table A in S5 Text). Violations of the topological condition not only make the results from fMRI inconsistent with those from neurophysiology but also make it difficult to perform quantitative analysis on retinotopic maps. Typically, fMRI signals are analyzed on a voxel-by-voxel basis. The most prominent method is based on the population receptive field (pRF) model, which decodes perception parameters of the fMRI signals at each voxel. However, influenced by low SNR and low spatial resolution, the retinotopic maps from the analysis are typically not topological. In the past decade, significant efforts have been devoted to improving their accuracy and stability. For instance, the reliability of perception size was discussed [12,13], and methods that could decode perception parameters by model-free back projections with potentially more precise perception shape were proposed [21][22][23]. Although they have greatly improved the pRF model, none of the methods has explicitly investigated the topological condition.
In addition, a number of methods have been proposed to correct or reduce topological violations in fMRI-based retinotopic maps from the pRF analysis. Yet none of them has fully solved the problem. For instance, traditional spatial smoothing methods, including Gaussian smoothing [14,24], Laplacian smoothing [25], mesh-spectrum-based smoothing [26], cannot guarantee the topological condition because retinotopic maps are defined on two-dimensional surfaces, yet these methods only smooth the maps along one dimension (e.g., polar angle) at a time. Another promising approach registers the noisy "raw" retinotopic maps to an ideal template and assigns coordinates based on the registration function/morphing [8]. Although it can generate smooth and topological retinotopic maps, the method relies on diffeomorphic registration of a noisy parametrized surface to a predefined template, which is not easy to optimize, especially when the data are noisy. Furthermore, even if the registration issue can be solved, defining the right template remains a challenging problem. Other methods based on fitting predefined algebraic models (such as Schwartz's complex "log" model [27] or Schira's "Double-Sech" model [9]) can guarantee topological condition. However, the models introduce a lot of assumptions on the data [8,28]. The topological condition was considered in our previous work [29], but it was only applied to V1 with limited validation. And, because the manually drawn boundaries were not accurate, the results were not very precise near the boundary.
We propose a topological smoothing method to generate smooth and topological retinotopic maps for multiple visual regions without assuming algebraic models. A geometric concept, Beltrami coefficient (BC) [30], is adopted to formulate the optimization problem to define, quantify, and ensure the topological condition [31]. The Beltrami coefficient is a complex number that can quantify triangle-wise geometric conformality (i.e., the angle preserving property) and topological condition for surface-to-surface mapping. When we consider the topological condition in the mapping from the polygon in Fig 1A to 1B, the overall topological condition can be checked for each triangle attached to point P i . Namely, one may check whether its mapping result, f i , is moved out of any of the polygon edges, and � f n f j . If f i is moved out of any of the edges, e.g., � f j f k , the mapping is not topological. Without losing generality, when considering the topological condition for polygon edge � f j f k , we can construct a coordinate system (Fig 2A and 2B). In Fig 2B, we illustrate three mapping results for P i : Case (1) if f i coincides with the green dot such that triangles ΔP i P j P k and Δf i f j f k are similar (ΔP i P j P k~Δ f i f j f k ), the magnitude of the Beltrami coefficient associated with triangle ΔP i P j P k is zero (|μ| = 0). This is the case of conformal mapping. Case (2) if f i is above the f (1) axis (f i (2) >0), shown as the yellow dot, the magnitude of the Beltrami coefficient associated with triangle ΔP i P j P k is greater than zero but less than one (0<|μ|<1). This is the case of quasiconformal mapping with respect to edge � f j f k . And Case (3) if f i is below the f (1) axis, shown as the red dot, the magnitude of the Beltrami coefficient associated with triangle ΔP i P j P k is greater than one (|μ|>1) and the mapping is no longer topological. In summary, the complex-valued Beltrami coefficient encodes relative mapping information for each triangle and can be used to quantify the topological condition [32,33].
To ensure the topological condition, we defined an objective function such that, after optimally adjusting the retinal visual coordinates of each point on the cortical surface, the maximum magnitude of the Beltrami coefficients on the new retinotopic map is less than one (i.e., topological). An efficient numerical method was then developed to solve the optimization problem.
Our proposed method can work on multiple visual areas simultaneously without assuming any prior model and allow us to draw more precise boundaries between visual areas. To validate the method, we conducted experiments on both synthetic and real retinotopic map data from the Human Connectome Project (HCP) [34] and compared the performance of our method with that of existing methods based on the quality of fits to the raw fMRI time series.

HCP retinotopy data
The Human connectome project (HCP) [34] provides a large publicly available retinotopy dataset collected on 7T MRI scanners. The data collection, conducted on 181 healthy young adults (22-35 years; 109 females, and 72 males) with normal or corrected-to-normal visual acuity, involved carefully designed retinotopy stimuli and resulted in a substantial amount of fMRI data (30 min, 1,800 time-points) acquired at very high spatial and temporal resolutions (1.6 mm isotropic voxels, 1-second temporal sampling). The dataset provides an exciting opportunity to evaluate our method.

Overview of the pipeline
The overall proposed processing pipeline is shown in Fig 3. In the HCP dataset, the visual stimuli consist of rotating wedges, expanding/shrinking rings, and moving bars [3,35]. The carriers of the stimuli are made of dynamic color textures that can better activate visual neurons ( Fig 3A). We denote a point in the visual field by v = (v (1) ,v (2) )2R 2 , where v (1) is the eccentricity (distance to the fovea in degrees of visual angle), and v (2) is the polar angle relative to the positive horizontal line (Fig 3B). Both high-quality structural MRI and fMRI scans were acquired by the HCP group (Fig 3C and 3D). The processing pipeline consists of fMRI preprocessing [36], compressive population receptive field decoding [3], and topological smoothing. The raw retinotopic maps from pRF decoding (RM; Fig 3G) contain many topological violations ( Fig 3H). Topological smoothing is the proposed method in this paper to fix topological violations and smooth the retinotopic maps ( Fig 3I). In the following subsections, we briefly describe fMRI preprocessing and pRF decoding (although we used the publicly available decoded (raw) retinotopic maps from [20]), and then introduce the topological smoothing method. We list the key definitions and symbols in Table A in S1 Text.

fMRI preprocessing
The goal of fMRI preprocessing is to detect, in a robust, sensitive, and valid way, the time series of brain activations of the voxels on the visual cortical surface that are associated with the visual stimuli. The HCP preprocessing pipeline [36] is standardized. First, a cortical surface is extracted from structural MRI using Freesurfer [37] and resampled ( Fig 3F). Then, the raw fMRI data from each imaging session are co-registered across time to reduce the influence of head movements and other motion artifacts ( Fig 3E). Finally, the co-registered fMRI data are projected onto the cortical surface. During the projection, spatial smoothing is applied along the cortical surface to improve the quality of the fMRI signal. The preprocessing eventually generates a resampled cortical surface as well as the fMRI activation time series on each point of the cortical surface.

pRF decoder
We briefly describe the population receptive field model (pRF) [3,38] to introduce some necessary notations. For readers familiar with pRF model, the visual coordinates we used here are related to the x-y coordinates by x = v (1) cos v (2) ,y = v (1) sin v (2) . For each voxel P = (X,Y,Z)2R 3 on the visual cortical surface, pRF is used to determine its receptive field, including its center location v and size σ in the visual field. Namely, the collection of the mapping between P and (v,σ) constitutes a raw retinotopic map. Assuming that the voxel's population receptive field model is r(v 0 ;v,σ) and the hemodynamic function is h(t), the predicted fMRI signal of the voxel can be written as:ŷ ðv; s; nÞ ¼ bð R rðv 0 ; v; sÞsðt; v 0 Þdv 0 Þ n � hðtÞ; ð1Þ where β is the activation level, n is the exponent of the power function. We used the standard population receptive field model, i.e., r v 0 ; v; s ð Þ ¼ exp À ðx 0 À xÞ 2 þðy 0 À yÞ 2 (1) cos v (2) ,y = v (1) sin v (2) (similarly for x 0 ,y 0 ). The center of the receptive field v and population receptive field size σ can be estimated by minimizing the Least Square difference between the measured and predicted fMRI activation levels: where y(P) is the fMRI signal at voxel P. Iterations of this procedure across all the voxels on the visual cortical surface generate the raw retinotopic map. The quality of fit of the pRF model for voxel can be evaluated by

Topological smoothing
Quasiconformal mapping and Beltrami coefficient. Most surface mappings can be modeled by quasi-conformal mapping, a generalization of conformal mapping. Specifically, the local deformation of a quasiconformal map can be characterized by its associated Beltrami coefficient [39]. Mathematically, f:C!C is quasiconformal if it satisfies the Beltrami equation, for some complex-valued Lebesgue measurable function μ f satisfying (2) and � u ¼ u ð1Þ À iu ð2Þ . Intuitively speaking, if kμ f k 1 <1, the map preserves the orientations of all triangles and therefore preserves the topological condition. However, if kμ f k 1 >1 the topological condition is violated because the orientations of some triangles are flipped.
Surface flattening. Let f r : v i 7 ! P i be the retinotopic map between the retina and V1 cortical surface. To use the tool of the Beltrami coefficient, we need to flatten the cortical surface from 3D to a 2D unit disk. Because topological proprieties are transitive (that is if mappings f r and c are both topological, the composed mapping f r �c is also topological) and mutual (if mappings f r is topological, f r −1 is also topological), we can reformulate the topological condition by mapping the 3-dimensional cortical surface to a 2D parametric space with a topological transformation [37,40]. Here, we performed a conformal mapping (angle preserving) of the visual cortical surface to a parametric space c: Although it is angle preserving and does not introduce any angle distortions, conformal mapping may introduce big metric distortions. To reduce metric distortions, we cut a geodesic disk patch that contains the region of interest from the cortical surface by: (1) picking a point on the cortical surface that roughly corresponds to the fovea as the center of the disk, (2) calculating the geodesic distances from all points on the cortical surface to the center of the disk [41], and (3) keeping the points of the cortical surface whose geodesic distances to the center of the disk are within a certain value. The geodesic disk patch (the region-of-interest-ROI-shown in gray color in Fig 4B) is mapped to a unit disk ( Fig 4C). We refer the readers to [31,42] for conformal mapping and [43] for the computation of geodesic distance.
In this work, we focus on the 2D retinotopic mapping f = f r where w is the weight function, r = (@/@u (1) ,@/@u (2) ) is the gradient operator, s is a positive scalar emphasizing the smoothness and m^f is the Beltrami coefficient off written as (see S2 Text): In Eq 5, we have rewritten a point in the 2D parametric domain as a complex number u = u (1) +iu (2) and the corresponding visual coordinate asf ¼f ð1Þ þ if ð2Þ . Because we will only discuss the problem in 2D in the subsequent sections of the paper, we interchangeably regard each complex number as a position vector or vice versa, that is, u = u (1) +iu (2) is regarded as a complex number or as a position vector u = (u (1) ,u (2) ). We denote a position vector with a bold u and a complex number with an italic u. The same notation applies to v. The Beltrami off can be used to quantify the Beltrami of f r −1 because the Beltrami coefficient does not change with the conformal mapping, that is, m f ¼ m f r À 1 if c is a conformal mapping (see S2 Text for mathematical details). Intuitively, the Beltrami coefficient measures the angle distortion in a mapping (see S4 Text), incorporating a conformal mapping c does not introduce any angle distortion and therefore does not change the original Beltrami coefficient. We can directly estimate the Beltrami coefficient of the retinotopic mapping in the 2D parametric domain on the unit disk based on this property. Model solver. To solve Eq (4) efficiently, we divide the problem into two subproblems, namely smoothing and topological projection. Specifically, during smoothing, we ignore the topological condition and find a smoothed map by, which can be solved by Laplacian smoothing [44,45] efficiently. For the second subproblem, we ignore the smoothness term and try to make the results topological.
Since we require the results to be both topological and smooth, we solve the subproblems iteratively. Namely, the output of topological projection is fed to smoothing, and so on so forth.
The key challenge is to solve Eq (7) efficiently. Because the topological condition is formulated with the Beltrami coefficient (kmv k 1 < 1), we need tools to compute the Beltrami coefficient and recover the function after processing the Beltrami. In the following, we describe the Beltrami coefficient computation when a function is given and the function recovery when the Beltrami coefficients are given in the discrete setting. Then we elaborate on the details about solving each subproblem.
Model discretization. We now convert the problem to the discrete domain because a cortical surface, S, is typically described by vertices V S and triangular faces F S : S = (F S ,V S ). The visual field is also discretized. Retinotopic mapping can be reformulated as finding the corresponding visual coordinate v i in the visual field for every vertex P i 2V S . We denote W i as w(v i ) to emphasize/deemphasize voxels with high and low quality of pRF fits.
Discrete Beltrami coefficient computation. Given the analytical form of a mapping function f, we can directly compute the complex-valued Beltrami coefficient μ f according to Eq (5). However, in the discrete case, the function is only defined on vertices, i.e., we only know the mapping function f that maps the vertices on the unit disk to the vertices on the visual field ( Fig 5A), that is, v 1 = f(u 1 ), v 2 = f(u 2 ), and v 3 = f(u 3 ). To approximate the derivatives, f is interpolated linearly on each triangle, i.e., for any point u within the triangle, 3 . The coefficients, B 1 , B 2 , and B 3 are called the barycentric coefficients. Specifically, B 1 (similarly for B 2 and B 3 ) is the ratio between the areas of triangles Δuu 2 u 3 and Δu 1 u 2 u 3 : . Now we can compute f's partial derivatives to u (1) and u (2) and therefore the Beltrami coefficient μ f according to Eq (5) (see S3 Text for the explicit expression). Clearly μ f is a constant within each triangle, as we approximate the mapping function linearly within each triangle.
In the discrete case, the mapping function is interpolated on each triangle. The gradient operator rf (1) (u) can be written out and is a constant vector on each triangle. The divergence r�G is approximated on the dual of each vertex (a convex cell consisted of the circumcenters of the neighbor triangles). Specifically, consider a vertex with its neighbors (Fig 5B). We assume that r�G is a constant on the dual of the center vertex, i.e., polygon D, enclosed by circumcenters of the neighbor triangles. Namely, we approximate the term with According to the divergence theorem [46], the divergence R D r�Gdσ can be written as the integral on the boundary, i.e.
where n is the dual cell unit normal. Since G is a constant vector on each triangle (in the discrete case), the integration can be written as: where G T j is the value of G within triangle T j , and s i = (u k −u j ) ? denotes a vector perpendicular to u k −u j with the same norm as u k −u j . Substituting Eq (12) into Eq (11), we have a set of linear equations for each f (1) (u i ) and its neighbors. Finally, we can write r�Arf (1) =0 in a matrix form (see S3 Text for the explicit form) and solve f (1) efficiently with the Dirichlet boundary condition. Similarly, f (2) can be solved efficiently. Next, we describe how we can optimize Eq (7) by solving two subproblems iteratively.

Solving the first subproblem: Laplacian smoothing
The first subproblem can be solved by Laplacian smoothing of f (1) and f (2) separately. To get a smooth scalarf ð1Þ from a noisy f (1) , we solve the following Laplacian smoothing equation [44,45]:f where the r is defined in the parametric domain, i.e., r = (@/@u (1) ,@/@u (2) ). To find the solution, we set the partial derivatives of Eq (13) to zero. It leads to the following equation ðÀ r � r þ 2s IÞf ð1Þ ¼ 2sf ð1Þ . It can be solved efficiently with linear algebra because r 2 = r�r can be written in a matrix form (the details are provided in S3 Text), and σ is optimally chosen to minimize the generalized cross-validation score [47]. Similarly, f (2) can be smoothed tof ð2Þ . Since the Beltrami coefficients quantify local deformations, Laplacian smoothing on visual coordinates will smooth the Beltrami coefficients. If most triangles are correctly orientated, the smoothing will also reduce the number of flipped triangles.

Solving the second subproblem: Topological projection
The second subproblem is to make f topological. It means that the associated Beltrami coefficients must satisfy the condition kμ f k 1 <1 on the entire cortical surface. To ensure topological condition, for any mapping whose |μ| is greater than 1, we scale its magnitude by μ 0 = μ/(|μ| +�). The procedure is called "topological projection" [48,49] or "chopping" [50]. The goal is to find a topological mapping that is close to the original non-topological mapping. The new μ 0 shares the same argument of μ, and|μ 0 | is less than 1 when we choose a small positive �. Because the Beltrami coefficient uniquely encodes the mapping [31,51], the projection process changes the mapping when we solve equation Eq (11). The process of making the mapping topological can be illustrated with the example in Fig 2B, in which the target f i is below the horizontal axis, and the magnitude of the associated Beltrami coefficient is greater than one (jm Df i f j f k j > 1). If we shrink the Beltrami coefficient for Δf i f j f k and retain the Beltrami coefficients for all the other triangles, the new Beltrami coefficient jm 0 Df i f j f k j is less than 1. When we solve the mapping function according to the new Beltrami coefficient (Eq (11)), f i is moved and placed above the horizontal axis, making the triangle Δf i f j f k topological. Note that other target points, f l , f m , and f n , also move to maintain the Beltrami coefficients of their associated triangles, e.g., when f i moves up, f l moves accordingly so that the Beltrami coefficient of triangle Δf i f k f l stays the same. In summary, given an initial non-topological mapping, we can project the Beltrami coefficient and solve Eq (11) to obtain a topological mapping that is close to it.
Removing phase jumping. We adopt the widely used polar coordinate system for the visual field, i.e., eccentricity r = v (1) and polar angle θ = v (2) 2[0,2π). The polar coordinate system has many advantages, but the polar angle θ is not continuous numerically near the x axis even when the mapping itself is. As shown in Fig 6A, points v 1 and v 2 share the same eccentricity and are both very close to the x axis, but their polar angles are very different because the polar angle jumps from 2π to 0 near the positive x axis. We call the phenomena "phase jumping". It makes the coordinates around the x axis discontinuous.
Phase jumping is not an intrinsic property of retinotopic mapping. The issue can be addressed with a transformed coordinate system: If we rotate the polar coordinate system, the phase jumping axis can be shifted. To avoid phase jumping, we shift the zero polar angle to the negative x axis for the left hemisphere ( Fig 6B). Now phase jumping occurs near the negative x axis. Because retinotopic mapping of the left hemisphere is in Quadrants I and IV, the shift

PLOS COMPUTATIONAL BIOLOGY
Topological smoothing of retinotopic maps allows us to avoid phase jumping. For the right hemisphere, retinotopic mapping is in Quadrants II and III; we do not need any shift. By doing this, we can avoid unnecessary visual coordinate discontinuities.
Removing phase reversal. We apply another linear transformation on the polar angle to remove phase reversal across visual areas. Phase reversal is an important intrinsic property of retinotopic mapping. As shown in Fig 6C, if one walks from V3d to V3v along the yellow arc, the corresponding polar coordinate v (2) on the visual field first increases and then decreases. We assign a visual area with a positive (negative) visual field sign if the polar coordinate increases (decreases) when we move clockwise on an equal eccentricity curve, e.g., the yellow arc in Fig 6C. A phase reversal across the boundary of two cortical areas causes a change of visual field sign and therefore breaks the topological condition for the points near the boundary between visual areas. To make our proposed smoothing method applicable to multiple visual areas simultaneously, we apply another linear transformation to remove the visual field sign changes across visual area boundaries. The combination of two transformations, one for removing phase jumping and the other for eliminating visual field sign changes, is summarized in Table 1. We call the transformed polar angle "extended polar angle". Because polar angles and visual area labels have a one-to-one correspondence with the extended polar angles, one can recover the original polar angles and visual area labels from the extended polar angles. With accurate extended polar angles, one can draw accurate interior boundaries between different visual areas (V1 and V2, V2 and V3, etc).
In addition, these transformations make the smoothing more precise, especially near the boundaries of visual areas, because they remove the sharp polar angle changes across the boundaries.
Algorithm. We now summarize the overall topological smoothing algorithm in Alg.  Table 1 to get adjusted visual coordinates v 0 v i v 0 i ; assign the initial retinotopic coordinates for each vertex on the unit disk repeat Compute Beltrami coefficient μ for the mappingf : u !v using Eq (5).

Left Hemisphere Right Hemisphere
V1v

Synthetic data
We first evaluated the performance of the algorithm on synthetic data generated with the following function, We chose ln function (k = 0.5 in this paper) as it is a good approximation of the retinotopic mapping from the visual field to the flattened cortical surface [54]. Note that in the original study, Schwartz used the model to map from the visual field to the parametric coordinates of a flattened cortical surface. It is the inverse function that we studied in this paper. We followed the convention only in this subsection, i.e., generating a smoothedûðvÞ from a noisy u(v), instead of a smoothedvðuÞ from a noisy v(u). Note that, because the topological condition is reciprocal, the conclusion based on this approach still applies to the inverse function. More specifically, the visual field grid mesh is first created (Fig 7A) with the following configurations: eccentricity spans from 0˚to 4.5˚and polar angle ranges from −π/2 to π/2. We then sampled 12 points in each direction and formed the quadrilateral grid mesh (144 points in total). Based on the quadrilateral mesh, we created the triangular faces by connecting diagonal points for each quadrilateral. Then we mapped the triangular mesh to the parametric space (Fig 7B), according to Eq (14). Ideally, the mapping points should be located on the grid (Fig 7B). We then added a small amount of noise with a peak signal-noise ratio (PSNR) of 10 (Fig 7C), and a large amount of noise with a PSNR = 5 (Fig 7D). We

PLOS COMPUTATIONAL BIOLOGY
Topological smoothing of retinotopic maps neighborhood. All existing smoothing procedures treat the u (1) and u (2) coordinates separately. In this section, we fixed the smooth parameter s = 0.001 (Eq (7)) and the boundary values were initially inferred from the Average Smoothing results and boundary change tolerance is � = 0.01. We list the performance metric, i.e., the deviation from the ground truth of all smoothers in Tables 2 and 3. Value deviation is the average difference relative to the ground truth. Angle distortion (the angle spanned by the tangent vector of u (1) 's contour curves and the tangent vector of u (2) 's contour curves, at the same location; see S4 Text for details) quantifies the local angle changes (in degree unit) after the mapping, and the number of flipped trianglesT measures violations of the topological condition. All the experiments were repeated 50 times to reduce fluctuation caused by random noise. The number of flipped trianglesT is the median of the 50 experiments. In addition, we ran a non-parametric permutation test [55] to conduct pairwise comparisons of the results from the proposed method and each of the other methods (including no smoothing). The null hypothesis is that the results from the proposed smoother are not significantly different from those of the other methods. We report the p value, i.e., the probability that the null hypothesis was true. In brief, if the p value is less than 0.01, the results from the proposed smoothing method are statistically different from those of the other methods.
We found that (1) all smoothing methods reduced value deviation, angle distortion, and number of flipped triangles (Tables 2 and 3); (2) the proposed method achieved the best topological result and the smallest angle deviation (all p<0.00001); (3) when the noise was relatively small ( Table 2), violations of the topological condition were largely corrected by most of the methods, and the average smoothing method achieved very good topological results; (4) when the noise was relatively high (Table 3), the average smoothing method lost its ability to correct topological violations; and (5) the proposed method always performed best in angle preserving, which was not a surprise as the magnitude of Beltrami coefficient |μ| is related to angle distortion of the mapping (see S4 Text), and we penalized big |μ| 0 s in the proposed method (Eq Table 2

Method
Value

Synthetic Data II
We performed another set of synthetic data experiment in which noise was directly added to the fMRI signals. We started with the ground-truth parameters: receptive field centers spanning the right visual space [0:0.08:6.4]˚×[0,0.08: 4.8]˚; receptive field size = 0.4˚; gain = 1; and exponent of the power function n = 0.5. Then, we followed Kay's pRF model (AnalyzePRF [38,56] with stimulus from the HCP group, 200×200 in spatial resolution with 1800 seconds) to generate the ideal fMRI signal. We further normalized the fMRI signal to unit variance, added a specific amount of white noise, decoded the pRF model for the fMRI signal (with noise), and computed the number of flipped triangles. We varied the noise level until the percentage of flipped triangle was about 20%. (Specifically, in our experiments, noise with a 2.5 standard deviation induced about 22% flipped triangles). With pRF decoding results at this noise level as input, we applied the topological smoothing method and generated new visual coordinates. Then we compared our work with the pRF model and several post-processing methods (Table 4). We found that (1) the proposed method can fully ensure the topological condition, with 0 flipped triangles, while other smoothing methods had at least 6% of flipped triangles; (2) the proposed method had minimal angle distortion; and (3) all smoothing methods improved the receptive field center estimates, with the median smoothing method achieving minimal center

V1 retinotopic map
Although 7-Tesla MRI systems have dramatically improved the signal-to-noise ratio, the retinotopic maps in the HCP dataset are still quite noisy and non-topological. We first applied our proposed smoother on the V1 retinotopic maps of all observers in the HCP dataset.
Since there is no ground truth, we evaluated the performance of the method by the goodness of fit of the predicted fMRI signals with the smoothed visual coordinates. Specifically, we first decoded the visual coordinates v, population perception field size σ, BOLD response level g for the 181 observers by the pRF method [3]. Then we applied several smoothing methods on their V1 retinotopic maps. Based on the smoothed coordinates, we computed fMRI time series predictions at every vertex of the cortical surface. To evaluate the impact of changing the receptive field centers, we set the response level g and the exponent of the power function n the same as the original pRF in computing fMRI time series predictions.
We measured the goodness of fit of the smoothing methods to the fMRI time series (with the code from AnalyzePRF [38,56]). The results are listed in Table 5, with the following metrics: (1) the average number of flipped triangles across subjects, (2) the average percentage of the non-topological area Ar n , (3) the average Normalized-Root-Mean-Square-Error (NRMSE) of the fits, and (3) the average variance � R 2 accounted for across all the vertices in V1. We observe that: (1) all the smoothing methods reduced the number of flipped triangles, while only the proposed method eliminated topology violations completely; (2) the percentage of the non-topological area Ar n was about 8~20% of the total area with all but the proposed smoothing method; (3) The proposed method decreased R 2 , mainly due to the constrain of the topological condition. The number of flipped triangles for individual observers by all methods can be found in the supplementary data (see Table A in S5 Text). The results also show that the proposed method was always better in preserving topology while the goodness of fit to the fMRI time series was reduced from the original result.
For a more intuitive comparison, Fig 9 shows the original retinotopic map of the left hemisphere of the first observer in the HCP dataset and the smoothed versions by different methods. The original retinotopic map (Fig 9A) is not topological because there are many flipped triangles. Fig 9B-9E show the results from the Average, Median, Laplacian, and proposed smoothing methods, respectively. Although all smoothing methods improved the smoothness of the maps, topological violation still occurred in the results of all but the proposed smoothing methods, especially near the fovea. Only our proposed method can generate topology-preserving results (Fig 9E). Fig 9F shows  mesh for each smoothing method in Fig 9H-9K. The green curves from the left to right represent 0.5˚to 6˚eccentricities with a 0.5˚interval, and the blue curves, from the bottom to up, represent -90˚to 90˚polar angles with a 15˚interval. We can also clearly see topological violations in Fig 9G-9J, especially in the fovea region. Only the proposed method can generate topological and smooth results (Fig 9K). The same experiment was conducted on the data of all observers and the results are available in S5 Text. These results clearly show that the proposed method can generate smooth and topological retinotopic maps in V1.

Retinotopic maps of the V1/V2/V3 complex
We took the initial visual area labels from the multimodal surface matching (MSM) results of the HCP group [57,58], removed phase-jumping and phase-reversal for the left and right hemispheres in multiple visual areas after polar angle transformations, and achieved topological and smooth results across them (Fig 10). Fig 10A-10E show the visual coordinates in the eccentricity-extended-polar coordinate space, in which the x axis is the eccentricity, and the y axis is the extended polar angle (Table 1), for the raw data, average smoothed data, median smoothed data, Laplacian smoothed data, and topological smoothed data, respectively. All smoothing methods reduced topology violations, compared to the raw data. However, only our proposed smoothing method generated topological mapping such that the contour curves (the green and blue curves in Fig 10K) were smooth and had no self-intersection. We performed the same comparison across all observers (Table A in S5 Text) and showed that our proposed topological smoothing was the best in topology preserving.

Boundary delineation
In the eccentricity-extended-polar coordinate space, the boundaries between different visual areas in the V1/V2/V3 complex are defined by specific polar angles. Fig 11A shows the boundaries directly inferred from the raw data of the first observer. Fig 11B-11E show the boundaries generated from the smoothed data of the same observer, after average, median, Laplacian, and topological smoothing. The new boundaries largely coincided with manually labeled visual areas from the average data, represented by the colors of the cortical surface Fig 11). The green curves (from the left to the right) are eccentricity contours from 1˚to 6˚with a 1˚interval, and the purple curves represent boundaries of visual areas (V3d/V2d/V1d/V1v/V2v/V3v, from up to down). As in the retinotopic maps, topological violations occurred in boundary delineation from the raw data and after smoothing with all but our proposed method. We present all observers' area boundaries after topological smoothing in S5 Text. They were all smooth and continuous.

Discussion
In this study, we modeled the topological condition and generated topological and smooth retinotopic maps. The Beltrami coefficient, a metric of quasiconformal mapping, was used to define the topological condition. We developed a mathematical model to quantify topological smoothing as a constrained optimization problem and elaborated an efficient numerical method to solve it. The method was applied to V1, V2, and V3 simultaneously, and was robust to inaccurate boundary labeling. Experiments with both simulated and real retinotopy data demonstrated that the proposed method could generate topological and smooth retinotopic maps. As a result, we were able to improve boundary delineation. To our best knowledge, conventional methods have not fully considered topological constraints for multiple regions in retinotopic smoothing. Our novel algorithm made the retinotopic maps from BOLD fMRI topological and consistent with results from neurophysiology. Our work improved the quality of retinotopic mapping and built a solid foundation for future retinotopy related research.

Quality assessment
In this study, we enforced the topological condition in retinotopic mapping. This enforcement could potentially introduce large visual coordinate changes to the raw data compared to other smoothing methods. It is important to monitor whether the smoothing results are acceptable or have relatively big deviations. Such quality assessment can provide the right level of confidence when drawing conclusions based on the smoothing results. In Fig 12, the visual coordinate changes of two typical observers are shown on the cortical and parametric surfaces respectively. Fig 12A shows the studied visual areas on the 3D cortical surface. Fig 12B and 12C shows the  (1) , and the y axis is the extended polar angle v (2) , in (H)-(K), the green and blue curves are levels sets, i.e., the contours of eccentricity v (1) and extended polar angle v (2) , respectively. distortion in visual space. The colors indicate the amount of visual coordinate changes for the proposed smoothing process. The first observer's results (Fig 12A) suggested that the smoothing was generally good, with on average an about 1.12˚change of visual coordinates. In contrast, the other observer's results (Fig 12C) suggested that some sub-regions have relatively big visual coordinate changes, with on average an about 1.75˚change of visual coordinates. We also show the mean visual coordination distortion distribution of all 181 subjects in Fig 12D. Together with Table 5, the quality assessment results ensure our retinotopic mappings well respects the original fMRI signals.

Out-most boundary condition
In this work, visual areas, V1/V2/V3, and the initial boundaries between them were determined by surface registration. Following the HCP protocol [36], MSM [58] was adopted to register each observer's cortical surface to the retinotopy atlas obtained from the HCP [34]. After the registration, the boundaries for each observer were determined along equal polar angle curves, and our extended polar angles were initialized with these boundaries. The outmost eccentricity boundaries were manually determined based on the eccentricity signal range (in our case, approximately from 0.2˚to 7˚). The initial boundaries might not be accurate since they were determined by surface registration. With extended polar angles, our algorithm was robust enough to fix the interior visual area boundaries. Our method requires that the visual coordinates of the out-most boundaries are specified when solving the LBS. Even so, if the visual coordinates of the out-most boundary are not properly given, we can still update the out-most boundaries and recover interior visual coordinates. We tested two approaches to set the out-most boundary values. One approach is to set a fixed boundary for all observers to simulate the case in which a good boundary configuration is given, and the other approach is to adjust the out-most boundary values by fitting a smooth spline of the interior visual values to simulate the case in which the boundary coordinates are given but not precise. Each approach has its own pros and cons. The first approach is numerically stable but ignores individual differences, while the second approach adjusts the boundary values based on retinotopic data, which has higher accuracy but may fail to set reasonable boundary values. We compared both approaches in the experiments. Fig 13A  and 13C are the results from the first approach, and Fig 13B and 13D are the results from the second approach. We found a relatively large difference on the out-most boundaries (red boundaries in Fig 13A and 13B), while the interior regions were less affected. The mean difference between the results of the two approaches was 0.23˚±0.15˚, which is relatively small. In this paper, we took a combined approach: the boundary values were set within a small tolerance range (0.5˚) of the initial average boundary values. Nevertheless, we found that, although the out-most boundaries may not be precise, the interior results were not much impacted by the boundary values. Therefore, we can use the pre-defined boundary configuration for all the subjects.
It is worth noting that our method is applicable to the registration results from Freesurfer [37] or MSM [58] or Bayesian register [8]. Here we applied it to successfully process data from all 181 HCP observers. In our code repository, we have published all the necessary manuals and atlas for others to reproduce our work. We also provided a simple software tool to manually draw the out-most boundaries on visual cortical surfaces and generate corresponding visual coordinates.

Robustness with inaccurate interior boundaries
When smoothing the retinotopic maps of the V1/V2/V3 complex, we applied linear transformations to the polar angles based on the visual region labels derived from the anatomy. Here, we explain why our proposed smoothing method was robust even with inaccurate boundary specifications.
We used the first observer's data as an example to illustrate that the proposed method was robust to inaccurate interior boundaries (Fig 14). In Fig 14A, we show the raw visual polar angle data on the parametric disk for the first observer, with a white curve at 0.3-degree eccentricity. Then we applied a linear transformation before smoothing, according to Table 1. For a vertex with an inaccurate boundary label, the polar angle would be transformed with a wrong offset. Since the transformation adds or subtracts multiples of π as the offset according to its visual region label in Table 1, a wrong label makes the transformed value dramatically different from its neighboring value. Therefore, an incorrectly transformed vertex looks like a spike and destroys its topological condition (red points in Fig 14B).
In the next step, we applied the topological smoothing (smoothed visual polar angle is the slope-shaped surface in Fig 14B). Since it enforces topology conditions, our topological

PLOS COMPUTATIONAL BIOLOGY
smoothing method replaces the spike with neighboring values. Namely, the proposed method can fix the wrong offset and is therefore robust to inaccurate interior boundaries.
We further verified that the robustness was not due to Laplacian smoothing. In Fig 14C, we plot the raw data and topologically smoothed results of the white curve. The white curve is formed by points with eccentricity v (1) = 0.3˚with tolerance 0.05˚. We now consider the data near the white curve. If the eccentricity measurement is accurate, the 2D topological condition can be considered as a monotonicity requirement. In Fig 14C, the blue dots are the raw extended polar angle data, and the green curve is the topological smoothing result. Even though there were outliers (the red circled points, mainly due to the noisy eccentricity data), the main portion of the trend was monotonic. In Fig 14D, we plot the results with only Laplacian smoothing without topology constraints as the green curve, which is crooked in order to fit the data but not monotonic (topological). This means Laplacian smoothing can smooth but cannot fix topological violations. In contrast, our topological smoother ensures topological condition, which essentially makes the mapping "monotonic". It is the topological constraint rather than Laplacian smoothing that makes the smoothing more robust to inaccurate interior boundaries.
The new topology-projection module in the proposed method makes it robust to inaccurate interior boundaries. The smoothed extended polar angles from topological smoothing can be used to infer visual area labels.
We computed visual coordinate changes following an alteration of the V1/V2 boundary. If the visual coordinates did not change too much when the boundary was altered, the results were robust to interior boundaries. We started from the first subject's data (Fig 15A). We expanded V2v by moving points in V1v into V2v if they were adjacent to V2v. The expansion was repeated to simulate different boundary misplacements (three rounds in Fig 15B and five rounds in Fig 15C). Then we applied our smoothing method and compared the results with no perturbation. The boxplot in Fig 15D shows that the visual coordinate differences were relatively small (the average is less than 0.5˚) with V2v boundary expansion. The results suggest that the proposed method is robust to interior boundary placement.

The trade-off between smoothness and quality of fit
The smoothing weight s in Eq (7) controls the trade-off between smoothness and quality of fit under topological condition. Although other studies have attempted to find the optimal s for linear smoothers [59], the topology constraints made it challenging. We used the following way to set the proper s: find the largest weight while ensure that the average difference between the smoothed visual coordinates and the raw data is within a certain range (1˚in this work). To empirically determine the s value, we computed the percentage of the non-topological area and the Normalized-Root-Mean-Square-Error (NRMSE) of the fit to the fMRI time series for the raw data, when s = 0, and when s was set properly for the first observer. For the raw data: NRMSE = 0.07944×10 −3 , 24.7% of the area is non-topological; when s = 0, NRMSE = 0.08546×10 −3 , 0% of the area is non-topological; and when s = 0.0121, NRMSE = 0.8551×10 −3 , 0% of the area is non-topological. The NRMSE values for s = 0 and s = 0.0121 were quite close. Therefore, we prefer the latter result in this work.

Beltrami coefficient vs. Jacobian
In this work, the topological condition is ensured by enforcing kμk 1 <1. Another geometric concept, the determinant of the Jacobian matrix J can also be adopted to enforce the topological condition. For example, it has been adopted by previous work to monitor the local visual field signs of retinotopic maps [14], which indicate the trend of visual polar angle changes (increasing or decreasing) on an iso-level eccentricity curve. Constraining the determinant of the Jacobian to be positive can also ensure the topological condition and visual field sign within a neighborhood. Mathematically, the Beltrami coefficient and the determinant of the Jacobian matrix are related via the following equation: We prefer the Beltrami coefficient because it is more suitable for retinotopic mapping. First, with the Beltrami coefficient, we can define the topological constraint as a penalty term to smoothing energy (Eq 7) because the Beltrami coefficient quantifies angle distortions and retinotopic mapping is considerably angle preserving [54,61]. In contrast, J cannot be added as a term in smoothing energy because it quantifies area changes, and retinotopic mapping is not area-preserving. Numerically, one must separately treat the constraint and adjust visual coordinates to ensure that J is postive. Second, as indicated by Theorem 1 in S2 Text, the Beltrami coefficient map uniquely determines a diffeomorphic mapping between two 2D surfaces, making it possible to recover/solve the mapping if we know the Beltrami coefficients. Therefore, we can divide the smoothing process into two subproblems and solve each separately and efficiently. To our knowledge, there are no efficient and stable methods to utilize Jacobian J to reconstruct diffeomorphic mappings. In summary, although it may be possible to use J to quantify the topological condition, the Beltrami coefficient is more suitable for modeling retinotopic maps.

Smoothness and topological projection
The proposed topological smoother produces topological and smooth retinotopic maps. Although these two conditions are related, the topological condition does not ensure smoothness, nor does smoothness ensure the topological condition. In fact, the topological condition can be ensured without smoothing. Empirically, no one expects retinotopic maps to be topological but extremely irregular and rough.
Smoothing does to some degree improve the topological condition. Although we do not know the exact relationship between visual coordinate smoothing and the topological condition (monitored by Beltrami coefficient), the Beltrami coefficient quantifies the local deformations based on its definition. Smoothing can reduce the outliers in Beltrami coefficients. Because 80% of the triangles in the HCP dataset are orientation preserving (See Table A in S5 Text), smoothing the visual coordinates may reduce the norm of the Beltrami coefficients in average and improve the topological condition. If the data have more flipped triangles (e.g., 50% triangles' are flipped), smoothing may introduce more non-topological triangles and make the problem worse.
The topological condition cannot be ensured by reiterative smoothing. As shown in Fig 15D, repeated Laplacian smoothing did not generate fully topological results. The topological condition is mainly fulfilled by topological projection. However, after fixing the topological condition, a major concern is the results could be too irregular or deviate too much from the original pRF solutions. The proposed Topological smoother can produce results that are both topological and smooth. Roughly speaking, if the fMRI signal to noise ratio is very high (say SNR >20), the pRF solutions would be smooth and topological. For example, the pRF solution of the averaged fMRI signals from many subjects (https://osf.io/bw9ec/wiki/home/) has very few flipping triangles in V1. This is why the energy function in the topological smoother has a term that penalizes non-smoothness and therefore may improve the fit to the data and improve the topological condition. We believe the approach can balance the topological condition and smoothness.

Extension to higher visual regions
Although we focused on the V1/V2/V3 complex in this study, our method can be in principle extended to higher visual areas. However, there are some additional challenges. As we move to higher visual areas, the signal-to-noise ratio of the retinotopic map is further reduced. It is extremely challenging for all smoothing methods to balance the topological condition and goodness of fit. As shown in Fig 16, we applied our proposed smoothing method on the V3B retinotopic maps of the first four observers in the HCP dataset. Instead of smoothing V3B with V1-V3 together, we performed separate smoothing on V3B and concatenated the results with those on V1-V3. Although topological smoothing was feasible in higher visual regions, extra care must be taken because of the low signal-to-noise ratio of BOLD fMRI activation.

Benefits of the extended polar angles
Our previous work on topology correction [29] was only applicable to V1; it could not be used to improve boundary delineation. In this work, by introducing extended polar angles, we extended the framework to the V1/V2/V3 complex and achieved topology-preserving smoothing in multiple visual areas simultaneously. Although our previous work [29] might be potentially appliable to V2 and V3, accurate V1/V2 and V2/V3 boundaries were necessary. In the current work, given the out-most boundaries of the V1/V2/V3 complex, our work not only generated more precise smoothing within the interior of each region, but also automatically delineated the boundaries between them. The new development on boundary delineation may significantly improve the robustness and accuracy of retinotopic mapping. In addition, our prior work was mainly evaluated with visual coordinate the changes. How much impact it had on the goodness of fit of the fMRI time series was not evaluated. Here we added fits to the fMRI time series in performance evaluation, which may better justify our work and alleviate the concern that the proposed method may over-smooth the retinotopic maps. With extensive evaluation based on both synthetic data and real retinotopic map data from 181 HCP observers, the current study demonstrated that our proposed method may improve the accuracy and stability of retinotopic mapping, especially the interior boundary conditions. Elaboration on the reduced R 2 in Table 5 The topological smoother can completely remove topology violations in retinotopic maps. However, as shown in Table 4, it reduced the measure of explained variance, R 2 . In general, the pRF model should have a better R 2 , since R 2 values are computed at each point on the cortical surface. There are no constraints, nor penalty for non-topological results in R 2 . However, neurophysiology studies strongly suggested that retinotopic mapping is topological for normal subjects. If the pRF results are not topological, the problem is very likely caused by the low spatial resolution and SNR of the fMRI data, an issue that cannot be addressed with current MRI technologies. In this work, we assume that the neurophysiology results are correct and attempted to fix topological violations in retinotopic maps by slightly updating the visual coordinates. We used pRF results as input data without 1) fine tuning pRF model parameters, such as the gain or the exponent of the power function [62], and 2) iterating with the pRF model to search for an optimal solution with balanced fitting quality and topological condition. Our ongoing work [62] tries to integrate the pRF model and topological conditions by defining a combined energy term. We expect such a balanced approach may improve topological condition without sacrificing data fitting quality.

Relation to registration
The Beltrami coefficient (BC) encodes 2D-to-2D mapping up to a normalization [60]. Manipulating the Beltrami coefficients of the registration function between a subject and a template can also achieve diffeomorphic registration and generate topological results [63]. The benefits of this approach was discussed in [8]. However, one main issue is that the template itself might be biased. One potentially promising solution is to combine registration and smoothing: while registration provides good boundaries, topological smoothing provides good interior smoothing [62].

Conclusion and future work
We adopted the Beltrami coefficients to quantify the topological condition in retinotopic mapping, formed a concise but fundamental model, and provided an efficient numerical solver to generate smooth retinotopic maps. To our knowledge, the proposed topological smoother is the first method that guarantees the topological condition in retinotopic mapping. Our work improved the quality of the retinotopic maps and built a solid foundation for future research based on retinotopy.