# Hybrid-CSR: Coupling Explicit and Implicit Shape Representation
for Cortical Surface Reconstruction

Shanlin Sun, Thanh-Tung Le
University of California, Irvine
{shanlins<EMAIL_ADDRESS>Chenyu You
Yale University
<EMAIL_ADDRESS>Hao Tang
Meta
<EMAIL_ADDRESS>Kun Han, Haoyu Ma, Deying Kong, Xiangyi Yan, Xiaohui Xie
University of California, Irvine
{khan7, haoyum3, deyingk, xiangyy4<EMAIL_ADDRESS>
###### Abstract

We present Hybrid-CSR, a geometric deep-learning model that combines explicit
and implicit shape representations for cortical surface reconstruction.
Specifically, Hybrid-CSR begins with explicit deformations of template meshes
to obtain coarsely reconstructed cortical surfaces, based on which the
oriented point clouds are estimated for the subsequent differentiable poisson
surface reconstruction. By doing so, our method unifies explicit (oriented
point clouds) and implicit (indicator function) cortical surface
reconstruction. Compared to explicit representation-based methods, our hybrid
approach is more friendly to capture detailed structures, and when compared
with implicit representation-based methods, our method can be topology aware
because of end-to-end training with a mesh-based deformation module. In order
to address topology defects, we propose a new topology correction pipeline
that relies on optimization-based diffeomorphic surface registration.
Experimental results on three brain datasets show that our approach surpasses
existing implicit and explicit cortical surface reconstruction methods in
numeric metrics in terms of accuracy, regularity, and consistency.

##  1 Introduction

Figure 1: Contour Representation Toy Example. A polygon target in (a) can be
represented by explicitly deforming a circle via optimizing chamfer distance,
which gives the orange contour in (d). In contrast, our hybrid method can
generate more accurate reconstructed contour in (g) by coupling explicit
(oriented point clouds) and implicit shape representation (indicator grid).

Cortical surface reconstruction is the task to extract both inner and outer
surfaces of the cerebral cortex from brain MRI scans, with the inner surface
situated between the cortical gray matter (GM) and white matter (WM) while the
outer surface between the cerebrospinal fluid and the cortical gray matter.
Accurate and detailed reconstruction of the cortical surface can be used to
facilitate brain mapping [23, 27, 68, 26], identify biomarkers for
neurological disorders[30, 22, 2], and enable pre-surgical planning [21].
However, extracting accurate cortical surface with genus of 0 is still very
challenging due to structural complexity of the cerebral cortex and partial
volume effect (PVE) [7] in medical imaging.

Traditional methods [17, 24, 25, 37] typically first segment the volumetric
structures and then extract the cortical meshes via Marching Cubes [46, 44] or
level set-based approaches [52, 45]. Specifically, FreeSurfer [24], the
current standard for cortical surface reconstruction, obtains WM surfaces by
applying mesh tessellation to segmented WM. With the help of convolutional
neural network and a novel spherical embedding, FastSurfer [32] achieves
faster and better brain segmentation. While deep learning-based volumetric
segmentation models have demonstrated remarkable performance in medical
imaging [60, 34, 31, 12, 73], the reconstructed cortical meshes may not
accurately delineate tissue boundaries due to PVE problem. Furthermore, the
predicted brain segmentation may contain topological defects, necessitating
time-consuming topology correction algorithms to be applied for genus-0
reconstructed meshes.

One way to overcome the partial volume effect is to utilize deep implicit
functions, which can represent complex shapes with fine-grained details by
representing the surface implicitly as the zero level-set of continuous
functions [53, 61, 50, 64, 67]. Recently, DeepCSR has been proposed to
reconstruct the cortical surface leveraging the deep implicit functions in
[16]. As deep implicit functions can be trained efficiently using simple L1 or
L2 loss functions, the optimization process is straightforward and robust in
terms of convergence. However, DeepCSR still requires a time-consuming
spherical level-set evolution algorithm [56] to remove holes and handles from
the reconstructed meshes, due to the absence of topology awareness.

Alternatively, there is another branch of work that exploit explicit topology
to resolve topology defects and partial volume effects simultaneously [8, 42,
62, 33]. By progressively deforming explicit meshes from template meshes,
explicit cortical reconstruction methods can guarantee reconstructed meshes
inherit the desired topology, thus avoiding time-consuming post-processing. At
the same time, sub-voxel variations of cerebral cortex surfaces can also be
captured. Nonetheless, these methods suffer from two main problems. Firstly,
the often used Chamfer distance loss to align predicted and target meshes
tends to get trapped in local minima easily, failing to distinguish bad
samples from the true one [1, 51]. Although a weighted Chamfer distance has
been proposed to prioritize fitting local regions with high curvatures in
Vox2Cortex [8], the issue is only alleviated but not resolved completely,
still resulting in suboptimal assignments between two sets of points [57, 39].
Secondly, while explicit regularizations can reduce self-intersections, they
also lead to lower geometric accuracy. Neural Mesh Flow [29] tries to avoid
explicit regularization by deforming template meshes through continuous
diffeomorphic flow, however, it does not perform well in modeling sharp and
large deformations.

In this work, we propose a novel approach for cortical surface reconstruction
called Hybrid-CSR, which integrates both implicit and explicit shape
representations. The method involves an initial step to deform template meshes
into a coarsely reconstructed cortical surface, based on which oriented point
clouds are estimated for the subsequent implicit cortical surface
reconstruction. Specifically, we apply differentiable Poisson surface
reconstruction [55] to bridge oriented point clouds (explicit) to indicator
grids (implicit), from which watertight cortical meshes can be extracted via
Marching Cubes. In addition, we propose to apply optimization-based
diffeomorphic surface registration to realize topology correction. Our Hybrid-
CSR offers several advantages over existing techniques. Unlike voxel-based
methods, our approach is not susceptible to partial volume effects. Compared
to deep implicit function-based method, Hybrid-CSR produces fewer topology
defects and demonstrates superior efficiency in inference. Finally, compared
to mesh-based methods, our hybrid approach offers a higher level of
expressiveness and flexibility.

The efficacy of our proposed method can be demonstrated by using a simple
contour representation toy example, as shown in Fig. 1. The target contour
(Fig. 1(a)) is a polygon that can be either reconstructed through the
deformation of a source contour (circle) or the extraction of the zero level
set from an indicator field. As can be seen in Fig. 1(d), the results obtained
from explicit contour deformations optimization are not optimal due to the
“regularizer’s dilemma” [29]. Instead, our proposed hybrid method optimizes
the positions and normals of the oriented points (initialized by Fig. 1(e))
and minimizes the difference between the indicator map reconstructed from the
ground truth and the optimized oriented point clouds (Fig. 1(c)(i)). Hybrid
shape representation can produce a well-captured contour from the predicted
indicator field, as demonstrated in Fig. 1(g). More details about this toy
example can be found in the supplementary.

In summary, our main contributions are:

  * •

We propose Hybrid-CSR, the first cortical surface reconstruction framework
coupling the explicit and implicit shape representation based on
differentiable Poisson surface reconstruction.

  * •

We propose a new topology correction pipeline based on optimization-based
diffeomorphic surface registration.

  * •

We demonstrate on multiple brain datasets that Hybrid-CSR surpasses implicit
and explicit reconstruction methods in terms of accuracy, regularity as well
as consistency.

##  2 Related Work

Traditional cortical surface reconstruction was accomplished through a
sequence of image-processing steps, with FreeSurfer [24], being a widely used
approach. While accurate, these methods are constrained by time limitations,
with each case taking up to 7-8 hours to complete, thus hindering practical
application in clinical settings. Deep learning methods, therefore, are
proposed to address the time limitation of traditional approaches and show
potential improvement for cortical surface reconstruction task.

##### Implicit representation

FastSurfer[32] utilizes eigenfunctions of the Laplace-Beltrami operator to
parametrize the surface and generate the final spherical map by scaling the 3D
spectral embedding vector to unit length. Furthermore, recent research in 3D
computer vision has focused on deep implicit representations [50, 64, 63, 53],
which have shown great potential in improving the accuracy and efficiency of
surface reconstruction. For instance, SegRecon [28] employs a 3D CNN for
simultaneously learning segmentation and surface reconstruction by utilizing
3D signed distance function. Additionally, DeepCSR[16] and CortexODE [47]
leverage deep implicit functions to represent the surface, which is then
extracted using Marching Cubes [46, 44] to produce the output mesh.

##### Explicit representation

Explicit representation methods provide an alternative approach to cortical
surface reconstruction, which often learns deformation networks to directly
transform a source mesh into a target mesh. For example, CorticalFlow [43]
implements a flow Ordinary Differential Equation (ODE) framework to learn to
deform a reference template towards a targeted object. Other methods, such as
Voxel2Mesh [70] and Vox2Cortex [8], learn deformable mesh models that take as
input a template mesh or sphere initialization and iteratively deform the mesh
by learning deformation field of the vertices.

##### Diffeomorphic Transformation

A diffeomorphism is an invertible mapping where the forward and backward
transformations are smooth. It is widely used in medical registration
problems. Many works [18, 6, 40, 19, 5] usually assume the velocity field is
stationary and defined in the grid space [4] so that they can apply scaling
and squaring method [3] to do fast integration. CorticalFlow and
CorticalFlow++ [62] learn a discrete stationary velocity field and integrate
the deformations associated with the template meshes using interpolation and
traditional ODE solvers. Recently, with the power of neural ordinary
differential equation solver [14, 13], optimizing a neural diffeomorphic flow
efficiently became possible. NDF [65] addresses organ shape representation and
registration simultaneously by decomposing the implicit shape representation
into continuous diffeomorphic transformations and template shape
representation. Neural Mesh Flow[29] focuses on generating manifold mesh from
images or point clouds via conditional continuous diffeomorphic flow.

##  3 Method

Figure 2: The overview of Hybrid-CSR. The architecture takes as input a 3D MRI
and template meshes and predicts a voxel-wise segmentation map and cortical
surface meshes. This illustration presents the first two steps of our whole
pipeline: coarse mesh deformation and oriented point cloud estimation. The
coarse mesh deformation module offsets the vertices of template meshes twice
using GCN based on image and shape descriptive features. The oriented point
cloud estimation module predicts the positions and normals sequentially with
GCN-based feature encoder and GLU-based final layer.

In this section, we will provide a comprehensive description of our proposed
method, Hybrid-CSR, which aims to couple explicit and implicit surface
representations to enable accurate cortical surface reconstruction from 3D
MRI. The pipeline of Hybrid-CSR consists of three main parts, including coarse
mesh deformation (Sec. 3.2), oriented point cloud estimation (Sec. 3.3),
topology correction (Sec. 3.4) and surface refinement (Sec. 3.5). Prior to
delving into the specifics of our method, we will first describe the methods
of reconstructing surfaces from oriented point clouds using differentiable
Poisson surface reconstruction (Sec. 3.1.1) and modelling continuous
deformations using diffeomorphic flow (Sec. 3.1.2).

###  3.1 Preliminary Knowledge

####  3.1.1 Poisson Surface Reconstruction

Poisson surface reconstruction (PSR) [36] aims to recover an indicator
function $\chi\in\mathbb{R}^{n}$ from K sampled points
$\mathcal{P}=\left\\{p_{i}\in\mathbb{R}^{3}\right\\}_{i=1}^{K}$ with normals
$\mathcal{N}=\left\\{n_{i}\in\mathbb{R}^{3}\right\\}_{i=1}^{K}$, by satisfying
that $\chi$ changes sharply between positive and negative values at the
surface boundary along the direction orthogonal to the surface. We consider
the case where $n:=r\times r\times r$, and $d=3$, where $r$ is the resolution
of the indicator grid. In practice, PSR first constructs a point normal field
$\mathbf{q}\in\mathbb{R}^{n\times d}$ from $\mathcal{P}$ and $\mathcal{N}$.
Then, it formulates the reconstruction of $\chi$ as a Poisson equation:
$\nabla^{2}\chi:=\nabla\cdot\nabla\chi=\nabla\cdot\mathbf{q}$, where the
Laplacian of $\chi$ is equal to the divergence of the normal vector field
$\mathbf{v}$, subject to the boundary condition that $\chi$ is zero at
infinity. This is equivalent to minimizing a quadratic energy function such
that $\min_{\chi}\|\nabla\chi-\mathbf{q}\|_{2}^{2}$.

Unlike [36], which encodes the indicator function $\chi$ as a linear
combination of sparse basis functions and solves the partial differential
equation (PDE) using a finite element solver on an octree, Differentiable
Poisson Surface Reconstruction (DPSR) [55] represents $\chi$ on a 3D grid in a
discrete Fourier basis and employs a spectral solver [11]. The unnormalized
indicator function $\chi^{\prime}$ is given by

$\tilde{\chi}=\tilde{g}_{\sigma,r}(\mathbf{u})\odot\frac{i\mathbf{u}\cdot\tilde{\mathbf{q}}}{-2\pi\|\mathbf{u}\|^{2}},\quad\chi^{\prime}=\operatorname{IFFT}(\tilde{\chi})$
(1)

where spectral domain signal is denoted as tilde symbol, i.e.,
$\tilde{\chi}=FFT(\chi)$, $\mathbf{u}\in\mathbb{R}^{n\times d}$ denotes the
spectral frequencies, $\tilde{\mathbf{q}}$ represents the fast Fourier
transform (FFT) of $\mathbf{q}$, IFFT($\tilde{\chi}$) represents the inverse
FFT of $\tilde{\chi}$, and $\tilde{g}_{\sigma,r}(\mathbf{u})$ is a Gaussian
smoothing kernel of bandwidth $\sigma$ at grid resolution $r$ in the spectral
domain. We denote the element-wise product as
$\odot:\mathbb{R}^{n}\times\mathbb{R}^{n}\mapsto\mathbb{R}^{n}$, and the dot
product $(\cdot):\mathbb{R}^{n\times d}\times\mathbb{R}^{n\times
d}\mapsto\mathbb{R}^{n}$, and the L2-norm as
$\|\cdot\|^{2}:\mathbb{R}^{n\times d}\mapsto\mathbb{R}^{n}$ . Finally, the
normalized indicator function is obtained by subtracting the mean of the
unnormalized indicator function at $\mathcal{P}_{up}$ and re-scaling it,
written as

$\chi=\frac{m}{\operatorname{abs}\left(\left.\chi^{\prime}\right|_{\mathbf{x}=0}\right)}\left(\chi^{\prime}-\frac{1}{|\\{\mathcal{P}\\}|}\sum_{\mathbf{c}\in\\{\mathcal{P}\\}}{\chi^{\prime}|_{\mathbf{x}=\mathbf{c}}}\right).$
(2)

####  3.1.2 Diffeomorphic Flow

Diffeomorphic flows can establish dense point correspondences between source
and target surfaces, and preserve the desired geometric topology. In our
methods, it is used for mesh registration in topology correction pipeline. Let
$\Phi(\boldsymbol{p},t):\Omega\subset\mathbb{R}^{3}\times[0,1]\mapsto\Omega\subset\mathbb{R}^{3}$
define a continuous, invertible trajectory from the initial position
$\boldsymbol{p}=\Phi(\boldsymbol{p},0)$ to the final position
$\boldsymbol{p}^{\prime}=\Phi(\boldsymbol{p},1)$, satisfying such ordinary
differential equation (ODE) and the initial condition:

$\frac{\partial\Phi(\boldsymbol{p},t)}{\partial
t}=\boldsymbol{v}(\Phi(\boldsymbol{p},t),t)\quad\text{ s.t.
}\quad\Phi(\boldsymbol{p},0)=\boldsymbol{p},$ (3)

where $\boldsymbol{v}(\boldsymbol{p},t):\Omega\times[0,1]\mapsto\Omega$
indicates the velocity vector of coordinate $\boldsymbol{p}$ at time t. If
$\boldsymbol{v}$ is Lipschitz continuous, a solution to Eq. 3 exists and is
unique in the interval $[0,1]$, which ensures that any two deformation
trajectories do not cross each other [15].

###  3.2 Coarse Mesh Deformation

As is shown in Fig. 2, the coarse mesh deformation module takes template
meshes $\mathcal{M}_{T}=\left(\mathcal{V}_{T},\mathcal{E}_{T}\right)$ as input
and deforms them into
$\mathcal{M}_{c}=\left(\mathcal{V}_{c},\mathcal{E}_{T}\right)$ given the image
features extracted from the voxel-based neural network. Our coarse mesh
deformation module follows the design from Vox2Cortex [8], using the same mesh
templates, residual Unet-based segmentation branch [72, 34], GCN-based mesh
deformation branch [59, 9], feature sampling, as well as training losses.
However, we deform the template meshes only twice, in contrast to the four
steps in Vox2Cortex. Coarse mesh deformation enables a good initialization for
the subsequent oriented point cloud estimation, resulting in fewer outlier
points. In our experiments, we utilize smaller templates ($\approx 42000$
vertices per surface) VoxCortex provides with their official implementation
for both training and inference.

###  3.3 Oriented Point Cloud Estimation

The oriented point cloud estimation module takes the coarsely deformed meshes
$\mathcal{M}_{c}$ as input and predicts the coordinates and normals of an
upsampled oriented point cloud $\mathcal{O}_{up}$ for Poisson surface
reconstruction. This is achieved through encoding features on the graph
(mesh), followed by the estimation of point positions and normals in sequence.

##### Feature Encoding

Same as Vox2Cortex, we applied a GCN-based residual block to encode point
features based on the image features. The residual block is composed of three
GCNs with subsequent batch normalization layers and ReLU activations. The
input residuum is added before the last ReLU output and reshaped with nearest-
neighbor interpolation if needed. Let $\mathbf{f}_{i}\in\mathbb{R}^{d_{out}}$
define the encoded feature of a vertex $\boldsymbol{v_{i}}\in\mathcal{V}_{c}$.

##### Position Estimation

Our method predicts multiple offsets associated with each vertex to upsample
the point cloud and improve shape representation accuracy. Instead of using
GCNs, which tend to generate over-smoothed offsets, we employ the gated linear
unit (GLU) [20] to estimate the displacements
$\boldsymbol{d_{i}}\in\mathbb{R}^{(S\times 3)}$ with an upsample scale of $S$
for each vertex. This can be expressed as:

$\boldsymbol{d_{i}}=(\mathbf{W}_{0}\mathbf{f}_{i}+\mathbf{b}_{0})\odot\sigma(\mathbf{W}_{1}\mathbf{f}_{i}+\mathbf{b}_{1})$
(4)

where $\mathbf{W}_{0},\mathbf{W}_{1}\in\mathbb{R}^{(S\times 3)\times d_{out}}$
together with $\mathbf{b}_{0},\mathbf{b}_{1}\in\mathbb{R}^{(S\times 3)}$
represent linear projections and $\sigma$ represents sigmoid function. The
upsampled displacements are then added to the vertex $\boldsymbol{v}_{i}$ to
obtain the point cloud position $\boldsymbol{p}_{i}^{up}\in\mathbb{R}^{S\times
3}$. In our experiments, $S$ is set to be 7.

##### Normal Estimation

We begin by extracting multi-scale image features for each point in
$\mathcal{P}_{up}$. Features from points that have been displaced from the
same vertex in $\mathcal{M}_{c}$ are stacked together, and point features
$\mathbf{f}_{i}^{up}$ are learned through a residual GCN layer. Based on
$\mathbf{f}_{i}^{up}$ and $\boldsymbol{p}_{i}^{up}$, the normals
$\boldsymbol{n}_{i}^{up}\in\mathbb{R}^{S\times 3}$ are predicted by a GLU,
formulated as

$\begin{split}\boldsymbol{n}_{i}^{up}=(\mathbf{W}_{0}^{up}[\mathbf{f}_{i}^{up},\boldsymbol{p}_{i}^{up}]+\mathbf{b}_{0}^{up})\odot\\\
\sigma(\mathbf{W}_{1}^{up}[\mathbf{f}_{i}^{up},\boldsymbol{p}_{i}^{up}]+\mathbf{b}_{1}^{up})\end{split}$
(5)

where $\mathbf{W}_{0}^{up},\mathbf{W}_{1}^{up}\in\mathbb{R}^{(S\times
3)\times(d_{out}+S\times 3)}$ together with
$\mathbf{b}_{0},\mathbf{b}_{1}\in\mathbb{R}^{(S\times 3)}$ represent linear
projections and $\sigma$ represents sigmoid function. Now, the upsampled
oriented point cloud
$\mathcal{O}_{up}=(\mathcal{P}_{up}=\\{\boldsymbol{p}_{i}^{up}\\}_{i=1}^{K},~{}\mathcal{N}_{up}=\\{\boldsymbol{n}_{i}^{up}\\}_{i=1}^{K})$
is achieved.

##### Surface Reconstruction

In practice, the Differentiable Poisson Surface Reconstruction (DPSR) method
acquires a uniformly discretized point normal field $\mathbf{v}$ by
rasterizing the predicted oriented point cloud normals $\mathcal{N}_{up}$ onto
four uniformly sampled 3D grids. Utilizing Eq. 2, an $r^{3}$ indicator grid
$\hat{\chi}$ can be obtained for each cortical structure, and subsequently
used to reconstruct cortical meshes
$\hat{\mathcal{M}}=\left(\hat{\mathcal{V}},\hat{\mathcal{E}}\right)$ via
Marching Cubes. The smoothness and resolution of the indicator grid
$\hat{\chi}$ can be modulated by the values of $\sigma$ and $r$. In our
experiments, we set $\sigma$ to 2 and $r$ to 256 for both training and
inference.

##### Training

For each cortical surface, we acquire the ground truth indicator grid $\chi$
by running DPSR on densely sampled oriented point clouds sampled from the
(pseudo) ground truth meshes. Next, we use a Sobel filter [41] to detect edges
in $\chi$. The resulting edge map is then smoothed using a Gaussian filter
with a kernel size of 7 and a standard deviation of 1. Finally, we apply a
weighted mean square error (wMSE) to measure the difference between predicted
and ground truth indicator grids, especially along the boundary regions, and
train Hybrid-CSR end-to-end, given as

$\mathcal{L}_{\text{DPSR}}=\|w_{edge}\odot\left(\hat{\chi}-\chi\right)\|^{2}$
(6)

where $w_{edge}\in\mathbb{R}^{n}$ is the smoothed edge map of the ground truth
indicator grid. More implementation details are included in the supplementary.

###  3.4 Topology correction

We propose a pipeline to fix the topological defects present in
$\hat{\mathcal{M}}$. The WM meshes are prone to having holes while the pial
meshes tend to have handles, which prevent the reconstructed meshes from being
genus-0. Four cortical meshes in $\hat{\mathcal{M}}$ can be post-processed in
parallel. For simplicity, we will note $\hat{\mathcal{M}}$ as a single mesh
(either WM or pial) in this section.

(a) WM

(b) Pial

Figure 3: Illustration of Topology Correction and Mesh Registration. Blue
contours represent $\hat{\mathcal{M}}$, red contours represent
$\hat{\mathcal{M}}_{tc^{\prime}}$, and dotted red contour represent
$\hat{\mathcal{M}_{tc}}$.

We first binarize the predicted indicator grids $\hat{\chi}$ into masks,
selecting the largest connected component that contains the majority of
voxels. The binary masks are then transformed into signed distance grids [10]
and smoothed with a Gaussian filter having a standard deviation of 1. In Fig.
3, we demonstrate how we transform $\hat{\mathcal{M}}$ into topologically
correct
$\hat{\mathcal{M}}_{tc^{\prime}}=(\hat{\mathcal{V}}_{tc^{\prime}},\hat{\mathcal{E}}_{tc})$.
To fill the holes in WM meshes (Fig. 3(a)), we extract the $\tau_{wm}$-level
set of the smoothed signed distance grid, where $\tau_{wm}=0.5$. Similarly, to
address the handles in pial meshes (Fig. 3(b)), we extract the
$\tau_{pial}$-level set of the smoothed singed distance grid, where
$\tau_{pial}=-1.8$.

Then, we apply optimization-based diffeomorphic registration to map
$\hat{\mathcal{M}}_{tc^{\prime}}$ to $\hat{\mathcal{M}}$, obtaining the
accurate, topologically correct cortical reconstruction results, denoted as
$\hat{\mathcal{M}}_{tc}=(\hat{\mathcal{V}}_{tc},\hat{\mathcal{E}}_{tc})$, as
indicated by the dotted red contours in Fig. 3. To do so, we assume
${\boldsymbol{v}}$ in Eq. 3 is stationary and can be modeled via a neural
field [71, 66], such that
$\mathcal{F}_{\theta}(\boldsymbol{p})=[\boldsymbol{v}_{p_{x}},\boldsymbol{v}_{p_{y}},\boldsymbol{v}_{p_{z}}]^{T}$,
where $\boldsymbol{p}$ and $\theta$ denote the network 3D coordinates and
parameters. Then, the initial value problem (IVP) in Eq. 3 can be solved with
a Differentiable ODE Solver (NODE) [14] whose dynamic function is set to be
$\mathcal{F}_{\theta}$.

In the forward pass, the destination position
$\hat{\boldsymbol{p}}_{tc}\in\hat{\mathcal{V}}_{tc}$ starting from
$\hat{\boldsymbol{p}}_{tc^{\prime}}\in\hat{\mathcal{V}}_{tc^{\prime}}$ is
estimated by integrating $\mathcal{F}_{\theta}(\boldsymbol{p})$ from $t=0$ to
$t=1$ via 4th-order Runge–Kutta methods with step size being 0.2. For
backpropagation, NODE adopts the adjoint sensitivity method [58], which
retrieves the gradient by solving the adjoint ODE backwards in time and allows
solving with O(1) memory usage no matter how many steps the ODE solver takes.
The network parameters $\theta$ are optimized with chamfer distance using the
Adam optimizer for 75 iterations with a step size of 3e-4. To compute the
chamfer distance, we sample 150,000 points with normals from both
$\hat{\mathcal{M}}_{tc}$ and $\hat{\mathcal{M}}_{tc^{\prime}}$ per iteration.

###  3.5 Learning-based Surface Refinement

A learning-based diffeomorphic surface registration model, similar to
CortexODE [47], is utilized to generate refined cortical meshes
$\hat{\mathcal{M}}_{f}=(\hat{\mathcal{V}}_{f},\hat{\mathcal{E}}_{tc})$ given
$\hat{M}_{tc}$. It will remove some artifacts introduced during the topology
correction procedures and improve the alignment of the cortical surfaces with
the original 3D MRI data. Different from the CortexODE, our refined pial
surfaces $\hat{\mathcal{M}}_{f}^{pial}$ are transformed from
$\hat{\mathcal{M}}_{tc}^{pial}$ instead of $\hat{\mathcal{M}}_{tc}^{wm}$.

##  4 Experiments

###  4.1 Datasets

To evaluate the performance of our method on the reconstruction of cortical
surfaces from MRI images, we use three publicly available datasets: the
Alzheimer’s Disease Neuroimaging Initiative (ADNI) dataset [35], the OASIS-1
dataset [49], and the test-retest (TRT) dataset [48]. We obtain the pseudo-
ground truth surfaces generated from Freesurfer v5.3 [24] for all three
datasets. We strictly follow pre-processing pipeline from [8]. Specifically,
we first register the MRIs to the MNI152 scan. After padding the input images
to have shape $192\times 208\times 192$, we resize them to $128\times
144\times 128$. The intensity values are min-max-normalized to the range
$[0,1]$.

##### ADNI

We use a subset of the ADNI dataset [35] containing a total of $419$
T1-weighted (T1w) brain MRI from subjects aged from $55$ to $90$ years old. We
stratify the dataset into $299$ scans for training ($\approx 70\%$), $40$
scans for validation($\approx 10\%$), and $80$ scans for testing ($\approx
20\%$). We report all of our experiment results on the test set.

##### OASIS

For the OASIS dataset [49], we use all of $416$ T1-weighted (T1w) brain MRI
images. We stratify the dataset into $292$ scans for training ($\approx
70\%$), $44$ scans for validation ($\approx 10\%$), and $80$ scans for testing
($\approx 20\%$). We report all of our experiment results on the test set.

##### Test-retest

To analyze the consistency of our approach, we evaluate all $120$ scans from
three different subjects, where each subject is scanned twice in $20$ days.

###  4.2 Performance Comparison

####  4.2.1 Competing Methods

For all competing methods, we train their models with their official
implementations using their suggested experimental settings and pick the best
checkpoints on the validation set for comparisons.

##### DeepCSR [16]

is an implicit surface-based that directly predicts implicit surface
representations for the coordinates in the MRI images. The surface can be
reconstructed using either occupancy field or signed distance function. We
reproduce the method in both ways and observe that signed distance function
yields better performance.

##### Vox2Cortex [8]

is a deformation-based model proposed to retrieve cortical surfaces by
deforming a generic template and employing a joint graph neural network and
U-Net. The results shown below are generated based on the higher-resolution
templates with $\approx 168,000$ vertices for each structure. Our framework is
built upon Vox2Cortex [8], and therefore it serves as our baseline model.

##### CorticalFlow++ [43]

is a diffeomorphic-based method that successively deforms template meshes by
integrating a discrete stationary velocity grid using traditional ODE solvers.

##### CortexODE [47]

is a multi-stage deformation approach that obtains the white matter initial
surface from a volumetric segmentation model, then leverages neural ordinary
differential equations to deform an initial WM surface into refined WM and
pial surfaces by learning a diffeomorphic flow.

####  4.2.2 Evaluation Metrics

We evaluate our method as well as other approaches using three commonly used
metrics including average symmetric surface distance (ASSD), normal
consistency (NC), and self-intersection faces ratio (SI). To calculate ASSD
and NC, we sample $100$K points uniformly from both predicted and target
meshes. For measuring regularity, we determine SI faces using PyMesh [74]
library.

####  4.2.3 Results Discussion

##### Accuracy

Table 1:  Cortical Surface Reconstruction Performance Comparison in terms of
average symmetric surface distance (ASSD), normal consistency (NC), and self-
intersection face ratio (SI) on ADNI and OASIS datasets. Best values are
highlighted. ASSD results are in mm. All results are listed in the format
“mean value $\pm$ standard deviation”. Hybrid-CSR represents the cortical
reconstruction results $\hat{M}$ without any topology correction operation.
+TC and +TC+Refine indicate the performance of reconstructed cortical surfaces
$\hat{M}_{tc}$ and $\hat{M}_{f}$ separately. While $\downarrow$ means smaller
metric value is better, $\uparrow$ indicates larger metric value is better.

(a) ANDI dataset

| Left Pial | Left WM | Right Pial | Right WM
---|---|---|---|---
Method | ASSD (mm) $\downarrow$ | NC $\uparrow$ | SI (%) $\downarrow$ | ASSD (mm) $\downarrow$ | NC $\uparrow$ | SI (%) $\downarrow$ | ASSD (mm) $\downarrow$ | NC $\uparrow$ | SI (%) $\downarrow$ | ASSD (mm) $\downarrow$ | NC $\uparrow$ | SI (%) $\downarrow$
DeepCSR [16] | .368 $\pm.082$ | .908 $\pm.015$ | 0 | .390 $\pm.162$ | .934 $\pm.016$ | 0 | .394 $\pm.083$ | .914 $\pm.012$ | 0 | .388 $\pm.172$ | .936 $\pm.014$ | 0
Vox2Cortex [8] | .339 $\pm.055$ | .918 $\pm.010$ | .741 $\pm.221$ | .346 $\pm.073$ | .926 $\pm.011$ | .719 $\pm.214$ | .350 $\pm.037$ | .915 $\pm.009$ | 1.025 $\pm.237$ | .335 $\pm.061$ | .927 $\pm.010$ | .745 $\pm.199$
CorticalFlow++ [42] | .296 $\pm.079$ | .925 $\pm.011$ | .164 $\pm.093$ | .271 $\pm.071$ | .936 $\pm.009$ | .058 $\pm.032$ | .270 $\pm.044$ | .924 $\pm.010$ | .187 $\pm.104$ | .268 $\pm.073$ | .933 $\pm.009$ | .067 $\pm.032$
CortexODE [47] | .258 $\pm.073$ | .929 $\pm.010$ | .112 $\pm.072$ |  .234 $\pm.064$ | .938 $\pm.010$ | .013 $\pm.011$ | .214 $\pm.035$ | .927 $\pm.009$ | .173 $\pm.091$ |  .231 $\pm.052$ | .939 $\pm.009$ | .004 $\pm.005$
Hybrid-CSR (ours) | .254 $\pm.054$ | .886 $\pm.012$ | 0 | .264 $\pm.055$ | .899 $\pm.010$ | 0 | .250 $\pm.041$ | .886 $\pm.012$ | 0 | .257 $\pm.045$ | .901 $\pm.010$ | 0
\+ TC | .267 $\pm.056$ | .926 $\pm.010$ | 0 | .268 $\pm.056$ |  .939 $\pm.010$ | 0 | .262 $\pm.042$ | .926 $\pm.010$ | 0 | .260 $\pm.046$ | .940 $\pm.010$ | 0
\+ TC + Refine |  .203 $\pm.049$ |  .935 $\pm.009$ | .090 $\pm.075$ | .244 $\pm.056$ |  .939 $\pm.010$ | .042 $\pm.025$ |  .200 $\pm.034$ |  .935 $\pm.008$ | .080 $\pm.071$ | .240 $\pm.043$ |  .941 $\pm.010$ | .011 $\pm.010$

(b) OASIS dataset

| Left Pial | Left WM | Right Pial | Right WM
---|---|---|---|---
Method | ASSD (mm) $\downarrow$ | NC $\uparrow$ | SI (%) $\downarrow$ | ASSD (mm) $\downarrow$ | NC $\uparrow$ | SI (%) $\downarrow$ | ASSD (mm) $\downarrow$ | NC $\uparrow$ | SI (%) $\downarrow$ | ASSD (mm) $\downarrow$ | NC $\uparrow$ | SI (%) $\downarrow$
DeepCSR | .424 $\pm.075$ | .898 $\pm.016$ | 0 | .312 $\pm.124$ | .941 $\pm.010$ | 0 | .444 $\pm.087$ | .895 $\pm.018$ | 0 | .344 $\pm.158$ | .941 $\pm.011$ | 0
Vox2Cortex | .401 $\pm.040$ | .900 $\pm.012$ | 1.110 $\pm.270$ | .302 $\pm.037$ | .928 $\pm.008$ | .994 $\pm.193$ | .405 $\pm.044$ | .898 $\pm.012$ | 1.321 $\pm.252$ | .303 $\pm.042$ | .929 $\pm.009$ | 1.022 $\pm.186$
CorticalFlow++ | .326 $\pm.058$ | .913 $\pm.011$ | .147 $\pm.100$ | .225 $\pm.038$ | .937 $\pm.007$ | .054 $\pm.060$ | .318 $\pm.057$ | .913 $\pm.011$ | .192 $\pm.123$ | .227 $\pm.046$ | .935 $\pm.008$ | .076 $\pm.068$
CortexODE | .279 $\pm.052$ | .919 $\pm.009$ | .277 $\pm.096$ |  .183 $\pm.036$ |  .943 $\pm.007$ | .032 $\pm.025$ | .280 $\pm.052$ | .918 $\pm.010$ | .151 $\pm.060$ |  .182 $\pm.052$ |  .943 $\pm.008$ | .022 $\pm.020$
Hybrid-CSR (ours) | .298 $\pm.045$ | .879 $\pm.013$ | 0 | .220 $\pm.039$ | .905 $\pm.008$ | 0 | .301 $\pm.049$ | .880 $\pm.013$ | 0 | .218 $\pm.047$ | .907 $\pm.009$ | 0
\+ TC | .311 $\pm.045$ | .915 $\pm.011$ | 0 | .220 $\pm.039$ | .941 $\pm.008$ | 0 | .313 $\pm.049$ | .915 $\pm.012$ | 0 | .219 $\pm.046$ | .942 $\pm.009$ | 0
\+ TC + Refine |  .274 $\pm.043$ |  .921 $\pm.010$ | .034 $\pm.022$ | .198 $\pm.035$ |  .943 $\pm.008$ | .040 $\pm.022$ |  .275 $\pm.052$ |  .920 $\pm.010$ | .029 $\pm.018$ | .199 $\pm.033$ |  .943 $\pm.008$ | .037 $\pm.026$

As shown in Tab. 1, our proposed method surpasses other competing methods in
terms of surface reconstruction accuracy (ASSD and NC). Compared with any
implicit-based and explicit-based methods including DeepCSR, Vox2Cortex, and
CorticalFlow++, our hybrid approach without refinement can generate
significantly more accurate results on all cortical surfaces. Furthermore, we
got comparable performance with CortexODE on both ADNI and OASIS datasets.
While CortexODE performs slightly better than us on in terms of ASSD on WM
surfaces, we excel them in terms of ASSD on pial surfaces and NC on all
surfaces. We observe that WM surfaces in CortexODE have a higher performance
than our Hybrid-CSR because it is based on volumetric segmentation models
which are typically suitable for extracting structures with details but few
ambiguities in topology, i.e., white matter. Our method, on the other hand,
performs better in more diverse shapes.

##### Regularity

To assess the regularity of Hybrid-CSR, we calculate the percentage of self-
intersecting faces and present them in Tab. 1. Thanks to the non-self-
intersection property of Marching Cubes, we guarantee to achieve $0$ self-
intersection, same as DeepCSR [16]. Although a small number of self-
intersecting faces are introduced after surface refinement, we still
significantly surpass Vox2Cortex, CorticalFlow++, and CortexODE on all
surfaces, except for CortexODE [47] on WM surfaces.

##### Consistency

Table 2:  Cortical Surface Reconstruction Consistency Comparison in terms of
ASSD on TRT dataset. Hybrid-CSR here indicates the result going through the
whole post processing procedures.

Method | ASSD (mm) |  $>1$mm |  $>2$mm
---|---|---|---
Vox2Cortex | .301 $\pm.171$ | 3.12% | .69%
CortexODE | .280 $\pm.164$ | 2.21% | .38%
FreeSurfer | .301 $\pm.176$ | 3.35% | .94%
Hybrid-CSR + TC + Refine |  .258 $\pm.051$ |  1.61% |  .18%

We evaluate the consistency of Hybrid-CSR (with topology correction and
surface refinement), Vox2Cortex [8], and CortexODE [47] (which are all trained
on OASIS), and Freesurfer on TRT dataset. We generate cortical surfaces from
MRI images of the same subject on the same day and measure the ASSD of the
resulting reconstructions. The brain morphology of two consecutive scans taken
on the same day should be similar to each other, except for the variations
caused by the imaging process. The result from Table 2 shows that we
outperform Vox2Cortex [8], CortexODE [47], and Freesurfer on the consistency
aspect.

###  4.3 Results Visualization

Figure 4: Hybrid-CSR Results Visualization. $\mathcal{M}_{gt}$ is the ground
truth. $\hat{\mathcal{M}}$ refers to the meshes reconstructed from poisson
surface reconstruction. $\hat{\mathcal{M}}_{tc}$ represents the meshes going
through topology correction. And $\hat{\mathcal{M}}_{f}$ indicates the refined
meshes. Colors indicate point-to-surface distances, with darker colors
indicating larger prediction errors.

The cortical surface reconstruction results generated by Hybrid-CSR are shown
in Fig. 4, where $\hat{\mathcal{M}}$ denotes the surfaces reconstructed from
PSR that are already well-aligned to the ground truth. After topology
correction, $\hat{\mathcal{M}}_{tc}$ is smoother and topologically correct,
and the surface refinement procedure further completes the details in the
cortical reconstruction, as shown in $\hat{\mathcal{M}}_{f}$. More
visualization of Hybrid-CSR with competing methods would be presented in
supplementary material.

###  4.4 Ablation Study

We conduct two ablation experiments of WM surface reconstruction on OASIS
validation set. For fair comparisons, all experiments below are based on the
same coarse mesh deformation module pre-trained for 50 iterations, then the
oriented point cloud module will be finetuned for another 50 iterations using
different settings. The evaluation metrics are ASSD and 90-percentile
Hausdorff distance (HD90).

Table 3:  Ablation Study on Upsample Ratio $\boldsymbol{S}$.

|  | Left WM | Right WM
---|---|---|---
Method | # params | ASSD (mm) | HD90 (mm) | ASSD (mm) | HD90 (mm)
$S=1$ | 6.49M | .357 $\pm.034$ | .881 $\pm.109$ | .352 $\pm.039$ | .867 $\pm.114$
$S=4$ | 6.54M | .323 $\pm.036$ | .727 $\pm.105$ | .317 $\pm.042$ | .712 $\pm.110$
$S=7$ | 6.59M | .314 $\pm.036$ | .692 $\pm.107$ | .309 $\pm.044$ | .678 $\pm.116$
$S=10$ | 6.64M |  .312 $\pm.037$ |  .688 $\pm.107$ |  .308 $\pm.043$ |  .669 $\pm.113$

Table 4:  Ablation Study on Network Type

Network | Left WM | Right WM
---|---|---
Type | ASSD | HD90 | ASSD | HD90
GCN | .314 $\pm.036$ | .692 $\pm.107$ | .309 $\pm.044$ | .678 $\pm.116$
GLU |  .304 $\pm.036$ |  .661 $\pm.108$ |  .298 $\pm.044$ |  .649 $\pm.114$

##### Hybrid Representation and Upsample Ratio

In this study, we employed GCN to estimate point positions and normals for
Hybrid-CSR with varying upsample ratios denoted by $S$. As shown in Table 3,
increasing the upsample ratio results in more accurate reconstruction of the
white matter surface, while not significantly increasing the number of
parameters that need optimization. We have chosen the value of $S=7$ in our
experiments to strike a balance between training efficiency and accuracy of
the reconstructed surfaces.

##### GLU Estimation

In this study, we apply two different networks, i.e., GCN and GLU, to estimate
the oriented point clouds. As can be observed from Tab 4, GLU can generate
significantly more accurate reconstruction results on WM surfaces. GCN tends
to generate over-smoothing predictions, while GLU can produce more acute
deformations. Since PSR is robust to outlier points, GLU might be an
advantageous choice. More discussions on GLU-based point cloud estimation are
included in the supplementary.

##  5 Limitations and Further Directions

Our proposed method for cortical surface reconstruction has some limitations
that need to be considered. One issue is the running time, as our method uses
neural fields-based diffeomorphic surface registration to correct topology,
which takes about 2 minutes for 75-iteration optimizations. This makes it
slower than methods [42, 62, 8, 33] relying on explicit surface
reconstruction. To address this issue, a persistent homology prior can be
imposed on the reconstructed indicator field to further eliminate the need for
topology correction procedures in our proposed method. Additionally, memory
consumption can be reduced by handling anatomical structures separately or
sequentially. CorticalFlow++ has shown that reconstructing the gray matter
surface starting from the white matter surface can be beneficial. These
improvements will help us achieve even more accurate and efficient cortical
surface reconstruction.

Another limitation lies in the data sources, where we used pseudo labels
generated from Freesurfer without manually removing faulty labels from our
training and test set. This may lead to inaccurate evaluation results and
potentially overfitting to some noisy labels. Additionally, our method is not
designed to handle MRI scans with significant domain shifts. That is, Hybrid-
CSR may generate undesired results when applied to brain MRI scans with tumors
or other morphological changes.

##  6 Conclusion

This paper introduces Hybrid-CSR, a novel cortical surface reconstruction
framework that leverages both explicit and implicit shape representations
through differentiable Poisson surface reconstruction. The proposed method
includes a new topology correction pipeline that employs continuous
diffeomorphic flow optimization. The effectiveness of Hybrid-CSR is evaluated
on multiple brain datasets, and it shows competitive performance in terms of
accuracy, regularity and consistency. Moreover, the versatility of our method
extends beyond the cortical surface reconstruction problem and has the
potential for wide application in a variety of surface reconstruction
scenarios.

## References

  * [1] Panos Achlioptas, Olga Diamanti, Ioannis Mitliagkas, and Leonidas Guibas.  Learning representations and generative models for 3d point clouds.  In International conference on machine learning, pages 40–49. PMLR, 2018.
  * [2] Liana G Apostolova, Paul M Thompson, Steve A Rogers, Ivo D Dinov, Charleen Zoumalan, Calen A Steiner, Erin Siu, Amity E Green, Gary W Small, Arthur W Toga, et al.  Surface feature-guided mapping of cerebral metabolic changes in cognitively normal and mildly impaired elderly.  Molecular imaging and biology, 12:218–224, 2010.
  * [3] Vincent Arsigny, Olivier Commowick, Xavier Pennec, and Nicholas Ayache.  A log-euclidean framework for statistics on diffeomorphisms.  In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2006: 9th International Conference, Copenhagen, Denmark, October 1-6, 2006. Proceedings, Part I 9, pages 924–931. Springer, 2006.
  * [4] John Ashburner.  A fast diffeomorphic image registration algorithm.  Neuroimage, 38(1):95–113, 2007.
  * [5] Brian B Avants, Charles L Epstein, Murray Grossman, and James C Gee.  Symmetric diffeomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain.  Medical image analysis, 12(1):26–41, 2008.
  * [6] Guha Balakrishnan, Amy Zhao, Mert R Sabuncu, John Guttag, and Adrian V Dalca.  Voxelmorph: a learning framework for deformable medical image registration.  IEEE transactions on medical imaging, 38(8):1788–1800, 2019.
  * [7] Miguel Ángel González Ballester, Andrew P Zisserman, and Michael Brady.  Estimation of the partial volume effect in mri.  Medical image analysis, 6(4):389–405, 2002.
  * [8] Fabian Bongratz, Anne-Marie Rickmann, Sebastian Pölsterl, and Christian Wachinger.  Vox2cortex: fast explicit reconstruction of cortical surfaces from 3d mri scans with geometric deep neural networks.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 20773–20783, 2022.
  * [9] Michael M Bronstein, Joan Bruna, Taco Cohen, and Petar Veličković.  Geometric deep learning: Grids, groups, graphs, geodesics, and gauges.  arXiv preprint arXiv:2104.13478, 2021.
  * [10] M Akmal Butt and Petros Maragos.  Optimum design of chamfer distance transforms.  IEEE Transactions on Image Processing, 7(10):1477–1484, 1998.
  * [11] Claudio Canuto, M Yousuff Hussaini, Alfio Quarteroni, and Thomas A Zang.  Spectral methods: fundamentals in single domains.  Springer Science & Business Media, 2007.
  * [12] Hu Cao, Yueyue Wang, Joy Chen, Dongsheng Jiang, Xiaopeng Zhang, Qi Tian, and Manning Wang.  Swin-unet: Unet-like pure transformer for medical image segmentation.  In Computer Vision–ECCV 2022 Workshops: Tel Aviv, Israel, October 23–27, 2022, Proceedings, Part III, pages 205–218. Springer, 2023.
  * [13] Ricky TQ Chen, Brandon Amos, and Maximilian Nickel.  Learning neural event functions for ordinary differential equations.  arXiv preprint arXiv:2011.03902, 2020.
  * [14] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud.  Neural ordinary differential equations.  Advances in neural information processing systems, 31, 2018.
  * [15] Earl A Coddington and Norman Levinson.  Theory of ordinary differential equations.  Pure & Applied Mathematics S. McGraw-Hill Education, London, England, Mar. 1984.
  * [16] Rodrigo Santa Cruz, Leo Lebrat, Pierrick Bourgeat, Clinton Fookes, Jurgen Fripp, and Olivier Salvado.  Deepcsr: A 3d deep learning approach for cortical surface reconstruction.  In Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, pages 806–815, 2021.
  * [17] Robert Dahnke, Rachel Aine Yotter, and Christian Gaser.  Cortical thickness and central surface estimation.  Neuroimage, 65:336–348, 2013.
  * [18] Adrian V Dalca, Guha Balakrishnan, John Guttag, and Mert R Sabuncu.  Unsupervised learning for fast probabilistic diffeomorphic registration.  In Medical Image Computing and Computer Assisted Intervention–MICCAI 2018: 21st International Conference, Granada, Spain, September 16-20, 2018, Proceedings, Part I, pages 729–738. Springer, 2018.
  * [19] Adrian V Dalca, Guha Balakrishnan, John Guttag, and Mert R Sabuncu.  Unsupervised learning of probabilistic diffeomorphic registration for images and surfaces.  Medical image analysis, 57:226–236, 2019.
  * [20] Yann N Dauphin, Angela Fan, Michael Auli, and David Grangier.  Language modeling with gated convolutional networks.  In International conference on machine learning, pages 933–941. PMLR, 2017.
  * [21] Alessandro De Benedictis and Hugues Duffau.  Brain hodotopy: from esoteric concept to practical surgical applications.  Neurosurgery, 68(6):1703–1723, 2011.
  * [22] Rahul S Desikan, Howard J Cabral, Fabio Settecase, Christopher P Hess, William P Dillon, Christine M Glastonbury, Michael W Weiner, Nicholas J Schmansky, David H Salat, Bruce Fischl, et al.  Automated mri measures predict progression to alzheimer’s disease.  Neurobiology of aging, 31(8):1364–1374, 2010.
  * [23] Christophe Destrieux, Bruce Fischl, Anders Dale, and Eric Halgren.  Automatic parcellation of human cortical gyri and sulci using standard anatomical nomenclature.  Neuroimage, 53(1):1–15, 2010.
  * [24] Bruce Fischl.  Freesurfer.  Neuroimage, 62(2):774–781, 2012.
  * [25] Bruce Fischl, Martin I Sereno, and Anders M Dale.  Cortical surface-based analysis: Ii: inflation, flattening, and a surface-based coordinate system.  Neuroimage, 9(2):195–207, 1999.
  * [26] Bruce Fischl, Martin I Sereno, Roger BH Tootell, and Anders M Dale.  High-resolution intersubject averaging and a coordinate system for the cortical surface.  Human brain mapping, 8(4):272–284, 1999.
  * [27] Matthew F Glasser, Timothy S Coalson, Emma C Robinson, Carl D Hacker, John Harwell, Essa Yacoub, Kamil Ugurbil, Jesper Andersson, Christian F Beckmann, Mark Jenkinson, et al.  A multi-modal parcellation of human cerebral cortex.  Nature, 536(7615):171–178, 2016.
  * [28] Karthik Gopinath, Christian Desrosiers, and Herve Lombaert.  Segrecon: Learning joint brain surface reconstruction and segmentation from images.  In Medical Image Computing and Computer Assisted Intervention–MICCAI 2021: 24th International Conference, Strasbourg, France, September 27–October 1, 2021, Proceedings, Part VII 24, pages 650–659. Springer, 2021.
  * [29] Kunal Gupta.  Neural mesh flow: 3d manifold mesh generation via diffeomorphic flows.  University of California, San Diego, 2020.
  * [30] Xiao Han, Jorge Jovicich, David Salat, Andre van der Kouwe, Brian Quinn, Silvester Czanner, Evelina Busa, Jenni Pacheco, Marilyn Albert, Ronald Killiany, et al.  Reliability of mri-derived measurements of human cerebral cortical thickness: the effects of field strength, scanner upgrade and manufacturer.  Neuroimage, 32(1):180–194, 2006.
  * [31] Ali Hatamizadeh, Yucheng Tang, Vishwesh Nath, Dong Yang, Andriy Myronenko, Bennett Landman, Holger R Roth, and Daguang Xu.  Unetr: Transformers for 3d medical image segmentation.  In Proceedings of the IEEE/CVF winter conference on applications of computer vision, pages 574–584, 2022.
  * [32] Leonie Henschel, Sailesh Conjeti, Santiago Estrada, Kersten Diers, Bruce Fischl, and Martin Reuter.  Fastsurfer-a fast and accurate deep learning based neuroimaging pipeline.  NeuroImage, 219:117012, 2020.
  * [33] Andrew Hoopes, Juan Eugenio Iglesias, Bruce Fischl, Douglas Greve, and Adrian V Dalca.  Topofit: Rapid reconstruction of topologically-correct cortical surfaces.  In Medical Imaging with Deep Learning, 2021.
  * [34] Fabian Isensee, Jens Petersen, Andre Klein, David Zimmerer, Paul F Jaeger, Simon Kohl, Jakob Wasserthal, Gregor Koehler, Tobias Norajitra, Sebastian Wirkert, et al.  nnu-net: Self-adapting framework for u-net-based medical image segmentation.  arXiv preprint arXiv:1809.10486, 2018.
  * [35] Clifford R Jack Jr, Matt A Bernstein, Nick C Fox, Paul Thompson, Gene Alexander, Danielle Harvey, Bret Borowski, Paula J Britson, Jennifer L. Whitwell, Chadwick Ward, et al.  The alzheimer’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–691, 2008.
  * [36] Michael Kazhdan, Matthew Bolitho, and Hugues Hoppe.  Poisson surface reconstruction.  In Proceedings of the fourth Eurographics symposium on Geometry processing, volume 7, page 0, 2006.
  * [37] June Sic Kim, Vivek Singh, Jun Ki Lee, Jason Lerch, Yasser Ad-Dab’bagh, David MacDonald, Jong Min Lee, Sun I Kim, and Alan C Evans.  Automated 3-d extraction and evaluation of the inner and outer cortical surfaces using a laplacian map and partial volume effect classification.  Neuroimage, 27(1):210–221, 2005.
  * [38] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [39] Soheil Kolouri, Gustavo K Rohde, and Heiko Hoffmann.  Sliced wasserstein distance for learning gaussian mixture models.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 3427–3436, 2018.
  * [40] Julian Krebs, Hervé Delingette, Boris Mailhé, Nicholas Ayache, and Tommaso Mansi.  Learning a probabilistic model for diffeomorphic registration.  IEEE transactions on medical imaging, 38(9):2165–2176, 2019.
  * [41] D Kroon.  Numerical optimization of kernel based image derivatives.  Short Paper University Twente, 3, 2009.
  * [42] Leo Lebrat, Rodrigo Santa Cruz, Frederic de Gournay, Darren Fu, Pierrick Bourgeat, Jurgen Fripp, Clinton Fookes, and Olivier Salvado.  Corticalflow: a diffeomorphic mesh transformer network for cortical surface reconstruction.  Advances in Neural Information Processing Systems, 34:29491–29505, 2021.
  * [43] Léo Lebrat, Rodrigo Santa Cruz, Frédéric de Gournay, Darren Fu, Pierrick Bourgeat, Jurgen Fripp, Clinton Fookes, and Olivier Salvado.  Corticalflow: A diffeomorphic mesh deformation module for cortical surface reconstruction.  arXiv preprint arXiv:2206.02374, 2022.
  * [44] Thomas Lewiner, Hélio Lopes, Antônio Wilson Vieira, and Geovan Tavares.  Efficient implementation of marching cubes’ cases with topological guarantees.  Journal of graphics tools, 8(2):1–15, 2003.
  * [45] Chunming Li, Chenyang Xu, Changfeng Gui, and Martin D Fox.  Distance regularized level set evolution and its application to image segmentation.  IEEE transactions on image processing, 19(12):3243–3254, 2010.
  * [46] William E Lorensen and Harvey E Cline.  Marching cubes: A high resolution 3d surface construction algorithm.  ACM siggraph computer graphics, 21(4):163–169, 1987.
  * [47] Qiang Ma, Liu Li, Emma C Robinson, Bernhard Kainz, Daniel Rueckert, and Amir Alansary.  Cortexode: Learning cortical surface reconstruction by neural odes.  IEEE Transactions on Medical Imaging, 2022.
  * [48] Julian Maclaren, Zhaoying Han, Sjoerd B Vos, Nancy Fischbein, and Roland Bammer.  Reliability of brain volume measurements: a test-retest dataset.  Scientific data, 1(1):1–9, 2014.
  * [49] Daniel S Marcus, Tracy H Wang, Jamie Parker, John G Csernansky, John C Morris, and Randy 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, 19(9):1498–1507, 2007.
  * [50] Lars Mescheder, Michael Oechsle, Michael Niemeyer, Sebastian Nowozin, and Andreas 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.
  * [51] Trung Nguyen, Quang-Hieu Pham, Tam Le, Tung Pham, Nhat Ho, and Binh-Son Hua.  Point-set distances for learning representations of 3d point clouds.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 10478–10487, 2021.
  * [52] Stanley Osher, Ronald Fedkiw, and K Piechor.  Level set methods and dynamic implicit surfaces.  Appl. Mech. Rev., 57(3):B15–B15, 2004.
  * [53] Jeong Joon Park, Peter Florence, Julian Straub, Richard Newcombe, and Steven Lovegrove.  Deepsdf: Learning continuous signed distance functions for shape representation.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 165–174, 2019.
  * [54] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  Advances in neural information processing systems, 32, 2019.
  * [55] Songyou Peng, Chiyu ”Max” Jiang, Yiyi Liao, Michael Niemeyer, Marc Pollefeys, and Andreas Geiger.  Shape as points: A differentiable poisson solver.  In Advances in Neural Information Processing Systems (NeurIPS), 2021\.
  * [56] Dzung L Pham, Pierre-Louis Bazin, and Jerry L Prince.  Digital topology in brain imaging.  IEEE Signal Processing Magazine, 27(4):51–59, 2010.
  * [57] François Pomerleau, Francis Colas, Roland Siegwart, et al.  A review of point cloud registration algorithms for mobile robotics.  Foundations and Trends® in Robotics, 4(1):1–104, 2015.
  * [58] Lev Semenovich Pontryagin.  Mathematical theory of optimal processes.  CRC press, 1987.
  * [59] Nikhila Ravi, Jeremy Reizenstein, David Novotny, Taylor Gordon, Wan-Yen Lo, Justin Johnson, and Georgia Gkioxari.  Accelerating 3d deep learning with pytorch3d.  arXiv preprint arXiv:2007.08501, 2020.
  * [60] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18, pages 234–241. Springer, 2015.
  * [61] Shunsuke Saito, Tomas Simon, Jason Saragih, and Hanbyul Joo.  Pifuhd: Multi-level pixel-aligned implicit function for high-resolution 3d human digitization.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 84–93, 2020.
  * [62] Rodrigo Santa Cruz, Léo Lebrat, Darren Fu, Pierrick Bourgeat, Jurgen Fripp, Clinton Fookes, and Olivier Salvado.  Corticalflow++: Boosting cortical surface reconstruction accuracy, regularity, and interoperability.  In Medical Image Computing and Computer Assisted Intervention–MICCAI 2022: 25th International Conference, Singapore, September 18–22, 2022, Proceedings, Part V, pages 496–505. Springer, 2022.
  * [63] Stewart Shipp.  Structure and function of the cerebral cortex.  Current Biology, 17(12):R443–R449, 2007.
  * [64] Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein.  Implicit neural representations with periodic activation functions.  Advances in Neural Information Processing Systems, 33:7462–7473, 2020.
  * [65] Shanlin Sun, Kun Han, Deying Kong, Hao Tang, Xiangyi Yan, and Xiaohui Xie.  Topology-preserving shape reconstruction and registration via neural diffeomorphic flow.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 20845–20855, 2022.
  * [66] Shanlin Sun, Kun Han, Deying Kong, Chenyu You, and Xiaohui Xie.  Mirnf: medical image registration via neural fields.  arXiv preprint arXiv:2206.03111, 2022.
  * [67] Matthew Tancik, Pratul Srinivasan, Ben Mildenhall, Sara Fridovich-Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan Barron, and Ren Ng.  Fourier features let networks learn high frequency functions in low dimensional domains.  Advances in Neural Information Processing Systems, 33:7537–7547, 2020.
  * [68] Arthur W Toga and Paul M Thompson.  Mapping brain asymmetry.  Nature Reviews Neuroscience, 4(1):37–48, 2003.
  * [69] Nanyang Wang, Yinda Zhang, Zhuwen Li, Yanwei Fu, Wei Liu, and Yu-Gang Jiang.  Pixel2mesh: Generating 3d mesh models from single rgb images.  In Proceedings of the European conference on computer vision (ECCV), pages 52–67, 2018.
  * [70] Udaranga Wickramasinghe, Edoardo Remelli, Graham Knott, and Pascal Fua.  Voxel2mesh: 3d mesh model generation from volumetric data.  In Medical Image Computing and Computer Assisted Intervention–MICCAI 2020: 23rd International Conference, Lima, Peru, October 4–8, 2020, Proceedings, Part IV 23, pages 299–308. Springer, 2020.
  * [71] Yiheng Xie, Towaki Takikawa, Shunsuke Saito, Or Litany, Shiqin Yan, Numair Khan, Federico Tombari, James Tompkin, Vincent Sitzmann, and Srinath Sridhar.  Neural fields in visual computing and beyond.  In Computer Graphics Forum, volume 41, pages 641–676. Wiley Online Library, 2022.
  * [72] Zhengxin Zhang, Qingjie Liu, and Yunhong Wang.  Road extraction by deep residual u-net.  IEEE Geoscience and Remote Sensing Letters, 15(5):749–753, 2018\.
  * [73] Hong-Yu Zhou, Jiansen Guo, Yinghao Zhang, Lequan Yu, Liansheng Wang, and Yizhou Yu.  nnformer: Interleaved transformer for volumetric segmentation.  arXiv preprint arXiv:2109.03201, 2021.
  * [74] Q Zhou.  Pymesh—geometry processing library for python.  Software available for download at https://github. com/PyMesh/PyMesh, 2019.

Supplementary Material

Figure 5: The pipeline of Hybrid-CSR.

##  Appendix A Pipeline Review

Fig. 5 depicts the whole pipeline of Hybrid-CSR. From the template meshes
$\mathcal{M}_{T}$, Hybrid-CSR first obtains coarsely deformed cortical meshes
$\mathcal{M}_{c}$, given which, we estimate the positions and normals of
upsampled oriented point cloud $\mathcal{O}_{up}$. Then the cortical surfaces
$\hat{\mathcal{M}}$ can be reconstructed via poisson surface reconstruction
from $\mathcal{O}_{up}$. To fix the topology defects in $\hat{\mathcal{M}}$,
we extract non-zero level set from the signed distance grids, obtaining
topologically correct meshes $\hat{\mathcal{M}}_{tc^{\prime}}$, and apply
optimization-based diffeomorphic registration to recover the accurate and
smooth genus-0 cortical surfaces $\hat{\mathcal{M}}_{tc}$. Lastly, we refine
$\hat{\mathcal{M}}_{tc}$ using a learning-based diffeomorphic transformation
model, to achieve our final reconstruction results $\hat{\mathcal{M}}_{f}$.

##  Appendix B Implementation Details

Hybrid-CSR framework is implemented using PyTorch [54] and executed on a
system equipped with an NVIDIA RTX A6000 GPU and an Intel i7-7700K CPU.

###  B.1 Toy Example

The target contour is controlled by 40 pivot points and the source circle
contour includes 200 pivot points. The mesh deformation is modeled via neural
fields [71]. The positions and normals of oriented point cloud are optimized
directly.

####  B.1.1 Network Architecture

To model contour displacements, we encoder positions with random Fourier
mapping [67] and approximate the deformation field with a SIREN [64] model.
The gaussian scale and embedding length of position encoding are 5 and 128.
The hidden features size and hidden layers number are 256 and 2. Our hybrid
method directly optimizes the positions and normals of oriented points.

####  B.1.2 Optimization

To optimize neural fields for contour deformation, 1000 points with normals
are separately sampled from ground truth contour and deformed coutour in each
iteration. The loss function consists of geometry-consistency loss and
regularization loss. For the geometry-consistency loss, we add up the chamfer
distance $\mathcal{L}_{cd}$ and normal distance $\mathcal{L}_{nd}$ between two
sets of 1000 sampled oriented points. The regularization loss consists of edge
length $\mathcal{L}_{edge}$ [69], as well as normal consistency regularization
$\mathcal{L}_{nc}$. The total mesh loss is as below:

$\mathcal{L}=\mathcal{L}_{cd}+0.02*\mathcal{L}_{nd}+0.005*\mathcal{L}_{edge}+0.005*\mathcal{L}_{nc}$
(7)

We apply Adam optimizer [38] with a learning rate of $1e^{-4}$ for 3000
iterations to update the parameters of neural fields.

For the hybrid method, we first uniformly sample 1000 points with normals from
the deformed contour obtained above, as initializations. To optimize the
positions and normals of oriented points, we minimize the L2 loss between the
indicator map reconstructed from the ground truth and the optimized oriented
points. We apply Adam optimizer with a learning rate of $3e^{-3}$ for 1000
iterations.

####  B.1.3 Results of diffeomorphic transformation

Figure 6: Explicit contour representation by diffeomorphic transformtion from
circle source contour.

In NMF [29], they present that diffeomorphic transformation can avoid the
“regularizer‘s dilemma”, but in our toy experiment, we found neural ode (NODE)
[14] is not suitable to model large and sharp deformations, as is shown in
Fig. 6. We model the dynamic function of NODE, i.e., neural velocity fields,
using the same neural fields as that for contour deformation and optimize with
only chamfer distance and normal distance without any regularization loss,
optimized via Adam optimizer with a learning rate of $1e^{-4}$ for 3000
iterations.

###  B.2 Coarse Mesh Deformation

####  B.2.1 Network Architecture

We apply Vox2Cortex [8] to deform template meshes. Same as Vox2Cortex,

  * •

we train the volumetric segmentation branch and mesh deformation branch end-
to-end;

  * •

volumetric segmentation branch is based on the Res-Unet;

  * •

we take four cortical surfaces template as one and use residual GCN-based
modules to encode features on graphs;

  * •

we use the same feature extraction strategy. That is, in the first step of
mesh deformation, the volumetric feature associated with meshes vertices are
from the 4th, 5th, 6th and 7th layers of Res-Unet. And in the second step of
mesh deformation, the meshes vertices are extracted from the 3rd, 4th, 7th and
8th layers of Res-Unet;

Different from Vox2Cortex, in the coarse mesh deformation module of Hybrid-CSR

  * •

we deform template meshes in two steps, instead of four steps. In other words,
our superiority in performance doesn’t come from more steps of surface
reconstruction;

  * •

In both training and inference, we use smaller templates ($\approx 42000$
vertices per surface).

####  B.2.2 Training

Same as Vox2Cortex, we apply the loss function composed of voxel loss
$\mathcal{L}_{v}ox$, curvature-weighted chamfer loss $\mathcal{L}$, normal
distance, laplacian smoothing, normal consistency as well as edge length
regularizations. The coarse mesh deformation module as well as segmentation
branch are first optimized for 50 epochs and then they will be optimized
together with the oriented point cloud estimation module for another 100
epochs. The other implementation details can be found in the supplementary
material of Vox2Cortex.

###  B.3 Oriented Point Cloud Estimation

####  B.3.1 Gated Linear Unit (GLU) for Point Estimation

Let $S$ denote the upsample ratio,
$\boldsymbol{p}_{i}^{up}\in\mathbb{R}^{S\times 3}$ denote the position of
upsampled oriented point clouds,
$\boldsymbol{v}_{i}^{up}\in\mathbb{R}^{S\times 3}$ denote the vertex of
deformed meshes repeated by $S$ times and $\boldsymbol{p}_{i}^{up}$ denote the
upsampled displacements of vertex $\boldsymbol{v}_{i}$. Let’s also define the
intermediate position of upsampled oriented point clouds as
$\boldsymbol{p^{\prime}}_{i}^{up}\in\mathbb{R}^{S\times 3}$, such that

$\boldsymbol{p}_{i}^{up}=(1-\boldsymbol{m}_{i})\odot\boldsymbol{v}_{i}^{up}+\boldsymbol{m}_{i}\odot\boldsymbol{p^{\prime}}_{i}^{up}$
(8)

where each dimension of $\boldsymbol{m}_{i}$ is between 0 and 1, controlling
the “confidence” of the intermediate results. We can represent
$\boldsymbol{p^{\prime}}_{i}^{up}$ as
$\boldsymbol{v}_{i}^{up}+\boldsymbol{d^{\prime}}_{i}$, where
$\boldsymbol{d^{\prime}}_{i}^{up}$ is intermediate displacement associated
with upsampled vertex $\boldsymbol{v}_{i}^{up}$. Therefore, the Eq. 8 can be
rewritten as:

$\boldsymbol{p}_{i}^{up}=\boldsymbol{v}_{i}^{up}+\boldsymbol{m}_{i}\odot\boldsymbol{d^{\prime}}_{i}^{up}$
(9)

As we have demonstrated in the main paper, using GLU, the displacements can be
written as
$\boldsymbol{d_{i}}=(\mathbf{W}_{0}\mathbf{f}_{i}+\mathbf{b}_{0})\odot\sigma(\mathbf{W}_{1}\mathbf{f}_{i}+\mathbf{b}_{1})$.
Then, we have

$(\mathbf{W}_{0}\mathbf{f}_{i}+\mathbf{b}_{0})\odot\sigma(\mathbf{W}_{1}\mathbf{f}_{i}+\mathbf{b}_{1})=\boldsymbol{m}_{i}\odot\boldsymbol{d^{\prime}}_{i}^{up}$
(10)

Thus, we suppose $\mathbf{W}_{0}\in\mathbb{R}^{(S\times
3)\times(d_{out}+S\times 3)}$ and $\mathbf{b}_{0}\in\mathbb{R}^{(S\times 3)}$
modulate the displacement vector, while $\mathbf{W}_{1}\in\mathbb{R}^{(S\times
3)\times(d_{out}+S\times 3)}$ and $\mathbf{b}_{1}\in\mathbb{R}^{(S\times 3)}$
modulate the “confidence ” of the predicted displacement vector. In our
experiments, $S=7$ and $d_{out}=64$.

####  B.3.2 Network Architecture for Normal Estimation

Figure 7: Network Architecture of Normal Estimation Module. We apply GCN to
encode features on graph and GLU to estimate normals of upsampled positions.

As is shown in Fig. 7, GCN is used to encode features on graph and GLU is used
to estimate normals of upsampled positions. The input of GCN is the
concatenation of grouped positions, grouped image features associated with
points and graph feature generated by the previous GCN layer. The positions
and associated features will be grouped into the same node if they are
displaced from the same vertices. We sample the image features from the 2nd,
3rd, 8th and 9th layers of Res-Unet given the point positions via linear
interpolation, and the channel number of these volumetric features are 32, 64,
16 and 8. The graph feature from the previous layer is in length of 64. Thus,
the input feature channel number of GCN is $S*(3+(32+64+16+8))+64=925$ and the
output feature channel number is 64.

The output of GCN concatenated with the grouped positions is taken as the
input of GLU, whose input feature channel number is therefore $64+S*3=85$. The
output of GLU is the normals of upsampled point clouds, thus the output
channel number is $S*3=21$.

####  B.3.3 Training

We apply weighted mean square error $\mathcal{L}_{DPSR}$ to measure the
difference between predicted and ground truth indicator grids. The weight map
is the smoothed edge map of the ground truth indicator grid. Together with
mesh-based loss proposed in Vox2Cortex, $\mathcal{L}_{DPSR}$ is used to
optimize Hybrid-CSR in an end-to-end manner for additional 100 epochs. The
parameters of oriented point cloud estimation module are optimized via an Adam
optimizer of a learning rate of $5e^{-5}$.

###  B.4 Topology Correction

We model the neural velocity fields $\mathcal{F}_{\theta}$, using the same
network architecture as described in Sec. B.1.1 and Sec. B.1.3. Other
implementation details are included in the main paper.

Different from the toy example, in the procedure of topology correction, the
source surface $\hat{\mathcal{M}}_{tc^{\prime}}$ and target surface
$\hat{\mathcal{M}}_{tc}$ have already been well-aligned, so that diffeomorphic
transformation optimized by chamfer distance is able to provide accurate
surface registration performance.

###  B.5 Surface Refinement

The network architecture of surface refinement is the same as CortexODE [47].
But there are some differences in training pial surface refinement model. In
the original CortexODE, they learn to map the ground truth WM surface to the
ground truth pial surface. Since the ground truth WM and pial surface share
the same topology, they can train with the L2 distance. However, during
inference, the source surface is the predicted WM surface obtained from
Marching Cubes so there exists a discrepancy between inference and training.
Instead, we learn to map the topological correct pial surface to the ground
truth pial surface, supervised by chamfer distance. And during inference, the
initial pial surface is also generated by the topology correction procedure.
In terms of other implementation details, we follow the original CortexODE.

##  Appendix C Visual Comparisons with Competing Methods

Figure 8: Visual Comparison between Hybrid-CSR and Other Competing Methods

In this section, we provide the comprehensive visual comparison between our
Hybrid-CSR and other competing methods, including Vox2Cortex [8],
CorticalFlow++ [62], CortexODE [47] as well as DeepCSR [16]. From Fig. 8, we
can see our proposed Hybrid-CSR in general can generate surfaces in lighter
colors, i.e., smaller point-to-surface distances, compared with other methods.


# Note on Using Singular Value Decomposition to Solve for Single-Phase Flow
and Identify Isolated Clusters in Heterogeneous Pore Networks

Ilan Ben-Noah  Juan J. Hidalgo  Marco Dentz

###### Abstract

Diluted pore networks or networks that are obtained directly from image
analysis are commonly used to evaluate flow through complex porous media.
These networks are often subjected to singular points and regions, that is,
isolated sub-networks that are not connecting the inlet and outlet. These
regions impair the invertibility of the modified Laplacian matrix and affect
the accuracy and robustness of pore network models. Searching and eliminating
these singularities is commonly done by preconditioning or searching
algorithms like breadth-first search or depth-first search. Here, we propose
using singular value decomposition (SVD) to simultaneously solve for single-
phase flow and locate isolated regions. We demonstrate the performance of the
method for networks obtained from images of a partially saturated medium, and
for diluted networks close to the percolation threshold.

††journal: Advances in Water Resources

[inst1]Institute of Environmental Assessment and Water Research (IDAEA),
Spanish National Research Council (CSIC), Barcelona, Spain. [inst2]Department
of Environmental Physics and Irrigation, Institute of Soil, Water and
Environmental Sciences, The Volcani Institute, Agricultural Research
Organization, Rishon LeZion, Israel

Singular value decomposition (SVD) is an efficient and robust method to solve
for single phase flow in heterogeneous networks with isolated clusters

SVD can be used to locate singular points and clusters

Singular value decomposition is robust, and allows for the determination of
the percolation threshold in heterogeneous networks.

##  1 Introduction

Conceptualization of the media as a pore network (rather than a bundle of
capillaries) goes back to the work of Fatt [1]. In its simplest form, a porous
medium is conceptualized as a lattice of conductors, disregarding the pore
shape. More complex (and realistic) pore network models relate the bond
geometry to their conductance [2] and map out the pore structure obtained from
pore space partitioning [3].

Depending on the flow conditions and complexity of the porous medium, the
resulting pore network of the fluid-filled pore space may be characterized by
subnetworks or clusters that are hydraulically isolated, that is, clusters
that are not connecting the inlet and outlet of the pore network [4]. For
example, imbibition or drainage in a porous medium can be modeled as the
activation or deactivation of conductors in an equivalent pore network model
(PNM) [5]. The probability of a cluster connecting inlet and outlet, can be
determined by percolation theory based on the occupation probability ($p$) of
either lattice bonds or lattice sites and their conductance [6, 7, 5, 8].
Similarly, during dilution of a pore network by randomly eliminating bonds,
i.e., assigning zero conductances to them, the resulting network may contain
isolated sites. These features cause numerical problems since they lead to
singular or ill-conditioned graph Laplacians. This issue is not restricted to
diluted regular networks but rather a general issue in any system of linear
equations with singular values, such as models derived directly from image
analysis of porous and fractured media. In order to address this singularity,
different types of iterative search algorithms have been used to find and
remove isolated sites and clusters [4, 9].

In this note, we proposed an alternative algorithm based on singular value
decomposition (SVD) of the graph Laplacian to identify isolated clusters in
the network. SVD decomposes a matrix into the product of a diagonal and two
orthonormal matrices [10]. It is used, for example, to determine the
pseudoinverse of a singular matrix and calculate low-rank approximations of a
matrix. SVD finds broad applications in data-driven science and engineering
[11, 12]. The idea here is to use SVD to identify nodes belonging to isolated
clusters by solving for single-phase flow in the PNM using the pseudoinverse
of the Laplacian matrix. The internal nodes of zero pressure belong to
isolated clusters. The proposed method is general and applies to the
identification of isolated clusters in networks.

##  2 Methods

###  2.1 Graph Laplacian

The graph Laplacian or Kirchoff matrix is a representation of the graph
underlying the PNM. The Laplacian here is weighted by the conductances of its
bonds of the graph. The conductance or weight matrix $\mathbf{G}$ contains the
conductances between $G_{ij}$ between sites $i$ and $j$. If all non-zero
conductances are unity, $\mathbf{G}$ is equal to the adjacency matrix of the
graph. The weighted degree matrix is a diagonal matrix that contains the sum
of the conductances of the bonds adjacent to a node, that is,
$D_{ij}=\delta_{ij}\sum_{[ki]}G_{ik}$, where $\sum_{[ki]}$ denotes summation
over all sites $k$ that are connected to site $i$. With these definitions, the
weighted graph Laplacian is given by $\mathbf{L}=\mathbf{D}-\mathbf{G}$. The
continuity equation for single-phase flow through the PNM can then be written
as

$\displaystyle\mathbf{L}\cdot\mathbf{P}=\mathbf{0},$ (1)

where $\mathbf{P}$ is the vector of pressures at all sites of the network.
This can be seen by defining the flow rate between sites $i$ and $j$ as

$\displaystyle Q_{ij}=G_{ij}(P_{i}-P_{j}),$ (2)

which is a statement of momentum conservation. Note that $Q_{ij}=-Q_{ji}$.
Continuity, or mass conservation, is expressed by

$\displaystyle\sum_{[ji]}Q_{ij}=0.$ (3)

The latter is equivalent to Eq. (1). The flow rates at sites are defined by

$\displaystyle Q_{i}=\sum_{[ji]}Q_{ij}$ (4)

From the mass conservation statement (3), it follows that the graph Laplacian
has the eigenvector $(1,\dots,1)$ with the eigenvalue zero, which is the
trivial solution of equal pressure at all sites. Thus, the determinant of the
graph Laplacian is zero, i.e., $\mathbf{L}$ is not invertible. The
multiplicity of the zero eigenvalue is equal to the number of isolated
clusters in the graph, or equivalently, the rank of the graph Laplacian is
equal to the number of sites minus the number of connected components[13].

We seek a solution to Eq. (1) for prescribed pressure at the inlet and outlet
nodes. To this end, we define the modified graph Laplacian
$\mathbf{L}^{\prime}$ such that $L^{\prime}_{ij}=\delta_{ij}$ for all sites
$i,j$ that are in the inlet or outlet boundaries and $L^{\prime}_{ij}=L_{ij}$
for all other sites. With this definition, the boundary value problem can be
written as

${\mathbf{L}}^{\prime}\cdot\mathbf{P}=\mathbf{b},$ (5)

where the vector $\mathbf{b}$ on the right side of Eq. (5) contains the
boundary conditions with $b_{i}=p_{0}$ for the sites $j$ at the inlet boundary
and $b_{i}=0$ else. For networks that consist of a single connected cluster,
the modified Laplacian $\mathbf{L}^{\prime}$ is an invertible (non-singular)
matrix. The trivial zero eigenvalue proper to the graph Laplacian is
eliminated by specifying non-trivial boundary conditions. Only in the presence
of isolated clusters has the modified graph Laplacian $\mathbf{L}^{\prime}$ a
zero eigenvalue and thus is not invertible. This can be illustrated with the
example of a single isolate internal site of the network. In this case, a row
in the modified Laplacian is zero and thus, its determinant is zero as well.
In order to solve the flow equation (5), and to determine the sites belonging
to isolated regions, we use the singular value decomposition.

###  2.2 Singular Value Decomposition

The singular value decomposition is the decomposition of a diagonal matrix
$\mathbf{L}$ into the product of two orthonormal matrices $\mathbf{U}$ and
$\mathbf{V}$ and a diagonal matrix $\boldsymbol{\Sigma}$ with positive real
entries $\sigma_{i}$ as

$\mathbf{L}=\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{T},$ (6)

The entries $\sigma_{i}$ are called the singular values of $\mathbf{L}$. The
number of non-zero $\sigma_{i}$ is defined as the rank of the matrix
$\mathbf{L}$. The columns of $\mathbf{V}$ (termed the right singular vectors
of $\mathbf{L}$), and the columns of $\mathbf{U}$ (the left singular vectors),
always form an orthogonal set with no assumptions on $\mathbf{L}$. If
$\mathbf{L}$ is an invertible matrix, its inverse is given by

$\mathbf{L}^{-1}=\mathbf{V}\boldsymbol{\Sigma}^{-1}\mathbf{U}^{T}.$ (7)

If $L$ is a singular matrix (i.e., at least one $\sigma_{i}=0$), its pseudo-
inverse $\mathbf{L}^{+}$ is determined as follows [10]. First, the pseudo-
inverse $\boldsymbol{\Sigma}^{+}$ of the diagonal matrix $\boldsymbol{\Sigma}$
is determined by inverting every non-zero element on the diagonal and leaving
the zero elements unchanged. In the numerical calculations, elements are
inverted if they are larger than a given threshold and set equal to zero if
they are smaller. The pseudo-inverse $\mathbf{L}^{+}$ is then given by

$\displaystyle\mathbf{L}^{+}=\mathbf{V}\mathbf{\Sigma}^{+}\mathbf{U}^{T}$ (8)

Once the pseudo-inverse is obtained, the pressure distribution is determined
by a matrix multiplication as

$\mathbf{P}=\mathbf{L}^{+}\mathbf{b}.$ (9)

The sites $i$ belonging to isolated clusters have zero pressure $P_{i}=0$
because they are not connected to inlet and outlet. This approach does not
require matrix preconditioning and can be used for singular and invertible
matrices. At the same time, it is robust and accurate.

##  3 Demonstration and Validation

In this section, we demonstrate and validate the proposed SVD approach for two
types of media characterized by isolated clusters. First, we consider a pore
network obtained by pore space partitioning of an image of the fluid phase of
a two-dimensional millifluidic device under partially saturated flow
conditions. We use the SVD-based method laid out in Section 2 on the original
network to solve the single-phase flow problem (5). The results are compared
to the direct inversion of the modified Laplacian in the pore network obtained
from pre-conditioned images, in which isolated clusters are removed. The
results are compared in terms of relative conductivity, flow rate, and
pressure statistics to assess the accuracy of the proposed method. The total
flow rate ($Q_{m}$) is defined by the sum of the flow rates at the network’s
inlet sites.

$\displaystyle Q_{m}={\sum_{i}Q_{i}},$ (10)

We determine the probability density function (PDF) of the flow rate at the
sites through areal sampling from the network flow simulations and image pore
space partitioning

$\displaystyle f(Q_{k})=\frac{\sum_{i}\mathbb{I}(Q_{i},Q_{k})A_{i}}{\Delta
Q_{k}\sum_{i}A_{i}},$ (11)

where $A_{i}$ [L2] is the area of pore $i$, and $Q_{i}$ is the flow rate at
the $i$-th site. The indicator function $\mathbb{I}$ is one if the flow rate
$Q_{i}$ is within the interval $(Q_{k},Q_{k}+\Delta Q_{k}]$ and zero
otherwise.

The distribution of the pressure at the sites is determined in analogy to the
distribution of flow rates

$\displaystyle f(P_{k})=\frac{\sum_{i}\mathbb{I}(P_{i},P_{k})A_{i}}{\Delta
P_{k}\sum_{i}A_{i}},$ (12)

Second, we consider synthetic media obtained by dilution of regular networks.
We evaluate the robustness, that is, the ability of the SVD method to solve
the flow equation of a media with many singular points and isolated clusters.
To this end, we use the SVD method to solve the steady flow equation on two
types of diluted networks based on two-dimensional hexagonal and square
lattices at different degrees of dilution. We compare the percolation
properties determined from the SVD-based flow solution to theoretical values
from percolation theory. For this media, it is convenient to evaluate it in
terms of relative network conductivity, which is defined in terms of the total
flow rate relative to the total flow rate in a corresponding regular network
as

$\displaystyle k_{r}=\frac{Q_{m}}{Q_{m}^{(0)}},$ (13)

where the superscript $0$ denotes to the flow rate for the corresponding
regular network.

Figure 1: a) Binary image of the liquid phase (white) configuration during
multi-phase flow through a milifluidic device (image provided by Yves
Méheust). b) Binary image of the processed image after removing the
disconnected and non-percolated liquid phase. In the small plot in b), the
partitioned liquid-filled pore space in a zoom-in section (blue square), pore
throats, and stagnant matrix are in black. c) Pressure distribution in the
pores using SVD to solve the non-processed image (a), singular regions are
assigned with a value of 0, P values smaller than 0.1%P(0), corresponding to
the singular sites, are colored in black. d) Pressure distribution (not using
SVD) in the pores of the processed image (b).

###  3.1 Pore network from a partially saturated medium

Figure 1 shows the fluid phase (in white) distribution in a milifluidic device
before (Fig. 1a) and after (Fig. 1b) processing the image to exclude the non-
percolated disconnected fluid clusters. The device is the same used for the
experimental work described in [14]. The original image (Fig. 1a) contains
many disconnected regions of different sizes. The small regions are formed by
wetting liquid phase films surrounding the grains, while the larger clusters
are formed by the entrapment of the liquid phase by the non-wetting gas phase.
The processed image (Fig. 1b) does not include isolated fluid clusters. The
pre-processing of the pixelized image was conducted using the MATLAB image
analysis toolbox [15] using the 8-connectivity criterion to identify isolated
clusters.

The pore-scale partitioning both in the original and pre-processed images is
performed following the methodology presented in [16]. This method uses
information about the curvatures of the distance map (obtained from the binary
image) to separate distinct pores (small plot in Fig. 1b) and to evaluate the
pore body and throat sizes. The construction of the detailed network follows
the methodology presented in [17]. A parabolic parameterization of the pore
geometry is assumed to evaluate the bond conductivities in the resulting pore
network models. For laminar flow, the bond conductance is evaluated by the
harmonic sum of the local conductivities along the pore channel with variable
width by using Boussinesq [18] equation for rectangular apertures. The pore
sizes and conductance distributions and additional information about the
milifluidic device and pore partitioning are provided in Figure SF1 in the
supplementary material.

As outlined above, for the network constructed from the original image, the
flow equation (5) is solved using the pseudoinverse (8) obtained from SVD,
which gives the pressures at all sites through Eq.(9). For the network
obtained from the pre-processed image, the flow equation is solved by direct
inversion of Eq. 5 and by using the pseudoinverse from SVD.

Figure 2: Top panel: Flow rates at the network sites, for a) the network
obtained from the original image, evaluated using SVD, and b) pre-processed
image, using direct matrix inverse. Flow rates lower than 1% of the mean flow
rate are colored in black. Bottom panel: Area weighted probability density
functions of c) the flow rate at the sites and d) the pressure at the sites
for the networks obtained from the original and pre-processed images. The flow
rates are normalized by the pre-processed image’s mean flow rate ($\langle
Q\rangle$).

Figure 2a and b show the spatial distribution of flow rates in the pore
networks obtained from the original and pre-processed images. The singular
disconnected regions attain a zero pressure value. In the network resulting
from the original image, the solution using SVD assigns zero pressure to the
isolated clusters. The clusters obtained from the SVD analysis on the original
image coincide with the clusters identified from the pre-processing step of
the binary image. The spatial flow rate distribution in the two networks is
essentially identical. This is confirmed by the probability density function
of bond flow rates and pressures for both networks, shown in the bottom panel
of Figure 2. Moreover, the device flow rate $Q_{m}$ defined in Eq. (10) for
the networks from the original and pre-processed images differ by about
$0.2\%$ further confirming the performance and accuracy of the proposed SVD-
based method. In order to assess whether this difference is a result of the
solution of the flow equation by different methods (SVD versus direct
inversion) or due to the pre-processing, we use the SVD method to solve for
the flow on the pre-processed image. The difference in the total flow rates
obtained by the two methods is smaller than $10^{-11}$, that is, they are
essentially identical. This implies that there is a slight difference in the
identification of isolated clusters by the SVD method on the pore network
extracted from the direct image and by the connectivity-based method on the
pixel image. In other words, the connected percolating clusters in the pore
networks extracted from the direct and pre-processed images are not exactly
the same. Thus, the SVD-based method provides an accurate method for the
identification of isolated clusters and solution of single phase flow in the
network obtained from the direct image without the need for pre-processing,
which may introduce some inaccuracies.

###  3.2 Diluted pore networks

Here, we evaluate the ability of the SVD-based method to solve the flow
equation and identify isolated clusters at bond occupation probabilities close
to the percolation threshold. To this end, we consider two types of diluted
lattices based on two-dimensional Hexagonal 2D lattices with coordination
number $Z_{s}=3$ and a two-dimensional square 2D lattice with coordination
number $Z_{s}=4$. For simplicity, we use a uniform conductance field
($G_{ij}=1$ for all $i,j$), unit bond length of $\ell=1$, $N=100$ sites along
each face of the lattices, and a unit pressure gradient $\Delta P/L$. The
domain length is $L=N\ell$ for the square lattice, and $L={3N\ell}/{4}$ for
the hexagonal lattice, which is oriented with the long diagonal in the main
flow direction. The total number of sites in the regular network is denoted by
$N_{s}$.

Starting from a regular lattice with coordination number $Z_{s}$, the network
is diluted to obtain a target occupation probability

$\displaystyle p_{b}=\frac{\langle Z\rangle}{Z_{s}},$ (14)

where $\langle Z\rangle$ is the mean coordination number of the diluted
network. The dilution method applied here is a simplification of the protocol
described in [4]. Each bond is randomly assigned a dilution or elimination
number between zero and one. Then, all bonds with an elimination number larger
than $p_{b}$ are assigned zero conductivity. In this way, $N_{R}=50$ network
realizations are generated for each $p_{b}$ value. Then, for each realization,
the flow equation (5) is solved using the pseudoinverse given by Eq. 8.

We determine for each target $p_{b}$ value the number of realizations for
which a percolating cluster prevails. A percolating cluster exists if the
permeability defined by Eq. (13) is larger than 0. The percolation frequency
$p_{p}$ for a given $p_{b}$ is defined as the ratio of the number $N_{p}$ of
realizations with percolating clusters and the total number $N_{R}$ of
realizations,

$\displaystyle p_{p}=\frac{N_{p}}{N_{R}}.$ (15)

For $p_{b}$ smaller than a critical value $p_{c}$, the bond percolation
threshold, there is no percolating cluster, that is, $p_{p}=0$. The
theoretical bond percolation threshold for an infinite hexagonal lattice is
$p_{c}=1-2\sin({\pi}/{18})\approx 0.6527$, and for and infinite square lattice
it is $p_{c}=\frac{1}{2}$ [5]. From the SVD-based numerical flow simulations
in the diluted networks, the percolation threshold $p_{c}$ is estimated as the
bond occupation probability $p_{b}$ for which $50\%$ of the lattice
realizations percolate [19].

As shown in Figure 3a, the numerical values obtained for $p_{c}$ are in very
good agreement with the theoretical values both for the hexagonal and square
lattices. Slight deviation of the numerical from the theoretical values can be
attributed to the finite size of the lattices with $N=100$, and the finite
number of realizations with $N_{R}=50$. Figure 3b shows the dependence of
permeability on $p_{b}$. For $p_{b}<p_{c}$, permeability tends to zero, then
for $p_{b}>p_{c}$ it increases monotonically toward $k_{r}=1$, which is the
value for the non-diluted lattices. For $p_{b}$ far from $p_{c}$, the
evolution is approximately linear. These behaviors are in accordance with the
literature [20, 21]. The variability of the $k_{r}$ values between
realizations, illustrated by the width of the error bars, is not linear, with
zero variance for $p_{b}=1$ and $p_{b}\ll p_{c}$ and a large variance for
intermediate values. These results demonstrate the robustness of the proposed
SVD method.

Figure 3:  The effect of the occupation probability ($p_{b}$) on: a) the
percolation probability ($p_{p}$) and b) the relative lattice permeability
($k_{r}$) for square (blue) and hexagonal (red) lattices. Lattices with $100$
sites in each main direction and $50$ realizations for each $p_{b}$ value are
used for this analysis. The symbols in a) denote the theoretical values of
$p_{c}$ for the square and hexagonal lattices.

##  4 Conclusion

We have shown that SVD is a robust, easy-to-implement, and accurate method to
identify sites belonging to isolated clusters in heterogeneous networks.
Isolated clusters render the modified graph Laplacian non-invertible. Thus, we
use the pseudoinverse determined from SVD to solve for single-phase flow in
the heterogeneous networks, which assigns zero pressure to sites belonging to
isolated clusters. The SVD method is demonstrated and validated for pore
network models of different complexity, which are characterized by singular
points and isolated clusters. This shows that SVD provides an alternative to
iterative search algorithms for the identification of isolated clusters and a
powerful tool to solve for single-phase flow on heterogeneous networks without
pre-conditioning.

## Acknowledgments

The authors thank Yves Méheust for the image of the milifluidic device. I.B.N,
M.D. acknowledge funding from the European Union’s Horizon 2020 research and
innovation programme under the Marie Sklodowska-Curie grant agreement No.
HORIZON-MSCA-2021-PF-01 (USFT). I.B.N. , J.J.H. and M.D. acknowledge the
support of the Spanish Research Agency (10.13039/501100011033), Spanish
Ministry of Science, Innovation and Universities through grant HydroPore II
PID2022-137652NB-C41. J.J.H. acknowledges the support of the Spanish Research
Agency (10.13039/501100011033), the Spanish Ministry of Science and Innovation
and the European Social Fund “Investing in your future” through the “Ramón y
Cajal” fellowship (RYC-2017-22300).

##  5 Open Research

No data was used in this research.

## References

  * Fatt [1956] I. Fatt,  The network model of porous media,  Transactions of the AIME 207 (1956) 144–177.
  * Sahimi [2011] M. Sahimi, Flow and transport in porous media and fractured rock: from classical methods to modern approaches, John Wiley & Sons, 2011.
  * Blunt [2017] M. J. Blunt, Multiphase flow in permeable media: A pore-scale perspective, Cambridge university press, 2017.
  * Raoof and Hassanizadeh [2010] A. Raoof, S. M. Hassanizadeh,  A new method for generating pore-network models of porous media,  Transport in porous media 81 (2010) 391–407.
  * Hunt and Sahimi [2017] A. G. Hunt, M. Sahimi,  Flow, transport, and reaction in porous media: Percolation scaling, critical-path analysis, and effective medium approximation,  Reviews of Geophysics 55 (2017) 993–1078.
  * Berkowitz and Ewing [1998] B. Berkowitz, R. P. Ewing,  Percolation theory and network modeling applications in soil physics,  Surveys in Geophysics 19 (1998) 23–72.
  * Broadbent and Hammersley [1957] S. R. Broadbent, J. M. Hammersley,  Percolation processes: I. crystals and mazes,  in: Mathematical proceedings of the Cambridge philosophical society, volume 53, Cambridge University Press, 1957, pp. 629–641.
  * Stauffer and Aharony [2018] D. Stauffer, A. Aharony, Introduction to percolation theory, CRC press, 2018.
  * Guo et al. [2021] X. Guo, K. Yang, H. Jia, Z. Tao, M. Xu, B. Dong, L. Liu,  A new method of central axis extracting for pore network modeling in rock engineering,  Geofluids 2021 (2021) 1–20.
  * Golub and Van Loan [2013] G. H. Golub, C. F. Van Loan, Matrix computations, JHU press, 2013.
  * Brunton and Kutz [2019] S. L. Brunton, J. N. Kutz, Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control, Cambridge University Press, 2019.
  * Bisgard [2020] J. Bisgard, Analysis and linear algebra: the singular value decomposition and applications, volume 94, American Mathematical Soc., 2020.
  * Bapat [2010] R. B. Bapat, Graphs and matrices, volume 27, Springer, 2010.
  * Jiménez-Martínez et al. [2017] J. Jiménez-Martínez, T. Le Borgne, H. Tabuteau, Y. Méheust,  Impact of saturation on dispersion and mixing in porous media: Photobleaching pulse injection experiments and shear-enhanced mixing model,  Water Resources Research 53 (2017) 1457–1472.
  * The MathWorks Inc. [2022] The MathWorks Inc., Matlab version: 9.13.0 (r2022b), 2022. URL: https://www.mathworks.com.
  * Ben-Noah et al. [2024a] I. Ben-Noah, J. J. Hidalgo, M. Dentz, Efficient pore space characterization based on the curvature of the distance map, 2024a. arXiv:2403.12591.
  * Ben-Noah et al. [2024b] I. Ben-Noah, J. J. Hidalgo, M. Dentz, Evaluation of different network models in variably saturated media, 2024b. arXiv:2403.13519.
  * Boussinesq [1868] J. Boussinesq,  Mémoire sur l’influence des frottements dans les mouvements réguliers des fluids,  Journal de mathématiques pures et appliquées 13 (1868) 377–424.
  * Friedman and Seaton [1998] S. Friedman, N. Seaton,  Percolation thresholds and conductivities of a uniaxial anisotropic simple-cubic lattice,  Transport in porous media 30 (1998) 241–250.
  * Kirkpatrick [1973] S. Kirkpatrick,  Percolation and conduction,  Reviews of modern physics 45 (1973) 574.
  * Friedman et al. [1995] S. P. Friedman, L. Zhang, N. A. Seaton,  Gas and solute diffusion coefficients in pore networks and its description by a simple capillary model,  Transport in porous media 19 (1995) 281–301.