# Volumetric Parameterization of the Placenta to a Flattened Template

S. Mazdak Abulnaga, Esra Abaci Turk, Mikhail Bessmeltsev, P. Ellen Grant,
Justin Solomon, and Polina Golland This work was supported in part by NIH
R01HD100009, NIH U01HD087211, NIH R01EB017337, NIH NIBIB NAC P41EB015902,
Wistron Corporation, ARO W911NF2010168, AFOSR FA9550-19-1-031, NSF
IIS-1838071, NSF CHS-1955697, NSERC RGPIN-2019-05097, FRQNT 2020-NC-270087,
NSF GRFP, NSERC PGS D, and MathWorks Fellowship.S. M. Abulnaga, J. Solomon,
and P. Golland are with the Computer Science and Artificial Intelligence Lab,
Massachusetts Institute of Technology, Cambridge, 02139, USA (e-mail:
abulnaga@mit.edu).E. Abaci Turk and P. E. Grant are with the Fetal-Neonatal
Neuroimaging and Developmental Science Center, Boston Children’s Hospital,
Harvard Medical School, Boston, MA, 02115, USA.M. Bessmeltsev is with the
Department of Computer Science and Operations Research, Université de
Montréal, Montréal, QC, H3C 3J7, Canada.

###### Abstract

We present a volumetric mesh-based algorithm for parameterizing the placenta
to a flattened template to enable effective visualization of local anatomy and
function. MRI shows potential as a research tool as it provides signals
directly related to placental function. However, due to the curved and highly
variable in vivo shape of the placenta, interpreting and visualizing these
images is difficult. We address interpretation challenges by mapping the
placenta so that it resembles the familiar ex vivo shape. We formulate the
parameterization as an optimization problem for mapping the placental shape
represented by a volumetric mesh to a flattened template. We employ the
symmetric Dirichlet energy to control local distortion throughout the volume.
Local injectivity in the mapping is enforced by a constrained line search
during the gradient descent optimization. We validate our method using a
research study of 111 placental shapes extracted from BOLD MRI images. Our
mapping achieves sub-voxel accuracy in matching the template while maintaining
low distortion throughout the volume. We demonstrate how the resulting
flattening of the placenta improves visualization of anatomy and function. Our
code is freely available at https://github.com/mabulnaga/placenta-flattening.

###### Index Terms:

Anatomy visualization, injective map, fetal MRI, placenta, shape, volumetric
mesh parameterization.

##  I Introduction

We present a volumetric method for mapping the placenta to a flattened
template to enable examination of local placental anatomy and function. The
placenta is an organ that connects the fetus to the maternal blood system to
provide oxygen and nutrients to enable fetal growth [1]. Placental
insufficiency can affect the mother’s health and cause long-term problems for
a child’s health and development [1, 2]. Blood oxygen level dependent (BOLD)
MRI has been used as a research tool to quantify oxygen transport within the
placenta [3, 4, 5, 2, 6, 7], providing initial evidence for clinical utility
of MRI for functional assessment. Further, T2-weighted MRI captures anatomical
information such as the vascular structure of the placenta [8]. Since MRI
provides direct measurements of placental function, it is a promising research
tool for studying the placenta [3, 4, 2, 7, 9, 10].

MRI has been used to study placental oxygen transport [2, 11], fetal-maternal
circulation [12], uterine contractions [7], and to characterize the
relationship between placental age and its appearance in MRI [13]. A promising
research direction is to use MRI to develop biomarkers of pathology. BOLD and
T2* MRI show promise for identifying intrauterine growth restriction (IUGR)
[14], for quantifying differences in discordant twins with IUGR [4], and for
characterizing preeclampsia [10]. While emerging as a promising tool to study
placental function and health, effective visualization of MRI signals in the
placenta is an open problem.

Interpreting MRI scans of the placenta is challenging due to its in vivo
shape. The placental shape is determined by the curved shape of the uterine
wall, making it difficult to visualize the interior of the organ using
standard techniques, for example, using cutting planes. As the location of the
placenta’s attachment to the uterine wall, its shape, and its appearance vary
greatly across subjects, a standard visualization of the functional and
anatomical images of the placenta is yet to be demonstrated.

As part of the post-delivery examination process, the placenta is placed on a
flat surface, and measurements of its disk-like shape and notes about its
anatomical features are recorded [15]. Motivated by this framework for
examination of the placenta, we present an algorithm for mapping the placental
shape observed in an MRI scan to a flattened template that resembles the
organ’s well-studied ex vivo shape.

We represent the placental volume as a tetrahedral mesh and find the
deformation to a flattened template as the solution to an optimization
problem. We provide a general formulation that readily accepts a broad set of
templates and shape distortion measures. We compute an injective map by
constraining the minimization of the objective function, to ensure that
anatomical fidelity is preserved. In this paper, we select two parallel planes
as our flattened template surface and measure shape distortion by the
symmetric Dirichlet energy to penalize local deformations [16]. We propose and
evaluate alternative templates including an ellipsoid and a relaxation of the
parallel planes template. Moreover, we provide guidelines for template
selection based on the desired visualization or computation task.

We evaluate the proposed method by applying it to segmented BOLD images of a
research study. We demonstrate effective mapping of the highly variable
placental shapes in our dataset. We highlight the utility of our method for
research by presenting enhanced visualization of local anatomical and
functional features, and demonstrate the possibility of quantifying MRI
signals within functional regions of the placenta.

###  I-A Related Work

We build on state-of-the-art mesh parameterization methods to estimate the
optimal deformation of the placenta onto a template. Parameterization is the
problem of mapping from a shape to a parameter domain, for example, from a
surface mesh to a plane [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
30] or to a sphere [31, 32, 33, 34, 35, 36]. In geometry processing,
parameterization methods are useful for applications such as texture mapping,
smoothing, and remeshing [37]. In medical imaging, parameterization has been
used to improve visualization of anatomy and to facilitate population studies
[38]. Surface-based parameterization has been used in many applications
including cortical surface mapping [32, 33, 34, 17, 36, 35], heart and vessel
flattening [39, 40], colon unfolding [28, 29], bone unfolding [41, 42], and
placental flattening [27]. For a more detailed treatment of parameterization
methods in geometry processing and in medical imaging, we refer the reader to
surveys [43, 44, 37, 38]. In contrast to most algorithms that focus on mapping
2D surface meshes, we address the parameterizion of volumetric meshes to a 3D
canonical template volume.

While surface parameterization methods are not directly applicable to
volumetric meshes, they provide a framework for computing the desired
bijective map with minimal distortion. Surface parameterization approaches use
Tutte’s embedding to guarantee a bijective planar parameterization through a
linear map to a closed, convex boundary [20, 45]. A bijective spherical
parameterization can be achieved using Tutte’s embedding with an inverse
stereographic projection [32, 36] or by using one of several proposed
extensions of the embedding to the sphere [31, 46]. As this initial
parameterization is likely to have high distortion, the map is refined by
minimizing a cost function that includes a distortion metric to achieve a
parameterization with desirable properties such as bijectivity and minimal
distortion. Existing methods optimize distortion metrics to compute an angle-
preserving (conformal) map [24, 25, 32, 33, 36, 28, 29, 42], an area-
preserving (authalic) map [30, 17], a distance-preserving (isometric) map [21,
22, 23, 34, 17], or a map that minimizes a combination of undesired
distortions [21, 26, 16, 47]. To achieve a bijective map, a barrier functional
is used to avoid flipped triangles and boundary intersections [18].

In the volumetric case, the existence of an exact bijective parameterization
is not guaranteed [48, 49]. The challenge is to enforce injectivity when
mapping to a defined parameterized domain. Many of the distortion measures
defined in 2D are easily extendable to 3D. Aigerman et al. [48] develop an
algorithm that projects an existing map to the space of injective and bounded-
distortion maps. They demonstrate successful parameterization to a polycube
domain. In contrast, we enforce injectivity at every step of the optimization
to explicitly control distortion while maintaining an injective map. We use
gradient descent with line search to prevent tetrahedra from folding. This
line search method was originally proposed for 2D parameterization [18] and
has been effective for 3D deformation of tetrahedral meshes of cubes [21, 23]
and animal shapes in computer graphics [23]. Here, we minimize local
distortion using the 3D extension of the symmetric Dirichlet energy [16]. We
propose a problem-specific unique volumetric template to serve as our
parameterized domain during the optimization. We demonstrate successful
parameterization of a dataset of highly variable placental shapes.

Most closely related to our work is the method by Miao et al. [27] that
represents the placenta volume as a stack of parallel surfaces spanning the
thickness of the organ, spaced at discrete intervals. The surfaces are
estimated based on the medial axis of the placenta. Each surface is flattened
independently using mean value coordinates [19] to iteratively move each
interior vertex to the average of its neighbors and to map the surface to a
disk. There are two limitations to this approach. First, due to the highly
variable shape of the placenta, parameterizing to a fixed boundary leads to
high distortion, affecting important image features. Second, because the 2D
surfaces are each parameterized independently, no correspondence or alignment
across layers is enforced, which distorts important depth-wise image features.
In contrast, we propose a continuous volumetric parameterization without a
fixed boundary, ensuring consistency and minimal distortion throughout the
placental volume. Our work presents a standardized method for visualizing
signals within the placenta, and is the first step towards developing a common
coordinate system for placental analysis.

A preliminary version of our method was presented at the International
Conference on Medical Image Computing and Computer-Assisted Intervention
(MICCAI) 2019 [50]. Here, we derive the parameterization in greater detail by
expanding on the definition of our template and the optimization procedure. We
provide a template formulation that is robust to irregularities in the meshing
procedure and extensions to closely match the ex vivo shape. We demonstrate
our algorithm on a significantly larger research dataset of 111 placental
shapes. Finally, we conduct extensive experiments and demonstrate that the map
is robust to placental shape, produces uniform distortion, and improves
visualization of local anatomy and function.

This paper is organized as follows. In the next section, we formulate the
optimization problem for parameterizing placental volumes to our proposed
flattened template, outline the optimization procedure to compute an injective
map, and provide implementation details. In Section III, we introduce the BOLD
MRI data and report the experimental results and comparisons with the baseline
method. In Section IV, we discuss potential applications of our method,
outline limitations in using the proposed template as a common space, and
present future directions to incorporate anatomical information to create a
common coordinate system. Our conclusions follow in Section V.

##  II Methods

We model the placental shape as a tetrahedral mesh that contains $N$ vertices
and $K$ tetrahedra. In our experiments, we extract the mesh from the placenta
segmentation in an MRI scan as described in Section III-A. We parameterize the
mapping via mesh vertex locations and interpolate the deformation to the
interior of each tetrahedron using a locally affine (piecewise-linear) model.
As the mapping is geometry-based, it is independent of the imaging modality
used to obtain the input mesh.

Figure 1: Visual depiction of the mapping procedure. Left: original placental
mesh, with boundary triangles shown. The solid black line indicates an area
that will be cut to illustrate the interior tetrahedra. Center: the original
placental mesh after cutting the solid lined-area, with the interior
tetrahedra shown in dark gray; right: the resulting flattened placental mesh,
depicting the interior tetrahedra. We use Z and X to denote the variables in
the original image space and the template space respectively.

###  II-A Problem Formulation

Let $x_{i}\in\mathbb{R}^{3}$ define the 3D coordinates of mesh vertex $i$ in
the template coordinate system, and let $X=\\{x_{1}\dots
x_{N}\\}\in\mathbb{R}^{3\times N}$ be the coordinates of all vertices in the
mesh. Tetrahedron $k$ is represented by a matrix $X_{k}\in\mathbb{R}^{3\times
4}$ whose columns are the 3D coordinates of the four corner vertices of
tetrahedron $k\in\\{1,\ldots,K\\}$ in the template space. Similarly, we use
$Z\in\mathbb{R}^{3\times N}$ to denote the coordinates of mesh vertices
$\\{z_{i}\\}_{i=1}^{N}$ in the original image space, and
$Z_{k}\in\mathbb{R}^{3\times 4}$ to define tetrahedron $k$ in the original
image space. We use $\partial$ to denote the boundary of a mesh, where
$\partial Z$ defines the set of vertex indices of $Z$ that lie on the mesh
boundary. The elements of $\partial Z$ remain constant throughout the
parameterization. Fig. 1 provides a visualization of the flattening procedure
and the notation.

We formulate placental parameterization as an optimization problem over the
mesh vertices that seeks to map the vertices to the template space while
minimizing a measure of local distortion:

$\phi(X,\theta)=\mathop{\mathchoice{\vtop{\halign{#\cr$\hfil\displaystyle\sum_{m\in\partial
Z}A_{m}\>T\left(x_{m},\theta\right)\hfil$\crcr\kern
2.0pt\nointerlineskip\cr\@mathmeasure\displaystyle{\sum_{m \in\partial Z }
A_{m} \> T\left(x_{m},
\theta\right)}\@mathmeasure\displaystyle{\upbrace}\@mathmeasure\displaystyle{\upbraceg}\@mathmeasure\displaystyle{\upbracegg}\@mathmeasure\displaystyle{\upbraceggg}\@mathmeasure\displaystyle{\upbracegggg}$\displaystyle\bracelu\leaders\hbox{$\bracemid$}\hfill\bracemu\leaders\hbox{$\bracemid$}\hfill\braceru$\crcr}}}{\vtop{\halign{#\cr$\hfil\textstyle\sum_{m\in\partial
Z}A_{m}\>T\left(x_{m},\theta\right)\hfil$\crcr\kern
2.0pt\nointerlineskip\cr\@mathmeasure\textstyle{\sum_{m \in\partial Z } A_{m}
\> T\left(x_{m},
\theta\right)}\@mathmeasure\textstyle{\upbrace}\@mathmeasure\textstyle{\upbraceg}\@mathmeasure\textstyle{\upbracegg}\@mathmeasure\textstyle{\upbraceggg}\@mathmeasure\textstyle{\upbracegggg}$\textstyle\bracelu\leaders\hbox{$\bracemid$}\hfill\bracemu\leaders\hbox{$\bracemid$}\hfill\braceru$\crcr}}}{\vtop{\halign{#\cr$\hfil\scriptstyle\sum_{m\in\partial
Z}A_{m}\>T\left(x_{m},\theta\right)\hfil$\crcr\kern
2.0pt\nointerlineskip\cr\@mathmeasure\scriptstyle{\sum_{m \in\partial Z }
A_{m} \> T\left(x_{m},
\theta\right)}\@mathmeasure\scriptstyle{\upbrace}\@mathmeasure\scriptstyle{\upbraceg}\@mathmeasure\scriptstyle{\upbracegg}\@mathmeasure\scriptstyle{\upbraceggg}\@mathmeasure\scriptstyle{\upbracegggg}$\scriptstyle\bracelu\leaders\hbox{$\bracemid$}\hfill\bracemu\leaders\hbox{$\bracemid$}\hfill\braceru$\crcr}}}{\vtop{\halign{#\cr$\hfil\scriptscriptstyle\sum_{m\in\partial
Z}A_{m}\>T\left(x_{m},\theta\right)\hfil$\crcr\kern
2.0pt\nointerlineskip\cr\@mathmeasure\scriptscriptstyle{\sum_{m \in\partial Z
} A_{m} \> T\left(x_{m},
\theta\right)}\@mathmeasure\scriptscriptstyle{\upbrace}\@mathmeasure\scriptscriptstyle{\upbraceg}\@mathmeasure\scriptscriptstyle{\upbracegg}\@mathmeasure\scriptscriptstyle{\upbraceggg}\@mathmeasure\scriptscriptstyle{\upbracegggg}$\scriptscriptstyle\bracelu\leaders\hbox{$\bracemid$}\hfill\bracemu\leaders\hbox{$\bracemid$}\hfill\braceru$\crcr}}}}\limits_{\text{Template
match}}+\mathop{\mathchoice{\vtop{\halign{#\cr$\hfil\displaystyle\lambda\sum_{k=1}^{K}V_{k}\>\mathcal{D}\left(X_{k}\right)\hfil$\crcr\kern
2.0pt\nointerlineskip\cr\@mathmeasure\displaystyle{\lambda\sum_{k=1}^{K} V_{k}
\>
\mathcal{D}\left(X_{k}\right)}\@mathmeasure\displaystyle{\upbrace}\@mathmeasure\displaystyle{\upbraceg}\@mathmeasure\displaystyle{\upbracegg}\@mathmeasure\displaystyle{\upbraceggg}\@mathmeasure\displaystyle{\upbracegggg}$\displaystyle\bracelu\leaders\hbox{$\bracemid$}\hfill\bracemu\leaders\hbox{$\bracemid$}\hfill\braceru$\crcr}}}{\vtop{\halign{#\cr$\hfil\textstyle\lambda\sum_{k=1}^{K}V_{k}\>\mathcal{D}\left(X_{k}\right)\hfil$\crcr\kern
2.0pt\nointerlineskip\cr\@mathmeasure\textstyle{\lambda\sum_{k=1}^{K} V_{k} \>
\mathcal{D}\left(X_{k}\right)}\@mathmeasure\textstyle{\upbrace}\@mathmeasure\textstyle{\upbraceg}\@mathmeasure\textstyle{\upbracegg}\@mathmeasure\textstyle{\upbraceggg}\@mathmeasure\textstyle{\upbracegggg}$\textstyle\bracelu\leaders\hbox{$\bracemid$}\hfill\bracemu\leaders\hbox{$\bracemid$}\hfill\braceru$\crcr}}}{\vtop{\halign{#\cr$\hfil\scriptstyle\lambda\sum_{k=1}^{K}V_{k}\>\mathcal{D}\left(X_{k}\right)\hfil$\crcr\kern
2.0pt\nointerlineskip\cr\@mathmeasure\scriptstyle{\lambda\sum_{k=1}^{K} V_{k}
\>
\mathcal{D}\left(X_{k}\right)}\@mathmeasure\scriptstyle{\upbrace}\@mathmeasure\scriptstyle{\upbraceg}\@mathmeasure\scriptstyle{\upbracegg}\@mathmeasure\scriptstyle{\upbraceggg}\@mathmeasure\scriptstyle{\upbracegggg}$\scriptstyle\bracelu\leaders\hbox{$\bracemid$}\hfill\bracemu\leaders\hbox{$\bracemid$}\hfill\braceru$\crcr}}}{\vtop{\halign{#\cr$\hfil\scriptscriptstyle\lambda\sum_{k=1}^{K}V_{k}\>\mathcal{D}\left(X_{k}\right)\hfil$\crcr\kern
2.0pt\nointerlineskip\cr\@mathmeasure\scriptscriptstyle{\lambda\sum_{k=1}^{K}
V_{k} \>
\mathcal{D}\left(X_{k}\right)}\@mathmeasure\scriptscriptstyle{\upbrace}\@mathmeasure\scriptscriptstyle{\upbraceg}\@mathmeasure\scriptscriptstyle{\upbracegg}\@mathmeasure\scriptscriptstyle{\upbraceggg}\@mathmeasure\scriptscriptstyle{\upbracegggg}$\scriptscriptstyle\bracelu\leaders\hbox{$\bracemid$}\hfill\bracemu\leaders\hbox{$\bracemid$}\hfill\braceru$\crcr}}}}\limits_{\text{Distortion}},$
(1)

where $x_{m}$ denotes the coordinates of boundary vertex $m$, function
$T\left(\cdot\right)$ measures the distance of the boundary vertices to the
template shape, $\theta$ represents the parameters defining the template,
$A_{m}$ is a constant defined to be the normalized barycentric area associated
with vertex $m$ on the boundary of the mesh (i.e., $A_{m}$ is proportional to
the sum of the areas of the triangle faces that include vertex $m$ and
$\sum_{m\in\partial Z}A_{m}=1$), $V_{k}$ is a constant defined to be the
normalized volume of tetrahedron $k$ in the original shape space
($\sum_{k=1}^{K}V_{k}=1$), and $\lambda$ is a parameter that governs the
trade-off between the template match and the shape distortion. We define the
specific terms of the cost function below.

###  II-B Shape Distortion

The formulation of the distortion term in Eq. (1) naturally accepts a range of
distortion metrics provided they are differentiable in vertex positions $X$.
The distortion term regularizes the mapping.

We model the deformation of tetrahedron $Z_{k}$ in the original image space to
its image $X_{k}$ in the template space as an affine transformation. The
Jacobian matrix

$J(X_{k})=\left(X_{k}B\right)\left(Z_{k}B\right)^{-1}$ (2)

defines the deformation of tetrahedron $k$ based on the new vertex coordinates
$X_{k}$ while ignoring the shared translation component. The constant matrix
$B\in\mathbb{R}^{4\times 3}$ extracts three basis vectors defining the
tetrahedron; the vectors originate from the first corner vertex of the
tetrahedron and point to its other three corner vertices (see Appendix A-A for
the definition of $B$).

We measure tetrahedron-specific distortion using the symmetric Dirichlet
energy density

$\displaystyle\mathcal{D}(X_{k})$
$\displaystyle=\mathcal{D}\left(J\left(X_{k}\right)\right)$
$\displaystyle=\|J\left(X_{k}\right)\|_{F}^{2}+\|\left(J\left(X_{k}\right)\right)^{-1}\|_{F}^{2}=\sum_{i=1}^{3}\left(\sigma_{k,i}^{2}+\sigma_{k,i}^{-2}\right)$
(3)

where $\|\cdot\|_{F}$ is the Frobenius norm and
$\\{\sigma_{k,1},\sigma_{k,2},\sigma_{k,3}\\}$ are the singular values of
matrix $J(X_{k})$ [21, 18, 16]. We chose the symmetric Dirichlet energy as it
penalizes expansion and shrinking equally and prevents tetrahedra from
expanding or shrinking dramatically as the corresponding singular values or
their reciprocals become unbounded. This distortion energy favors a locally
injective map by preventing tetrahedra from collapsing to a point.

Figure 2: Placenta parcellation procedure: a) placenta mesh with boundary
vertex normals shown; b) spectral embedding of the second smallest eigenvector
of the Laplacian $L$; c) placenta boundary segmented into the maternal
(yellow) and fetal (blue) sides. The placenta margin is the set of vertices on
the edge of the two clusters; d) the expanded placenta margin is shown in
green.

###  II-C Template Definition

We propose two template parameterizations to resemble the ex vivo placental
shape. The first contains two parallel planes to enforce a uniform thickness,
and the second models the placental shape as an ellipsoid.

####  II-C1 Parallel Planes Template

We construct the parallel planes template to mimic the post-delivery
examination of the placenta where the maternal side of the organ is placed
flat on an examination table. We define the template as two parallel planes,
one for the maternal side and and one for the fetal side of the organ. The two
planes ensure a constant thickness of the placenta, which facilitates 2D
visualization as we demonstrate later in the paper. Function $T(\cdot)$
measures the distance to the appropriate plane:

$T(x_{m},h)=\begin{cases}(x_{m}^{(3)}-h)^{2}&\text{if
}m\in\mathcal{F}(\partial Z),\\\ (x_{m}^{(3)}+h)^{2}&\text{if
}m\in\mathcal{M}(\partial Z),\\\ 0&\text{otherwise},\end{cases}$ (4)

where $h$ is the template half-height ($\theta=h$ in Eq. (1)), $x_{m}^{(3)}$
is the third coordinate of vertex $m$ in the template coordinate system, and
$\mathcal{F}(\partial Z)$, $\mathcal{M}(\partial Z)$ denote the subset of
vertex indices of the fetal and maternal sides of the placental boundary
surface $\partial Z$, respectively. The template half-height $h$ can be
treated as either a fixed parameter or a variable. We use the term “placenta
margin” to denote the region of high curvature that separates the fetal and
maternal sides of the placenta. The placenta margin is not mapped to either
plane; the margin vertices are driven only by the distortion term.

The template term in Eq. (4) requires identifying the maternal and fetal sides
of the placenta boundary. As we expect the boundary normals to cluster into
two distinct groups that correspond to the two sides of the placenta, we
construct a similarity measure based on the angle between boundary vertex
normals and project the vertices to a 1D space using spectral embedding [51,
52], described below.

For boundary vertex $i$, we determine the unit normal $\hat{n}_{i}$ as a
weighted average of normals over the mesh boundary triangles that contain
vertex $i$, normalized to have unit length. The weighting is proportional to
the triangle areas. We construct an affinity matrix $W\in\mathbb{R}^{|\partial
Z|\times|\partial Z|}_{\geq 0}$ so that

$w_{i,j}=\exp\left\\{\gamma\left\langle\hat{n}_{i},\hat{n}_{j}\right\rangle
l\left(x_{i},x_{j}\right)\right\\}$ (5)

for any two boundary vertices $i$ and $j$ that are connected by any path
containing three or fewer edges (the three-ring neighborhood), and set
$w_{i,j}=0$ otherwise. The parameter $\gamma$ penalizes the local effect of
variation in the orientation of the normals, $\langle\cdot,\cdot\rangle$ is a
dot product, and function $l(x_{i},x_{j})$ measures the approximate geodesic
distance between vertex positions $x_{i},x_{j}$, normalized by the maximum
distance across all boundary vertices connected within a three-ring
neighborhood. Smaller values of parameter $\gamma$ cause sensitivity to fast
changes in the orientation of normal vectors, and would likely cluster a small
region of the mesh with high variation of normals, such as a corner-like
region. We approximate geodesic distance using Dijkstra’s shortest path
algorithm [53] along boundary mesh edges, with edge weights given by the
Euclidean distance. The weighting by geodesic distance corrects for meshing
irregularities along the boundary. This prevents clusters forming in small but
dense areas within the fetal or maternal sides of the boundary and ensures the
appropriate clusters form to capture the large regions of interest.

The normalized Laplacian of $W$ is defined as

$L=I-D^{-\frac{1}{2}}WD^{-\frac{1}{2}},$ (6)

where $D$ is a diagonal matrix with $d_{i,j}=\sum_{j}w_{i,j}$ and $I$ is the
identity matrix. The second smallest eigenvector of $L$ defines a spectral
embedding of the boundary vertices that can be thresholded to segment the
boundary into two components with consistent orientation of the normals [51].
After segmenting the placenta boundary, we identify the fetal and maternal
sides. We use the fact that the maternal side of the placenta boundary has a
more convex shape since it is attached to the curved uterine wall. We
construct the convex hull of the boundary to assign the maternal label to the
cluster with the larger number of vertices on the hull.

We define the margin of the placenta as the connected region separating the
fetal and maternal sides. We identify vertices on the edge of the two
clusters, i.e., those with neighbors in the other cluster, and assign them to
the margin. We use the approximated geodesic distance to define this
neighborhood, and assign all vertices within a specified distance. This
distance, or the margin width, is independent of the template half-height h.
Using the geodesic distance accounts for irregularities in the mesh and
ensures the margin separates the two sides of the placenta with consistent
width. Fig. 2 provides a visualization of the placenta parcellation procedure.
The vertex labeling as belonging to the fetal side, maternal side, or the
placenta margin is defined at initialization on the original shape and remains
fixed throughout the optimization.

####  II-C2 Ellipsoid Template

An ellipsoid template is proposed as a simple parameterization to model the
unfolded placenta. Function $T(\cdot)$ measures the distance to an ellipsoid,

$T(x_{m},R)=(x_{m}^{T}Rx_{m}-1)^{2},$ (7)

where $R$ are the ellipsoidal parameters,
$R=\text{diag}(r_{x}^{-2},r_{y}^{-2},r_{z}^{-2})$ ($\theta=R$ in Eq. (1)).

###  II-D Optimization

We minimize the objective function $\phi(\cdot)$ in Eq. (1) using gradient
descent over the vertex locations and the template half-height. We initialize
the mapping by the identity transformation by setting $X=Z$. The gradient of
$\phi(\cdot)$ is given in Appendix A-B.

We optimize using line search to compute an injective map by preventing
tetrahedra from “flipping” at every iteration [18]. A tetrahedron “flips”
after crossing the singularity point of zero volume, causing tetrahedra to
fold onto one another. At every iteration, we determine the largest step size
value $\eta$ such that moving the current vertex locations $X$ by
$-\eta\nabla\phi(X)$ avoids singularities for all tetrahedra [18]. The
(signed) volume $V_{k}$ of tetrahedron $k$

$V_{k}\left(\eta\right)=\frac{1}{6}\det\Big{[}\big{(}X_{k}-\eta\nabla\phi(X_{k})\big{)}B\Big{]}$
(8)

is a cubic polynomial of $\eta$ whose real and positive roots are the step
sizes that lead to the tetrahedron flipping. The upper limit $\eta_{max}$ in a
line search is found by finding the smallest positive real root of Eq. (8)
over all tetrahedra [18]. We optimize using backtracking line search. At every
iteration, $\eta$ is set to a fraction $\beta\in(0,1)$ below $\eta_{max}$,
$\eta=\beta\cdot\eta_{max}$, and is decreased by a factor $\rho\in(0,1)$ until
the new set of vertex locations decreases the objective function [54]. Since
the original value of $\eta_{max}$ can result in an increase of the objective
function, precluding convergence, backtracking line search is needed to ensure
convergence by guaranteeing that the objective function decreases in every
iteration. Algorithm 1 describes the optimization procedure.

Algorithm 1 Volumetric Map with Local Injectivity

Input: $Z,\theta,\lambda$

Output: $X$

1:$X^{(0)}\leftarrow Z$

2:$\theta^{(0)}\leftarrow{\theta}$

3:$\phi^{(0)}\leftarrow\phi(X^{(0)},\theta^{(0)},\lambda$) // initial
objective function

4:$\nabla_{X}\phi^{(0)}\leftarrow\nabla_{X}\phi(X^{(0)},\theta^{(0)},\lambda)$
// initial gradient

5:$l\leftarrow 1$ // iteration counter

6:while $\|\nabla_{X}\phi^{(l-1)}\|_{F}>\epsilon$ do

7: $\eta_{max}\leftarrow\min\left(\text{Eq.}\eqref{eqn-tetvolume}\right)$ //
over all tetrahedra

8: $\eta\leftarrow\beta\cdot\eta_{max}$

9: $\textit{X}^{(l)}\leftarrow X^{(l-1)}-\eta\nabla_{X}\phi^{(l-1)}$

10: $\theta^{(l)}\leftarrow\theta^{(l-1)}-\eta\nabla_{\theta}\phi^{(l-1)}$

11: $\phi^{(l)}\leftarrow\phi(X^{(l)},\theta^{(l)})$

12: while $\phi^{(l)}\geq\phi^{(l-1)}\ \textbf{do}$ // backtracking line
search

13: $\eta\leftarrow\eta\cdot\rho$

14: $\textit{X}^{(l)}\leftarrow X^{(l-1)}-\eta\nabla_{X}\phi^{(l-1)}$

15: $\phi^{(l)}\leftarrow{\phi(X^{(l)},\theta^{(l)}})$

16: $l\leftarrow l+1$

17:$X\leftarrow X^{(l)}$

###  II-E Intensity Mapping

Once the optimal vertex coordinates are found, we map image intensity inside
the placenta segmentation to the template coordinate system. Since our mapping
is injective and piecewise affine, we use barycentric coordinates to determine
the unique transformation for any point inside the placenta. Specifically,
point $x$ inside tetrahedron $k$ can be represented as a convex combination of
the tetrahedron’s corner vertices, i.e., $x=X_{k}\alpha$, through the point’s
barycentric coordinates $\alpha$ as follows:

$\left[\alpha_{2}\;\;\alpha_{3}\;\;\alpha_{4}\right]^{T}=\left(X_{k}B\right)^{-T}\left(x-X_{k}^{1}\right),$
(9)

where $X_{k}^{1}$ is the first column of matrix $X_{k}$. Setting
$\alpha_{1}=1-(\alpha_{2}+\alpha_{3}+\alpha_{4})$ completes the vector to
ensure $\sum_{i=1}^{4}\alpha_{i}=1$ [55]. The point’s coordinates
$z=Z_{k}\alpha$ in the original volume are used to pull the image intensity
into the template coordinate system using linear interpolation.

###  II-F Implementation Details

We generate tetrahedral meshes from segmentation labelmaps using iso2mesh
[56], a MATLAB toolbox that provides an interface to TetGen [57], a commonly
used tetrahedral meshing library. Prior to applying the mapping, we center the
mesh and rotate it to align its principal axes to the axes of the template. We
assume that the algorithm convergences when the Frobenius norm of the gradient
is lower than $\epsilon=10^{-4}$.

We implemented the algorithm in MATLAB using GPU functionality to parallelize
the computation of the gradient and the line search. We ran our experiments on
an NVIDIA Titan Xp (12GB) GPU. The parameterization ran for an average of 5097
iterations for the parallel planes template and 7203 iterations for the
ellipsoid template on our dataset of meshes containing a mean of 6708
tetrahedra. The mean and 90 percentile wall clock times were $17$ and $31$
minutes respectively for the parallel planes template and $26$ and $49$
minutes respectively for the ellipsoid template.

We employ grid search to set the values of the hyper-parameters. We set the
shape distortion parameter $\lambda=1$ as it was in the optimal trade-off
range between the template fit and shape distortion, while producing a minimal
final objective value. We observe that $\lambda\in[0.1,5]$ produced similar
results (Fig. 3), thus the method is not very sensitive to the parameter
$\lambda$. We initialize the template half-height $h$ and the ellipsoid radius
$r_{z}$ as half of the placenta thickness, estimated from the 95 percentile of
the histogram of the distance transform values inside the segmentation
boundary. We initialize the ellipsoid radius $r_{y}$ as half of the placenta
width, estimated from the 95 percentile of the histogram of shortest path
geodesic distances along the placenta surface. We initialize the ellipsoid
radius $r_{x}$ to maintain the placenta volume. We set the spectral clustering
parameter $\gamma=20$ and use a boundary geodesic distance of five voxels (15
mm in our application) as the half-width of the placenta margin. We set the
line search parameters to $\beta=0.9$, $\rho=0.5$. Our implementation is
freely available at https://github.com/mabulnaga/placenta-flattening.

Figure 3: The trade-off between template match and Dirichlet distortion for
$\lambda\in[10^{-3},10^{2}]$, approaching optimality for $\lambda\in[0.1,5]$.

##  III Experimental Results

In this section, we present our experimental results in data from a clinical
research study to assess the quality of the transformation by visualizing the
mapped BOLD MRI signals in the flattened space and quantify the resulting
levels of distortion. We compare our approach to a 2D baseline
parameterization method.

###  III-A Data

We illustrate the utility of the proposed approach on a set of $111$ MRI scans
acquired in 78 pregnant women. MRI BOLD scans were acquired on a 3T Skyra
Siemens scanner (single-shot GRE-EPI, 3 mm isotropic voxels, interleaved slice
acquisition, TR=5.8–8s, TE=32–36ms, FA= 90∘) using 18-channel body and
12-channel spine receive arrays. There were a total of 18 twin and 60
singleton pregnancies, with gestational ages ranging from $26$ weeks, $6$ days
to $37$ weeks, $6$ days. Scans of $33$ of the $60$ singleton pregnancy
subjects were acquired in the supine and lateral positions, yielding $66$ of
the total $111$ segmentations as the shape of the placenta changed between the
two maternal positions. The twin pregnancies are monochornionic, meaning one
placenta is shared by the twin fetuses.

The images were pre-processed to correct for signal nonuniformity and fetal
motion during acquisition [58]. Following preprocessing, the placenta was
manually segmented by a trained expert. The resulting segmentation label maps
were used as input to the meshing software. The resulting meshes had (mean
$\pm$ std.) $6708\pm 1813$ tetrahedra and $2801\pm 731$ surface triangles.

###  III-B Method Variants and Baseline

We evaluate our method on three proposed templates: ellipsoid, parallel
planes, and single plane. The single plane template is a variant of the
parallel planes that allows the fetal vertices to relax while keeping the
maternal vertices fixed to the bottom plane. This template is a relaxation of
the parallel planes, and more closely mimics the ex vivo examination process
in which the maternal side is placed on an examination table. To implement
this template, we first map the placenta to the parallel planes template and
then run the optimization with a zero weight assigned to the fetal plane in
the template term in of the optimization function.

We compare our methods with a 2D parameterization approach. We emulate the
prior placenta parameterization method of Miao et al. [27] in which 2D curved
surfaces spanning the placenta volume are independently flattened to a disk in
$\mathbb{R}^{2}$ in a complex multi-step computational pipeline. These
surfaces are constructed by cutting Euclidean distance level sets from the
estimated medial axis of the placenta, with the surface boundaries defined to
minimize local curvature change. We derive these 2D surfaces by intersecting
the flattened volume with horizontal planes spaced one voxel apart. The
intersection yields a set of parameterized surfaces. We map the surfaces back
to the original placenta volume using the known barycentric coordinates,
yielding a set of surfaces spanning the original placenta volume. We use
harmonic parameterization [59] to map each surface to a disk. In order to
align the disks, we chose the highest point in each surface to map to the
north point of the disk. The radius of the disk can be selected arbitrarily,
and we chose one that leads to zero-mean distributions of area and metric
distortion.

###  III-C Evaluation

We evaluate the quality of the transformation by visualizing the mapped MRI
signals and by measuring the quality of the template match and distortion. We
visually assess the quality of our transformation by mapping the MRI intensity
patterns to the template coordinate system. We examine the flattened mesh and
anatomical landmarks of the mapped MRI signal in the flattened and original
spaces. We compare our method with the 2D baseline by quantifying differences
in areal distortion and visualizing the consistency in the mapped BOLD signal.

We define the average template matching error as the square-root of the
template matching term in Eq. (1), $(\sum_{m\in\partial
Z}A_{m}T(x_{m},\theta))^{\frac{1}{2}}$. We quantify shape distortion using the
final Dirichlet distortion term in Eq. (1), and compute the fraction over the
minimum of the symmetric Dirichlet energy for the identity transformation. We
also examine the distribution of the log-determinant of the Jacobian matrix
$\log_{2}\det\left(J\left(X_{k}\right)\right)$ of tetrahedron $k$ to quantify
local volumetric distortion [60]. For comparisons with the 2D baseline, we
compute area distortion using the ratio of triangle areas
$\log_{2}(\hat{A}_{m}/A_{m})$, where $A_{m}$ is the area of triangle $m$ in
the original mesh and $\hat{A}_{m}$ is the area of the same triangle in the
flattened mesh. We compute metric distortion using the ratio of edge length
$\log_{2}\left(x_{ij}/z_{ij}\right)$, where $z_{ij}$ is the length of edge
$(i,j)$ in the original mesh and $x_{ij}$ is the length of the same edge in
the flattened mesh. Local metric distortion is used to measure how far the
resulting map is from a local isometry and therefore how much it distorts the
local signal patterns. Similarly, a volume (3D) or area (2D) preserving map is
desirable to map the MRI signal intensities to the parameterized domain as
described in Section II-E.

###  III-D Results

####  III-D1 Visual Results

Fig. 4 illustrates the mapping results for four subjects. The algorithm
enables robust mapping of the highly variable placental shapes encountered in
the dataset. The ellipsoid template successfully unfolds the placenta,
although curved regions remain throughout the volume. The parallel planes
successfully unfolds the placentae and produces flat maternal and fetal sides.
Mapping to the single plane template maintains the curvature in the fetal side
of the organ while producing a flat maternal side. Mapping the placenta from
the same subject acquired in two scanning positions is visually more similar
in the canonical domain. The estimated placenta margin is connected and
encapsulates the curved region that separates the fetal and maternal sides,
highlighting the effectiveness of spectral clustering in parcellating the
placenta boundary. The parallel planes template facilitates visualization of
the maternal and fetal sides using the standard 2D cutting planes. Mapping
using the parallel planes template results in meshes with a flat geometry
which enables computation in the Euclidean domain. The remainder of our
experiments focus on the parallel planes template. We perform quantitative
comparisons of the three templates in Section III-D2. We provide guidelines
for template selection in Section IV.

Figure 4: Flattening results of tetrahedral meshes from three subjects. In
each subfigure, the top row corresponds to the mesh in the original volume;
the bottom row is the resulting flattened mesh. For each mesh using the
parallel planes or single plane template, the region in red represents the
placenta margin. The fetal side is facing upwards: (a) parallel planes; (b)
ellipsoid; and (c) single plane flattening in a singleton subject; (d)
parallel planes mapping in a twin subject. Sub-figures (e) and (f) demonstrate
mapping the placenta as acquired in the supine and left lateral positions in a
singleton subject using the parallel planes template. The ellipsoid template
results in an unfolded shape though multiple curved areas in the fetal and
maternal sides remain. The relaxation to the single plane template maintains
the curvature of the fetal side (c). The same placenta acquired in two scan
positions (e-f) appears visually more similar in the canonical template than
in the original volume.

Fig. 5 illustrates how the mapping enhances the visualization of landmarks.
The mapped BOLD intensity patterns clearly visualize local anatomy as is
apparent in the honeycomb structure of the cotyledons, which are the circular
structures that support the exchange of oxygen and nutrients between the
surrounding maternal blood and the fetal blood [1]. The cotyledons receive
highly-oxygenated maternal blood through the spiral arteries which appear
hyperintense in BOLD MRI, and receive less-oxygenated fetal blood through the
chorionic villi, which appear darker. In the flattened space, we observe this
effect as a darkening of intensity patterns closer to the fetal side. Our
mapping brings forth the contextual information about the cotyledon
distribution across the placenta, with consistent cotyledon structure
maintained as we move from the fetal to the maternal side of the organ. The
spatial relationship across cotyledons is difficult to observe in the original
volume while the flattened view provides spatial information relevant for
assessing placental function [4].

Figure 5:  Visualizing anatomical landmarks in one twin pregnancy subject,
moving depthwise from the fetal side (top) to the maternal side (bottom) of
the placenta. Left column: flattened slices with mapped MRI signals. Right
column: three orthogonal cross-sections of the original volume. Segmentation
of the placenta is outlined in yellow, and the red marker outlines the
approximate center of one cotyledon. The cotyledons, characterized by a
honeycomb structure of small hyperintense regions, are immediately apparent in
the flattened images. The spatial relationships between the marked cotyledon
and others are lost in the original images.

Fig. 6 demonstrates improved visualization of spatial relationships between
landmarks in a twin subject. The flattened view provides contextual
information surrounding the two marked cotyledons, e.g., the approximate
distance between them and a clear delineation of two regions. The curved
geometry of the in vivo placenta makes it difficult to visualize the
cotyledons together, losing the spatial relationships that are distinctly seen
in the flattened view. As shown in Fig. 5 and Fig. 6, several cross-sections
of the original image are required to identify the landmarks and their
surrounding context.

####  III-D2 Quantitative Results

Table I reports the mean template matching error and Dirichlet distortion for
each template. We observe sub-voxel accuracy of matching the template close to
the minimum of the symmetric Dirichlet energy. We observe similar shape
distortion between the ellipsoid and parallel planes templates, though the
ellipsoid template has a larger template matching error with greater variance
across subjects.

TABLE I: Template mismatch and Dirichlet distortion per template.

Cost Term (mean$\pm$std. dev.) | Ellipsoid | Parallel Planes | Single Plane
---|---|---|---
Dirichlet* (%) | $3.4\pm 2.6$ | $4.2\pm 2.0$ | $2.0\pm 1.0$
Template (voxels) | $0.45\pm 0.32$ | $0.26\pm 0.05$ | $0.26\pm 0.05$

*The Dirichlet distortion is reported as a percentage over the minimum Dirichlet energy as given by the identity transformation.

Fig. 7 reports the distributions of local volumetric and metric distortion
measures for the three templates evaluated. We report distributions of the
mean as the mean is sensitive to outliers and thus tight distributions of the
mean values indicate well behaved measures across all subjects. Both the
single plane and parallel planes templates yield narrow distributions of
distortion centered at $0$ with few outliers. The ellipsoid template exhibits
slight shrinking across the subjects and a wider interquartile range,
indicating difficulty in some placentae conforming to an ellipsoid template.
Given the large variability in placental shape, we observe that the parallel
planes template more accurately models the desired flattened shape.

The optimization results in a template that more accurately models the
placental thickness. We measure the absolute difference between the final and
initial value of the template half-height $h,|\Delta
h|=|h_{final}-h_{initial}|=1.5\pm 1.0\text{mm}$ for the parallel planes
template. We observe $h_{final}=13.1\pm 2.6\text{mm}$, which is in line with
the reported average placenta thickness after delivery of $28\pm 5\text{mm}$
[61].

We evaluate the robustness of our method to differences between singleton and
twin pregnancy cases. For the parallel planes template, we compare
distributions of local volumetric and metric distortion across the singleton
and twin subjects using Welch’s t-test to test if the distributions have the
same means. The hypothesis is not rejected for the distributions of volumetric
distortion ($p=0.16$) or metric distortion ($p=0.40$).

To illustrate the benefits of the flattened representation, we examine
placental shape differences between singleton ($n=60$) and twin ($n=18$)
pregnancies. Specifically, we compute the major and minor radii and the area
of the central cross-section, and the placental thickness in the flattened
representation. The major and minor radii of the placenta are found by fitting
the mesh vertices to an ellipse by minimizing the algebraic distance. We
perform an ANCOVA analysis while accounting for gestational age as a
covariate. Fig. 8 shows scatter plots of these quantities as a function of
gestational age. We found significant differences ($p<5\cdot\ 10^{-4}$) in the
major and minor radii and the area, while differences in placental thickness
were not significant. The same analysis using the fetal and maternal cross-
sections yielded similar results. This study suggests that the placenta has
larger extent but not substantially different thickness in twin pregnancies
when compared to singleton pregnancies.

We also evaluate the spatial distributions of distortion. Fig. 9 reports
volumetric distortion as a function of the vertex distance from the origin
(x-y plane) and vertex height from the $x^{(3)}=0$ plane, indicating
distortion is evenly distributed throughout the placenta volumes.

Figure 6:  Visual evaluation of the spatial relationships between two
cotyledons in a twin pregnancy subject. Left to right: (a) the original and
flattened tetrahedral meshes; (b) one slice of the flattened volume; (c) three
orthogonal cross-sections of the original image. The yellow outline marks the
placenta boundary. The red and green ticks mark the approximate centers of two
cotyledons. The views in the original image correspond to the center of the
red cotyledon. We observe the spatial relationship between the two cotyledons
in the flattened view clearly. Contextual information and spatial
relationships are difficult to visualize in the original image due to the
curvature of the uterine wall that determines the in vivo shape of the
placenta. Figure 7:  Box-and-whisker plots of (a) local volumetric
$\left(\log_{2}\det\left(J\right)\right)$ and (b) metric
$\left(\log_{2}\frac{x_{ij}}{z_{ij}}\right)$ distortions for the three
templates. We report distributions of the mean distortion across the 111
cases. Black horizontal lines inside the box indicate the median, the boxes
extend to the 25 and 75 percentiles, and the whiskers reach to the most
extreme values not considered outliers (black bars). The outliers, shown as
red crosses, are points farther than 1.5 times the interquartile range. Since
both distortion measures are defined as $log_{2}$ of ratios, $1$ corresponds
to doubling, and $-1$ corresponds to halving. The mean distortions are
centered near zero, with narrow interquartile ranges and few outliers for the
parallel planes and single plane templates. The ellipsoid template
demonstrates shrinking in the mean distortion with a wider interquartile
range. Figure 8: Comparison of placental shape measurements in twin and
singleton pregnancy scans: a)-c) major radius, minor radius, and the logarithm
of the area of the central cross-section respectively; d) placenta thickness.
Blue circles represent singleton scans and red triangles represent twin
pregnancies. Dashed blue lines indicate linear regression fit for the
singleton population. We conduct an ANCOVA analysis and find significant
differences across the singleton and twin scans for the major radius, minor
radius, and area, but not for the placenta thickness. Figure 9:  Local
volumetric distortion as a function of vertex distance: mean volumetric
distortion as a function of (a) vertex’s planar distance from the origin
$((x^{(1)})^{2}+(x^{(2)})^{2})^{\frac{1}{2}}$ and (b) absolute vertex height
from the $x^{(3)}=0$ plane. Solid line indicates mean distortion and shaded
areas extend to $\pm 1$ standard deviation across all vertices in the 111
cases. There are small variations in distortion as a function of distance,
indicating distortion is evenly distributed throughout the volumes. In (a), we
observe shrinking near the origin, which is expected, given the positively
curved placental shape.

####  III-D3 Comparison to 2D Baseline

Fig. 10 compares the distributions of areal distortion and the mapped
intensity patterns for four subjects. The baseline method of parameterizing to
a disk produces high distortion near the boundary, and introduces misalignment
of the flattened surfaces, as is visible by the through-plane images. The
baseline 2D mapping parameterizes the surfaces independently of each other,
lacks coupling across surfaces, and produces strong intensity artifacts. These
results support the need for a volumetric parameterization method.

Fig. 11 compares our method’s distributions of areal and metric distortion
with the baseline 2D approach. We achieve significantly lower distortion
across all cases. We conduct a two-sample Kolmogorov-Smirnov test to assess
the signficance. Our method produces significantly lower areal and metric
distortion per placenta ($p<10^{-24}$).

Figure 10:  Areal distortion and image mapping for our method and the baseline
2D method, demonstrated for four subjects: (a) distributions of areal
distortion on a central plane; (b) through-plane view of the mapped MRI values
in the parameterized space by our method (top) and the baseline 2D method
(bottom). Our method results in more homogeneous distributions of distortion.
The baseline severely distorts the placenta at the boundary. The through-image
view of the mapped intensities demonstrates intensity artifacts due to
independent parameterization of different planes.  Figure 11:  Areal and
metric distortion, for our method and the 2D baseline parameterization method.
The distributions are across the triangles (areal) and edges (metric) of all
parameterized surfaces.

##  IV Discussion

Our volumetric parameterization approach successfully maps the placenta to a
canonical template that closely resembles the ex vivo shape, thereby enhancing
visualization of anatomy and function. Experiments comparing cotyledon
distributions in the flattened and original views demonstrate that spatial
relations of anatomical landmarks are more clearly seen, and fewer slices are
required to span the placenta volume. Both local and global distortions are
successfully controlled while ensuring good template fit. Our algorithm is
robust to the greatly varying placental shapes encountered in the study,
successfully handling both twin and singleton cases. Our volumetric approach
results in significantly less distortion and ensures consistent
parameterization throughout the volume to a greater extent than that of a
baseline 2D parameterization approach. The significant differences in
distortion are due to the baseline approach fixing the placental boundary to a
disk. The natural placental boundary is often not disk-like, as can be seen in
Fig. 4, resulting in large distortion. Our approach parameterizes the placenta
to a free boundary, thereby maintaining the natural shape, and is therefore
better suited for placental parameterization.

We evaluate three templates, the ellipsoid, parallel planes, and single plane.
The proposed parallel planes template results in a parameterization with less
distortion compared to the ellipsoid. The single plane, which is a
modification of the parallel planes, results in slightly lower distortion. The
single plane template matches closely the ex vivo examination process as the
curvature in the fetal side is maintained. The parallel planes template is
more natural for visualization using standard 2D cutting planes. The
flexibility in the proposed approach enables the user to select a desired
template. We recommend the parallel planes template for most uses. In addition
to its advantages for visualization, the parallel planes template enables
computation in the Euclidean domain by producing flattened geometry within the
placenta. The single plane template minimizes local distortion and is
recommended for regional analysis within the volume of the placenta.

Our results demonstrate improved visualization of the BOLD MRI signals. Since
our model is geometry-based, it can also be used with other modalities,
provided the placenta segmentation is given. Additional modalities, such as
diffusion MRI [62] and T2-weighted MRI [8] are used to study placental
function and anatomy. This provides an opportunity to map and register several
signals to the common template for visualization and analysis.

An alternative and natural formulation of this work would be to use a physical
modeling approach, where the parameterization is modeled as a physical
deformation under gravity. Unfortunately, to the best of our knowledge, there
is minimal published work on the tissue properties of the placenta. The
symmetric Dirichlet energy is an elastic energy and is an appropriate proxy
for a physical deformation. The formulation of our template as two parallel
planes was inspired by a physical system, where the bottom plane (maternal
side) is the post-delivery examination table, and the top plane forces the
placenta to flatten. By optimizing for the template height, the optimization
results in a template that minimizes compressive and expansive forces on the
placenta.

###  IV-A Limitations and Future Work

A limitation of our template formulation is the invariance to rotation.
Rotations about the $x^{(3)}-$axis do not affect the value of the objective
function, which presents difficulty for aligning placentae, or the placenta of
the same subject over time. Additional registration may be required for proper
alignment. Furthermore, our template does not provide detailed correspondence.
For example, although we map the maternal side to the bottom plane, we are
unable to distinguish the left or right side of the placenta. This is a
challenging problem that will motivate the future methodological development.

Going forward, we plan to improve the estimation of the placenta margin.
Currently, we overestimate the margin to prevent distortion caused by
incorrectly assigning vertices to the maternal or fetal side. An improvement
is to refine the margin size throughout the optimization, converging on the
margin that estimates the placental thickness.

We plan to improve the template using anatomical landmarks such as the
umbilical cord insertion. The cord cannot be assumed to be inserted in the
center of the placenta, especially in twin cases, and must be identified in
the image. We are currently collecting anatomical images to identify landmarks
as they are not easily seen in the lower-resolution BOLD (functional) images.
The use of anatomical landmarks is necessary for the development of a common
coordinate system of the placenta. A potential approach to refining
correspondences between the placenta and the template would be to co-register
the anatomic and functional images prior to flattening, then compute spatial
statistics using mapped landmarks in the flattened space to enable comparative
studies.

We envision several potential clinical research studies enabled by this work.
The flattened view can be used to visualize regional MRI signal differences,
for example within cotyledons, the cord insertion, and the placental regions
responsible for each twin pair. Visualization of regional differences will
support quantification of signal differences, for example across discordant
twins with IUGR. Further, the canonical domain enables studies comparing shape
differences across populations as in Fig. 8. In Fig. 5 and Fig. 6, we
demonstrate visualization of the spatial relationships across cotyledons.
Assessing the signal variations across cotyledons is a promising direction to
study the possible conditions during pregnancy and after birth. For example,
the number of cotyledons at birth has been shown to be a predictor of blood
pressure in childhood [61].

In future work, we plan to quantify the similarity between the proposed
parameterization and the ex vivo placental shape using, for example, an ex
vivo system for placental perfusion and MR imaging [63]. By imaging a pregnant
mother near gestation and comparing the flattened placenta with the ex vivo
images, we can quantify the similarities between the two. A strong correlation
between the in vivo and ex vivo shapes will enable pathology sampling of
potential biomarkers identified during pregnancy. Furthermore, characteristics
of placental shape have been shown to correlate with pathology including
preeclampsia [64], and conditions that affect the fetus later in life such as
hypertension [65] and coronary heart disease [66].

##  V Conclusion

We developed a volumetric mesh-based mapping of the placenta to a flattened
template represented by two parallel planes, resembling its ex vivo flattened
shape. Our method successfully parameterizes the placenta with minimal
distortion and guarantees local injectivity. In a clinical study of placental
BOLD MRI images, we demonstrated that the mapping improves visualizing local
anatomy and function, and its potential for becoming a standardized method for
representing and visualizing placental images.

## References

  * [1] K. Benirschke and S. G. Driscoll, “The pathology of the human placenta,” in Placenta, pp. 97–571, Springer Berlin Heidelberg, 1967.
  * [2] E. Abaci Turk et al., “Placental MRI: Developing accurate quantitative measures of oxygenation,” Topics in Magnetic Resonance Imaging, vol. 28, no. 5, pp. 285–297, 2019.
  * [3] A. Sørensen et al., “Changes in human fetal oxygenation during maternal hyperoxia as estimated by BOLD MRI,” Prenatal Diagnosis, vol. 33, no. 2, pp. 141–145, 2013.
  * [4] J. Luo et al., “In vivo quantification of placental insufficiency by BOLD MRI: a human study,” Scientific Reports, vol. 7, no. 1, p. 3713, 2017.
  * [5] R. Pratt, J. Deprest, T. Vercauteren, S. Ourselin, and A. L. David, “Computer-assisted surgical planning and intraoperative guidance in fetal surgery: a systematic review,” Prenatal Diagnosis, vol. 35, no. 12, pp. 1159–1166, 2015.
  * [6] M. C. Schabel et al., “Functional imaging of the nonhuman primate placenta with endogenous blood oxygen level–dependent contrast,” Magnetic resonance in medicine, vol. 76, no. 5, pp. 1551–1562, 2016.
  * [7] E. Abaci Turk et al., “Placental MRI: Effect of maternal position and uterine contractions on placental bold MRI measurements,” Placenta, vol. 95, pp. 69 – 77, 2020.
  * [8] J. Torrents-barrena et al., “Fully automatic 3d reconstruction of the placenta and its peripheral vasculature in intrauterine fetal MRI,” Medical image analysis, vol. 54, pp. 263–279, 2019.
  * [9] A. Sørensen, J. Hutter, M. Seed, P. E. Grant, and P. Gowland, “T2*-weighted placental MRI: basic research tool or emerging clinical test for placental dysfunction?,” Ultrasound in Obstetrics & Gynecology, vol. 55, no. 3, pp. 293–302, 2020.
  * [10] A. E. Ho et al., “T2* placental magnetic resonance imaging in preterm preeclampsia: An observational cohort study,” Hypertension, vol. 75, no. 6, pp. 1523–1531, 2020.
  * [11] J. Hutter, P. J. Slator, L. Jackson, A. D. S. Gomes, A. Ho, L. Story, J. O’Muircheartaigh, R. P. A. G. Teixeira, L. C. Chappell, D. C. Alexander, M. A. Rutherford, and J. V. Hajnal, “Multi-modal functional MRI to explore placental function over gestation,” Magnetic Resonance in Medicine, vol. 81, no. 2, pp. 1191–1204, 2019.
  * [12] A. Melbourne, R. Aughwane, M. Sokolska, D. Owen, G. Kendall, D. Flouri, A. Bainbridge, D. Atkinson, J. Deprest, T. Vercauteren, A. David, and S. Ourselin, “Separating fetal and maternal placenta circulations using multiparametric MRI,” Magnetic Resonance in Medicine, vol. 81, no. 1, pp. 350–361, 2019.
  * [13] M. Pietsch, A. Ho, A. Bardanzellu, A. M. A. Zeidan, L. C. Chappell, J. V. Hajnal, M. Rutherford, and J. Hutter, “Applause: Automatic prediction of placental health via u-net segmentation and statistical evaluation,” Medical Image Analysis, vol. 72, p. 102145, 2021.
  * [14] M. Sinding, D. A. Peters, J. B. Frøkjaer, O. B. Christiansen, A. Petersen, N. Uldbjerg, and A. Sørensen, “Placental magnetic resonance imaging T2* measurements in normal pregnancies and in those complicated by fetal growth restriction,” Ultrasound in Obstetrics & Gynecology, vol. 47, no. 6, pp. 748–754, 2016.
  * [15] D. J. Gersell, “ASCP survey on placental examination,” American journal of clinical pathology, vol. 109, no. 2, pp. 127–143, 1998.
  * [16] J. Schreiner, A. Asirvatham, E. Praun, and H. Hoppe, “Inter-surface mapping,” ACM Trans. Graph., vol. 23, no. 3, pp. 870–877, 2004.
  * [17] B. Timsari and R. M. Leahy, “Optimization method for creating semi-isometric flat maps of the cerebral cortex,” in Proc.SPIE, Medical Imaging: Image Processing, vol. 3979, pp. 698–709, 2000.
  * [18] J. Smith and S. Schaefer, “Bijective parameterization with free boundaries,” ACM Trans. Graph., vol. 34, no. 4, pp. 70:1–70:9, 2015.
  * [19] M. S. Floater, “Mean value coordinates,” Computer Aided Geometric Design, vol. 20, no. 1, pp. 19–27, 2003.
  * [20] W. T. Tutte, “How to draw a graph,” Proceedings of the London Mathematical Society, vol. s3-13, no. 1, pp. 743–767, 1963.
  * [21] M. Rabinovich, R. Poranne, D. Panozzo, and O. Sorkine-Hornung, “Scalable locally injective mappings,” ACM Trans. Graph., vol. 36, no. 4, 2017.
  * [22] K. Hormann and G. Greiner, “MIPS: An efficient global parametrization method,” Curve and Surface Design: Saint-Malo 1999, pp. 153–162, 2000\.
  * [23] X.-M. Fu, Y. Liu, and B. Guo, “Computing locally injective mappings by advanced MIPS,” ACM Transactions on Graphics (TOG), vol. 34, no. 4, p. 71, 2015.
  * [24] P. Mullen, Y. Tong, P. Alliez, and M. Desbrun, “Spectral conformal parameterization,” vol. 27, pp. 1487–1494, 2008.
  * [25] B. Lévy, S. Petitjean, N. Ray, and J. Maillot, “Least squares conformal maps for automatic texture atlas generation,” ACM Trans. Graph., vol. 21, no. 3, pp. 362–371, 2002.
  * [26] M. Desbrun, M. Meyer, P. Schröder, and A. H. Barr, “Implicit fairing of irregular meshes using diffusion and curvature flow,” in Proceedings of the 26th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH ’99, pp. 317–324, ACM, 1999.
  * [27] H. Miao et al., “Placenta maps: in utero placental health assessment of the human fetus,” IEEE Transactions on Visualization and Computer Graphics, vol. 23, no. 6, pp. 1612–1623, 2017.
  * [28] S. Halier, S. Angenent, A. Tannenbaurn, and R. Kikinis, “Nondistorting flattening maps and the 3-d visualization of colon CT images,” IEEE Transactions on Medical Imaging, vol. 19, no. 7, pp. 665–670, 2000.
  * [29] S. Haker, S. Angenent, A. Tannenbaum, and R. Kikinis, “Non-distorting flattening for virtual colonoscopy,” in Medical Image Computing and Computer-Assisted Intervention – MICCAI 2000, pp. 358–366, 2000.
  * [30] M. Desbrun, M. Mark, and A. Pierre, “Intrinsic parameterizations of surface meshes,” Computer Graphics Forum, vol. 21, no. 3, pp. 209–218, 2002.
  * [31] C. Gotsman, X. Gu, and A. Sheffer, “Fundamentals of spherical parameterization for 3D meshes,” ACM Trans. Graph., vol. 22, no. 3, pp. 358–363, 2003\.
  * [32] S. Haker, S. Angenent, A. Tannenbaum, R. Kikinis, G. Sapiro, and M. Halle, “Conformal surface parameterization for texture mapping,” IEEE Transactions on Visualization and Computer Graphics, vol. 6, no. 2, pp. 181–189, 2000.
  * [33] L. M. Lui, Y. Wang, T. F. Chan, and P. Thompson, “Landmark constrained genus zero surface conformal mapping and its application to brain mapping research,” Applied Numerical Mathematics, vol. 57, no. 5-7, pp. 847–858, 2007.
  * [34] B. Fischl, M. I. Sereno, and A. M. Dale, “Cortical surface-based analysis: II: inflation, flattening, and a surface-based coordinate system,” Neuroimage, vol. 9, no. 2, pp. 195–207, 1999.
  * [35] B. Fischl, “Freesurfer,” Neuroimage, vol. 62, no. 2, pp. 774–781, 2012\.
  * [36] D. Tosun and J. L. Prince, “Hemispherical map for the human brain cortex,” in Proc.SPIE, Medical Imaging: Image Processing, vol. 4322, pp. 290–301, 2001\.
  * [37] A. Sheffer, E. Praun, and K. Rose, “Mesh parameterization methods and their applications,” Foundations and Trends in Computer Graphics and Vision, vol. 2, no. 2, pp. 105–171, 2007.
  * [38] J. Kreiser, M. Meuschke, G. Mistelbauer, B. Preim, and T. Ropinski, “A survey of flattening-based medical visualization techniques,” in Computer Graphics Forum, vol. 37, pp. 597–624, 2018.
  * [39] S. Oeltze, A. Kuß, F. Grothues, A. Hennemuth, and B. Preim, “Integrated visualization of morphologic and perfusion data for the analysis of coronary artery disease.,” in EuroVis, pp. 131–138, 2006.
  * [40] M. Termeer, J. O. Bescós, M. Breeuwer, A. Vilanova, F. Gerritsen, and E. Gröller, “Covicad: Comprehensive visualization of coronary artery disease,” IEEE Transactions on Visualization and Computer Graphics, vol. 13, no. 6, pp. 1632–1639, 2007.
  * [41] H. Martinke, C. Petry, S. Großkopf, M. Suehling, G. Soza, B. Preim, and G. Mistelbauer, “Bone Fracture and Lesion Assessment using Shape-Adaptive Unfolding,” in Eurographics Workshop on Visual Computing for Biology and Medicine (S. Bruckner, A. Hennemuth, B. Kainz, I. Hotz, D. Merhof, and C. Rieder, eds.), 2017.
  * [42] T. Vrtovec, S. Ourselin, L. Gomes, B. Likar, and F. Pernuš, “Automated generation of curved planar reformations from mr images of the spine,” Physics in Medicine & Biology, vol. 52, no. 10, p. 2865, 2007.
  * [43] K. Hormann, K. Polthier, and A. Sheffer, “Mesh parameterization: Theory and practice,” in ACM SIGGRAPH ASIA 2008 Courses, SIGGRAPH Asia ’08, 2008.
  * [44] M. S. Floater and K. Hormann, “Surface parameterization: a tutorial and survey,” in Advances in Multiresolution for Geometric Modelling, pp. 157–186, Springer Berlin Heidelberg, 2005.
  * [45] M. Floater, “One-to-one piecewise linear mappings over triangulations,” Mathematics of Computation, vol. 72, no. 242, pp. 685–696, 2003.
  * [46] N. Aigerman, S. Z. Kovalsky, and Y. Lipman, “Spherical orbifold tutte embeddings,” ACM Transactions on Graphics (TOG), vol. 36, July 2017.
  * [47] D. Ezuz, J. Solomon, and M. Ben-Chen, “Reversible harmonic maps between discrete surfaces,” ACM Transactions on Graphics (TOG), vol. 38, no. 2, pp. 1–12, 2019.
  * [48] N. Aigerman and Y. Lipman, “Injective and bounded distortion mappings in 3d,” ACM Transactions on Graphics (TOG), vol. 32, no. 4, p. 106, 2013.
  * [49] M. S. Floater and V. Pham-Trong, “Convex combination maps over triangulations, tilings, and tetrahedral meshes,” Advances in Computational Mathematics, vol. 25, no. 4, pp. 347–356, 2006.
  * [50] S. M. Abulnaga, E. Abaci Turk, M. Bessmeltsev, P. E. Grant, J. Solomon, and P. Golland, “Placental flattening via volumetric parameterization,” in Medical Image Computing and Computer Assisted Intervention – MICCAI 2019, pp. 39–47, 2019.
  * [51] A. Y. Ng, M. I. Jordan, and Y. Weiss, “On spectral clustering: Analysis and an algorithm,” in Advances in Neural Information Processing Systems, pp. 849–856, 2002.
  * [52] U. Von Luxburg, “A tutorial on spectral clustering,” Statistics and Computing, vol. 17, no. 4, pp. 395–416, 2007.
  * [53] E. W. Dijkstra, “A note on two problems in connexion with graphs,” Numerische mathematik, vol. 1, no. 1, pp. 269–271, 1959.
  * [54] J. Nocedal and S. Wright, Numerical optimization.  Springer, 2006.
  * [55] A. F. Möbius, “Der barycentrische calcul,” Verlag von Johann Ambrosius Barth, 1827.
  * [56] Q. Fang and D. A. Boas, “Tetrahedral mesh generation from volumetric binary and grayscale images,” in 2009 IEEE International Symposium on Biomedical Imaging: From Nano to Macro, pp. 1142–1145, 2009.
  * [57] H. Si, “Tetgen, a Delaunay-based quality tetrahedral mesh generator,” ACM Transactions on Mathematical Software (TOMS), vol. 41, no. 2, pp. 11:1–11:36, 2015.
  * [58] E. A. Turk et al., “Spatiotemporal alignment of in utero BOLD-MRI series,” Journal of Magnetic Resonance Imaging, vol. 46, no. 2, pp. 403–412, 2017.
  * [59] P. Joshi, M. Meyer, T. DeRose, B. Green, and T. Sanocki, “Harmonic coordinates for character articulation,” ACM Trans. Graph., vol. 26, no. 3, 2007.
  * [60] A. D. Leow et al., “Statistical properties of Jacobian maps and the realization of unbiased large-deformation nonlinear image registration,” IEEE Transactions on Medical Imaging, vol. 26, no. 6, pp. 822–832, 2007\.
  * [61] D. Barker, C. Osmond, S. Grant, K. L. Thornburg, C. Cooper, S. Ring, and G. Davey-Smith, “Maternal cotyledons at birth predict blood pressure in childhood,” Placenta, vol. 34, no. 8, pp. 672–675, 2013.
  * [62] P. J. Slator et al., “Combined diffusion-relaxometry MRI to identify dysfunction in the human placenta,” Magnetic resonance in medicine, vol. 82, no. 1, pp. 95–106, 2019.
  * [63] J. N. Stout, S. Rouhani, E. A. Turk, C. G. Ha, J. Luo, K. Rich, L. L. Wald, E. Adalsteinsson, W. H. Barth Jr, P. E. Grant, et al., “Placental MRI: Development of an MRI compatible ex vivo system for whole placenta dual perfusion,” Placenta, vol. 101, pp. 4–12, 2020.
  * [64] E. Kajantie, K. Thornburg, J. G. Eriksson, L. Osmond, and D. J. Barker, “In preeclampsia, the placenta grows slowly along its minor axis,” International Journal of Developmental Biology, vol. 54, no. 2-3, pp. 469–473, 2009.
  * [65] D. J. Barker, K. L. Thornburg, C. Osmond, E. Kajantie, and J. G. Eriksson, “The surface area of the placenta and hypertension in the offspring in later life,” The International journal of developmental biology, vol. 54, p. 525, 2010.
  * [66] J. G. Eriksson, E. Kajantie, K. L. Thornburg, C. Osmond, and D. J. Barker, “Mother’s body size and placental size predict coronary heart disease in men,” European heart journal, vol. 32, no. 18, pp. 2297–2303, 2011.

## Appendix A

###  A-A Jacobian Formulation

The Jacobian defined in Eq. (2) defines the transformation of tetrahedron $k$
using the three basis vectors $X_{k}B$ defining the tetrahedron. The matrix
$B\in\mathbb{R}^{3\times 4}$

$B=\begin{bmatrix}\ -1&-1&-1\\\ 1&0&0\\\ 0&1&0\\\ 0&0&1\end{bmatrix}$

is a linear operator that extracts the basis vectors. Matrix $X_{k}$
represents the vertices of tetrahedron $k$,
$\\{x_{i},x_{i+1},x_{i+2},x_{i+3}\\}$ as column vectors. The basis used to
compute the Jacobian is then

$X_{k}B=\begin{bmatrix}|&|&|\\\ x_{i+1}-x_{i}&x_{i+2}-x_{i}&x_{i+3}-x_{i}\\\
|&|&|\end{bmatrix}.$

###  A-B Gradient of the Objective Function

The gradient of the objective function in Eq. (1) per vertex is computed as
the sum of the gradient of the template and the distortion terms. The gradient
of the parallel planes template term with respect to the boundary vertex $m$
is linear,

$\frac{\partial T(x_{m},h)}{\partial
x_{m}}=\begin{cases}2(x_{m}^{(3)}-h)&\text{if }m\in\mathcal{F}(\partial Z),\\\
2(x_{m}^{(3)}+h)&\text{if }m\in\mathcal{M}(\partial Z),\\\
0&\text{otherwise}.\end{cases}$

The gradient of the ellipsoid template term with respect to boundary vertex
$m$ is

$\frac{\partial T(x_{m},R)}{\partial x_{m}}=4A_{m}(x_{m}^{T}Rx_{m}-1)Rx_{m}.$

The gradient of the distortion term for tetrahedron $k$ is:

$\displaystyle\frac{d\,\|J(X_{k})\|_{F}^{2}}{d\,X_{k}}$
$\displaystyle=X_{k}\left[2B\left(Z_{k}B\right)^{-1}\left(Z_{k}B\right)^{-T}B^{T}\right],$
$\displaystyle\begin{split}\frac{d\,\|J(X_{k})^{-1}\|_{F}^{2}}{d\,X_{k}}&=-\left[2X_{k}B\left(B^{T}X_{k}^{T}X_{k}B\right)^{-1}B^{T}Z_{k}^{T}Z_{k}B\right]\times\\\
&\quad\quad\quad\quad\left[\left(B^{T}X_{k}^{T}X_{k}B\right)^{-1}B^{T}\right].\end{split}$

We also optimize over the template height. For the two parallel planes
template, the corresponding derivative is

$\frac{\partial\phi(X,h)}{\partial h}=-2\left[\sum_{j\in\mathcal{F}(\partial
Z)}A_{j}\left(x_{j}^{(3)}-h\right)-\sum_{j\in\mathcal{M}(\partial
Z)}A_{j}\left(x_{j}^{(3)}+h\right)\right].\\\ $

For the ellipsoid template, the corresponding derivative with respect to the
parameter $r_{x}$ is

$\frac{\partial\phi(X,R)}{\partial
r_{x}}=-\frac{4}{r_{x}^{3}}\sum_{m\in\partial
Z}A_{m}\left(x_{m}^{T}Rx_{m}-1\right).\\\ $

An equivalent form is found for $r_{y}$ and $r_{z}$.