Coupled Reconstruction of Cortical Surfaces by Diffeomorphic Mesh Deformation
Authors: Hao Zheng
Hongming Li
Yong Fan
Venue: NeurIPS 2023 poster (2023)
Keywords: Brain MRIs
cortical surface reconstruction
deep learning
TL;DR: We develop a new DL-based framework to jointly reconstruct the inner, outer, and midthickness surfaces with topology constrains and estimate cortical thickness directly from 3D MRIs.

Abstract: Accurate reconstruction of cortical surfaces from brain magnetic resonance images (MRIs) remains a challenging task due to the notorious partial volume effect in brain MRIs and the cerebral cortex's thin and highly folded patterns. Although many promising deep learning-based cortical surface reconstruction methods have been developed, they typically fail to model the interdependence between inner (white matter) and outer (pial) cortical surfaces, which can help generate cortical surfaces with spherical topology. To robustly reconstruct the cortical surfaces with topological correctness, we develop a new deep learning framework to jointly reconstruct the inner, outer, and their in-between (midthickness) surfaces and estimate cortical thickness directly from 3D MRIs. Our method first estimates the midthickness surface and then learns three diffeomorphic flows jointly to optimize the midthickness surface and deform it inward and outward to the inner and outer cortical surfaces respectively, regularized by topological correctness. Our method also outputs a cortex thickness value for each surface vertex, estimated from its diffeomorphic deformation trajectory. Our method has been evaluated on two large-scale neuroimaging datasets, including ADNI and OASIS, achieving state-of-the-art cortical surface reconstruction performance in terms of accuracy, surface regularity, and computation efficiency.

Coupled Reconstruction of Cortical Surfaces by Diffeomorphic Mesh Deformation
Hao Zheng
Hongming Li Yong Fan
University of Pennsylvania
hzheng1 @upenn. edu, (hongming... li,yong faang..fnn@@eeennn openmmediiiin upenn. edu
Abstract
Accurate reconstruction of cortical surfaces from brain magnetic resonance images (MRIs) remains challenging task due the notorious partial volume effect in brain MRIs and the cerebral cortex's thin and highly folded patterns. Although many promising deep leamiin-bbsst cortical surface reconstruction methods have been developed, they typically fail to model the interdependence between inner (white matter) and outer (pial) cortical surfaces, which can help generate cortical surfaces with spherical topology. To robustly reconstruct the cortical surfaces with topological correctness, we develop new deep learning framework jointly reconstruct the inner, outer, and their in-between (midthickness) surfaces and estimate cortical thickness directly from 3D MRIs. Our method first estimates the midthickness surface and then learns three diffeomorphic flows jointly to optimize the midthickness surface and deform it inward and outward the inner and outer cortical surfaces respectively, regularized by topological correctness. Our method also outputs cortex thickness value for each surface vertex, estimated from its diffeomorphic deformation trajectory. Our method has been evaluated on two large-scale neuroimaging datasets, including ADNI and OASIS, achieving state-of- the-art cortical surface reconstruction performance in terms of accuracy, surface regularity, and computation efficiency.
Introduction
The analysis of the cerebral cortex using magnetic resonance imaging (MRI) is crucial in under- standing neurodegenerative diseases 99 43) and psychological disorders 442 Since the cerebral cortex is thin (a few millimeters thick) and highly folded sheet between the inner (white matter: WM) and outer (pial) surfaces, voxel-based segmentation methods cannot accurately capture its complicated morphology [17]. Instead, triangular meshes have been widely used for cortical surface reconstruction (CSR) [4] in order to accurately measure cortical thickness, volume, and gyrification. Although well-established methods for CSR (e.g., FreeSurfer [17], BrainSuite (48) produce promising results, they often require significant computational resources (e.g., 6h/subject [17] and may necessitate manual editing in order to attain sub-voxel precision.
Recently, deep learning (DL) methods have achieved significant improvement in CSR terms of both accuracy (sub-voxel error) and efficiency (orders of magnitude faster) [10] 13 19
20
26 30
31
These DL methods can broadly classified into two categories according to representation of the output surfaces.
(I) Implicit surface representation, such as signed distance function 13 [9]. occupancy filed [13], and level set [41, can predicted by neural networks,
and cortical surface is then generated by marching cubes algorithm [27]. (II) Explicit surface reconstruction methods take a coarse or fine initial mesh as input and directly predict a target
mesh [00] 20 26 30 31 59
A detailed comparison of the existing DL-based CSR methods is
summarized in Table
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Table I: A comparison of existing DL-based CSR methods and ours. For input, "volume": whole 3D volume; "cube": a sub-volume; FS: fine surface mesh; CS: coarse surface template; IS: implicit surface representation. For output, explicit surface: triangular mesh; implicit surface: signed distance function or occupancy field. Network
Input
Qurput
Joint surf.
Direct thickness estimation
Methods
Img.
Surf.
GNN/MPP
Explicit surf.
reconstruction
DR
18
volume
vox2surf
17
vol.
FastCSR
11
volume
PiaINN(3]
cube
FS
/
CorticalFlow 26
volume
vox2eortex
[M
volume
CS
TopoFit 20
volume
CS
Co 10
cube
FS
Ours
volume
FS+18
However, the existing DL-based CSR methods are subject to limitations. First, the interdependence between the inner and outer cortical surfaces is generally ignored, and therefore separate or multi- stage DL models typically trained to reconstruct both the inner and outer cortical surfaces and intersections between them may occur. Even if both the inner and outer surfaces can be reconstructed simultaneously [10], they loosely combined with topological constraints. Second, complex DL architectures are commonly used, with separate learning of image and surface/vertex features using both convolutional neural networks (CNNs) and graph neural networks (GNNs)/multi-layer perceptrons (MLPs). Moreover, graph convolutions become less scalable as the number of vertices in the mesh grows accommodate complex shapes, and often fail learn diffeomorphic mappings to produce genus-zero regular meshes [10]. Third, using a coarse mesh template incurs difficulty of learning large deformations for highly folded cortical regions and may lead non-smooth deformation and undesirable artifacts. Last but not least, the cortical thickness estimation is neglected in all existing DL-bbased CSR methods (i.e., needs separate step to compute cortical attributes as in conventional pipelines [17, (56)). It may serve as an anatomical constraint to couple the inner and outer cortical surfaces, help generate cortical surfaces with topological correctness, and facilitate quantitative analysis of the cortical thickness.
In order robustly reconstruct the cortical surfaces with topological correctness, we develop DL-bbased approach to simultaneously reconstruet both the inner and outer cortical surfaces and estimate the cortical thickness by optimizing and deforming an initialized midthickness surface. First. our method explicitly couples the inner and outer surfaces by jointly learning three diffcomorphic flows optimize the initialization midthickness surface to lie halfway between the inner and outer cortical surfaces and deform it to inner and outer cortical surfaces, respectively. Second, instead of designing a complex mixed architecture of CNNs and GNNs/MLPPs, our method employs a single model of 3D CNNs to predict the diffeomorphic flows from multi multi-channel input, consisting of 3D brain MRI, a ribbon segmentation map that encodes structural information of the cerebral cortex, and signed distance function that implicitly encodes the initialization surface. Third, our method calculates the diffeomorphic deformation trajectories in continues coordinate space rather than on a 3D voxelwise grid, achieving higher sub-voxel accuracy from fine- rre-grained velocity fields while maintaining reasonable computational efficiency even as the number of mesh vertices increases. Fourth, we devise an efficient and reliable approach to initialize a fine midthickness surface from cortical ribbon segmentation result, followed by topology correction to ensure genus zero. Finally, a vertex-wise thickness estimation can be obtained by tracing the geodesic trajectory of each vertex during the mesh deformation process. In summary, our new DL framework differs from the existing DL-based CSR approaches in its coupled reconstruction of multiple surfaces and the simultaneous cortical thickness estimation, which facilitates robust reconstruction of the cortical surfaces with spherical topology. Ablation studies and comparison experiments on large public datasets (ADNI (22 and OASIS (33) have demonstrated that our method attains superior performance over the state-of-the-ard methods [10 13 26 30
2 Related Works
2.1 Learning-based Cortical Surface Reconstruction
Recent years have witnessed a surge of interest in geometric DI.-based methods for general computer
vision tasks [15 29
49 52 54 57 58) and biomedical object reconstruction (24 25 37 38, 55).

---

However, their applications to biomedical tasks are limited to organs with simple shapes such as the liver and heart. The cerebral cortex, on the other hand, is a highly folded, thin structure with a significantly complex shape, necessitating more advanced approaches.
Table [ summarizes two categories of CSR methods, i.e., implicit and explicit surface reconstruction. The implicit surface reconstruction methods typically learn function that maps 3D coordinates to a continuous implicit representation of the surface, such as signed distance function (39] and occupancy
field [34. Compared to well-established pipelines 177, these methods have significantly improved inference efficiency (reducing time cost from hours to minutes). Despite their ability to generate surfaces at any desired resolutions [13] 18 19 41), these methods often rely on a marching cubes algorithm [27] to obtain a triangular mesh and post- prosessoess topology correction methods 7 for ensuring correct spherical topology of the reconstructed surfaces.
The explicit surface reconstruction methods typically learn a mapping from an initialization mesh to
a target surface and the mapping can be modeled as a deformation model [52 53 555 Particularly, Pixe12Mesh [52] utilizes GNNs to learn vertex-wise deformations of an ellipsoid and increase the mesh resolution in a coarse-to-fine fine manner. Such a strategy has also been adopted in Voxe12Mesh (55) to reconstruct simple human organ surfaces (e.g.. liver, hippocampus) from CT/MR images. Similarly, PiaINN 31 learns to deform an input inner cortical (i.e.. WM) surface to a target outer (pial) surface by a sequence of deformation blocks; TopoFit 200 warps a topologically-cont correct template surface to fit the WM surface by series of graph convolutional blocks; and Vox2Cortex [10] deforms a brain template surface to target cortical surfaces by leveraging combined CNNs and GNNs. Another line of research leverages ODE to parameterize the deformation of vertices as diffeomorphic flows, such as CoticalFlow methods [26 47 and CortexODE (30.. However, all these methods may produce intersecting inner and outer cortical surfaces because they predict the inner and outer cortical surfaces either using separate DL models or using a joint DL model without explicit constraints to penalize the generation of intersecting surfaces.
2.2 Diffeomorphic Deformation
diffeomorphism theoretically differentiable and invertible and guarantees smooth and one-to-one
mapping (46). Diffcomorphic deformation has been widely used in image registration 4 35 51
can be generated from a velocity field v by integrating an ordinary differential equation (ODE) ..
dø(x, 1) di
= v(o(x,)) t), and thus o(x, t) = o(xx 0) + J. v t), t)dt,
(I)
where oxx, 0) x. The velocity field can be stationary 4 or time-varying (8. Standard numerical integration techniques, such as the Euler method and the Runge- Kutge-Kuttaaan method M. can be used to perform the integration. CoticalFlow methods 26 47 and CorteXODE (30] parameterize the ODE by a neural network [1]. CoticalFlow requires the deformation field to be Lipschitz each time step to ensure bijective mapping with Lipschitz inverse and derives stability condition for the numeric approximation of Φb As a chain of deformation modules trained in sequential stages to predict a series of deformation fields that deform the initial mesh template to the target surface in a coarse- to-fine manner, there is a risk of generating self-intersections due to the low-resolution template in the coarsest level model and the whole process is prone to error accumulation in multiple steps. CorteXODE (30 is built upon NODE [12] and operates on a topologically -orreeeeel nigh-resolution initialization mesh The sufficient condition of diffeomorphism can be satisfied if the deformation network is Lipschitz continuous and a sufficiently small step size is used in numerical approximation. Most learning- beaseng-bassed based image registration approaches 66] 28 36] adopt stationary velocity field (SVF), and the integral (displacement field) is computed on a voxel-gri -eel-grrd using the scaling and squaring method (3). yielding comparable registration accuracy to conventional methods while significantly improving efficiency. Our work lies at the intersection of these methods in that our method uses CNNs parameterize multiple SVFs based on a multi-channel input and train the SVFs jointly to optimize cortical surfaces under prrpeeyy-rrssrving and inverse-consistent transformation regularizations.
3 Methodology
As illustrated in Fig. D. our framework consists of two parts: a pipeline that estimates a midthickness surface, represented both explicitly as 3D mesh and implicitly as 3D distance map (Sect. 3.1) and
3

---

Image
Segnsentation
iaput
VEL
FCN
DDM
DDM
3D DNN
DUM
DCM
VF2
Unifermuephic Deformation Medale
Surt. Extraction
Roter
Topologgy
Inital Surfhce
VF
Trassfomation
Midthickaess Surface
Mid sUrt
Velacity Field
Surface initialization
(b) Cortical surface reconstruction network
Figure 1: Overview. (a) Surface initialization. An FCN generates cortical ribbon segmentation maps from the raw image. An SDF derived, followed by topology correction and marching cubes algorithms, and then used for initializing the midthickness surface. (b) Cortical surface reconstruction network. takes a multi-channel input and learns three parallel velocity fields (VFs) in parallel. VFi used to generate deformation field for optimizing the midthickness surface Sar- VF2 (or VF3) is used to deform SM to the WM surface Sw (or pial surface Sç) and then VF3 (or VF2) deforms Sw (or Sa) to midthickness surface Sh (or S to constrain topology between cortical surfaces.
an end-to-end fully convolutional network (FCN) that reconstructs multiple surfaces and estimate cortical thickness simultancously (Sect. 3.2). Loss functions are presented in (Sect. 3.3
3.1 Midthickness Surface Initialization
The human cerebral cortex is 2D sheet with an average thickness of 2.5mm and has highly folded geometry with peaks (i.e., gyri) and grooves (i.e., sulci) [17]. The existing CSR methods typically use a smoothed template surface estimated from group of subjects [TO. 47 or WM surface [30 to initialize the surface reconstruction. However, our preliminary experiment shows that the closer the initial surface is to its target surface, the higher the reconstruction accuracy is (see Fig. 2 (()). Thus, we propose to extract the midthickness layer which lies halfway between the inner and outer cortical surfaces as initialization surface, which brings three advantages: (1) Deforming a surface from the midthickness surface reduces the learning difficulty and improves accuracy by avoiding learning "large" deformations. It also strikes balance between optimizing the inner and outer surfaces, making it less challenging. (2) The distance between the midthickness surface and the pial surface can help prevent topology errors in
of deep sulci, this layer still provides clear separation between them. (3) Deforming the midthickness surface inward and outward to the inner and outer surfaces establishes a one-to-one mapping that explicitly encodes the correspondence between surfaces, facilitating coupled surface learning and improving CSR accuracy.
A key challenge is how to obtain a surface close to the midthickness surface as possible. A straightforward method to inflate the WM surface along the normal direction by pre-determined distance 330 However, this method cannot accurately estimate the midthickness surface in that the cortical thickness varies across the cortex and the normal merely an approximation using neighboring faces. We propose a more accurate method by leveraging the cortex ribbon segmentation map (i.e., the filled interior area of WM and pial surfaces). Given an input brain MRI volume E RD, DD utilize 3D U-Net 44 to generate WM segmentation map Mw E R, D
and a GM segmentation map Mg E TR. DixDxx (see Fig. aa)). The network trained on a large- scale public neuroimaging dataset [22] by minimizing the cross-entropy loss between the prediction and ground truth which can be obtained using existing pipelines [17 45). Based on the predicted segmentation map Mw. we generate signed distance function (SDF), Kw E R xD, D3, using distance transform algorithm: d(v;)) SD F (1;) ) is the minimal Euclidean distance of voxel E the boundary voxels. Voxels with values equal zero represent the surface boundaries and voxels with negative or positive values encode their distances to the surface boundaries inward or outward, respectively. Similarly, we generate an SDF, KG, for the pial (gray matter) surface. A new SDF can be obtained by averaging the WM and GM SDFs: KM (Ka Kw)) )/2, whose 0-level defines the midthickness surface implicitly. To ensure the midthickness surface maintains a spherical topology, fast topology check and correction algorithm 7 30 is then applied to the implicit surface KM-

---

ASSD (pea)
Paal Surf
Var Pond
Pial Stett
Velocity filed
Urd
Vi
Pard
Pase
Paue
Vur Velscity filed
Mad
Mid
Valdy
gring mnvedd
Saat
Injedtyy
MM
DN Plerzz
Verto's trnjestory
(8)
(
and Var resp.
(c)
Zavrpuattion
Figure 2: (a) The relationship between the initialization surface position (X-axis) and the surface reconstruction accuracy (Y-axis) using PialNN BT. Orange and blue lines represent WM and pial surfaces reconstruction resp. (b) Illustration of a symmetric deformation trajectory and cortical thickness estimation. (c) Illustration of DDM. Each vertex is deformed by sampled velocities.
Finally, the initialization midthickness surface is extracted by the marching cubes algorithm 27 from the O-level of Kas and parameterized by a triangular mesh Sg-
3.2 Coupled Reconstruction of Cortical Surfaces
Based on the formulation in Eq. [ an accurate parameterization of the SVFs is crucial to modeling the diffcomorphic trajectory of each vertex and reconstructing the cortical surfaces. In this section, we will show how our method achieves the goals.
Feature Extraction from a Mutli-channel Input of Brain MRI, Cortical Ribbon, and Cortical Surface. The existing DL-based CSR studies have demonstrated that is critical to fuse both image features extracted from the MRIs using CNNs and geometry features extracted from surface meshes using GNNs/MLPs. However, besides employing two different network architectures, these methods typically learn the image and geometry features separately before fusion, which does not adequately utilize image texture and surface geometric information. To better model the SVFs, we propose to learn features from multi-channel input consisting of brain MRI, its cortical ribbon segmentation maps, and its midthickness surface represented as an SDF: Leoms = I@M eg @ M. where @ is channel-wise concatenation and MWeG is a multiclass mask (BG = O. WM = 0.5, GM = l).
Such a feature learning procedure brings two key benefits: (1) Utilizing heterogeneous features enables mutual knowledge distillation. The brain MRI contains detailed texture and semantic information but may include noise and irrelevant regions far from the target surfaces. The cortical ribbon segmentation maps contain structural/semantic information about the cortical sheet and can act as an attention guide for extracting informative features around its boundaries. The SDF implicitly embeds the surface location information and relative relation between all voxels. Together, the multi-channel input provides richer and complementary information for our model to reconstruct the surfaces. (2) Only single CNN is needed, which simplifies the model design and improves efficiency. Features can be extracted in a single forward pass for all coordinates. When scaling up the number of vertices in mesh we can interpolate in the feature space efficiently.
Coupled Learning of Cortical Surfaces. The goal is learn diffeomorphic deformations that deform the initialization midthickness surface Sy C R to its target WM and pial surfaces, Sw and Se. Taking into account the discrepancy between the initialization S, and the true midthickness surfaces SM, our method also learns diffeomorphic deformation to optimize the initialization midthickness surface. In total, our method learns function to model three diffeomorphic defor- mations Jo(Teomb: Sn) (#M, Φw, ΦG). using an FCN and several diffeomorphic deformation modules (DDMs). Specifically, the FCN has a similar architecture as U-Net [4), consisting of a 5-level hicrarchical decoder with skip connections as shown in Fig. bb) (see Supplementary Materials for details). By taking the multi-channel input, the FCN learns to estimate three dense SVFs jointly, denoted by Vas. Vw, and Vg. We then use V to compute ΦM that deforms So to the true midthickness surface SM, Vw to compute Φw that deforms SM inward to the WM surface Sw. and ve to compute dG that deforms SM outward the pial surface Sç. By doing so, establish one-to-one mapping across Sw, SM, and Sq-
However, since dw and dे are computed using different SVFs, they may cause non-invertible transformation around the midthickness surface. To address this issue, we utilize the property of diffcomorphic mapping to compute symmetric deformation trajectory of each vertex and devise
5

---

symmetric cycle function Leye (Eq. 2) for training. Fig. 2(b) illustrates the deformation of a vertex PMid outward to PGM using @G (Va) followed by the deformation inward Paraa using @w (Vw). which should be as close to PMid as possible. Similarly, PMrid should also be as close to PMid as possible when deformed inward by @w followed by the outward deformation by ΦG. The symmetric cycle loss is formulated as:
Leye
Σ Poww PESM
||
2
where Poice, represents deforming vertex p using velocity fields Va and Vo sequentially. ensures si- multaneous alignment of Páwo = = Pegcew for
Luye
wi Agys
Trojectey from Sana to Sqe
each vertex E SM. Fig. F demonsttraees the effective- ness of Loye in the coupled reconstruction of multiple surfaces in the network optimization. Moreover, we can trace the geodesic trajectory of each vertex during the mesh deformation process for estimating vertex-wise cortical thickness (Fig. (())..
Trajectory from Sigu Syne
Figure 3: Example vertex deformation tra- jectories with and without Leye-
Diffeomorphic Deformation Module (DDM). To numerically solve the ODE in Eq. I scaling and squaring (SS) method
3 can be applied on the image grid, followed by trilinear interpolation (Lint()) for vertices in contimuous coordinates. However, there are three limitations: (1) Interpolation cannot guarantee the invertibility of diffeomorphic mapping; (2) Numerical errors are amplified in the two sequential steps of SS+Lint(-): ), (3) Computation on grid voxels including no vertices is unnecessary. Hence, we utilize the DDM to directly compute the vertex-wise integral. As shown in Fig, (c), obtain the velocity vector for vertex with coordinate x by interpolating its neighboring velocity vectors (:= (xx)),, i.e., V Lint( 7 \(x)). ). The vertex then moves to a new coordinate V x T where T' is the total time steps. We can obtain overall deformation using this procedure in T steps.
3.3 Loss Functions
We design multiple loss functions to optimize the geometric precision of the reconstructed surfaces and regularize the SVFs for diffeomorphic deformation.
Mesh loss. It aims to minimize distances of the vertices between the predicted surface meshes Sw (and Sa) and their corresponding ground truth (GT) meshes S, by the bidirectional Chamfer distance [26]:
Lahw
Σ
min |p PP-p.1.
Σ
min |p. PES
-pll
P.ESw.
PESW
P. P.ESW,
where p and p. are the coordinates of vertices on meshes. We can compute LenG analogously. The mesh loss is L.h- = Leew + L.iG-
Trajectory loss. Starting from the midthickness surface, the trajectory length of the vertex moving to the WM and the pial surfaces should be equal. We propose to compute the mean square difference of the vertex's trajectories:
Laist
1 N
Σ (LMMa-cMr(p) Lardawwa(())) PER
where LMGA-GMM(P Et (p) the accumulated Euclidean distance (i.e., trajectory length) of T steps of deformation. This term encourages the midthickness surface can be deformed to the inner and outer cortical surfaces with the same deformation path length. In other words, the midthickness surface should lie halfway between the inner and outer cortical surfaces.
Symmetric cycle loss. We formulate it as Eq. 2 to encourage dw and èG be invertible.
Symmetric similarity loss. To optimize the midthickness surface to lie halfway between the inner and outer cortical surfaces, magnitude difference constraint is adopted directly on the SVFs:
L. = l|vG Vw 11:
(5)
where Vw represents reverse-directional vw.

---

Normal consistency loss. We also incorporate a normal consistency regularization term to promote robust learning of the surfaces and ensure their smoothness:
Luo E (1- (1-cos(nfo nf2)),
(6)
where e an edge, So and /1 are e's two neighboring faces with their unit normals fo and nf.
In summary, we combine all the losses to jointly optimize our DL model: L = Å Leh + A2Laisst AgLeye + Xg + AgCner where (=1..., are weights balance the loss terms. We empirically set X = 1 1 1, ..... ,) and As = 0.001.
4 Experiments
We evaluated our method for reconstructing both white-matter (WM) and pial surfaces large- scale datasets, (ADNI) [22] and OASIS [33], and compared it with state-of-the-art (SOTA) DL-based CSR methods. We also tested its robustness and performed ablation analyses.
Datasets. The ADNI-1 [2 dataset consists of 817 subjects and we randomly split it into 654, 50, and 113 subjects for training. validation, and testing. respectively. The OASIS-1 (33 dataset consists of 413 subjects and we randomly split into 330, 25, and 58 for training, validation, and testing, respectively. The models were trained the training set until they reached loss plateau on the validation set, after which their performance was evaluated on the test set. We followed pre-processing protocols in previous works [10] 13 26 30 for fair comparison. The Weeeghhted MRI scans were aligned rigidly to the MNI152 template and clipped to the size of 192 x 224 x 192 mmm isotropic resolution. The pseudo ground-truth of ribbon segmentation and cortical surfaces were generated using FreeSurfer v7.2.0 17. The intensity values of MRI scans, ribbon segmentation maps, and SDFs were normalized to [0, and the coordinates of the vertices were normalized to 1]. The WM and GM in the cortical ribbon segmentation maps were assigned values of 0.5 and respectively.
Implementation details. Our framework was implemented in PyTorch 40 and trained on an NVIDIA P100 GPU of 16 GB memory. The 3D U-Net 44 for ribbon segmentation was trained for 200 epochs using Adam (3 optimization and achieved an average Dice index of 0.96 on the testing set. The CSR model was trained for 400 epochs using Adam (81 0.9, 99 0.999,
16-10 learning rate le to optimize the midthickness surface and reconstruct the WM and pial surfaces for each hemisphere. The surface meshes had 130K vertices.
Evaluation metries. We utilized three distance-based metrics to measure the CSR accuracy, including Chamfer distance (CD), average symmetric surface distance (ASSD), and Hausdorff distance (HD). In particular, CD measures the mean distance between two sets of vertices [16 52]; ASSD and HD measure the average and maximum distances between two surfaces [13 [0]] They computed bidirectionally over ~30K points uniformly sampled from the predicted and target surfaces. A lower distance indicates better result. We used the 90th percentile instead of the maximum because HD is sensitive to outliers [21]. We also utilized the ratio of self-intersection faces (SIF) to measure surface quality 13 30
4.1 Comparison with Related Works
From the two categories of existing DL-based cortical surface reconstruction methods described in Section [ we selected representative ones from each category for comparison. The experimental results are summarized in Table and illustrated in Figure 4
Main Results & Analysis. It is evident that our method achieved substantial improvement on both the WM and pial surface reconstruction over other approaches. Since DeepCSR [3 predicts an SDF-based implicit surface and requires post-processing to correct topology and extract mesh, its results may contain no SIFs but were less accurate compared with the explicit CSR methods. Starting from the WM surface, PialNN [1 can achieve sub-voxel accuracy but its SIF ratio was relatively high. Vox2Cortex [30] can generate multiple surfaces from different template meshes. employs complex CNN and GNN models to model the deformation for each vertex but has no diffeomorphism
guarantee. The promising results of 26 30 47 indicated that using neural networks to parameterize
the ODE can facilitate the diffcomorphic deformation, yielding better CSR accuracy. Our method
7

---

Table Quantitative analysis of cortical surface reconstruction on geometric accuracy and surface quality. The Chamfer distance (CD), average symmetric surface distance (ASSD), Hausdorff distance (HD), and the ratio of the self-intersecting faces (SIF) were measured for WM and pial surfaces on the two datasets. The mean value and standard deviation are reported. The best ones are in bold.
L-Pial Surface
L.WM Surface
Method
(mum)
ASSD
HD (m)m)
(1F)
CD
(SSD)))))
HD (mom)
SIP
DeepCRR
[k]
0.96520.078
0.59:20.005
0.14
0,93820.076
0.567
1137
PlaiNN
0621=0003
0.465:0:044
8.002=0:106
0.37::00003
ADNI
CorticaFlew
0601
0.497=0:049
1.10628.115
0.14::00087
064120.057
0.65200042
0.09
0002100073
Corticallooess
+
0.55=0003
0.400=0.035
0.86600009
0099::00007
0,542=0.054
0.41=
0.76
0009:00022
sortexODE
AV
0.076:20.017
0.214:100200
0,055
0.022:10.012
0.458:0.016
0.192
0,436:10.014
0001:11000111
Vaa2Coriex
0.582:00:228
0.370:0:025
0,746000057
0059:100039
0.977
0.35320.022
0.722:0.055
00431100023
Ours
0.410:0.016
0.136:00002
0.293
0.055:100221
0.213:00.008
0.01
0.155:0.012
0.007110.010
DespCSR
0996--0.085
0.617-00070
1.21--0212
\
0.75=-0.081
0.39::0.067
L151
PlaINN
0.35=00.052
0.460=00.038
0.993=0.082
0..41
OASIS
Corticalllow
26
0.87:10.010
0.095:000477
0.10
0.047:10.086
0.637:0.035
0.062:0.0000
0.992
0.1011000700
CortesiFFoou+
+
0.51=000335
0.39920:030
0812=0.057
0088:00045
0.529=0.033
0.39:20.030
0810:20.055
0086:2002
SOneRODE
SB
0.011
0.282:0:021
0.461=0.062
0.026:20.015
0.465=0.018
0.07
0.435
0018:001000
Voc2Coriex
1P
0.588=0.032 0.381
0.750=0.063
0061
0:581=0.026
0.375=0.027
0.059
004220077
Ours
0.442-00014
0.161=0.012
0.348 0.308-0.025
0.057:+00025
0.218=0.007
0.073
0.159
0.008
achieved the overall best performance due to its explicit regularizations on the deformation trajectory of vertices and its better initialization surface. On the ADNI dataset, our method achieved <488 8% improvement in mean ASSD (of WM and pial surfaces) compared to the second best CortexODE (i.e.. 0.104mm v.s 0.203mm)) with competitive ieff-intersectio ratio (average ~0.021% 1%). On the OASIS dataset, our method achieved similar performance improvement. More quantitative results on the WM surface are reported in the Supplementary Materials. As shown in Fig. the CSR results obtained by our method had uniformly smaller errors across the whole surfaces.
DSspCSR
PiaINN
Coww
2exreuxxxx
CorlesODE
GT
Ensrr
Pial surfhee
03
0.4
WM suffice
Figure 4: Visualization of the reconstructed surfaces with distance compared to ground truth.
Runtime Analysis. It took ~0.8s for our method to obtain a ribbon segmentation map and another ~2s for the topology correction and surface initialization. For the surface reconstruction, the inference time for simultaneously reconstructing three (i.e., WM, pial, and midthickness) surfaces was ~ 11.5.. In comparison, the SOTA explicit CSR method CorteXODE (30) needed ~2s reconstruct two surfaces sequentially. Overall, our method is computationally as efficient as the SOTA alternatives.
4.2 Ablation Studies
Input. Our method takes as input a multi-channel 3D images Tomb and an initialization midthickness surface Sg. We conducted two experiments to analyze the influence of the input components on the CSR performance (Table Top). First, by gradually removing components from Loomb, observed a significant drop in accuracy, indicating the contribution of both the SDF and segmentation maps to the final results. Second, we investigated the impact of the fine structure of the initialization surface by applying Laplacian smoothing on Sy to generate oversmoothed initialization surfaces. The results revealed that the accuracy of the model decreased with the increasing of smoothing steps which resulted in a coarser initialization surface Sg-
Loss functions. We evaluated the contribution of different losses of our method to the surface reconstruction performance in terms of both accuracy (CD, ASSD, HD) and topological correctness (SIF), with the results summarized in Table 3 Middle. Through the ablation studies, we observed
8

---

Table 3: Ablation studies with the surface reconstruction results quantified in terms of CD, ASSD,
HD, and SIF. The setting I0, SO, and T = 5 all refer to our complete setting (cf. Table 2 Top:
The impact of the input image and initialization surface. Middle: The impact of the loss functions.
Bottom: The impact of the deformation steps in DDM.
Input
Initial mesh # of Lap. Sm.
L-Pial Surface
L-WM Surface
Setting
SDF
Seg
CD (mm)
ASSD (mrmn))
HD (trem))
CD (xzamra)
ASSD (mom)
HD (vrr)))
0,110=0.016
0.13620.012
0.29340:006
0.213=000088
000000
0.155220012
0.426=000017
0.167=0.017
0.222==0:011 0.33840:038
0:005500:006
0.16
12
/
0
0.45310.021
0,201 20.026
0.0377
0.250:00.013
0.085-100008
0.184400.016
MI
10
0,416=0.016
0.141=0.013
0,313=000228
0.225=0:010
0.084:0.0007
0.184
M2
/
20
0.429=00018
0.163=0.017
0.30140.040
0.235::00012
0.091+00.009
0.190
Sentin
Zez Laist Lage Cus Lnn
Less
(Dimm))
(3SD)mmm)
P-Pial
Sunface
(mmm)
SIFTS)
CD
(mS))
L-WM
(1mm)
SIFTS)
50
/
0410+0016
+0012 0.277 40000
0035-50.122
0.21=:0.005
0071 +0.MS
01155
4-000
0:12200010
10012
0211500:000
001310000
0.00
52
0412::0016
0.138:1:02
0.300:1:0027
0.037:10.022
0:000
0003:10.007 0.158:10013
0.008:100111
53
/
0.009
0.125+-00012 0.300+-0:027 0.275+-0.100
0.209 0.073+0.007 0.156-+-0013 0.000+00.011
S4
0.01110015
0.129:10011 0.278:10:21 0.522
0.203:00.009 0.069:00.066 0.153:10013 0,0091100012
Number of steps
Surface
Surface
(T)in DDM
CD (mmm)
ASSD (170m))
HD (mm)
SIF(%)
(mm)
ASSD (mm)
(mm)
SIF(%)
0,410:100016
0.136=00012
0.293:0.026
0.035=0:021
0.213=00008
0.071+00000
0.155::0.012
0.007=0.010
0.411::0.015
0.137:0.012
0.026
0.038:0.0222
0,208
0.157:10.012 0.072:10.006
0.008:100011
3
0.411:10.016
0.108110.012
0.207.10.026
0.044_1..025
0.005
0.002.10.006
0.158.10.013
0.00910.0012
2
0.415::00016
0.143=00013
0.309:00028
0.432=0.110
0.204=00000
0.075:t0.005
0.166:+:0012
0.011=000013
that each component played its own role in a complementary way. The first row (referred to as SO) corresponds our complete network setting, while the last row (S4) represents using Chamfer distance alone. The results of setting S4 indicated that the model generated surfaces well matched to the ground truth data at the cost of high topologiaal errors, particularly in highly curved regions, reflected by the results that the SIF ratio was worse on the pial surface than on the WM surface. Enforcing equality of the trajectories from the midthickness surface to the WM and pial surfaces (S3, Ldist) helped optimize the midthickness surface, thereby preventing deformation in an arbitrary direction and reducing self ieffiintersectinn However, the geometric accuracy slightly decreased, which might be caused by the difficulty in accessing highly curved regions or deep sulci under such strong topology constraints. The proposed symmetric cycle loss (S2, Laye) promoted the invertibility of deformations, yielding a significant reduction of self-intersections on the meshes since our method jointly reconstructs both the inner and outer surfaces by deforming the midthickness surface inward and outward with two VFs. Such invertibility also facilitates accurate estimation of the cortical thickness from the trajectory, as illustrated by a sample vertex deformation trajectory in Fig. 2 obtained with settings of S2 and S3. Moreover, the inclusion of regularization terms on the smoothness of SVFs (S1, L..) and surfaces (SO, Lin)) contributed to enhancement in surface quality. Overall, our proposed method struck balance between geometric accuracy and topology quality.
Deformation steps in DDM. Table Bottom shows the evaluation results of different numbers of deformation steps (I') in DDM. As T increased, the performance first improved and then saturated, indicating that five steps were sufficient to deform the midthickness surface to the inner or outer surfaces.
4.3 Reproducibility
We carried out two experiments on two datasets: a paired ADNI 5&3T dataset [22) consisting of 1.5T and 3T MRIS of the same subjects, and the Test-Retest
Table 4: Reproducibility analysis.
L-WM Surface
Method
CD (mum)
ASSD (m7m))
HD (mrm)
dataset 32 comprising 40 MRIs collected within
ADNI-pai pair
Ours
0.520=20055
0.337
0.738=0.151
CortexODE
0.521
0.340=0.000
0,741=0.154
short period for each of the 3 subjects. In these sce-
DeepCSR
0.618:10.105 0.397:0.000
0.823=0.211
narios, the cortical surfaces of the same subject should be nearly identical. Following the experimental setup outlined in 10 13 30), we utilized the iterative closest- point algorithm (ICP) to align image pairs and com- puted the geometric distance between surfaces. The results for the left WM surfaces are presented in Table
FreeSurfer
0.556.10.049 0.36410.054
0.764_0.111
Ours 0.451:100019 0,236-000030
0.492-0.059
CoonexODE
0.457=0.021 0.238+0.031 0.504=00071
DeepCSR
0.5054400047 0,297:00.053 0.610=0.100
FreeSarfer
0.447610.0055 0.253 0.519
4 (more in Supplementary Materials),

---

demonstrating that our method obtained superior reproducibility compared with FreeSurfer and was comparable the SOTA DL methods.
4.4 Cortical Thickness
In contrast to the alternative methods that rely on the ICP algorithm for registering WM and pial faces prior to calculating the Euclidean distance [I0], our proposed method directly provides vertex-wise cortical thickness estimation. To validate the corti- cal thickness estimation, we compared our method with FreeSurfer for estimating the cortical thickness. We identified 200 subjects from the ADNI-2GO [22] dataset (100 are diagnosed with Alzheimer's dis- ease and 100 are normal controls) and computed the average cortical thickness across 35 cortical re- gions based on surface parcellation provided by FreeSurfer [17]. Fig. Fioows the correlation between ours and FreeSurfer's s results, showcasing the effec- tiveness of our proposed framework in accurately capturing cortical thickness.
Cortial
20-
Coefficionts: b-0.00
2/5
3:0
Figure 5: Correlation between prediction (Y- axis) and FreeSurfer's thickness (X-axis) on 35 cortical regions (mmm)
Limitations and Future Directions
Despite achieving improved CSR accuracy and low SIF ratio, our method can be further improved by adopting post-processing methods and new loss functions in order to minimize the SIF ratio and improve surface quality. While focusing on cortical thickness estimation in this paper, recognize the value of incorporating other cortical attributes like surface area and sulci depth into CSR and analysis tasks. These attributes could serve as complementary constraints, enhancing overall performance. It should be noted that further analysis is merited to thoroughly evaluate the proposed method on a large cohort subjects (e.g., subjects in different stages of AD) although we have demonstrated the correlation between the estimated cortical thickness and that of FreeSurfer on balanced dataset.
Conclusion
We introduce DI. framework for cortical surface reconstruction by generating a midthickness surface to initialize coupled reconstruction of both the WM and pial surfaces. Specifically, the midthickness surface is estimated from 3D distance map from each MRI by generating cortical ribbon segmentation map that encodes structural information of the cerebral cortex. The estimated midthickness surface is represented optimized to lie at the center of the inner and outer cortical surfaces and deformed to the inner and outer cortical surfaces by three diffeomorphic flows that are learned jointly with CNNs optimized with a multi-channel input consisting of the brain MRI, the 3D distance map of midthckness surface, and the cortical ribbon segmentation map. Our proposed symmetric cycle loss helps learn diffcomorphic deformation and the numerical solution of DDM improves CSR accuracy and computation efficiency. Experiments on two large-scale neuroimage datasets have demonstrated the superior performance of our method. Moreover, our method generates an estimation of cortical thickness, facilitating statistical analyses of brain atrophy.
triangular mesh with spherical topology, and the mesh is
Acknowledgements
This work was supported in part by the NIH grants AG066650 and EB022573.
References
[1] V. I. Amold. Ordinary differential equations. Springer Science & Business Media, 1992.
10

---

[2] V. Arsigny. Processing data in lie groups: An algebraic approach. Application to non-linear registration and diffusion tensor MRI. PhD thesis, Citeseer, 2004.
[3] V. Arsigny, O. Commowick, X. Pennec, and N. Ayache. A log-euclidean framework for statistics on diffeomorphisms. In International Conference on Medical Image Computing and Computer-Assisted, pages 931. Springer, 2006.
[4] Ashburner. A fast diffeomorphic image registration algorithm. Neuroimage, 38(1):95-113, 2007.
[5] B. B. Avants, C. Epstein. M. Grossman, and J. C. Gee. Symmetric diffeomorphic image reg istration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain. Medical Image Analysis, 12(1):26-41, 2008
[6] G. Balakrishnan, A. Zhao, M. R. Sabuncu, J. Guttag, and A. V. Dalca. Voxelmorph: a learning framework for deformable medical image registration. IEEE Transactions on Medical Imaging, 38(8):1788-1800, 2019.
[7] P.-L. Bazin and D. L. Pham. Topology correction of segmented medical images using a fast marching algorithm. Computer Methods and Programs in Biomedicine, 88(2):182-199 2007,
[8] M. Beg, M. I. Miller, Trouvé, and L. Younes. Computing large deformation metric mappings via geodesic flows of diffeomorphisms. International Journal of Computer Vision, 61:139-157 157, 2005.
[9] M. Bertoux, Lagarde, F. Corlier, I.. Hamelin, J.-F. Mangin, O. Colliot, M. Chupin, M. N. Braskie, P. M. Thompson, M. Bottlaender, et al. Sulcal morphology in alzheimer's s disease: an effective marker of diagnosis and cognition. Newrobiology of Aging, 84:41-48 2019.
[10] Bongratz, A.-M. Rickmann, Pölsterl, and C. Wachinger. Vox2Cortex: Fast explicit reconstruction of cortical surfaces from 3D MRI scans with geometric deep neural networks. In Proceedings of the IEEF/CVF Conference on Computer Vision and Pattern Recognition, pages 20773-20783, 2022.
[11] R. L. Burden, J. D. Faires, and A. M. Burden. Numerical Analysis. Cengage learning, 2015.
[12] R. T. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud. Neural ordinary differential equations. Advances in Neural Information Processing Systems, 31, 2018.
[13] R. S. Cruz, L.. Lebrat, P. Bourgeat, C. Fookes, J. Fripp, and O. Salvado. DeepCSR: deep learning approach for cortical surface reconstruction. In Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, pages 806-815, 2021.
[14] A. M. Dale, B. Fischl, and M. I. Sereno. Cortical surface-based analysis: segmentation and surface reconstruction. Neuroimage, 172)1179-194, 1999.
[15] G. Fahim, K. Amin, and S. Zarif. Single-view 3D reconstruction: A survey of deep learning methods. Computers Graphics, 94:1 190- 2021.
[16] H. Fan, Su, and L. J. Guibas. point set generation network for 3d object reconstruction from single image. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 605-613, 2017.
[17] B. Fischl. Freesurfer. Neuroimage, 62(2):774-781, 2012.
[18] K. Gopinath, C. Desrosiers, and H. Lombaert. SegRecon: Learning joint brain surface re- construction and segmentation from images. In International Conference on Medical Image Computing and Computer Assisted Intervention, pages 650-659. Springer, 2021.
[19] Y. Hong, S. Ahmad, Y. Wu, S. Liu, and P.-T. Yap. Vox2Surf: Implicit surface reconstruction from volumetric data. In International Workshop Machine Learning in Medical Imaging in Conjunction with MICCAI 2021, pages 644-653, Springer, 2021.
[20] A. Hoopes, Iglesias, B. Fischl, Greve, and Dalca. Topofit: Rapid reconstruction of topologically-comen cortical surfaces. In Medical Imaging with Deep Learning, 2022.

---

[21] D. P. Huttenlocher. G. A. Klanderman, and W. J. Rucklidge. Comparing images using the Hausdorff distance. IEEE Transactions on Pattern Analysis and Machine Intelligence, 15(9):850- 863, 1993.
[22] C. R. Jack Jr, M. A. Bernstein, N. C. Fox, P. Thompson, G. Alexander, D. Harvey, B. Borowski, J. Britson, J. L. Whitwell, C. Ward, et al. The alzheimer's s disease neuroimaging initiative (ADNI): MRI methods. Journal of Magnetic Resonance Imaging: An Official Journal of the International Society for Magnetic Resonance in Medicine, 27(4):685-65 691, 2008.
[23] D. P. Kingma and Ba. Adam: A method for stochastic optimization. In Y. Bengio and Y. LeCun, editors, International Conference Learning Representations, 2015.
[24] F. Kong and S. C. Shadden. Learning whole heart mesh generation from patient images for computational simulations. IEEE Transactions on Medical Imaging, 2022.
[25] F. Kong, N. Wilson, and S. Shadden. A deep-learning approach for direct whole-heart mesh reconstruction. Medical Image Analysis, 74: 74:102222, 2021.
[26] L. Lebrat, R. Santa Cruz, F. de Gournay, D. Fu, Bourgeat. Fripp. C. Fookes, and O. Salvado. CorticalFlow: A diffeomorphic mesh transformer network for cortical surface reconstruction. Advances in Neural Information Processing Systems, 34, 2021.
[27] T. Lewiner, H. Lopes, A. W. Vieira, and G. Tavares. Efficient implementation of marching cubes' cases with topological guarantees. Journal of Graphics Tools, 8(2):1-15, 2003.
[28] H. Li, Y. Fan, and A. D. N. Initiative. MDReg-Net Multi-resolution diffeomorphic image registration using fully convolutional networks with deep self-supervision. Human Brain Mapping, 43(7):2218 20022211,,,,, 2022.
(29] B. Ma, Z. Han, Y.S. Liu, and M. Zwicker. Neural-pull: Learning signed distance function from point clouds by learning to pull space onto surface. In International Conference on Machine Learning, pages 7246-7257, 2021.
[30] Q. Ma, I. Li, F. C. Robinson, B. Kainz, D. Rueckert, and A. Alansary. CortexODE: Learning cortical surface reconstruction by neural ODEs. IEEE Transactions on Medical Imaging, 2022.
[31] Q. Ma, E. C. Robinson, B. Kainz, D. Rueckert, and A. Alansary. PialNN: A fast deep learning framework for cortical pial surface reconstruction. In International Workshop on Machine Learning in Clinical Neuroimaging, pages 73-81. Springer, 2021.
[32] J. Maclaren, Z. Han, S. B. Vos, N. Fischbein, and R. Bammer. Reliability of brain volume measurements: test-retest dataset. Scientific Data, 1(1):199, 2014.
[33] D. S. Marcus, T. H. Wang, Parker, J. G. Csernansky, C. Morris, and R. L. Buckner. Open access series of imaging studies (OASIS): cross-sectional mri data in young, middle aged, nondemented, and demented older adults. Journal of Cognitive Neuroscience, 1999):1988-1507, 2007.
[34] L. Mescheder, M. Oechsle, M. Niemeyer, S. Nowozin, and A. Geiger. Occupancy networks: Learning 3D reconstruction in function space. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 4460 4470, 2019.
[35] M. Modat, D. M. Cash, P. Daga, G. P. Winston, J. S. Duncan, and S. Ourselin. Global image registration using symmetric block-matching approach. Journal of Medical Imaging, 10):0240003-02003,, 2014.
[36] T. C. Mok and A. Chung. Fast symmetric diffeomorphic image registration with convolutional neural networks. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 4644-4653, 2020.
[37] D. H. Pak, M. Liu, S. S. Ahn, A. Caballero, J. A. Onofrey, L. Liang, W. Sun, and J. S. Duncan. Weakly supervised deep learning for aortic valve finite element mesh generation from 3D CT images. International Conference on Information Processing Medical Imaging, pages 637-648. Springer, 2021.
12

---

[38] D. H. Pak, M. Liu, T. Kim. L. Liang, R. McKay, W. Sun, and J. S. Duncan. Distortion energy for deep learning-bar based volumetric finite element mesh generation for aortic valves. In International Conference on Medical Image Computing and Computer Assisted Intervention, pages 485-494, Springer, 2021.
[39] J. Park, P. Florence, J. Straub, R. Newcombe, and S. Lovegrove. Deepsdf: Learning continuous the signed distance functions for shape representation. In Proceedings the IEEE/CVÉ Conference on Computer Vision and Pattern Recognition, pages 165-174, 2019.
[40] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al. Pytorch: An imperative style, high-performance deep learning library. Advances in Neural Information Processing Systems, 32, 2019.
[41] J. Ren, Q. Hu, W. Wang, W. Zhang, C. S. Hubbard, P. Zhang, N. An, Y. Zhou, L. Dahmani, D. Wang, et al. Fast cortical surface reconstruction from mri using deep learning. Brain Informatics, 9(1):1-16, 2022.
[42] L. M. Rimol, Nesvåg, D. J. Hagler Jr, Ø. Bergmann, C. Fennema-Notestine C. B. Hartberg, U. K. Haukvik, E. Lange, C. J. Pung, A. Server, et al. Cortical volume, surface area, and thickness schizophrenia and bipolar disorder. Biological Psychiatry, 71(6):552-660, 2012.
[43] J. M. Roe, D. Vidal-Pifeeiro, Ø. Sørensen, A. M. Brandmaier, S. Düzel, H. A. Gonzalez, R. A. Kievit, E. Knights, S. Kühn, U. Lindenberger, et al. Asymmetric thinning of the cerebral cortex across the adult lifespan is accelerated in alzheimer's disease. Nature Communications, 12(1):721, 2021.
[44] O. Ronneberger, P. Fischer, and T. Brox. U-net: Convolutional networks for biomedical image segmentation. In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 234-241,2011
[45] A. G. Roy, S. Conjeti, N. Navab, C. Wachinger, A. D. N. Initiative, et al. Quicknat: A fully convolutional network for quick and accurate segmentation of neuroanatomy. Neurolmage, 186:713-727, 2019.
[46] D. Ruelle and D. Sullivan. Currents, flows and diffeomorphisms. Topology, 14(4):319-377, 1975.
[47] R. Santa Cruz, L. Lebrat, D. Fu, P. Bourgeat, J. Fripp, C. Fookes, and O. Salvado. Corti- calfFow++: Boosting cortical surface reconstruction accuracy, regularity, and interoperability. In International Conference on Medical Image Computing and Computer Assisted Intervention, pages 496-505. Springer, 2022.
[48] D. W. Shattuck and R. M. Leahy. Brainsuite: an automated cortical surface identification tool. Medical Image Analysis, 6(2):129-142, 2002.
[49] A. Sinha, A. Unmesh, Q Huang, and K. Ramani. Surfnet: Generating 3d shape surfaces using deep residual networks. In Proceedings of the IEEE Conference on Computer Vision and Puttern Recognition, pages 6040-6049, 2017.
[50] A. A. Taha and A. Hanbury. Metries for evaluating 3D medical image segmentation: analysis, selection, and tool. BMC Medical Imaging, 15(1):1-28, 2015.
[51] Vercauteren, X. Pennec, A. Perchant, and N. Ayache. Diffeomorphic demons: Efficient non-parametric image registration. Neurolmage, 45(1) 45(1):S61-S72, 2009.
[52] N. Wang, Y. Zhang, Z. Li, Fu, H. Yu, W. Liu, X. Xue, and Y.-G. Jiang. Pixel2Mesh: 3D mesh model generation via image guided deformation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 43110)3300-3311 2020.
[53] C. Wen, Y. Zhang, C. Cao, Z. Li, X. Xue, and Y. Fu. Pixel2Mesh+4: 3D mesh generation and refinement from multi-view images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 21(02) 2180, 2023.
13

---

[54] X. Wen, J. Zhou, Y.S. Liu, H. Su, Z. Dong, and Z. Han. 3d shape reconstruction from 2d images with disentangled attribute flow. In Proceedings of the IEEE/CVF Conference on Computer Vision and Puttern Recognition, pages 3803-3813, 2022.
[55] U. Wickramasinghe, E. Remelli, G. Knott, and P. Fua. Voxel2Mesh: 3D mesh model generation from volumetric data. In International Conference on Medical Image Computing and Computer- Assisted Intervention, pages 299-308. Springer, 2020.
[56] Z. Wu, F. Zhao, F. Wang, W. Lin, L. Wang, D. Shen, and G. Li. Surface-basee analysis of the developing cerebral cortex. In Advances in Magnetic Resonance Technology and Applications, volume 2, pages 287-307. Elsevier, 2021.
[57] H. Xie, H. Yao, S. Zhang, S. Zhou, and W. Sun. Pix2Vox++: Multi-scale context-aware 3d object reconstruction from single and multiple images. International Journal of Computer Vision, 128(12):2919-29))
2020.
[58] S. Yang. M. Xu, Xie, Perry, and Xia. Single-view 3d object reconstruction from shape priors in memory. In Proceedings of the IEEE/CVF Conference Computer Vision and Pattern Recognition, pages 3152-3161, 2021.
[59] H. Zheng, H. Li, and Y. Fan. SuffN: Joint reconstruction of multiple cortical surfaces from magnetic resonance images. In International Symposium Biomedical Imaging, pages 1-4, 2023.
14


Reviews:
[{"review_id": "PvJuLvpx7u", "reviewer": "Reviewer_2ocD", "summary": "This paper develops a new deep-learning-based framework to jointly reconstruct the inner, outer, and their in-between (midthickness) surfaces and estimate cortical thickness. The proposed method has been evaluated on two neuroimaging datasets and shown superior cortical surface reconstruction performance. This work is potentially interesting to the computational neuroscience community and valuable for brain study.", "questions": "See limitations.", "limitations": "Experiments can be performed on more populations (e.g., babies) to evaluate the effectiveness. More explanations are needed to demonstrate the contributions of different losses in the ablation study (Table 3).", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The authors propose a novel deep-learning-based framework to robustly reconstruct the cortical surfaces with topological correctness. Different from existing methods, it explicitly couples the inner and outer surfaces by jointly learning three diffeomorphic flows. Instead of using mixed structures of CNN and GNNs/MLP, this work utilizes a single 3D CNN to leverage 3D brain MRI, a ribbon segmentation map, and a signed distance function altogether. Also, the diffeomorphic deformation trajectory is calculated by a new numerical computation strategy. Furthermore, the initialization of fine midthickness surfaces is effective, and thickness estimation is performed together with cortical reconstruction. Overall, this paper is well-written and easy to follow. The proposed method is novel and can potentially alleviate several limitations of existing methods. The experiments also demonstrate the effectiveness of the proposed method on two neuroimaging datasets.", "weaknesses": "The demographics for the two testing datasets are not illustrated, and experiments can be performed on more populations. \nI am confused about the effectiveness of different losses in Table 3.  \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your comments and suggestions. We will revise our paper accordingly. \n\n**C1: Lack demographics for the two testing datasets. And experiments can be performed on more populations for evaluation.**\n\n**A1**: We will provide detailed demographic information for all datasets used in our experiments. \n\nIn the present study, we used multiple datasets obtained from ADNI-1, ADNI-2&GO, OASIS-1, and a Test-retest dataset. The ADNI-1 dataset consisted of 817 T1-weighted (T1w) brain MRIs from subjects aged 55 to 90, including normal controls (NC), mild cognitive impairment (MCI), and Alzheimer’s disease (AD), the ADNI-Go&2 dataset consisted of 200 T1w scans from subjects aged 60 to 87, including 100 NC and 100 AD subjests, the OASIS-1 dataset consisted of 413 T1w scans from subjects aged 18 to 96 years, including NC and AD subjects, and the Test-retest dataset consisted of 120 T1w scans from three subjects aged 26 to 31. \n\nFurthermore, we have also evaluated our method on HCP dataset with T1w scans from young adult aged 22 to 35. Due to time limit, we evaluated our method on a subset of randomly sampled 300 scans with a splitting of (training: 230, validation: 20, test: 50). The ASSD (mm) measures on L-Pial and L-WM surface reconstruction were 0.156 and 0.073, respectively. \n\nIn summary, our method has been evaluated on MRI scans of subjects aged 18 to 96. \n\nWe will further evaluate our method on scans of younger subjects in our ongoing studies. Particularly, we will evaluate our method on the dHCP dataset [4] that consists of 886 T2w scans from neonatal subjects (20 to 45 weeks). \n\nMore details of the datasets used in our study can be found in their official websites. I will add the links to the supplementary materials.\n\n**C2: More explanations are needed to demonstrate the contributions of different losses in the ablation study (Table 3).**\n\n**A2**: Thanks for the instructive comments.\n\nWe adopted multiple evaluation metrics to evaluate how different losses/components of our method contribute to the surface reconstruction in our experiments. In addition to commonly used metrics, such as Chamfer distance (CD), average symmetric surface distance (ASSD), and the 90th percentile Hausdorff distance (HD) that assess surface **_geometric accuracy_** but are not equipped to quantify the (spherical) topology of reconstructed surfaces, we also adopted the ratio of self-intersection faces (SIF) to measure the surface **_topological quality_** (Line 285). Overall, these complementary metrics enabled to quantify the contributions of different losses in the ablation study in terms of both accuracy and topological correctness. \n\n+ The results of setting S4  (the model regularized with $L_{ch}$ alone) indicated that the model generated surfaces well matched to the ground truth data at the cost of high topological errors, particularly in highly curved regions, reflected by the result that the SIF was much worse on the pial surface than on the WM surface. \n\n+ The results of setting S3 (the model regularized with $L_{ch}+L_{dist}$) indicated that the trajectory loss ($L_{dist}$) helped enforce the midthickness surface to reside between WM and pial surfaces, prevent surface intersection, and mitigate challenges in learning large deformations (Lines 144-149), reflected by the result that SIF was significantly reduced, particularly on the pial surface. A caveat is that the topology constraint might also yield marginally decreased accuracy in highly curved or deep sulci regions (Lines 326-328).  \n\n+ Our method jointly reconstructs both the inner and outer surfaces by deforming the midthickness surface inward and outward with two velocity fields (VFs). The cycle loss ($L_{cyc}$) was effective for promoting invertibility of the VFs, as shown in Fig 3. Such the invertibility also enabled to estimate the cortical thickness from the trajectory. Unlike $L_{cyc}$ working on each vertex, the symmetric similarity loss ($L_{ss}$) directly operated on the two VFs at the voxel level to further ensure the invertibility of these two VFs. The results of S2 ($L_{ch}+L_{dist}+L_{cyc}$) and S1 ($L_{ch}+L_{dist}+L_{cyc}+L_{ss}$) settings showed that surface quality was further improved. Geometric accuracy changes were possibly attributed to the smoothing effect of these loss terms. \n\n+ Given that each vertex was deformed independently through sampled velocity from VFs, the normal consistency loss ($L_{nc}$) discouraged neighboring faces from forming sharp angles. These regularization terms improved both geometric accuracy and surface quality as indicated by the results of S0 setting using all losses above."}, {"review_id": "MZldv274Nc", "reviewer": "Reviewer_7qip", "summary": "This paper proposed a method for cortical surface reconstruction by deforming a middle layer inwards and outwards to be the gray matter and white matter surface respectively. First, segmentation network and topology correction are used to extract the initial surface. Then, three velocity fields are predicted to deform the initial surface to be the final white matter and gray matter. The experiments demonstrate the proposed method and the design. ", "questions": "Figure 1 is confusing to me, For the cyclic loss on the GM as the example, why the forward and backward deformation is from the same velocity field? In Eq. 2, they should use different velocities.  ", "limitations": "Please try to solve the questions and weakness above. ", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The idea of deforming the middle surface to generate white matter and gray matter is interesting. \n2. The experiments are sound and demonstrating, including the comparison between the SOTA methods and the ablation study for the components. \n\n", "weaknesses": "This paper seems to overstate its contribution. First, it claims that it devises an efficient and reliable approach for initialization. However, the initialization step is basically following CortexODE, where the CortexODE predicts the white segmentation, and initializes the WM surface, this paper initializes the midthickness surface. The pipeline of initialization is the same as proposed by CortexODE. Second, the diffeomorphic deformation. This paper claims that they design a new numerical computation strategy to calculate the diffeomorphic deformation trajectory, where the deformation is the accumulated velocity along the several time steps. The CortexODE also utilizes explicit fixed-step methods for ODE discretization instead of the scaling and squaring methods. \n\nThe idea of utilizing the middle state to generate the final result is also not new, a similar approach can be found in [1]. \n\n[1] Mok, Tony CW, and Albert Chung. \"Fast symmetric diffeomorphic image registration with convolutional neural networks", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your comments and suggestions. \n\n**C1: This paper seems to overstate its contribution.\\\nC1-1: First, the initialization pipeline is following/the same as CortexODE (seg + WM surf.).**\n\n**A1-1**: Our objectives and emphases are different from them.\n+ First, our method jointly reconstructs **_both_** the inner and outer cortical surfaces based on an initialization of the midthickness surface, different from CortexODE that relies on the WM segmentation for **_soely_** WM surface initialization and reconstruction. Our experiment results have demonstrated the effectiveness of our initialization method in terms of both surface location and quality (Lines 139-142; Fig. 2(a); Table 3, M1, M2), which has not been demonstrated by prior methods. \n\n+ Second, CortexODE needs to train a secondary model for pial surface reconstruction, entailing an additional step for pial surface initialization. Specifically, they deform the WM surface outward along the normal direction $\\vec{n}$ by a manually selected distance $\\rho$. However, this strategy is _not necessarily optimal_ for all cortical regions because the thickness distribution is not uniform across the entire surface. In contrast, our midthickness surface initialization leverages both WM and GM segmentation maps (Lines 167-169). Each vertex is determined by local WM and GM structures, resulting in a more precise approximation of the true midthickness surface.\n\n**C1-2: Second, the diffeomorphic deformation is the accumulated velocity along several time steps. CortexODE also utilizes explicit fixed-step methods for ODE discretization instead of the scaling and squaring methods.**\n\n**A1-2**: To ensure diffeomorphic deformation (i.e., stronger theoretical guarantees), a natural framework is to formulate it as a flow ODE, modeled as an initial value problem depicted in Eq. (1), and parameterize a set of learnable deformations. Although the problem formulation and high-level idea of numerical solutions (e.g., Euler method) are similar, our concrete solution distinctly differs from that of CortexODE. \n\nCortexODE employs the **_entire_** neural network to represent the continuous dynamics of the vertex trajectory, which is why it is called “neural ODE” [1] based method. The output of the network is computed using a black-box differential equation solver and provides deformation predictions at each time step. The memory cost is $O(1)$ but the time complexity is $O(T)$ (involving $T$ forward propagations), where $T$ is the number of steps. \n\nIn contrast, our approach utilizes a neural network to learn velocity fields (VFs) and bypasses the need for $T$ forward propagations of the network for velocity predictions by leveraging sampling the VFs $T$ times. Additionally, the integration operation is solely on sampled vertices (\\~130K). This is different from typical scaling and squaring methods used in image registration that apply the integration to the whole volume grid (\\~8M, i.e., 55 times larger), and therefore our method is more efficient. \n\n[1] Chen RT, Rubanova Y, Bettencourt J, Duvenaud DK. Neural ordinary differential equations. (NeurIPS 2018)\n\n**C2: Using the middle state to generate the final result is not new, a similar approach can be found in: [1] Mok, Tony CW, and Albert Chung. Fast symmetric diffeomorphic image registration with convolutional neural networks.**\n\n**A2**: We were aware of the referenced work, acknowledged it as [39], and discussed it in our main paper (Lines 119-128). Our method fundamentally distinguishes itself in three key aspects:\n \n+ [1] is designed for 2D image registration, and the deformation field is computed on a pixel grid utilizing scaling and squaring method. A direct extension to 3D images would incur significant computational costs, and the deformation field computed on a voxel grid is not necessarily optimal for cortical surface reconstruction that deforms _vertices_ (i.e., sub-voxel level) instead of _voxels_. \n+ [1] decomposes the transformation from image $Y$ to image $X$ into $\\phi_{YX}^{(1)}(Y)=\\phi_{YX}^{(0.5)}\\circ\\phi_{XY}^{(-0.5)}(Y)$ and estimates the “half” forward deformation $\\phi_{YX}^{(0.5)}$ and the \"half\" inverse deformation $\\phi_{XY}^{(-0.5)}$ by _integrating the velocity fields_ forward and backward respectively. That is, [1] deforms two known images to each other. In contrast, our method (i) computes _vertex-wise deformation_ from velocity fields directly; and (ii) deforms the midthickness surface to inner and outer surfaces and then deforms backwards to midthickness surface (Eq. 2, Fig. 2b). \n+ Their symmetric deformation requires $\\phi_{YX}^{(1)}(Y)\\simeq X$ and $\\phi_{XY}^{(1)}(X)\\simeq Y$ respectively. On the contrary, our cycle loss ensures simultaneous alignment of $p_{\\phi_W\\circ\\phi_G}\\simeq p\\simeq p_{\\phi_G\\circ\\phi_W}$ for each vertex $p$ on the midthickness surface $S_M$.  \n\n**C3: Figure 1 is confusing. Why is the forward and backward deformation from the same velocity field? In Eq. 2, they should use different velocities.**\n\n**A3**: Sorry for the confusion in Fig. 1. Your understanding regarding Eq. 2 was accurate. There was an error on the data flow arrows in Fig. 1. We realized this issue and rectified it in the supplementary materials (Fig.1). We will update Fig. 1 in the main paper.  \n\n"}, {"review_id": "2XnIRVresU", "reviewer": "Reviewer_9EVU", "summary": "This paper presents a deep-learning-based approach for cortical surface reconstruction (CSR) with topology constraints to ensure the physical meaning of reconstructed surfaces. Experiments on public benchmarks show that the proposed approach outperforms current state of the art networks.", "questions": "N/A", "limitations": "The authors did not discuss their method's limitations, or the potential future directions.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The overall presentation is clear and easy to follow.\n\nIn cortical surface reconstruction, the reservation of the reconstructed surface geometry is a crucial and challenging part. This paper aims at tackling the problem by generating a mid-thickness surface with topology constraints, which is intuitive and well-motivated.\n\nThe experimental results seem to be encouraging and the proposed method outperforms SOTAs by a clear margin.", "weaknesses": "The claimed contributions are over-stated, and several related works seem to be ignored and would potentially lead to misunderstanding for the readers.\n\"the cortical thickness estimation is neglected in all existing DL-based CSR methods\" is simply not true. There are plenty of related methods on cortical surface reconstruction/prediction/parcellation that consider cortical attributes (cortical thickness, surface area, cortical folding, diffusivity, structural and functional connectivity, etc) in their frameworks [1,2]. In [3], cortical thickness was already considered as a constraint during the network's supervised learning process, and their network was also designed to predict *both* the surfaces (inner and outer) and cortical thickness at the same time.\nSince all these works are closely related, the authors should at least be aware of and mention or discuss the relations and differences between these methods and their proposed one. \n\n[1] Wu et al. Surface-based analysis of the developing cerebral cortex. Advances in Magnetic Resonance Technology and Applications, 2021.\n[2] Zhao et al. Spherical Deformable U-Net: Application to Cortical Surface Parcellation and Development Prediction. IEEE Transactions on Medical Imaging, 2021.\n[3] Liu et al. Deep modeling of growth trajectories for longitudinal prediction of missing infant cortical surfaces. International Conference on Information Processing in Medical Imaging (IPMI), 2019.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**C1: Contributions over-stated, several related works ignored, and potential for reader misunderstanding.\\\nThe statement \"the cortical thickness estimation is neglected in all existing DL-based CSR methods\" is incorrect. Related methods on cortical surface reconstruction/prediction/parcellation consider cortical attributes (thickness, area, structural and functional connectivity, etc.) in [1,2].\\\nIn [3], cortical thickness was considered as a constraint in network learning, and it also predicts both inner and outer surfaces and thickness concurrently. The authors should acknowledge and discuss the relations and differences between these methods and their proposed one.**\n> [1] Wu et al. Surface-based analysis of the developing cerebral cortex. (AMRTA 2021)\\\n> [2] Zhao et al. Spherical Deformable U-Net: Application to Cortical Surface Parcellation and Development Prediction. (TMI 2021)\\\n> [3] Liu et al. Deep modeling of growth trajectories for longitudinal prediction of missing infant cortical surfaces. (IPMI 2019) \n\n**A1**: Thanks for your suggestions. We will revise the paper as suggested.\n\nFirst, we would like to clarify that the statement (\"the cortical thickness estimation is neglected in all existing **_DL-based CSR_** methods\") intended to depict the methods dedicated to “DL-based cortical surface reconstruction (CSR)” rather than the “entire pipeline of the cortical surface analysis”. We do not focus on non-DL methods or toolboxes/pipelines that enable various tasks in brain MRI analysis in separate steps (i.e., independent solutions for topology correction, cortical thickness/area estimation, and parcellation, etc.). This paper focuses on CSR with one important attribute, cortical thickness. We will further clarify our scope in the main paper. \n\nSecond, we want to elaborate the relations and differences between the mentioned methods and ours. \n\n1. [1] outlines the entire pipeline for analyzing developing cerebral cortex, encompassing the following four steps:\n+ (1-a) Cortical topology correction. A learning-based method is presented to adaptively correct detected topological errors, whereas we use a topology correction algorithm on a signed distance function for initializing the midthickness surface.\n+ (1-b) Cortical surface reconstruction. A two-step non-DL approach is presented to generate the inner and outer cortical surfaces _separately_. Specifically, a tessellation method is used to extract inner surface from the WM segmentation map, and then the inner surface is deformed to the outer surface driven by an external force derived from the Laplace's equation and an internal force to keep the surface tight and smooth, along with some constrains on the spherical topology. The middle surface is then reconstructed as the average of inner and outer surfaces. On the contrary, our DL-based method _jointly_ reconstructs multiple cortical surfaces in a united framework. \n+ (1-c) Cortical surface-based measurement on the reconstructed surfaces. Our main focus is on the cortical surface reconstruction from imaging data.\n+ (1-d) Cortical surface registration, parcellation, modeling, and atlas construction. The rest tasks are out of scope of this work and could be a future direction. \n\n2. [2] introduces a Spherical Deformable U-Net that operates on surface meshes, focusing on two tasks: cortical surface parcellation (a vertex-wise classification/segmentation task) and cortical attribute map development prediction (a vertex-wise dense regression task). In contrast, our method primarily focuses on cortical surface reconstruction and thickness estimation. \n\n3. [3] proposes a method for longitudinal prediction of cortical surfaces using a spatial graph convolutional neural network. Note that their method is not designed for CSR from brain MRIs but instead predicting longitudinal surface changes within the same individual. The initial cortical surfaces are generated by a toolbox as input to their framework. While they do incorporate a thickness constraint using the Euclidean distance between paired vertices on inner and outer surfaces, the predicted thickness is compared with the ground truth thickness by L2 norm. On the contrary, our approach traces the geodesic distance in the vertex deformation trajectory as thickness estimation, without utilizing ground truth thickness as supervision.\n\nIn summary, we will discuss our method’s relatedness to existing works in cortical surface analysis and clarify the focus and scope of our work. \n\n**C2: The authors did not discuss their method's limitations, or the potential future directions.**\n\n**A2**: Thanks. We will add a paragraph discussing our method’s limitations and future directions. For example, \n\n+ Despite achieving a low self-intersection ratio and improved cortical surface reconstruction accuracy, our method can be further improved by adopting and exploring post-processing methods and designing new loss functions to further minimize SIF and improve surface quality.\n\n+ While focusing on cortical thickness estimation in this paper, we recognize the value of incorporating other cortical attributes like surface area and sulci depth into CSR and analysis tasks. These attributes could serve as supplementary constraints, enhancing overall reconstruction quality.\n\n+ We demonstrate the correlation between the estimated cortical thickness and that of FreeSurfer in this paper. We should conduct more analyses on a broader range of subjects (e.g., health subjects and subject with dementia or Alzheimer’s disease) to evaluate the effectiveness of the proposed method and its limitations. \n\n\n"}, {"review_id": "ebJdxdoOdu", "reviewer": "Reviewer_YgxF", "summary": "The paper focuses on the cortical surface-mesh reconstruction in Magnetic Resonance Imaging (MRI) data. The goal is to reconstruct both the inner and outer surfaces of the cortical region. The proposed method begins by estimating an intermediate midthickness layer, which is positioned halfway between the inner and outer surfaces. This midthickness layer serves as the initialized surface mesh for the reconstruction process.\n\nTo evaluate the effectiveness of the proposed method, the authors present results obtained from two neuroimage MRI datasets. These datasets likely consist of MRI scans of the brain, specifically targeting the cortical region. The authors use these datasets to showcase the performance of their surface reconstruction technique.\n\nBy employing their proposed method, the authors aim to demonstrate the ability to accurately reconstruct the cortical surfaces from MRI data. This can be valuable for various applications in neuroscience and neuroimaging research, such as brain mapping, cortical thickness analysis, and studying cortical abnormalities or diseases.", "questions": "- In your framework, after computing the average of WM and GM SDFs, the Marching Cubes algorithm is used for mesh reconstruction. What are the advantages of incorporating ODE?\n\n- What is the average number of vertices in the initialized mesh, and how does it affect the output metrics?\n\n- The caption of Figure 2 (a) does not adequately explain the graphs in the Figure. Please provide a clearer description in the text.\n\n", "limitations": "The paper should thoroughly address the limitations associated with converting the segmentation map into an SDF volume, emphasizing how the resolution of the 3D segmentation volume affects accuracy. Additionally, discussing the potential implications or consequences of these limitations for downstream applications is crucial.\n\nTo provide a comprehensive analysis, it is essential to explicitly mention the failure cases and limitations of the proposed method. Discussing scenarios in which the method might not perform optimally or cases where it is not suitable would enhance the paper's credibility and assist readers in understanding the method's boundaries and potential weaknesses.\n", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper presents a significant advancement in cortical mesh reconstruction and thickness estimation. The key innovation lies in the incorporation of a constraint that ensures the model predicts both the inner and outer meshes while considering a plausible distance representing the cortical thickness between the white matter and gray matter. While previous studies have explored the use of ODE for cortical mesh reconstruction, the proposed method employs a distinct approach by leveraging this technique to obtain two mesh surfaces and incorporates several loss functions to optimize the results.\n\nOverall, the paper introduces a valuable contribution to the field by addressing the challenge of accurately reconstructing cortical meshes and estimating cortical thickness. The novel constraint and the unique application of ODE set this approach apart from previous works.", "weaknesses": "The choice of spatial resolution for applying Marching Cubes on the Signed Distance Function (SDF) volume has a direct impact on the granularity of the resulting mesh, specifically in terms of the number of vertices generated. Increasing the number of vertices can potentially reduce the geometry reconstruction error, but this comes at the expense of higher computational requirements and increased memory consumption. Surprisingly, the paper lacks any explanation or analysis regarding this crucial aspect.\n\nIt would be beneficial for the paper to provide a detailed discussion or analysis concerning the trade-off between spatial resolution, mesh granularity, and the associated computational and memory costs. This would help readers understand the implications of different spatial resolutions on the quality of the reconstructed geometry. Additionally, examining the potential impact on downstream applications or comparing the performance across different resolutions could provide valuable insights for future research and practical implementations.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your comments and suggestions. We will revise our paper accordingly. \n\n**C1: Surface granularity\\\nC1-1: The choice of spatial resolution for applying Marching Cubes on the Signed Distance Function (SDF) volume directly impacts resulting mesh granularity. More vertices can reduce the reconstruction error but increase computation cost. A detailed discussion/analysis on the trade-off between spatial resolution, mesh granularity, and the associated computational costs would clarify implications of different spatial resolutions on the geometry quality.**\n\n**A1-1**: The standard brain MRI image resolution is $256^3$. As a proof of concept, we apply Marching Cubes (MC) on $256^3$ and up-sampled $512^3$ WM SDF. The reconstructed WM surface ASSD is reduced from 0.265mm to 0.254mm, but the computation time increases from 0.23s to 2.45s for MC and from 2.8s to 19.3s for topology correction (i.e., \\~7X totally). The similar trend is observed on the pial surface with lower geometric accuracy. We also increase the mesh granularity by mesh subdivision (\\~1s for a mesh with 130K vertices), but the ASSD barely changes (0.256mm). However, reducing the resolution by down-sampling SDF substantially increases the reconstruction error (especially for the pial surface). In summary, simply applying Marching Cubes to SDF at higher resolutions does bring small improvement on geometric accuracy but incurs heavier computational cost; it may degrade significantly or fail if the surface initialization is obtained based on SDFs at lower resolutions. \n\n**C1-2: What is the avg. number of vertices in the initialized mesh? How does it affect the output metrics? Examining the potential impact on downstream applications could provide valuable insights for future research and practical implementations.**\n\n**A1-2**: Following previous works, we used meshes with \\~130K vertices in our experiments. \n\nTo test the impact of the number of vertices on the downstream tasks, we simplify the meshes by merging neighboring faces (\\~60K) and subdivide the meshes by dividing each face into four faces (\\~520K) and utilize them to train different networks. Due to time limit, we randomly choose a subset from ADNI-1 dataset (300, 20, 80 for training, validation, and test, respectively) to conduct the experiments. The ASSD results are as follows:\n||60K|130K|520K|\n|-|-|-|-|\n|L-WM|0.122|0.074|0.067|\n|L-pial|0.201| 0.156|0.154|\n\n**C2: In your framework, after computing the average of WM & GM SDFs, the Marching Cubes algorithm is used for mesh reconstruction. What are the advantages of incorporating ODE?**\n\n**A2**:\n**_First_**, higher geometric accuracy. Marching Cubes (MC) works on voxel-level but cortical surface reconstruction (CSR) requires sub-voxel accuracy. Simply applying MC yields worse ASSD (mm) (WM: \\~0.265, pial: \\~0.323) than our method (WM: \\~0.071, pial: \\~0.136). \n**_Second_**, better topological quality. Due to the partial volume effect and cerebral cortex's thin and highly-folded pattern (Lines 138-139), the contrast on the gray matter (GM)/cerebrospinal fluid (CSF) interface is very low. Directly applying MC to reconstruct pial surface inevitably generates topology errors (SIF: \\~8%). In contrast, the ODE formulation provides stronger theoretical guarantees for diffeomorphic deformation and we deform an initialization surface with good topology to the pial surface, resulting in better topological quality (SIF: \\~0.02%). \n**_Third_**, by coupling multiple CSR within a united framework, the geodesic trajectory of each vertex serves as an estimation of the cortical thickness. \n\n**C3: Fig. 2a caption does not adequately explain the graphs in the Figure. Please provide a clearer description in the text.**\n\n**A3**: We will add more description. This sub-figure illustrates the surface reconstruction accuracy (Y-axis) with initialization surfaces at different locations (X-axis) for WM and pial surfaces, respectively (using PialNN method). The X-axis represents the initial surface generated from WM SDF ($K_W$), GM SDF ($K_G$), or midthickness SDF ($K_M$). The Y-axis represents the ASSD metric of the reconstructed surface. \n\n\n**C4: Authors should mention the failure cases and limitations of the proposed method, e.g., scenarios where the method might not perform optimally or is not suitable. It would enhance the paper's credibility and help readers to understand the method's boundaries and potential weaknesses.**\n\n**A4**: The standard deviation is very small in our experimental results. We carefully checked the test results and did not find \"visually\" failed cases. Nevertheless, we will visualize results with the worst performance in terms of accuracy and topology correctness in the supplementary file. We postulated that the good generalization performance should be attributed to the fact that our model was trained on large datasets and both the training and testing datasets consisted of high-quality T1 MRI scans of adult subjects.\n\nThere are a few limitations of our method: \n+ Despite achieving a low self-intersection ratio and improved cortical surface reconstruction accuracy, our method can be further improved by adopting and exploring post-processing methods and designing new loss functions to further minimize SIF and improve surface quality.\n+ While focusing on cortical thickness estimation in this paper, we recognize the value of incorporating other cortical attributes like surface area and sulci depth into CSR and analysis tasks. These attributes could serve as supplementary constraints, enhancing overall reconstruction quality.\n+ We demonstrate the correlation between the estimated cortical thickness and that of FreeSurfer in this paper. We should conduct more analyses on a broader range of subjects (e.g., health subjects and subject with dementia or Alzheimer’s disease) to evaluate the effectiveness of the proposed method and its limitations.\n\nWe will add a section to discuss limitations and broader/societal impact.\n"}, {"review_id": "aPZ2nCsYCC", "reviewer": "Reviewer_DGUA", "summary": "This paper proposes a method for reconstructing cortical surface (CSR) (including white matter (WM) and gray matter (GM)) from MRI images. The method first estimates the “midthickness” surface, which is the surface between WM and GM, and learn jointly diffeomorphic flows to deform it to WM and GM, respectively. The proposed model is then evaluated on 2 brain dataset including ADNI and OASIS and compare with other CSR methods.\n", "questions": "- In terms of qualitative results, I check the visualization of every method that the paper compares with, and none of them has equally distributed errors across the surface as illustrated in the paper (Fig 4). I wonder how the visualization is constructed to demonstrate your result. \n- The paper seems to be an improvement work of CortexODE. Instead of using the initial surface that is extracted from WM segmentation, the authors extract the mid surface, which is in between the inner surface and the outer surface as the initial surface. Regarding loss functions, I wonder why the SIF of the L-Pial surface when using only Chamfer loss is so bad (2.522), but SIF of the L-WM surface is still good (0.009). Furthermore, it seems to me that adding auxiliary losses other than Chamfer loss makes the model even worse (Table 3b). Specifically, when using only Chamfer loss (S4), the model outperforms the complete network setting loss (S0) in almost every metric. Instead of using only Chamfer loss, why do the authors propose to use multiple losses but get worse performance?", "limitations": "The authors did not acknowledge any potential negative societal impact of their work.", "rating": 4, "confidence": 5, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- The submission is well-written, the authors adequately demonstrate their ideas in the paper. \n\n- In general, the proposed method provides a good quantitative result compared to other CSR methods.\n", "weaknesses": "- While I appreciate the reported improvement of the proposed method, I do not see a significant novel contribution of this work. Specifically, the idea of using Neural ODE to deform from the initial surface to target surface is widely explored in this task (CortexODE, CorticalFlow). The losses that the authors employed are also commonly applied in this task. Finally, the technical detail proposed in the paper does not seem to have a wide applicability in other areas.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your comments and suggestions.\n\n**C1: Novelty\\\nC1-1: Neural ODE is widely used for deformation in this task (CortexODE, CorticalFlow).**\\\n**A1-1**: While ODE based formation has been widely used for diffeomorphic (image/surface) deformation, our method is different from CortexODE and CorticalFlow in following aspects:\n\n1) Only CortexODE is based on “neural ODE” [1]. It uses the **_entire_** neural network (NN) to represent the continuous dynamics of the vertex trajectory. Deformation predictions are computed by a black-box differential equation solver, involving $T$ forward propagations (i.e., time complexity $O(T)$, T: # of steps). In contrast, our approach bypasses the need for $T$ forward propagations of the network by directly sampling VFs.\n\n2) CorticalFlow uses _several_ NNs to predict a chain of 3 deformations to reconstruct a surface in a coarse-to-fine manner. It requires stage-wise training for each NN and may incur errors in coarse mesh template. In contrast, our method trains _one_ deformation NN on a fine initialization surface.\n\n3) CortexODE & CorticalFlow do not consider physical relations between multiple cortical surfaces and need to train different models to reconstruct the WM and pial surfaces separately. In contrast, our method jointly reconstructs the WM, pial, and midthickness surfaces with cortical thickness estimation in an end-to-end learning framework.\n\n[1] Chen RT et. al. Neural ordinary differential equations (NeurIPS 2018)\n\n**C1-2: The losses are commonly used in this task**\\\n**A1-2**: SOTA CSR methods and their losses have been discussed in Table 1 & Sect. 2. They only used Chamfer, MSE, or edge loss. Our method incorporates multiple new regularization losses, including $L_{dist}$, $L_{cyc}$, and $L_{ss}$, which improved the surface reconstruction quality as demonstrated by the ablation studies (See **A3-2** below for explanation of each loss).\n\n**C1-3: Technical details in the paper seem not to have a wide applicability in other areas**\\\n**A1-3**: 1) The benchmark shows the effectiveness (ASSD \\~0.1mm, SIF \\~0.02%) and efficiency (\\~2s) of our method. It offers a fast alternative for large-scale neuroimage studies (e.g., neurodegenerative diseases). 2) It can be applied to reconstruct other structures from CT/MRI (Lines 76-77, [26,27,40]) and help downstream tasks, such as heart mesh reconstruction for fluid dynamics simulation [2].\n\n[2] Attar R, et. al. 3D cardiac shape prediction with deep neural networks: simultaneous use of images and patient metadata. (MICCAI 2019)\n\n**C1-4: The paper is like an improvement work of CortexODE (initialize midthickness surface)**\\\n**A1-4**: Our method is built upon many pioneering works in image registration and cortical surface reconstruction. In addition to the initialization of a midthickness surface, the joint reconstruction of both WM and GM surfaces under multiple regularization terms improved the robustness, accuracy, and computational efficiency of cortical surface reconstruction, as demonstrated by the extensive experimental results.\n\nPlease refer to my response \"_A1-1 to Reviewer 7qip_\" for detailed comparison on surface initialization, and \"_A1-1 above_\" for ODE solution.\n\n**C2: Other papers’ visualization does not have equally distributed errors. How is Fig. 4 constructed to demonstrate results?**\\\n**A2**: 1) The Nilearn package was used to map geometric errors onto surfaces. Others might use different tools with better shading and lighting. 2) We used more scales of color bar and thus more details were rendered. 3) When zooming in, we can observe that errors were not evenly distributed.\n\nIn all, the lighter color in the error map signifies our method's lower geometric error. We will improve the visualization to highlight differences of the results obtained by different methods.\n\n**C3: Loss function\\\nC3-1: S4, why is SIF of L-Pial bad but SIF of L-WM good?**\\\n**A3-1**: We postulated that this was caused by the fact that pial surfaces are more convoluted than WM surfaces and contrast on the GM/CSF interface is worse than that of WM/GM interface. When the model was trained with the Chamfer loss alone, the VF might not be smooth enough, incurring self-intersection.\n\n**C3-2: Table 3b, using auxiliary losses gets worse performance. Why use them?**\\\n**A3-2**: _Topology correctness and accuracy are equally important_ in the cortical surface reconstruction. To obtain topologically correct surfaces, we may need to sacrifice the accuracy.\n+ Chamfer loss $L_{ch}$ only considers the coordinates of vertices with no topology constrains, resulting in topological errors in highly curved regions (setting S4).\n+ Trajectory loss $L_{dist}$ helps enforce the midthickness surface to reside between WM and pial surfaces, avoid surface intersection, and mitigate challenges in learning large deformations (Lines 144-149). SIF was significantly reduced, particularly on the complex pial surface (setting S3).\n+ Cycle loss $L_{cyc}$ and symmetric similarity loss $L_{ss}$ promote invertibility of the VFs to couple multiple CSR and enable cortical thickness estimation from the trajectory. S2 & S1 settings showed that surface quality further improved. Geometric accuracy changes were possibly attributed to the smoothing effect of these loss terms.\n+ Normal consistency loss $L_{nc}$ discourages neighboring faces from forming sharp angles as vertex is deformed independently through sampled velocity from VFs. S0 setting showed that these regularization terms improved both geometric accuracy and surface quality. \n\n**C4: Any potential negative societal impact of the work?**\\\n**A4**: We use 2 renowned public datasets collected with standards that control sex, age, and health status etc. But there may be oversights in terms of ethnicity and socioeconomic status, potentially putting underrepresented groups at a disadvantage in the study. We will discuss the potential impact and the need for more diverse datasets for validation."}, {"review_id": "Zqgvwr2rKj", "reviewer": "Reviewer_vxoq", "summary": "This paper is solving a solid practical problem, handling the cortical thickness estimation problem that was historically ignored. Moreover, it improves the quality of surface reconstruction. Authors observed that deforming an over-smoothed template results in lower accuracy if the template is relatively further away from the target. Therefore, reconstructing the midthickness layer could help improve the accuracy.", "questions": "(1) Can you detail Figure1, so that it’s more informative? For example, the variables in the equations could appear in the figure at corresponding locations.\n\n(2) Can you split into shorter paragraphs, where each paragraph conveys one observation / insight?\n", "limitations": "Yes", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "(1) The problem to be solved is well-motivated, limitations of existing methods are analyzed.\n\n(2) Figures are nicely made.\n\n(3) Storytelling is good.\n", "weaknesses": "(1) Rather than integrating all the equations within the text, it would enhance clarity if you present them as numbered equations. Furthermore, providing explanations for each variable immediately after equations would help.\n\n(2) Paragraphs are too long. Most subsections are single-paragraphed, and all contents are presented within that one paragraph.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your comments and suggestions. We will revise our paper accordingly. \n\n**C1: Rather than integrating all the equations within the text, it will enhance clarity if you present them as numbered equations. Furthermore, providing explanations for each variable immediately after equations would help.**\n\n**A1**: Due to the space limit, we put equations in text in Section. 3.3. Camera-ready version allows one additional content page, we can present them as numbered equation and provide explanation accordingly if our paper gets accepted.  \n\n\n**C2: Paragraphs are too long. Most subsections are single-paragraphed, and all contents are presented within that one paragraph. Can you split into shorter paragraphs, where each paragraph conveys one observation/insight?**\n\n**A2**: Thanks. If more space allowed, we would definitely split long paragraphs into shorter ones. \n\n\n**C3: Can you detail Figure1, so that it’s more informative? For example, the variables in the equations could appear in the figure at corresponding locations.**\n\n**A3**: We will add more description in the caption and notation in the figure. Note that there was an error in the data flow in the main paper and we rectified it in the supplementary material (Fig. 1). We will update the Fig. 1 in the main paper. \n\nFig. 1(a) illustrates the midthickness initialization process. An FCN utilizes raw images to generate white matter (WM) and gray matter (GM) segmentation maps. From these, a signed distance function (SDF) is derived, followed by topology correction and Marching Cubes algorithms, and then used for initializing the midthickness surface. \n\nFig. 1(b) outlines the cortical thickness reconstruction network. It takes a composite input of raw images, SDF, and segmentation maps, simultaneously learning three parallel velocity fields (VFs): $\\symbf{v}_M$, $\\symbf{v}_W$, and $\\symbf{v}_G$. The $\\symbf{v}_M$ is used to generate deformation field for optimizing the midthickness surface $S_M$. The $\\symbf{v}_M$ (or $\\symbf{v}_G$) is used to deform $S_M$ to the WM surface $S_W$ (or pial surface $S_G$) and then deform $S_W$ (or $S_G$) to midthickness surface $S_M'$ (or $S_M''$) to constrain topology between WM and pial surfaces. \n"}]

Author rebuttal:

Meta-review: In this paper, the authors propose a new method for cortical surface reconstruction. Most reviewers find the proposed method interesting and novel, with convincing experimental results. After a careful evaluation of the reviews and rebuttal, the AC agrees with these assessment and recommends acceptance.
Decision: Accept (poster)