METHODS AND SYSTEMS FOR PRECISE QUANTIFICATION OF HUMAN SENSORY CORTICAL AREAS

Abstract

A sensory mapping method for a human brain is disclosed. The method includes the steps of flattening the cortical surface of the human brain, projecting functional imaging data onto the flattened surface, smoothing the functional imaging data, generating a sensory map, registering sensory maps across individuals and analyzing the maps in the common space. The flattening utilizes a conformal parametrization method. The smoothing utilizes a topological smoothing method that utilizes a diffeomorphic smoother. The registering is diffeomorphic. The sensory mapping method may further include a step of processing the functional imaging data to produce topological results.

Description

TECHNICAL FIELD

The present disclosure relates to imaging, and in particular to techniques for imaging brain structure and activity.

SUMMARY

Methods and systems for precise quantification of human sensory cortical areas are disclosed. In an exemplary embodiment, a sensory mapping method for a human brain comprises flattening the cortical surface of the human brain, projecting functional imaging data onto the flattened surface, smoothing the functional imaging data or processing the functional imaging data to produce topological results, generating the sensory map, registering sensory maps across individuals, and analyzing the maps in the common space. The sensory mapping method may further comprise rendering the sensory map with an image of stimuli that activate the cortical surface.

The flattening may comprise a conformal parametrization method. The flattening may comprise applying harmonic mapping to a unit disk and refining repetitively parametric coordinate for each point on a parameter domain to produce conformal mapping.

The smoothing may comprise a topological smoothing method that utilizes a diffeomorphic smoother. The smoothing may comprise applying Beltrami coefficient to quantify diffeomorphism and modeling the smoothing in an optimization framework.

The processing may comprise a topology-preserving cortical surface segmentation method. The processing may comprise segmenting the cortical surface of the human brain by preserving the prior connectivity pattern of the cortical surface areas, decoding the functional imaging data based on at least one topological condition within each cortical area, and repeating the segmentation and decoding.

The registering may be diffeomorphic. The registering may comprise a method that aligns a human brain sensory map of an individual to a predefined template.

In another exemplary embodiment, a sensory mapping method for a human brain comprises flattening the cortical surface of the human brain, projecting functional imaging data onto the flattened surface, processing the functional imaging data to produce topological results, and generating a sensory map. The processing may comprise a topology-preserving cortical surface segmentation method.

The contents of this section are intended as a simplified introduction to the disclosure, and are not intended to be utilized to limit the scope of any claim.

BRIEF DESCRIPTION OF THE DRAWINGS

With reference to the following description and accompanying drawings:

FIG. 1 illustrates a conformal parameterization process in accordance with various exemplary embodiments: (a) harmonic map of the region of interest; (b) cortical surface and region of interest; the checker-board texture inferred from the conformal map is applied on the region of interest; and (c) the resulting conformal map of the region of interest;

FIG. 2 illustrates a texture map between visual field image and visual cortex in accordance with various exemplary embodiments;

FIG. 3 illustrates a human visual system in accordance with various exemplary embodiments;

FIG. 4 illustrates several spaces and relations between them in accordance with various exemplary embodiments;

FIG. 5 illustrates (a) approximate a mapping in discrete and (b) a divergence approximation on the vertex ring in accordance with various exemplary embodiments;

FIG. 6 illustrates a function domain, mapping result, and mapping results with noise in accordance with various exemplary embodiments: (a) mapping domain on the visual field; (b) mapping result without noise; (c) mapping result with weak noise (Peak Signal-Noise Ratio (PSNR)=20 dB); and (d) mapping result with strong noise (PSNR=10 dB);

FIG. 7 illustrates smoothing results for map with big noise in accordance with various exemplary embodiments: (a) average smoothing; (b) median smoothing; (c) Laplacian smoothing; and (d) exemplary smoothing;

FIG. 8 illustrates a subject's left retinotopic mapping in accordance with various exemplary embodiments: (a) raw retinotopic mapping; (b)-(e) smoothing results of average, median, Laplacian, and proposed smoothing methods; (f) the raw mapping on cortical surface; (g)-(j) smoothing results on inflated surface of each method, in the same order of (b)-(e); and (k) the left inflated cortex;

FIG. 9 illustrates a retinotopic map experiment and the pipelines of population receptive field (pRF)/topological receptive field (tRF) models in accordance with various exemplary embodiments;

FIG. 10 illustrates a parametrization and segmentation process in accordance with various exemplary embodiments: (a) the visual field; (b) the region of interest (ROI); (c) conformal parametrization of ROI; and (d) the topology-prior (the label is defined for each point);

FIG. 11 illustrates comparison of various methods for the V1-V3 model in accordance with various exemplary embodiments: (a) Ground-truth; (b) pRF; (c) model-fitting; (d) smoothing; (e) registration; and (f) the proposed tRF;

FIG. 12 illustrates results of: (a) pRF; (b) model-fitting; (c) smoothing; (d) registration; and (e) tRF in accordance with various exemplary embodiments;

FIG. 13 illustrates the visual system in accordance with various exemplary embodiments;

FIG. 14 illustrates several spaces and registration in accordance with various exemplary embodiments;

FIG. 15 illustrates (a) approximate the mapping in discrete and (b) the divergence approximation on the vertex ring in accordance with various exemplary embodiments;

FIG. 16 illustrates template and subject in accordance with various exemplary embodiments: (a) template; (b) perfect subject; (c) subject data with a small added noise (Peak Signal-Noise Ratio is 20 dB); (d) Subject data with a big added noise (Peak Signal-Noise Ratio is 10 dB); (e) ground truth displacement (curves are for eccentricity/angle contour: landmarks are marked for (a)-(d));

FIG. 17 illustrates: (a)-(b) Before and (c)-(d) After Registration of the left hemisphere of the subject ‘102816’ in accordance with various exemplary embodiments;

FIG. 18 illustrates the template morphing process in accordance with various exemplary embodiments; and

FIG. 19 illustrates a flow chart of a sensory mapping method for a human brain in accordance with various exemplary embodiments.

DETAILED DESCRIPTION

The following description is of various exemplary embodiments only, and is not intended to limit the scope, applicability or configuration of the present disclosure in any way. Rather, the following description is intended to provide a convenient illustration for implementing various embodiments including the best mode. As will become apparent, various changes may be made in the function and arrangement of the elements described in these embodiments without departing from principles of the present disclosure.

For the sake of brevity, conventional techniques and components for mathematical processes, transforms, mapping, smoothing, and/or the like may not be described in detail herein. Furthermore, the connecting lines shown in various figures contained herein are intended to represent exemplary functional relationships and/or physical couplings between various elements. It should be noted that many alternative or additional functional relationships or physical connections may be present in exemplary methods and systems for imaging and/or components thereof.

Sensory maps of the human brain are of great scientific interest and could have important medical applications. fMRI makes it possible to define the sensory maps for an individual in vivo. However, due to the low signal-noise ratio and low spatial resolution of fMRI, the sensory maps are usually noisy and don't preserve topology.

In particular, in clinical settings, doctors require quantitative scores. The current mapping methods of the sensory areas of the human brain can only generate qualitative maps with high uncertainty. To address these shortcomings of prior approaches, exemplary embodiments provide smoothing tools, topology correction tools, quantification tools, and multi-subject alignment tools to generate a complete quantitative description of the maps. The quantitative scores are desirable, for example in connection with disease diagnosis and prognosis.

The sensory areas of the brain, including the auditory cortex, visual cortex, and so forth, contain topographical representations of the sensory space. Functional magnetic resonance imaging (fMRI) signals have produced subject-specific cortical maps. However, such maps are noisy, incomplete, and even contradict neurophysiological results. In contrast, the methods and systems in this disclosure can reduce the noise, fix topology violations, and quantify the properties of the maps across time and/or individuals. In particular, exemplary systems and methods can precisely smooth, register, and quantify sensory maps of the human brain. The methods and systems can be used to quantify properties of the sensory areas of the human brain that are associated with normal and abnormal development, aging, and diseases related to sensory processing.

An exemplary method for quantification of a sensory map comprises: (1) flattening the cortical surface; (2) projecting the functional imaging data onto the flattened surface; (3) smoothing the functional imaging data or processing the functional imaging data to produce topological results; (4) generating the sensory map; (5) registering sensory maps across individuals; (6) analyzing the maps in the common space. The sensory map may be a retinotopic map.

With reference now to FIGS. 1 through 18, in accordance with various exemplary embodiments, exemplary methods utilize concepts from computational geometry and provide techniques to reduce noise, fix topology violations, and precisely quantify the properties of the sensory maps. In an exemplary embodiment, the disclosed systems and methods may be applied to visual retinotopic maps with excellent results. However, all practical applications of principles of the present disclosure are contemplated herein.

FIGS. 1-2 illustrate methods and systems for flattening a cortical surface. FIGS. 3-8 illustrate methods and systems for smoothing the functional imaging data. FIGS. 9-13 illustrate methods and systems for processing the functional imaging data to produce topological results. FIGS. 13-18 illustrate methods and systems for registering sensory maps.

Exemplary systems and methods may use a conformal mapping method together with a geodesic disk to reduce angle distortion and geodesic distortion. Exemplary systems and methods may correct topology violations and perform topology-preserving smoothing. Exemplary systems and methods may reduce angle distortion as well as geodesic distortion and produce high quality mapping. Exemplary systems and methods may provide utilities to register or align sensory maps across subjects based on both structural and functional MRI data. The exemplary systems and methods disclosed herein may be integrated to allow the quantification of human brain sensory maps.

In an exemplary embodiment, in a sensory mapping method, visual cortical surfaces are conformally mapped to a topological disk on which local geometry structures are well preserved. Then retinotopy data are smoothed to ensure that topology is preserved on the maps. Specifically, the Beltrami coefficient is adopted as the metric of topology, and a topology-preserving smoothing formulation is solved efficiently. It is believed that the disclosed method may be the first method that can smooth the retinotopic maps and fix topology violations. Thereafter, quantitative properties of the maps, e.g., cortical magnification factor, can be measured precisely using differential geometry tools on the smoothed cortical maps.

In addition to use for analyzing the sensory maps of individual subjects, principles of the present disclosure also provide the capability to better align and analyze sensory maps across time and/or subjects. By considering the intrinsic structure in the sensory maps, exemplary methods can achieve better cortical registration across subjects.

Exemplary systems and methods make it possible to precisely quantify sensory maps of individuals and further improve the quality of sensory maps with data from multiple individuals. Although various exemplary embodiments disclosed herein are discussed in the context of visual retinotopic maps, those skilled in the art will appreciate that principles of the present disclosure can also be applied to all the other sensory areas, including the auditory cortex, somatosensory cortex, and olfactory cortex.

Methods and Systems for Flattening a Cortical Surface

As shown in FIG. 1, an exemplary method is disclosed herein for conformal disk parameterization, c:M→D, to flatten the cortical surface patch M to the disk domain D, the check-board region in and the disk in FIG. 1.

To achieve quality conformal mapping results, a new algorithm is introduced to achieve conformal mapping of the V1 retinotopic map on the unit disk D. First, the harmonic mapping is applied to the unit disk as the initial mapping. FIG. 1(a) illustrates harmonic mapping of a region of interest whose location in the cortical surface is labeled in FIG. 1(b). This produces the result that the boundary of M is on the boundary of the D. Second, the conformality of the patch is iteratively optimized until the BC within V1 is within the predefined error |μ|∞<ε, as shown in FIG. 1(c).

Harmonic Maps

The harmonic map h: M→D can be found by minimizing the following energy.

E(h)=∫_M|∇h|²dv_M

The harmonic map for disk-like surfaces h: M→D can be found by setting ∇E(h)=0, which yields the Laplace equation

Δh(u)_M=0

h/∂M=g

where Δ is the Laplacian operator and g: ∂M→∂G can be given by the arc length parameterization.

In the discrete case, the Laplace equation is a sparse linear system, written as L^hh=0. The matrix element

$L\frac{h}{i,j}$

is

$L_{i,j}^h=\{\begin{array}{c}\sum _{\left[u_i,u_j,u_k\right]\in N\left(i\right)}\frac{\left({\mathrm{As}}_j\right)·s_i}{\left[u_i,u_j,u_k\right]},\mathrm{if}i\ne j\\ -\sum _{k\ne i}{L}_{i,k},\mathrm{if}i=j\\ 0,\mathrm{otherwise}.\end{array}$

Here, the boundary vertices are set by their boundary length. Each vertex {v_i}_i=0ⁿ⁻¹ is mapped to the unit circle according to the ratio of the edge lengths along the boundary loop. The harmonic map h is efficiently obtained by solving the linear system.

Optimizing Conformality in the Primary Visual Cortex (V1)

After the harmonic mapping is achieved, the parametric coordinate is iteratively refined for each point on the parameter domain to achieve conformal mapping in the primary visual cortex (V1). Let r be the desired refinement map, such that the result r is conformal mapping to 3D cortical surface in V1. The mapping that this method is seeking r=c∘h⁻¹ is illustrated in FIG. 1 (a)-(c). As described above, if there is the Beltrami coefficient μ_r of r, r can be reconstructed. If c is conformal, the Beltrami coefficient of the composite map is μ_r=μh⁻¹. Since h is given, its inverse h⁻¹ can be computed. To compute the Beltrami coefficients for the inverse harmonic map in the discrete case, i.e., h⁻¹: D→M between the corresponding triangular faces, the first step is to use piece-wise rigid motion transformation R to translate each triangular face of M onto R². Then the next step is to compute the Beltrami coefficient, μ_R∘h⁻¹, of the map R∘h⁻¹. Since rigid motions are conformal, the resulting Beltrami coefficient μ_R∘h⁻¹ equals μ_h−1. With the defined μ_r the next step is to compute a new mapping r and achieve the desired conformality. Numerically, to achieve the best quality, the refinement is repeated until the mapping is conformal, i.e., μ_c falls below a predefined error bound ^εμ_r, ∥μ_r∥∞≤^εμ_r.

Because a goal is to achieve a conformal parameterization of the primary visual cortical surface (V1), μ_r=μ_h−1 is set for the faces in the V1 area. For other parts, μ_r is kept unchanged. The conformal procedure for the primary visual cortex (V1) is summarized in the following Algorithm 1. Finally, an observer-wise Mobius transformation is applied on the disk, a conformal transformation with analytical scaling and rotating, to eliminate conformal mapping ambiguity and roughly align the visual regions between observers.

Once the conformal parameterization is obtained, the next step is V1's retinotopic mapping from the parametric disk to the visual field. More specifically, the 2D map is ƒ_r: c(V1)→D, from the parametric domain of V1 cortex to the visual field. With the bijective conformal parameterization c, the retinotopic mapping can be evaluated through a planar-to-planar function ƒ:D_V1→D such that D_V1=c(V1) and ƒ=ƒ_r∘c⁻¹.

Algorithm 1: Conformal parametrization of the visual cortex.
   Input: Cut patch M = (PM, EM, FM), where PM = {Pi} is the mesh vertex set; vertex indices
   IP and face indices IF for the primary visual cortex respectively, and a tolerance εμr
   Output: Conformal parametric coordinates {ui | i ∈ IP}
1  begin
2  Compute the harmonic parametric h = {hi|Pi ∈ PM};
3  Initialize ui ← hi, δμ ← ∞;
4  Compute the Beltrami coefficient {μi|fi ∈ FM} for the inverse map u−1
5  while δμ > εμr do
6  Change BC to zero in visual cortex, i.e., μ0i = 0, ∀i ∈ IF, and μ0i = μi, ∀i ∈ IF;
7  Fix the boundary points and compute the refined coordinate {ui|Pi ∈ PM}
8  Compute the Beltrami coefficient {μi|fi ∈ FM} for the inverse map u−1
9  Update maximal BC in V1, by δμ ← max|μi|: i ∈ IF
10 end
11 return {ui}
12 end

Illustrative Visualization of the Retinotopic Map

Once a complete retinotopic map is processed by the method described herein, the following method provides an intuitive illustration of the perception centers by texture mapping, as shown in FIG. 2. This texture map becomes meaningful with the smooth and topological retinotopic map. Specifically, the visual coordinate for each visual cortex is first computed, as shown in FIG. 2a. Then the next step is to interpret visual coordinate for point within the triangular face with the utilities of OpenGL. The final step is to render the cortical surface with the visual image shown in FIG. 2b. This method can easily be used to check, invest, and verify discoveries relates to the retinotopic map.

Methods and Systems for Smoothing the Functional Imaging Data

SUMMARY

Retinotopic mapping, the mapping of visual input on the retina to cortical neurons, is an important topic in vision science. Typically, cortical neurons are related to visual input on the retina using functional magnetic resonance imaging (fMRI) of cortical responses to slowly moving visual stimuli on the retina. Although it is known from neurophysiology studies that retinotopic mapping is locally diffeomorphic (i.e., smooth, differentiable, and invertible) within each local area, the retinotopic maps from fMRI are not diffeomorphic, especially near the fovea, because of the low signal-noise ratio of fMRI. The present disclosure discloses a mathematical model that produces diffeomorphic retinotopic mapping from fMRI data. In an exemplary embodiment, the present disclosure adopts a geometry concept, the Beltrami coefficient, as the tool to define diffeomorphism, and models the problem in an optimization framework. Efficient numerical methods are applied to produce the model. Experimental results with both synthetic and real retinotopy data demonstrate that the method disclosed herein is superior to conventional smoothing methods.

Introduction

Almost half of the human cerebral cortex is dedicated to visual processing. Identifying and analyzing visual areas of the human brain is an important topic in vision science and neurobiology. Retinotopic mapping via functional magnetic resonance imaging (fMRI) provides a noninvasive way of defining the location and size of receptive filed of populations of visual neurons on the retina. Typically, cortical neurons are related to visual input on the retina using fMRI of cortical responses to slowly moving visual stimuli on the retina. The visual stimulus is carefully designed so that the stimulus pattern is unique with respect to each retinal location. By recording and analyzing the fMRI signals at each cortical location, one can decode the retinal visual input coordinate that generates such fMRI signals. Previous studies have revealed that much of the visual cortex is organized into visual regions (e.g., V1, V2) with a retinotopic map in each region. It is known from neurophysiology studies that retinotopic mapping is topology preserving (i.e., local neighborhood relationships are maintained) and diffeomorphic (i.e., smooth, differentiable, and invertible) within each visual area. However, the decoded visual coordinates from fMRI retinotopic mapping studies are not diffeomorphic or topological preserving. One reason is that the signal-noise ratio and spatial resolution of fMRI signals are low. Another reason is that the decoding process largely deciphers the fMRI signal point by point without considering the local smoothness, although Gaussian spatial smoother has been applied to fMRI signal processing.

Spatial smoothing was previously conducted on retinotopic maps to perform automatic delineation of visual regions. However, the smoothing is applied to the visual coordinates one dimension at a time, e.g., the x and y visual coordinates separately. There is no guarantee of diffeomorphism. Another approach is registering the noisy retinotopic maps to a diffeomorphic template and assigning the template retinotopic coordinates after registration. It relies on diffeomorphic registration of two surfaces, which is difficult to achieve. Even if the registration issue can be solved, how to define the right template will remain a big problem. Therefore, there is a need to solve the smoothing issue.

In an exemplary embodiment, a method is disclosed that applies the Beltrami coefficient to quantify the diffeomorphism as well as angle deviation. It is known that retinotopic mappings preserve angles to a considerable extent. The method further models the smoothing process as an optimization procedure with diffeomorphic constraints and elaborates the numerical steps to solve the optimization problem. The algorithm was tested on a synthetic dataset and a real dataset. The synthetic data is generated using an ideal retinotopic mapping model (ground truth) with added noise. The application of the method to the synthetic data generates diffeomorphic results without any large deviations, while other methods violate the condition. The real dataset is the Human Connectome Project 7T retinotopy (HCP) data (See Noah C. Benson and Others, “The Human Connectome Project 7 Tesla retinotopy dataset: Description and population receptive field analysis,” Journal of Vision, vol. 18, no. 13, pp. 1-22, December 2018.). Diffeomorphic retinotopic mapping was also achieved with HCP data. In addition, the Beltrami coefficients were obtained that quantify angle distortions for each subject. The method produces diffeomorphic results and minimizes the angle distortion, two important conditions for a valid retinotopic mapping research.

Problem Restatement

A problem involved here is restated in a biological intuitive and mathematically rigorous way. Assume there is a light spot on point v=(v⁽¹⁾, v⁽²⁾)∈² in the visual field, as shown in FIG. 3, E an illustration of visual system to the visual cortex. After the visual system's perception, propagation, and dispensation, the light spot will eventually active a population of neurons. The main purpose of retinotopic mapping is to find the center v and the extent σ∈⁺ of its receptive field on the retina, for each point P=(X,Y,Z)∈R³ on visual cortex.

fMRI provides a noninvasive way to determine v and σ for P, based on the following procedure. (1) Design a stimulus time sequence s(t; v), such that the stimulus sequence is unique for every visual coordinate, i.e., s(t; v₁)≈s(t; v₂), ∀v₁≈v₂; (2) Present the stimulus sequence to an individual and record the fMRI signals from the visual cortex; (3) For each fMRI time sequence on a cortical location, y(t; P), determine the corresponding receptive field, including its center location v and its size σ on the retina, that most-likely generated the fMRI signals. Specifically, given the neurons' spatial response r(v′; v, σ) (a predefined model depicts neural response around v) and the hemodynamic function h(t) (a model of the time course of neural activation to a stimulus), the predicted fMRI signal can be written as:

ŷ(v,σ)=β(∫r(v′;v,σ)s(t,v′)dv′)*h(t)  (1)

where, β is a coefficient that converts the units of response to the unit of fMRI activation. Then the parameters v and σ are given by minimizing the prediction error,

$\begin{array}{cc}\left(v,\sigma \right)=\underset{\left(v,\sigma \right)}{\arg \min }{\hat{y}\left(v,\sigma \right)-y\left(P\right)}^2.& \left(2\right)\end{array}$

The retinotopic mapping of the entire visual cortex is obtained when (v, σ) is solved for every point on the cortical surface. It is known that nearby visual cortex will share nearby visual coordinates, i.e., topology preserving. The notion retinotopic mapping preserves topology means that, if ƒ is the retinotopic mapping, ƒ:Pv, then a neighboring point P′ of P, will have visual coordinate v′ close to v. That is, ƒ is diffeomorphic in the domain of cortical surface.

Due to the low signal-noise ratio of fMRI, ƒ is not diffeomorphic. In an exemplary embodiment, a method is disclosed herein that achieves a diffeomorphic smoothing on ƒ, such that the new mapping ƒ:P{circumflex over (v)} is diffeomorphic and {circumflex over (v)} is close to v. The Beltrami coefficient is adopted to quantify the diffeomorphism,

$\begin{array}{cc}{\mu }_f=\frac{\left(E-G+2iF\right)}{\left(F+G+2\sqrt{EG-F^2}\right)},\mathrm{where}g(\hat{v}=\left(\begin{array}{cc}E\left(\hat{v}\right)& F\left(\hat{v}\right)\\ F\left(\hat{v}\right)& G\left(\hat{v}\right)\end{array}\right)& \left(3\right)\end{array}$

is the metric tensor. If |μ|<1 then ƒ is diffeomorphic. Furthermore, the bigger μ_ƒis, the larger angle deviation is. Therefore, the diffeomorphism defined in this way also provides a measure of angle deviation. The aim is to find {circumflex over (v)}={circumflex over (ƒ)}(P), such that,

$\begin{array}{cc}\hat{v}=\underset{\hat{v}}{\mathrm{argmin}}\int {v-\hat{v}}^2\mathrm{ds},s.t.{\mu }_f<1& \left(4\right)\end{array}$

Model and Method

Model in Discrete

In practice, cortical surface S is usually not an analytical function but consists of a bunch of vertices V_S and triangular faces F_S, i.e. S_S=(F_S, V_S). Similarly, visual space can also be discretized. Instead of taking another discretization, it is beneficial to take the same discretization as S. Then, the remaining problem is to find a visual coordinate v_i for every vertex P_i∈V_S.

In an exemplary embodiment, the problems are analyzed and solved on 2D instead of the original 3D cortical surface. Since diffeomorphism is transitive (If ƒ and c are both diffeomorphic function, then the composite function ƒ∘c is also diffeomorphic), a method in the present disclosure computes a conformal mapping (diffeomorphic and angle preserving) from the occipital region of the cortex to a parametric disk to simplify the problem to 2D. As shown in FIG. 4, an illustration of several spaces and relations between them, the gray colored region is mapped to the unit disk by c:Pu,u=(u⁽¹⁾, u⁽²⁾)∈².

Although conformal mapping is angle preserving, it introduces distance distortion. To reduce the distortion, the method cuts a geodesic disk patch: picks a point on the cortical surface that corresponds to the fovea as the center point, then calculates the geodesic distances from all cortical position to the center point and only keep the portion of cortical surface whose geodesic distance to the center point is within a certain value. If ƒ_c=ƒ∘c⁻¹:u∈²v∈² is solved with constraint |μ_ƒc|<1, then diffeomorphic retinotopic mapping is obtained by ƒ=ƒ_c ∘c. Accordingly, Eq. (4) is simplified to,

$\begin{array}{cc}\widehat{f_c}=\underset{f_c}{\mathrm{argmin}}\int {f_c-\widehat{f_c}}^2\mathrm{du},s.t.{\mu }_f<1& \left(5\right)\end{array}$

To solve Eq. (5) in discrete, the energy is formed as,

E=Σ _i |F _i −{circumflex over (F)} _i|²+λ₁|μ_i|²+μ₂|∇{circumflex over (F)}_i|²  (6)

where F_i=ƒ_c(μ_i) is the raw result from Eq. (2), {circumflex over (F)}_i is to be solved, μ_i=μ_ƒc(μ_i) is the Beltrami coefficient at μ_i, λ₁ and λ₂ are constants. Compared to Eq. (5), gradient of {circumflex over (F)},∇{circumflex over (F)}_i, is introduced to increase smoothness, and |μ_i|² is used to reduce non-diffeomorphism as well as angle distortion.

Numerical Method

It is difficult to minimize the energy in Eq. (6) directly, as it mixes the Beltrami coefficient and the mapping function together. A solution is to alternatively solve with respect to function and Beltrami coefficient in separate steps. Tools are utilized to convert the mapping function and Beltrami coefficient back and forth.

Beltrami Coefficient

Since there is a parametric coordinate u for the cortical surface S, the definition, Eq. (3), can be simplified to,

$\begin{array}{cc}{\mu }_f=\left(\frac{\partial f}{\partial u^{\left(1\right)}}+i\frac{\partial f}{\partial u^{\left(2\right)}}\right)/\left(\frac{\partial f}{\partial u^{\left(1\right)}}-i\frac{\partial f}{\partial u^{\left(2\right)}}\right).& \left(7\right)\end{array}$

Given the explicit form of function ƒ(u), μ_ƒis computed according to Eq. (7). In the discrete case, ƒ_c is interpreted linearly on each triangle. As shown in FIG. 5(a), for u within the triangle Δu_iu_ju_k, ƒ_c(u)=α_iv_i+α_jv_j+α_kv_k, where v_i=ƒ_c(u_i), v_j=ƒ_c(u_j), and v_k=f_c(u_k). FIG. 5(b) illustrates the divergence approximation on the vertex ring. The coefficients α_i, α_j, and α_k are called barycentric coefficients. Specifically, α_i (similarly for α_j and α_k) is the area portion of triangle Δuu_ju_k to Δu_iu_ju_k, i.e. α_i=Area(Δuu_ju_k)/Area(Δu_iu_ju_k). Now the Beltrami coefficient μ_ƒc is computed according to Eq. (7). It is easy to see that μ_ƒc is a constant for each triangle.

Linear Beltrami Solver (LBS)

LBS is introduced to recovery function ƒ=(ƒ⁽¹⁾), ƒ⁽²⁾) for the given Beltrami coefficient μ=ρ+iτ. According to the definition in Eq. (7), the following is the result,

$\begin{array}{cc}\left(\frac{\partial f}{\partial u^{\left(1\right)}}+i\frac{\partial f}{\partial u^{\left(2\right)}}\right)/\left(\frac{\partial f}{\partial u^{\left(1\right)}}-i\frac{\partial f}{\partial u^{\left(2\right)}}\right)=\rho +i\tau .& \left(8\right)\end{array}$

After re-organizing Eq. (8) and eliminating ƒ⁽²⁾, the following is derived,

$\begin{array}{cc}\nabla ·A\nabla f^{\left(1\right)}=0,\mathrm{where}A=\left(\begin{array}{cc}{\alpha }_1& {\alpha }_2\\ {\alpha }_2& {\alpha }_3\end{array}\right),{\alpha }_1=\frac{{\left(\rho -1\right)}^2+{\tau }^2}{1-{\rho }^2-{\tau }^2},{\alpha }_2=\frac{-2\tau }{1-{\rho }^2-{\tau }^2}\mathrm{and}{\alpha }_3=\frac{1+2\rho +{\rho }^2+{\tau }^2}{1-{\rho }^2-{\tau }^2},\nabla f^{\left(1\right)}=\left(\partial f^{\left(1\right)}/\partial u^{\left(1\right)}+\partial f^{\left(1\right)}/\partial u^{\left(2\right)}\right),& \left(9\right)\end{array}$

and the divergence ∇·on vector G=A∇ƒ⁽¹⁾=(G⁽¹⁾, G⁽²⁾) is defined as ∇·G=∂G⁽¹⁾/∂u⁽¹⁾+∂G⁽²⁾/∂u⁽²⁾. By solving the partial equation Eq. (9) with certain boundary conditions, f⁽¹⁾ is solved. Similarly, f⁽²⁾ is solved.

In the discrete case, the function is linearly interpreted on each triangle. The gradient of ƒ, ∇ƒ(u), can be written out. For the divergence on vector G, ∇·G, it is approximated on the dual cell of each vertex (a cell consisted of circumcenters of neighboring triangles). As shown in FIG. 5(b), consider the vertex u_i with its neighbors N(u_i), the divergence on vector G is approximated by,

$\begin{array}{cc}\nabla ·G=\frac{1}{D}{\int }_{\partial D}\mathrm{Gds}=\frac{1}{D}\sum _{T_j\in N\left(u_i\right)}{G}_{T_j}·\left(u_k-u_j\right)& \left(10\right)\end{array}$

FIG. 5(a) is an illustration of approximate the mapping in discrete. Applying the approximations, the result is a linear equation for each f_i together with its neighboring vertices, according to Eq. (9). Finally, f is solved efficiently by collecting all the equations and write them in a matrix form.

Laplacian Smoothing and Chopping

Now, the function and Beltrami coefficient are converted. The purpose is to make a smooth result. The next step is to apply the smooth operation on the function and Beltrami coefficient by the Laplacian smoothing. For instance, to get a smooth Beltrami coefficient v, the following is solved,

$\begin{array}{cc}v=\arg \underset{v}{\min }\int {\nabla v}^2+\lambda {v-\mu }^2\mathrm{ds}& \left(11\right)\end{array}$

where λ is a constant. Eq. (11) can be efficiently solved by its Euler-Lagrange equation (−∇·∇+2λI)υ=2λμ.

Although the non-diffeomorphism is penalized in Eq. (6), if there is a consistent point whose |μ|>1, it breaks the diffeomorphic result. To avoid this situation, one can shrink υ′=υ/|υ|, if|υ|>1. This process is called Chop.

Algorithm

An overall smoothing process is summarized in the following Algorithm 2

Algorithm 2: Diffeomorphic Smoother
   Input: Retinotopic coordinates {vi} of each vertex
   Result: Smoothed retinotopic coordinates {circumflex over (v)} = {{circumflex over (v)}i}
1  begin
2  | {circumflex over (v)}i ← vi, for each vortex
3  | repeat
4  | | Compute μ for {circumflex over (v)}, by Eq. (7):
5  | | Get ν by Laplacian smoothing μ, by Eq. (11);
6  | | Chop non-diffeomorphic according, to
   | | ν′ = ν/ |ν|;
7  | | Compute new mapping {circumflex over (v)} by LBS on ν′;
8  | until max |μ| < 1:
9  | return {circumflex over (v)}i'
10 end

Results

Synthetic Data

Performance of the algorithm may be evaluated on synthetic data. The ground truth mapping is given by,

u ⁽¹⁾ +iu ⁽²⁾=log(v⁽¹⁾+iv⁽²⁾)  (12)

Then log function is chosen as it is a good approximation of the retinotopic mapping from visual space to the flattened cortical surface. As shown in FIG. 6(a), the visual domain is mapped to the parametric space FIG. 6(b). Also, the results show the mapping with a small amount of noise in FIG. 6(c) and a large amount of noise in FIG. 6(d). Several smoothing methods are applied on the synthetic mapping data in both small and big noise cases, including Average Smoothing in FIG. 7(a), Median Smoothing in FIG. 7(b), and Laplacian Smoothing in FIG. 7(c) and their results were compared with the method in the present disclosure in FIG. 7(d). The results for the big noisy case are shown in FIG. 7. All of them are spatial smoothing procedures. As shown in FIG. 7(a), Average Smoothing takes the average around each neighborhood. As shown in FIG. 7(b), Median Smoothing takes the median value of each neighborhood. As shown in FIG. 7(c), Laplacian smoothing is the so-called “Perfect Smoothing”. Except for the method of the present disclosure as shown in FIG. 7(d), all other smoothing procedures process the coordinates u⁽¹⁾ and u⁽²⁾ separately. The following Table 1 lists the performance metric of all the methods.

TABLE 1
Comparing different smoothing methods by four metrics, mean value
deviation of function, the max of the function deviation, and
the number of overlapped triangles (non-diffeomorphic triangles).
Value deviation means error to the ground truth. For cells with
two values like “A/B”, “A” is the result for the
map with a small added noise (PSNR = 20 dB); while “B”
is the result for the map with a big added noise (PSNR = 10 dB).
	Value Deviation	# Flipped
	Method	Mean	Max	Triangles
	Before Smoothing	0.25/0.40	0.64/1.05	19/56
	Average Smoothing	0.18/0.24	0.62/0.77	3/8
	Median Smoothing	0.22/0.28	0.69/0.83	2/24
	Laplacian Smoothing	0.23/0.36	0.66/0.92	11/43
	The Present Smoother	0.19/0.22	0.65/0.75	0/0

It was found that (1) all the smoothing methods helped reduce the noise, and reduce flipped triangles (i.e. triangle violates the topology preserving condition); (2) with a small amount of added noise, the topology violation problem is not severe, average smoothing can also achieve almost diffeomorphic results; (3) however, with a large amount of added noise, only exemplary methods in the present disclosure can achieve a diffeomorphic result.

HCP Data

With 7-Tesla MRI systems, the signal-to-noise ratio of retinotopic mapping has been dramatically improved. Still, the mapping is noisy and non-diffeomorphic. A smoother in the present disclosure may be applied to data from a 7T MRI system. The details of the dataset are known (See Noah C. Benson and Others, “The Human Connectome Project 7 Tesla retinotopy dataset: Description and population receptive field analysis,” Journal of Vision, vol. 18, no. 13, pp. 1-22, December 2018.). Because smoothing is a primary interest here, the evaluation process used publicly available data of the decoded retinotopic coordinate for each cortical vertex of five subjects, and evaluated the performance of an exemplary smoother. The total number of flipped triangles for raw data is 1313, for the average, median, and Laplacian methods are 1289, 1312, and 1337, respectively. The smoother has no flipped triangles (i.e., zero).

For a persuasive and intuitive comparison, FIG. 8 shows the retinotopic mapping of the left hemisphere of the subject, and the smoothed results of different methods. As shown in FIGS. 8(a) and (f), the original mapping is not diffeomorphic. FIGS. 8(b)-(e) show smoothing results of average, median, Laplacian, and smoother methods respectively while FIG. 8 (g)-(j) show smoothing results on inflated surface of each of the methods in FIGS. 8(b)-(e) respectively. FIG. 8(k) is an illustration of the left inflated cortex. Although all smoothing methods improved the smoothness, only an exemplary smoother herein can generate diffeomorphic result, especially near the fovea.

Conclusion

The smoother of the present disclosure adopts the Beltrami coefficient to quantify diffeomorphism of retinotopic mapping and forms a mathematical model to generate diffeomorphic maps. The smoother produces results of diffeomorphism and minimizes angle distortion.

Methods and Systems for Processing Functional Imaging Data

In an exemplary embodiment, a method is disclosed that applies a topological receptive field model (tRF) for human retinotopic mapping, as shown in FIGS. 9-12.

SUMMARY

The mapping between visual inputs on the retina and neuronal activations in the visual cortex, i.e., retinotopic map, is an essential topic in vision science and neuroscience. Human retinotopic maps can be revealed by analyzing the functional magnetic resonance imaging (fMRI) signal responses to designed visual stimuli in vivo. Neurophysiology studies indicated that visual areas are topological (i.e., nearby neurons have receptive fields at nearby locations in the image). However, conventional fMRI-based analyses generate nontopological results because they process fMRI signals on a voxel-wise basis, without considering the neighbor relations on the surface. The present disclosure applies a topological receptive field (tRF) model which imposes the topological condition when decoding retinotopic fMRI signals. More specifically, the method parametrizes the cortical surface to a unit disk, characterizes the topological condition by tRF, and employs an efficient scheme to solve the tRF model. The framework was tested on both synthetic and human fMRI data. Experimental results showed that the tRF model removes the topological violations, improves model explaining power, and generates biologically plausible retinotopic maps. This framework is general and can be applied to other sensory maps.

Introduction

It is of great interest to quantify, simulate, and understand the relation between the visual inputs and the neuronal response with the retinotopic maps. A precise retinotopic map provides opportunities to understand or even simulate various aspects of the visual system. For instance, the complex-log model, which depicted a rough position map between the retina and the primal visual cortex (V1), explained the dynamics of spiral visual illusions. The retinotopic maps discovered some stable visual regions. Additionally, retinotopic map research holds great promise in understanding brain plasticity and improving rehabilitation's efficacy from visual impairments. Further, retinotopic maps have been clinically adopted to monitor the progress and recovery under amblyopia treatment, a disorder that affects about 2% of children and may cause significant visual impairment if untreated.

Typically, human retinotopic maps are obtained in vivo by analyzing cortical functional magnetic resonance imaging (fMRI) response signals to visual stimuli. Since the first fMRI work on retinotopic mapping, several analysis models were proposed to decode perception parameters from the noisy fMRI signals. Among them, the population receptive field (pRF) model generates state-of-the-art results and becomes a cornerstone in fMRI signal analysis of retinotopic maps.

Although the pRF model achieved great success, the pRF results are not topological. The topological condition is an essential requirement of retinotopic maps since neurophysiology studies have revealed nearby neurons have receptive fields at nearby locations in the image (the topological condition). The topological condition is also the requirement of the vision system's hierarchical organization: each visual area represents a unique map of a portion of retina. If there are duplicated representations in a visual area, this visual area should be further divided into more areas. Besides, it is challenging to infer accurate visual-related quantification without a topological retinotopic map, e.g., visual boundaries, cortical magnification factor, visual anisomery. Therefore, topological retinotopic map results are vital for neuroscience research.

A variety of methods have been developed to reduce topological violations in the postprocessing of pRF results. For instance, the model fitting is widely used to fit the decoded parameters with an algebraic model for V1-V3. Smoothing methods, e.g., Laplacian smoothing, process visual parameters one by one but cannot ensure the topological condition. Surface registration of retinotopic mesh is also developed for post-processing of pRF results. However, none of the existing methods have considered the topological condition when decoding fMRI signals. It is advantageous but challenging to impose the topological condition when decoding fMRI. First, the topological condition is different across visual areas, while the precise visual areas are delineated upon the decoding results. One may think the segmentation is available through the anatomical surface. However, all anatomical-based segmentation is not accurate, limited by its modality. Second, because the fMRI signal is noisy, it is challenging to segment visual areas from the noisy decoded results.

Various exemplary embodiments apply a topological receptive field (tRF) framework, which imposes the topological conditions by integrating the topology-preserving segmentation and topological fMRI decoding into a coherent system. The framework first segments the visual areas by preserving the prior connectivity pattern of visual areas based on initial decoding and then models the fMRI decoding with the topological condition within each visual area. Since the new fMRI decoding results can refine the segmentation, the framework repeats the segmentation and decoding until theoretically guaranteed convergence.

This exemplary framework has been validated on both synthetic data and human retinotopy data. The experiments showed the superiority of tRF over other methods, including the population receptive field (pRF) model, post-processing by model fitting, smoothing, or registration methods. In summary, (1) it is the first method that enforces the topological condition in decoding retinotopic fMRI signals; (2) it is also the first automatic visual area segmentation with the topology organization preserved. The framework is general and can be extended to other sensory maps because most of the human perception maps are spatially organized, e.g., auditory map is spatially related to sound frequency.

Retinotopic Mapping Experiment

Here is a brief introduction of an exemplary retinotopic mapping experiment. Such an experiment acquires both structural MM and fMRI signals, as illustrated in FIG. 9. The structural MRI (FIG. 9g) is used to reconstruct the cortical surface. The fMRI (FIG. 9b) is acquired to monitor neuron activities during the visual stimuli (FIG. 9a). The visual stimulation is carefully designed to encode the visual space (FIG. 9f) uniquely. Before the analysis, the fMRI volumes are preprocessed (FIG. 9c) and projected onto the cortical surface (FIG. 9h). Eventually, for each subject, there is a high-quality anatomical/structural cortical surface in the form of discrete manifold S=(V, E, F) (where V is the surface vertex set, E is the edge set, and F is the face set), together with fMRI time series y_i(t) for each vertex V_i∈V on the cortical surface.

pRF Model

The population receptive field (pRF) model decodes fMRI signals vertex by vertex. Namely, on the measurement y_i(t), the pRF predicts perception parameters, including the population perception center v_i=(v_i⁽¹⁾, v_i⁽²⁾) and population receptive field size σ_i. The model takes the polar coordinates for the visual fields, as shown in FIG. 9f. Namely, v⁽¹⁾ and v⁽²⁾ are the eccentricity and polar angle, respectively. Given the receptive model r(v′; v_i, σ), and the hemodynamic model h(t), the predicted fMRI is, ŷ(t; v_i, σ)=β(∫r(v′; v_i, σ)s(t, v′)dv′)*h(t), where β is the activation level (a time-independent scalar). The perception center v_i and population receptive field size σ_i are estimated by minimizing the error between the measured and predicted fMRI activation signal,

$\begin{array}{cc}\left(v_i,{\sigma }_i\right)=\arg \underset{\left(v_i,{\sigma }_i\right)}{\min }{E}_p=\int {\left(\hat{y}\left(t\text{;}{v}_i,{\sigma }_i\right)-y_i\left(t\right)\right)}^2\mathrm{dt}.& \left(13\right)\end{array}$

Solving Eq. 13 for all points generates the pRF retinotopic map. This model shows a typical pRF retinotopic map on the visual cortex in FIG. 9d. Still, a large portion of the result is not topological (the visual coordinates' contour curves are messed up in FIG. 9e). An exemplary framework applies the tRF (FIG. 9i & j) to improve the retinotopic maps.

Parametrization

First, a parametrization is established between the 3D visual cortex and the 2D planar disk. A geodesic patch containing the visual areas is cut. In specific, the framework picks a point p₀ on the cortex, then encloses a patch P comprising cortical points within a certain geodesic to P₀, as illustrated in FIG. 10b. The next step is to compute the conformal parametrization for patch P to planar disk D. If c:P→D is the conformal parameterization, then u_i=c(V_i) where u_i∈² is the parametric coordinate for V_i. The conformal mapping may be done via spherical conformal mapping of the double coverings of surface patches, followed by stereographic projection. With the parametrization c, the topology condition within the planar domains can be discussed. FIGS. 10(a), (c) and (d) illustrates the visual field, the conformal parametrization of the region of interest, and the topology-prior respectively.

Topology Conditions

The topology conditions have two inferences. First, the visual coordinates are topological with respect to the surface parametrization within each visual area (this condition is called the topological condition within each visual area). Second, because the human visual cortex is organized into several visual areas, and the visual areas hold the same organization (for instance, V1 is adjacent to V2d and V2v) across subjects (to distinguish, this condition is called topology-preserving condition across visual areas).

Both conditions can be quantified by the Beltrami coefficient of quasiconformal maps. Namely, the framework disclosed herein uses the Beltrami coefficient to enforce a topology-preserving surface segmentation and uses the concept again to ensure the topological condition for the retinotopic map within each visual area. The following will introduce the Beltrami coefficient's definition and then explain how it is adopted into these two conditions.

An orientation-preserving homeomorphism φ(z):→ is quasiconformal if it satisfies the Beltrami equation

$\frac{\partial \phi \left(z\right)}{\partial z}={\mu }_{\phi }\frac{\partial \phi \left(z\right)}{\partial z}$

where μ∂ is the Beltrami coefficient satisfying |μ_φ|_∞<1. In particular, if μ₁₀₀=0 then φ is conformal.

The topological condition requires the retinotopic map from each visual area to the retina to be orientation consistent. The orientation (biologically the visual field sign, VFS) for the visual area can be either positive or negative, depending on the parametrization and visual area. Once fixing a parametrization, the map (from parametric space to retina) orientation for a specific visual area is then fixed. Of course, one can treat the negative visual field area to be positive by flipping the parametrization vertically (u′⁽²⁾=−u⁽²⁾). Without losing generality, considering a positive visual area, it is required that the VFS is consistently positive. Mathematically, it requires the Beltrami coefficient μ_ƒ, associated with the retinotopic map from parametric space (within the unit disk) to visual field space ƒ:u_iv_i, satisfies |μ_ƒ<1.

A registration-akin method may be used for topology-preserving segmentation. The method diffeomorphically morphs the topology-prior (which has the right topology) to the subject's parametrization disk and then uses prior's label information for the subject. The segmentation is topology-preserving if the morphing is diffeomorphic. Let g: D→D be the morph function. If the Beltrami coefficient μ_g satisfies |μ_g|_∞<1, g is diffeomorphic. Compared to retinotopic registration, where visual coordinates are required when defining a template, the prior label is required in this exemplary framework, which is easier to give and reduces the model visual coordinates hypothesis.

tRF Model

Here is an introduction of the tRF model, which includes fMRI decoding T={(v_i, σ_i) V_i∈V} and visual areas segmentation L={L_i∈|V_i∈V}, (one value for an area) altogether.

Topology-Preserving Segmentation

A wrapping function g is searched that maximizes the label alignment between subject and wrapped-prior. However, there are no labels before segmentation. As an alternative, this framework maximizes the alignment between VFS of topology prior and VFS of a subject. Let J=(F_J, E_J, u_J, L_J) be the topology-prior mesh, where F_J is face set, E_J is edge set, u_J={u_j} is vertex set and L_J={L_j∈|u_j∈u_J} is label set, Ŝ_g:D→{−1,1} be the VFS for topology-prior after wrapping, and s:=sign(1 −|μ_ƒ|) be the VFS function of the subject, where

$\mu f=\frac{\partial v}{\partial \stackrel{\_}{u}}/\frac{\partial v}{\partial u}(v-v^{\left(1\right)}+{\mathrm{iv}}^{\left(2\right)}$

and u=u⁽¹⁾+iu⁽²⁾). This framework models the registration by E=∫_D|Ŝ_g−s|+λ_g∫|∇g|² du, such that |μ_g|<1, where |Ŝ_g−s| is the VFS mismatch cost, and λ_g is a positive constant that controls the smoothness of g. The mismatch cost is formulated to encourage the subject, and the wrapped topology-prior shares the same VFS at each position.Topological fMRI Decoding

Without losing generality, for a positive orientated visual area V⁺, the topological fMRI decoding model is to solve within the V⁺area, by T=arg min Σ_Vj∈V+∫|ŷ(t; v_j, σ_j)−y_j(t)|²dt+λ_v|∇v|₁² du, s. t. |μ_ƒ|<1, where |∇v|²=|∇v⁽¹⁾|²+|∇v⁽²⁾|^2. Similarly, the tRF model can be applied to the negative area by flipping the parametrization vertically. To avoid over smoothing, the values of λ_v are chosen based on the Generalized Cross-Validation.

tRF Model

The tRF finds the visual parameters for the visual cortex via two modules,

{circumflex over (L)}=L _J(g), where g=arg_g^min∫_D|Ŝ_g−s|+λ_g∫_D|∇g|²du,s.t.|μ_g|<1  (14)

{circumflex over (T)}=arg_{{circumflex over (T)}}^minΣ_ui∈D∫(ŷ−y_i)²dt+λ_v|∇v|₁²,s.t.s(u_i)=ŝ(u_i|{circumflex over (L)})  (15)

Numerical Method

A direct search for ({circumflex over (L)}, {circumflex over (T)}) is computationally infeasible because the number of parameters are typically high and ŷ is also computationally heavy. The problem may be divided into two subproblems, segmentation and fMRI decoding, and the two subproblems are updated iteratively. Since both subproblems constrain the norm of Beltrami, a technique is introduced in advance, which is called topology-projection.

Topology-Projection with Smoothing

Given a map ƒ, the topology-projection computes a smooth and topological map f′ with small changes of f According to Measurable Riemannian Mapping Theorem, the Beltrami coefficient uniquely encodes a map upon suitable normalization. Therefore, a map can be manipulated by its Beltrami coefficient. If a map whose |μƒ|>1 for some points, its norm is adjusted for those points by μ′_ƒ=αμ_ƒ=αμ_ƒ/|μ_ƒ/| (α=0.95 in this disclosure). Then μ′ is used to recover the topological map ƒ′. In specific, one can find ƒ′ by solving ∇·A∇ƒ′=0, with Dirichlet boundary condition. Here A is a 2-by-2 matrix upon μ′_ƒ[28], ∇ and ∇ are the gradient and divergence operators, respectively. ƒ′ can be solved efficiently by writing ∇ A∇ as a matrix. After the topology-projection, Laplacian smoothing is applied to ƒ′ to make it smooth.

Topology-Preserving Segmentation

To solve Eq. 14, it is divided into two subproblems: näive g searching and topology projection. Specifically, (1) for u∈u_J, the target of g(u) is updated to nearest point g′ such that s(u_i)=ŝ(g′), since it minimizes the VFS mismatch; (2) the topological projection is used to fix the topological condition for g′ (see SI for segmentation results).

Topological fMRI Decoding

Eq. 15 is also split into two subproblems: parameter searching and topology-projection. Specifically, (1) the topology-projection is used for each visual area respectively to get topological visual parameters T; (2) the gradient descent is used to update T′ parameters, T′←T′η

$\left\{\left(\frac{\partial \mathrm{Ep}}{\partial v\left(1\right)},\frac{\partial \mathrm{Ep}}{\partial v\left(2\right)},\frac{\partial \mathrm{Ep}}{\partial \sigma }\right)\right\}$

where η is a constant (updating step). The updated result T′ is further used to improve the segmentation, described in Eq. 14.

Algorithm

The overall procedures of the tRF framework is summarized in the following Algorithm 3. Note that the updating step is decayed by η^(t+1)=k_βη^(t), k_β<1, to ensure tRF converges. Its convergence proof is provided in SM.

Algorithm 3: The tRF framework
Input: Discrete region of interest S = (V, E, F); fMRI signal yi for each point; a tolerance
ϵ ∈   +; and a decay factor kβ < 1.
Results: The topology-preserving segmentation L and topological retinotopic map T = {(vi, σi)}
that sufficiently explains the fMRI signal.
Solve the pointwise pRF parameter (vi, σi) by Eq. 13
Compute the conformal disk parametrization ui = c(Vi), Vi ∈ V
Initialize T = {(vi, σi)}, and δ ← ∞
while δ > ϵ do
Compute the segmentation L, by solving Eq. 14
Update T′ on given segmentation L by solving Eq. 15
δ = max(T′ − T), T ← T′, and η(t+1) ← kβη(t)
end

Dataset

Synthetic Data

An exemplary framework has been evaluated on synthetic data. The evaluation process used double-sech model as ground truth and assigned a receptive field size σ=k₂v⁽¹⁾ (k₂=0.01, 0.02, 0.03 for V1, V2, and V3, respectively). The next step is to generate the normalized fMRI signal y₀(t) with a known specific stimulation pattern (see Kay, N., Benson, C., Jamison, K., Arcaro, M., Vu, A., Coalson, T., Essen, D. Van, Yacoub, E., Ugurbil, K., Winawer, J., Kay, K.: The HCP 7T Retinotopy Dataset, https://osf. io/bw9ec/.) and add two levels of noise to the fMRI signal y(t)=y₀(t)+r(t), where r˜N(0,γ)(γ=0.1 and γ=0.5 respectively).

Real Data with 7T MM System

An exemplary framework was also tested on the Human Connectome Project (HCP) dataset (See Benson, N.C., Jamison, K. W., Arcaro, M. J., Vu, A. T., Glasser, M. F., Coalson, T. S., Van Essen, D.C., Yacoub, E., Ugurbil, K., Winawer, J., Kay, K.: The Human Connectome Project 7 Tesla retinotopy dataset: Description and population receptive field analysis. J. Vis. 18, 1-22 (2018). https://doi.org/10.1167/18.13.23.). The HCP dataset is a publicly available retinotopy dataset on 7T fMRI scanners with high spatial and temporal resolutions.

Results

Results on Synthetic Data

The performance of an exemplary framework disclosed herein was compared with several methods, including the pRF model, the model fitting result, the Laplacian smoothing, and the registration method (which registers pRF to another parameterized Double-Sech template and uses the template's value). More specifically, given the noisy functional signal y(t) for each vertex, the errors were compared between methods' output and ground truth in Table 2, including the perception center error |Δv|, receptive field size error |Δσ|, and the numbers of triangles that violates the required VFS within its segmentation.

The results indicate that: (1) only the exemplary framework disclosed herein makes topological results (T_n=0), and the smoothing achieved the worst results, which means the smoothing does not contribute to the exemplary method's topological results in topology-projection; (2) on the other hand, the smoothing method also achieved decent precision (|Δv|˜2.), (3) the model-fitting method is promising in topological condition (T_n<5). However, it is only applicable to the V1-V3 complex but not to other visual areas (see Table 3; the flipping number is significant due to regions beyond V3); (4) the TPS registration method can ensure the topology in theory. However, because of the noisy measure of retinotopic coordinates, the anchors/landmarks are noisy and may destroy the topological condition. Those results suggest the framework of the present disclosure achieves the best accuracy under topological condition. An illustrative comparison is provided in FIG. 11 for various methods. FIGS. 11(a)-(f) illustrate the results for the V1-V3 model by ground-truth method in FIG. 11(a), pRF model in FIG. 11(b), model-fitting method in FIG. 11(c), smoothing method in FIG. 11(d), registration method in FIG. 11(e), and the tRF framework in FIG. 11(f).

TABLE 2
A performance comparison (the smaller, the better) on
synthetic data. Results for γ = 0.1 (small
noise) and γ = 0.5 (big noise) are separated by the “/” symbol.
	Method	|Δν|	|Δσ|	Tn
	pRF Method	2.924/2.313	0.902/1.100	393/483
	Smoothing	2.827/2.288	0.863/1.021	469/528
	Registration	3.549/3.518	0.974/0.971	440/444
	Model-fitting	3.559/3.590	0.902/1.100	3/5
	tRF	2.485/1.801	0.768/0.857	0/0

Results on the Real Data

An exemplary tRF was applied to subjects in the HCP dataset. Since no ground-truth is available for the real dataset, the methods were compared based on the fMRI fitting goodness (the RMSE of fitting error ē) and the number of non-topology triangles (T_n). The fitting error ē and topology violations T_n for the first three observers were listed in Table 3. The tRF achieved the best fMRI fitting and zero topology violations. The results also showed the intuitive comparison between the methods for the first observer in FIGS. 12a-e. The proposed tRF (FIG. 12e) fixed the topology violations (T_n=0) compared to the pRF model (FIG. 12a, T_n=870), model-fitting (FIG. 12b, T_n=489), smoothing (FIG. 12c, T_n=842), or registration (FIG. 12d, T_n=953). Besides, the tRF achieved the smallest fMRI fitting error (RMSE is 0.296), compared to pRF (0.300), model-fitting (0.335), smoothing (0.301), or registration method (0.305). The average results for the first twenty observers were also reported in the dataset. The tRF reduced the mean fitting error from the second-best method's 0.376 (pRF method) to 0.372, and the topology violations were reduced from 1234 (median value) to 0. Those results showed that the tRF not only reduced the topology violation (i.e., compatible with neurophysiology conclusions) but also improved/kept the fitting power.

TABLE 3
The ē/Tn (fitness/topology) of different methods for the first three HCP observers.
Observers	pRF	Model-Fitting	Smoothing	Registration	tRF
S1	0.300/870	0.335/489	0.301/842	0.305/953	0.296/0
S2	0.252/1119	0.288/934	0.253/1126	0.287/1141	0.250/0
S3	0.276/1194	0.296/808	0.276/1197	0.298/1082	0.273/0

Methods and Systems for Registering Sensory Maps

In an exemplary embodiment, a method is disclosed herein that registers sensory maps, as shown in FIGS. 13-18.

SUMMARY

It is an important topic in vision science and neurobiology to identify and analyze visual areas of the human brain (i.e. the visual atlas). Retinotopic mapping with functional magnetic resonance imaging (fMRI), provides a non-invasive way to delineate the boundaries of the visual areas. It is known in neurophysiology that retinotopic mapping is locally diffeomorphic (i.e. smooth, differentiable, and invertible) within each local area. However, the decoded visual coordinates from fMRI retinotopic mapping studies are not diffeomorphic because of the low signal-noise ratio of fMRI. The non-diffeomorphic problem is more severe near the most important fovea region because the size of the neuronal receptive field is much smaller and the retinotopic organization is more complicated. As a result, it is unreliable to directly delineate the visual atlas from a single subject's retinotopy.

To get a better visual atlas, especially for the fovea, it is preferable to use more subjects and take advantage of the group average. However, individuals' visual regions are different in size, shape and even location. Directly averaging several individuals' retinotopic maps cannot improve the result. Indeed, there are lots of methods and packages to register cortical surfaces, e.g. FreeSurfer, FSL, and Brain-Suite. These methods are designed for diffeomorphic cortical wrapping using anatomical scalar information (e.g. curvature, thickness) only.

There is an opportunity to use retinotopic maps to improve registration of cortical areas. First, unlike image registration where only scalar features are available, retinotopic mapping associates an estimated visual coordinate (although noisy) to each location of the visual cortex. Second, the quality of the estimated visual coordinates can be assessed with performance metrics, which can help emphasize high-quality locations and under-weight poor-quality locations. These features have been used to register retinotopic measurements to a template. The used algorithm morphs the subject's surface to fit template data and introduces several penalties to avoid overstretching. However, these penalties compromise registration accuracy.

There is a need to model retinotopic registration, especially the diffeomorphic constraint. Visual coordinate differences may be quantified after registration, e.g., with Euclidean distance. The question is how to define and quantify diffeomorphism. Diffeomorphism has been long discussed in geometry. The registration method is disclosed herein that applies the Beltrami coefficient, an important concept in quasiconformal geometry, to quantify the diffeomorphic condition and model the registration as a function optimization problem with constraints. The present registration method eliminates redundant registration limitations and ensures diffeomorphic result. Then the linear Beltrami solver is applied to solve the registration model. The present model may also adopt human-based landmarks as a rough/precise guideline to help registration. The method of the present disclosure improves the accuracy and robustness of retinotopic map registration

A purpose of exemplary systems and methods is to generate diffeomorphic retinotopic maps and improve the accuracy of the retinotopic atlas from fMRI measurements through the development of a specifically designed registration procedure. In consideration of retinotopic mapping features, a quasiconformal geometry-based registration model is formed and solved with efficient numerical methods. An exemplary registration method was compared with existing methods on synthetic data. The results demonstrate that it is superior to conventional methods for the registration of retinotopic maps. The application of an exemplary method to a real retinotopic mapping dataset also reduces registration errors.

Problem Restatement

Retinotopic Mapping

The idea of retinotopic mapping via fMRI is briefly introduced here. Assume a unit illumination light spot is on point v=(v⁽¹⁾, v⁽²⁾)∈² in the visual field, as shown in FIG. 13. After the visual system's processing, the light spot will eventually activate a population of neurons. The main purpose of retinotopic mapping is to find, the center v and the extent σ∈⁺ of its receptive field on the retina, for each point P=(X, Y, Z)∈³ on visual cortex.

fMRI provides a noninvasive way to determine v and a for P, based on the following procedure. (1) Design a stimulus time sequence s (t; v), such that the stimulus sequence is unique for every visual coordinate, i.e., s(t; v₁)≠s(t; v₂), ∀v₁≠v₂; (2) Present the stimulus sequence to an individual and record the fMRI signals from the visual cortex; (3) For each fMRI time sequence on a cortical location, y(t; P), determine the corresponding receptive field, including its central location v and its size σ on the retina, that most-likely generated the fMRI signals. Specifically, given the neurons' spatial response r(v′; v, σ) (a predefined model depicts of neural response around v) and the hemodynamic response function h(t) (a model of the time course of neural activation to a stimulus), the predicted fMRI signal can be written as:

{circumflex over (y)}(v,σ)=β(∫r(v′;v,σ)s(t;v_x,v_y)dv′)*h(t)  (16)

where, β is a coefficient that converts the units of response to the unit of fMRI activation, and * is the convolution operator. Then the perceptive center v and perceptive field size σ are estimated by minimizing the difference between predicted signal and measurement,

$\begin{array}{cc}\left(v,\sigma \right)=\arg \underset{\left(v,\sigma \right)}{\min }\int {\hat{y}\left(v,\sigma \right)-y\left(t\text{;}P\right)}^2\mathrm{dt},& \left(17\right)\end{array}$

The retinotopic mapping of the entire visual cortex is obtained when (v, σ) is solved for every point on the cortical surface. The goodness of the estimation for each location is inferred by the variance explained, R²=∫|ŷ−y|²dt/∫|y−y|²dt

Mathematical Model

In practice, cortical surface S is usually represented by a set of vertices P_S={P₁, P₂, . . . , P_n}={P_i} and triangular faces F_S={F_i}, denoted by S_S=(F_S, P_S). The fMRI signal is decoded for each location on cortical surface P_i∈P_S according to Eq. (17). The decoding process generates a raw visual coordinate v_i, perception size σ, and variance explained R_i² for each vertex P_i. S=(F_S, P_S, v_S, σ_s, R_s²) is used to denote the bundle data of cortical surface as well as the raw retinotopic measurements. Here the capitalized subscript is used to denote all data of a subject, e.g. P_S means all the point collection of the subject, and lowercase subscript to denote data of a point, e.g. P_i is for i-th point of the subject. If the visual coordinate v_S is accurate, then the visual atlas is of high quality. However, v_S is influenced by the fMRI signal-noise ratio, head movement, subject's behavior and so on.

The next question is how to assign a new visual coordinate {circumflex over (v)}, such that it is most likely inferred by the measurements from all the subjects' data. One approach is registering each subject to a presumed template or mathematical model and then assigning the registered template value to the subject. An exemplary approach is to find a registration method and a presumed template T=(F_T, P_T, v_T, σ_T, R_T²), such that T the overall registration cost from all subjects' data to the template is minimal, i.e., T=arg min_TΣ_JR(T, S_J), where R is the registration cost between the template T and J-th subject's data S_J.

A remaining problem may be to define the diffeomorphic registration function and its cost. It may be easy to understand the problems on 2D instead of the original 3D cortical surface. If a diffeomorphic function can be found from the cortical surface to a parametric domain, then the 3D registration can be simplified as a 2D problem, because the composite of two diffeomorphic functions is still diffeomorphic. A discrete conformal (angle preserving) mapping of the occipital region may be helpful. Specifically, as shown in FIG. 14 which is an illustration of several spaces and registration, a feature point is first defined that corresponds to the fovea as the center point. The geodesic distances from all vertices to the center point are calculated. Only the portion of cortical surface whose geodesic distance to the center point is within a certain value is considered, and this patch is mapped to a unit disk, c:Pu, u=(u⁽¹⁾, u⁽²⁾)∈².

Similarly, a similar portion of cortical surface is obtained and mapped to a parametric space, c′: P′u′,u′=(u′⁽¹⁾, u′⁽²⁾)∈² for template cortical surface (FIG. 14(c) gray color region). If a registration function ƒ:uu′ is found in the parametric space, the registration for the two cortical mesh as P′=c⁻¹∘ƒ∘c(P) can be obtained. Ideally, ƒtries to: (1) minimize the retinotopic visual coordinate differences between the subject and the template, i.e. ∫w|v_S(ƒ)−v_T|² du (w is proportional to the variance R²); (2) ensure diffeomorphism; (3) introduces some smoothness. In geometry, diffeomorphism can be quantified by Jacobian, Beltrami coefficient, and so forth.

Beltrami coefficient is adopted to quantify the diffeomorphism, considering that Beltrami coefficient also quantifies angle deviation. Angle deviation is considered here because retinotopic mapping preserves angle to a considerable extent. Since cortical surface is conformal to visual space to some extent and cortical surface is mapped to the parametric space by conformal mapping, the registration f should also preserve the angle to some extent. The Beltrami coefficient, associated with f is defined as,

$\begin{array}{cc}{\mu }_f=\left(\frac{\partial f}{\partial u^{\left(1\right)}}+i\frac{\partial f}{\partial u^{\left(2\right)}}\right)/\left(\frac{\partial f}{\partial u^{\left(1\right)}}-i\frac{\partial f}{\partial u^{\left(2\right)}}\right).& \left(18\right)\end{array}$

If |μ_ƒ|<1 then ƒ is diffeomorphic. Considering the diffeomorphic condition, angle minimization, and smoothness requirements, the next is to seek ƒminimizing the energy in Eq. (19),

E=∫w|v _S(ƒ)−v_T|²+λ₁|μ_ƒ|²+λ₂|∇ƒ|²du.  (19)

Numerical Method

It is costly to minimize energy in Eq. (19) directly as it mixes Beltrami coefficient μ_ƒ and mapping function ƒ. An exemplary solution is to alternatively solve with respect to the mapping function and Beltrami coefficient in separate steps. Tools may be utilized to convert the mapping function and Beltrami coefficient back and forth.

Beltrami Coefficient

Given the explicit form of function ƒ, it is possible to compute the Beltrami coefficient μ_ƒ according to Eq. (18). In the discrete case, ƒ is interpreted linearly on each triangle. As shown in FIG. 15(a), for u within the triangle, ƒ(u)=α_iv_i+α_j v_j+α_kv_k, where v_i=ƒ_c(u_i), v_j=f_c(u_j), and v_k=ƒ_c(u_k). The coefficients α_i, α_j, α_k are called the barycentric coefficients. Specifically, α_i (similarly for α_j and α_k) is the area portion of triangle Δuu_ju_k relative to Δu_iu_ju_k, i.e. α_i=Area(Δuu_ju_k)/Area(Δu_iu_ju_k). Now it is possible to compute the Beltrami coefficient μ_ƒc according to Eq. (18).

Linear Beltrami Solver (LBS)

The Linear Beltrami Solver (LBS) is introduced to recovery function ƒ=(f⁽¹⁾, ƒ⁽²⁾) for the given Beltrami coefficient μ=ρ+iτ. According to the definition, i.e. Eq. (19), the following is obtained,

$\begin{array}{cc}\left(\frac{\partial f}{\partial u^{\left(1\right)}}+i\frac{\partial f}{\partial u^{\left(2\right)}}\right)/\left(\frac{\partial f}{\partial u^{\left(1\right)}}-i\frac{\partial f}{\partial u^{\left(2\right)}}\right)=\rho +i\tau .& \left(20\right)\end{array}$

After reorganizing Eq. (20) and eliminating f⁽²⁾, the following is derived,

$\begin{array}{cc}\nabla ·A\nabla f^{\left(1\right)}=0,\mathrm{where}A=\left(\begin{array}{cc}{\alpha }_1& {\alpha }_2\\ {\alpha }_2& {\alpha }_3\end{array}\right),{\alpha }_1=\frac{{\left(\rho -1\right)}^2+{\tau }^2}{1-{\rho }^2-{\tau }^2},{\alpha }_2=\frac{-2\tau }{1-{\rho }^2-{\tau }^2}\mathrm{and}{\alpha }_3=\frac{1+2\rho +{\rho }^2+{\tau }^2}{1-{\rho }^2-{\tau }^2},\nabla f^{\left(1\right)}=(\partial f^{\left(1\right)}/& \left(21\right)\end{array}$

∂u⁽¹⁾+∂f⁽¹⁾/∂u⁽²⁾), and the divergence ∇·on vector G=A∇f⁽¹⁾=(G⁽¹⁾,G⁽²⁾) is defined as

$\nabla ·G=\partial G^{\left(1\right)}/\partial u^{\left(1\right)}+\frac{\partial G^{\left(2\right)}}{\partial u^{\left(2\right)}}.$

By solving the partial equation Eq. (21) with certain boundary conditions, ƒ⁽¹⁾ is solved. Similarly, ƒ⁽²⁾ is solved.

In the discrete case, gradient operator ∇ƒ⁽¹⁾(u) can be written out with the linear interpretation. For discrete divergence ∇·G operator, it is approximated on each vertex of its dual-cell (a cell includes neighboring triangles' circumcenters). In specific, as shown in FIG. 15(b), in consideration of the center vertex with its neighbors N (u_i), the method of the present disclosure approximates ∇·G (u_i) as the average of divergence on D, which is written as,

$\begin{array}{cc}\nabla ·G=\frac{1}{D}{\int }_{\partial D}\mathrm{Gds}=\frac{1}{D}\sum _{T_j\in N\left(z_i\right)}{G}_{T_j}·\left(u_k-u_j\right).& \left(22\right)\end{array}$

According to these approximations, there are a set of linear equations for ƒ_i and its neighbors. Finally, Eq. (21) is written in a matrix form and ƒ is solved efficiently.

Laplacian Smoothing and Chopping

Now it is possible to convert between the mapping function and Beltrami coefficient back and forth. During the registration, the smooth operation is applied to the function ƒ to make the registration smooth, namely the Laplacian smoothing. In specific, for a scalar function g, a smooth version of g′ is obtained by solving the following equation,

$\begin{array}{cc}{g}^{\prime }=\arg \underset{g}{\min }\int {\nabla g}^2+\lambda \int {g-g^{\prime }}^2,& \left(23\right)\end{array}$

where λ is a constant. Eq. (23) can be efficiently solved by its Euler-Lagrange equation, i.e. (−∇·∇+2λI)υ=2λμ. The smooth process is applied to ƒ⁽¹⁾ and ƒ⁽²⁾, separately. For diffeomorphic mapping, the Beltrami coefficient |μ|<1. Similarly, if there is a point whose |μ|>1, the mapping is non-diffeomorphic. To pursue a diffeomorphic mapping, the present method shrinks the magnitude of |mu| and keeps the argument arg μ (notice μ is a complex function), i.e. υ′=υ/|υ|, if |υ|>1. This is called Chop.

Simple Registration Regardless Diffeomorphism

To use the features from retinotopic mapping, instead of updating the registration function using scalar features, the present method searches (by brute-force) retinotopic visual coordinates within the nearby cortical surface region |f_i′−f_i|<r₀, and updates the registration by,

$\begin{array}{cc}{f}_i^{\prime }=\arg \underset{f_i^{\prime }}{\min }{\left(v_S\left(f_i^{\prime }\right)-v_T\right)}^2,f_i^{\prime }-f_i<r_0& \left(24\right)\end{array}$

Algorithm

The overall registration process is summarized in the following Algorithm 4.

Algorithm 4: Diffeomorphic Registration
   Input: Subject's data S = (FS, PS, vs, σs, Rs2),
   template T = (FT, PT, vT, σT, RT2) and ϵ
   Result: Registration function f
1  | begin
2  | Compute the conformal maps uS = cs (FS, PS)
   | and uT = cT (FT, PT);
3  | Take Idetical mapping as initial uT = f (uT),;
   | repeat
4  | | Update registration f′ according to Eq. (24);
5  | | Apply Laplacian smooth to f′
6  | | Compute Beltrami coefficient ν and Chop;
7  | | Compute new mapping f′ from ν′ by LBS;
8  | | until max |f′ − f| < ϵ;
9  | | return f
10 end

Results

Performance on Synthetic Data

Synthetic subject and template data in the scenario of retinotopic mapping are generated according to the double-sech model. The subject and template data are generated by two different sets of parameters, as shown FIGS. 16(a)-(e). The subject data is added with a small, as shown in FIG. 16(d), and big amount, as shown in FIG. 16(d), of uniform noise, respectively. The next step is to calculate the ground truth registration displacement from the subject to the template, as shown in FIG. 16(e). Then several methods are used to register the noisy subject data to the predefined template and get registration displacement. Eventually, the registration error (difference between method's registration displacement and ground truth displacement) is reported in the following Table 4.

The results indicate that: (1) an exemplary algorithm of the present disclosure achieves the smallest registration error and ensures diffeomorphism; (2) LDDMM and QCHR can also achieve almost diffeomorphic result, however, the method of the present disclosure is more precise; (3) TPS is fast, but its accuracy is poor; (4) LDDMM, D-Demos, and QCHR ensure diffeomorphism for image registration, however, they cannot guarantee overall diffeomorphism when the visual coordinates are considered as intensity of image; and (5) Because of an alternative scheme in an exemplary method of the present disclosure, it usually takes approximately 100 to 150 seconds to register.

TABLE 4
Comparing registration performance relative to the ground
truth. Results for small noise and big noise cases are reported
respectively (Running time is reported in average). Landmarks/anchors
are given for methods with a “*” symbol.
	Registration Error	#Overlap
Method	Mean	Max	Triangles	Time/(s)
TPS* [18]	0.75/4.47	1.87/14.43	2/18	0.36
Bayesian* [9]	3.89/3.98	6.59/9.45	2079/2469	110
D-Demos [19]	1.12/1.20	2.34/2.59	2377/2473	1.6
LDDMM [20]	0.59/0.64	1.27/1.49	1/1	86
QCHR* [11]	0.06/0.10	0.62/0.47	1/0	173
The Present	0.04/0.08	0.12/0.19	0/0	123
Method*

HCP Data

A registration method of the present disclosure has been applied to the first twenty subjects of the Human Connectome Project (HCP) retinotopic dataset. The original retinotopic result is known. The first step of test is to initialize σ template by averaging the subjects' data. Then, all subjects' data are registered to the initial template. FIG. 17 shows one subject's registration process. The curves are the boundaries of V1/V2/V3 inferred from the template. Before application of an exemplary registration method, as shown in FIGS. 17(a)-(b), the boundary is not aligned especially in the fovea region. FIGS. 17(a)-(b) result is from a registration result which used both the structural as well as some fMRI feature. FIGS. 17(c)-(d) illustrate results after application of an exemplary registration method. This means the alignment performed herein can further improve the result by incorporation of more retinotopic information.

Besides the subject registration, the overall registration cost is defined as the sum of visual coordinate error across subjects in the dataset. After registering all subjects' data, the average of the subjects is obtained and the same method to register the template is applied to the average data. This process is repeated until satisfactory results are obtained. FIG. 18(a)-(c) shows the first, second, and last iteration of the template, respectively. FIG. 18(d) shows the registration cost, overall visual error, is decreased during the process. The decreased overall visual error means the template is improved for human retinotopic maps. The registration method is also applicable in higher visual regions.

Sensory Mapping Method for a Human Brain

In various exemplary embodiments, a sensory mapping method 100 for a human brain is illustrated by the flow chart in FIG. 19. Method 100 comprises flattening the cortical surface of the human brain (step 102). In various embodiments, the flattening may comprise a conformal parametrization method. The flattening may comprise applying harmonic mapping to a unit disk and refining repetitively parametric coordinate for each point on a parameter domain to produce conformal mapping.

Method 100 further comprises projecting functional imaging data onto the flattened surface (step 104). Method 100 further comprises smoothing the functional imaging data or processing the functional imaging data to produce topological results (step 106). In various embodiments, the smoothing may comprise a topological smoothing method that utilizes a diffeomorphic smoother. In various embodiments, the smoothing may comprise applying Beltrami coefficient to quantify diffeomorphism and modeling the smoothing in an optimization framework.

In various embodiments, the processing may comprise a topology-preserving cortical surface segmentation method. In various embodiments, the processing may comprise segmenting the cortical surface of the human brain by preserving the prior connectivity pattern of the cortical surface areas, decoding the functional imaging data based on at least one topological condition within each cortical area, and repeating the segmentation and decoding.

In various embodiments, method 100 further comprises generating the sensory map (step 108). Method 100 may comprise registering sensory maps across individuals (step 110). In various embodiments, the registering may be diffeomorphic. In various embodiments, the registering may comprise a method that aligns a human brain sensory map of an individual to a predefined template. In some exemplary embodiments, method 100 comprises analyzing the maps in the common space (step 112).

While the principles of this disclosure have been shown in various embodiments, many modifications of structure, arrangements, proportions, the elements, materials and components, used in practice, which are particularly adapted for a specific environment and operating requirements may be used without departing from the principles and scope of this disclosure. These and other changes or modifications are intended to be included within the scope of the present disclosure.

The present disclosure has been described with reference to various embodiments. However, one of ordinary skill in the art appreciates that various modifications and changes can be made without departing from the scope of the present disclosure. Accordingly, the specification is to be regarded in an illustrative rather than a restrictive sense, and all such modifications are intended to be included within the scope of the present disclosure. Likewise, benefits, other advantages, and solutions to problems have been described above with regard to various embodiments. However, benefits, advantages, solutions to problems, and any element(s) that may cause any benefit, advantage, or solution to occur or become more pronounced are not to be construed as a critical, required, or essential feature or element.

As used herein, the terms “comprises,” “comprising,” or any other variation thereof, are intended to cover a non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements does not include only those elements but may include other elements not expressly listed or inherent to such process, method, article, or apparatus. Also, as used herein, the terms “coupled,” “coupling,” or any other variation thereof, are intended to cover a physical connection, an electrical connection, a magnetic connection, an optical connection, a communicative connection, a functional connection, and/or any other connection. When language similar to “at least one of A, B, or C” or “at least one of A, B, and C” is used in the specification or claims, the phrase is intended to mean any of the following: (1) at least one of A; (2) at least one of B; (3) at least one of C; (4) at least one of A and at least one of B; (5) at least one of B and at least one of C; (6) at least one of A and at least one of C; or (7) at least one of A, at least one of B, and at least one of C.

Claims

1. A sensory mapping method for a human brain, the method comprising: flattening the cortical surface of the human brain;projecting functional imaging data onto the flattened surface;smoothing the functional imaging data; andgenerating a sensory map.

2. The method of claim 1, wherein the smoothing comprises a topological smoothing method.

3. The method of claim 1, wherein the smoothing comprises: determining coordinates of each cortical area based on the functional imaging data;repeating a process until the maximum absolute value of the Beltrami coefficient for each coordinate is less than one.

4. The method of claim 1, wherein the smoothing comprises: applying Beltrami coefficient to quantify diffeomorphism; andmodeling the smoothing in an optimization framework.

5. The method of claim 1, wherein the flattening comprises a conformal parametrization method.

6. The method of claim 1, wherein the flattening comprises: applying harmonic mapping to a unit disk; andrefining repetitively parametric coordinate for each point on a parameter domain to produce conformal mapping.

7. The method of claim 1, further comprising rendering the sensory map with an image of stimuli that activate the cortical surface.

8. The method of claim 1, further comprising registering a plurality of sensory maps across different individuals.

9. The method of claim 8, further comprising analyzing the plurality of sensory maps in a common space.

10. The method of claim 8, wherein the registering is diffeomorphic.

11. A registration method that aligns a human brain sensory map of an individual to a predefined template, wherein the method comprises: registering a plurality of sensory maps, each sensory map corresponding to a human brain of a different individual; andanalyzing the plurality of sensory maps in a common space.

12. The method of claim 11, wherein the registering is diffeomorphic.

13. The method of claim 11, wherein the registering comprises: applying Beltrami coefficient to quantify diffeomorphism; andmodeling the registering in an optimization framework.

14. A sensory mapping method for a human brain, the method comprising: flattening the cortical surface of the human brain;projecting functional imaging data onto the flattened surface;processing the functional imaging data to produce topological results; andgenerating a sensory map.

15. The method of claim 14, wherein the processing comprises a topology-preserving cortical surface segmentation method.

16. The method of claim 14, wherein the processing comprises: segmenting the cortical surface of the human brain by preserving the prior connectivity pattern of the cortical surface areas;decoding the functional imaging data based on at least one topological condition within each cortical area; andrepeating the segmentation and decoding.

17. The method of claim 14, further comprising registering a plurality of sensory maps across different individuals.

18. The method of claim 17, wherein the registering is diffeomorphic

19. The method of claim 17, further comprising analyzing the plurality of sensory maps in a common space.

20. The method of claim 17, wherein the registering comprises: applying Beltrami coefficient to quantify diffeomorphism; andmodeling the registering in an optimization framework.