Placental Flattening via Volumetric Parameterization

TL;DR

This paper introduces a volumetric mesh-based algorithm to flatten placental shapes from MRI data into a canonical template, improving visualization of placental anatomy and function for better pregnancy assessment.

Contribution

The novel method maps in vivo placental shapes to a flattened template using symmetric Dirichlet energy minimization with local injectivity constraints.

Findings

Achieved sub-voxel accuracy in boundary mapping

Controlled distortion throughout the volume

Enhanced visualization of placental anatomy

Abstract

We present a volumetric mesh-based algorithm for flattening the placenta to a canonical template to enable effective visualization of local anatomy and function. Monitoring placental function in vivo promises to support pregnancy assessment and to improve care outcomes. We aim to alleviate visualization and interpretation challenges presented by the shape of the placenta when it is attached to the curved uterine wall. To do so, we flatten the volumetric mesh that captures placental shape to resemble the well-studied ex vivo shape. We formulate our method as a map from the in vivo shape to a flattened template that minimizes the symmetric Dirichlet energy to control distortion throughout the volume. Local injectivity is enforced via constrained line search during gradient descent. We evaluate the proposed method on 28 placenta shapes extracted from MRI images in a clinical study of placental function. We achieve sub-voxel accuracy in mapping the boundary of the placenta to the template while successfully controlling distortion throughout the volume. We illustrate how the resulting mapping of the placenta enhances visualization of placental anatomy and function. Our code is freely available at https://github.com/mabulnaga/placenta-flattening .

Equations

$$\phi(X)=\mathop{\mathchoice{\vtop{\halign{#\cr$\hfil\displaystyle\sum_{m=1}^{M}A_{m}T\left(x_{m}\right)\hfil$\crcr\kern 2.0pt\nointerlineskip\cr\@mathmeasure\displaystyle{\sum_{m=1}^{M} A_{m} T\left(x_{m}\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=1}^{M}A_{m}T\left(x_{m}\right)\hfil$\crcr\kern 2.0pt\nointerlineskip\cr\@mathmeasure\textstyle{\sum_{m=1}^{M} A_{m} T\left(x_{m}\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=1}^{M}A_{m}T\left(x_{m}\right)\hfil$\crcr\kern 2.0pt\nointerlineskip\cr\@mathmeasure\scriptstyle{\sum_{m=1}^{M} A_{m} T\left(x_{m}\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=1}^{M}A_{m}T\left(x_{m}\right)\hfil$\crcr\kern 2.0pt\nointerlineskip\cr\@mathmeasure\scriptscriptstyle{\sum_{m=1}^{M} A_{m} T\left(x_{m}\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}},$$

$$\mathcal{D}(J)=\Vert J{\Vert }_F^2+\Vert J^{-1}{\Vert }_F^2=\sum _{i=1}^3\left({\sigma }_i^2+{\sigma }_i^{-2}\right),$$

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

$$\frac{d\Vert J(X_k){\Vert }_F^2}{d{X}_k}$$

$$\frac{d\Vert J(X_k{)}^{-1}{\Vert }_F^2}{d{X}_k}$$

Full text

1 Introduction

The placenta is a critical organ that connects the fetus to the maternal blood system. Placental dysfunction increases the risk of pregnancy complications, with long-term effects on a child’s health and development and the mother’s health. It is therefore critical to monitor placental function and health in vivo. Ultrasound and MRI capture detailed information about the placental position, shape, and tissue properties [11]. Blood oxygen level dependent (BOLD) MRI has recently been demonstrated to assess oxygen transport within the placenta [11, 5], providing initial evidence for clinical utility of MRI for functional assessment of the placenta. Compared with ultrasound, MRI provides direct measurements of placental function, which provides signals necessary to study function and assess pathology [11, 5].

The in vivo shape of the placenta is determined by the curved surface of the uterine wall to which it is attached during pregnancy. This presents significant challenges for interpretation of the MRI scans. No common standard exists for visualizing the functional or anatomical images of the placenta whose in vivo shape and location of attachment to the uterine wall vary greatly across subjects. We present a novel algorithm for mapping the placental shape observed in an MRI scan to a flattened template that resembles the organ’s well-studied ex vivo flattened shape, to alleviate visualization challenges during in vivo examination and to facilitate clinical research and development of placental health biomarkers. Our work offers the first step toward developing a common coordinate system to enable statistical analysis.

We build on state-of-the-art mesh parameterization methods to represent and estimate the deformation of the placenta onto a template. Mesh parameterization is a topic of active research in geometry processing for mapping surfaces to canonical domains such as planes or spheres while guaranteeing desirable properties of the mapping such as injectivity [8, 10]. When applied to cortical mapping, this parameterization facilitates visualization and population studies [2, 12, 13]. The formulation may estimate the optimal map by minimizing a cost function that penalizes areal distortion [13], changes in geodesic distance [2], or more commonly a combination of different distortion measures [12]. In contrast to the inherently two-dimensional (2D) cortex however, the placenta is a fully three-dimensional (3D) organ. Image information along the depth direction from the maternal to the fetal side of the placenta is important for characterizing function. The parameterization must therefore map the entire volume.

In placenta imaging, a method for placental flattening has been recently demonstrated [6]. Their algorithm represents the placenta as a stack of parallel surfaces spanning the thickness of the organ. Each surface is flattened separately by mapping the boundary to a disk with each interior vertex moved to the average of its neighbors. This approach is limited in two ways. First, the lack of correspondence across surfaces throughout the volume results in through-plane artifacts, distorting important depth-wise information. Second, due to the variability in placental shape, constraining the mapping to a fixed boundary results in high distortion. In contrast, our method computes a continuous volumetric mapping with free boundary that ensures uniform consistency throughout the volume and enables explicit control of the resulting deformation.

To our knowledge, this is the first volumetric approach for mapping the placenta to a canonical domain. Our algorithm estimates the transformation of the volumetric mesh to a template as the solution of an optimization problem. This formulation readily accepts a broad family of templates and shape distortion functions. We choose to minimize the symmetric Dirichlet energy [9, 10] to penalize local deformations of the volumetric mesh. We evaluate our method on images from a clinical MRI study, demonstrating effective mapping of the highly variable placental shape to the template with minimal distortion. We demonstrate improved visualization of anatomical structures and their surrounding context, illustrating the promise of our algorithm to support clinical use of MRI in placental imaging.

2 Methods

We represent the placental shape as a tetrahedral mesh that contains $N$ vertices and $K$ tetrahedra. In our experiments, we extract such meshes from segmented MRI scans as described later. We parameterize the mapping via mesh vertex locations in the template coordinate system and interpolate the deformation to the interior of each tetrahedron using a locally affine (piecewise-linear) model. The mapping is geometry-based so it is independent of imaging modality used.

2.1 Problem Formulation

Let $X\in\mathbb{R}^{3\times N}$ be a matrix whose columns are the 3D coordinates of all mesh vertices in the template space with the $M$ boundary vertices forming the first $M$ columns. Let $X_{k}\in\mathbb{R}^{3\times 4}$ be a matrix whose columns are the 3D coordinates of the four corner vertices of tetrahedron $k$ ($k=1,\ldots,K$) in the template space. We formulate the mapping as an optimization problem over the mesh vertices that seeks to map to the template space while minimizing shape distortion. We formulate the general objective function as

[TABLE]

where $\{x_{m}\}_{m=1}^{M}$ are the boundary vertex coordinates in the template system, $T\left(\cdot\right)$ is a measure of distance from the template shape, $A_{m}$ is the normalized barycentric area of boundary vertex $m$ in the original space, $\mathcal{D}\left(\cdot\right)$ measures local distortion, $V_{k}$ is the normalized volume of tetrahedron $k$ in the original space, and $\lambda$ is a parameter that governs the trade-off between the template fit and the shape distortion. The distortion term regularizes the mapping.

We use a locally affine model to capture the deformation of tetrahedron $Z_{k}$ in the original image space to $X_{k}$ in the template space. The Jacobian matrix $J(X_{k})=\left(X_{k}B\right)\left(Z_{k}B\right)^{-1}$ captures the linear transformation of 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.

We measure local distortion using the symmetric Dirichlet energy density

[TABLE]

where $\|\cdot\|_{F}$ is the Frobenius norm and $\{\sigma_{1},\sigma_{2},\sigma_{3}\}$ are the singular values of matrix $J$ [8, 10]. We chose the symmetric Dirichlet energy since it penalizes expansion and shrinking equally and favors a locally-injective mapping.

The image intensities are mapped to the template space using barycentric coordinates (BC). The BC represent a voxel in a tetrahedron as a convex combination of the tetrahedon’s vertices. Since our mapping is injective and affine, the resulting voxel position is determined using its BC and the mapped vertices.

2.2 Template

After evaluating several volumetric templates (ellipsoid, sphere, cylinder; not shown), we find that the best results (i.e., clear anatomical structure and small local distortion) are achieved by encouraging constant height of the flattened placenta to mimic the post-delivery examination process, where the maternal side is placed on an examination table. The fetal side is flattened to the ease visualization. The function $T(\cdot)$ measures the distance to the appropriate plane:

[TABLE]

where $x^{(3)}$ refers to the third coordinate of point $x$ in the template coordinate system, $\partial Z$ denotes the mesh boundary in the original image space, and $\mathcal{F}(\partial Z)$, $\mathcal{M}(\partial Z)$ denote the fetal and maternal sides of $\partial Z$. We identify the maternal and fetal sides via spectral clustering [7] as described below.

We employ a similarity metric based on the angle between unit normals of boundary vertices ($\hat{n}$). We construct an affinity matrix $W\in\mathbb{R}^{M\times M}$ whose ($i,j$) element $w_{i,j}=\exp\left\{\gamma\left(\hat{n}_{i}^{\top}\hat{n}_{j}\right)\right\}$ for any two boundary vertices $i$ and $j$ that are connected by a path of at most $3$ edges (the $3$-ring neighborhood), and $w_{i,j}=0$ otherwise. The parameter $\gamma$ penalizes the variation in the orientation of the normals. We cluster the boundary vertices by thresholding the values of the second smallest eigenvector of the Laplacian $L=I-D^{-\frac{1}{2}}WD^{-\frac{1}{2}}$, where $D$ is a diagonal matrix with $d_{i,i}=\sum_{j}w_{i,j}$ and $I$ is the identity matrix [7]. Since the maternal side is more curved, we assign the corresponding label to the cluster with the larger number of vertices on the convex hull of the mesh.

We use the term rim to denote the highly curved region that separates the fetal and maternal sides. We first assign to the rim all vertices on the boundary of the two clusters, i.e., those with neighbors in the other cluster. The rim is then dilated to a set width based on the approximated geodesic distance along the mesh boundary. The geodesic distance accounts for mesh irregularities and ensures a consistent rim width.

2.3 Optimization

We minimize cost function $\phi(\cdot)$ in (1) using gradient descent. We initialize the mapping using the identity transformation. The gradient of the template term is linear in the vertices, $\frac{\partial T(x)}{\partial x^{(3)}}=2\left(x^{(3)}\pm h\right)$. We derive the gradient of the symmetric Dirichlet energy term defined in (2) using the chain rule for matrices:

[TABLE]

We employ line search to prevent tetrahedra from “flipping,” i.e., from crossing the singularity point of zero volume, thereby enforcing local injectivity [10]. In every iteration, we determine the largest value $\eta$ such that adjusting the current vertex locations $X$ by $-\eta\nabla\phi(X)$ avoids singularities for all tetrahedra. The (signed) volume of tetrahedron $k$ is computed as the determinant of matrix $\left(X_{k}-\eta\nabla\phi\left(X_{k}\right)B\right)$ and is a cubic polynomial of $\eta$. The smallest, positive real root provides the upper limit for $\eta$ in the line search [10].

2.4 Implementation Details

We generate tetrahedral meshes from segmentation labelmaps using iso2mesh[1]. Prior to mapping, we center the mesh and rotate it to align its principal axes with the template. We assume the algorithm convergences when the Frobenius norm of the gradient is lower than $1\times 10^{-4}$. We implemented the algorithm in MATLAB using GPU functionality to parallelize computation, and ran our experiments on an NVIDIA Titan Xp (12GB) GPU. The algorithm took an average of $4083$ iterations and took less than $20$ minutes to converge. Our implementation is freely available at https://github.com/mabulnaga/placenta-flattening.

3 Experiments

**Data: ** We validate the approach on a set of 28 MRI scans from 2 studies. The first is a twin study on 7 pregnant women (gestational age (GA): 28–34 weeks). All twin pregnancies had one placenta shared by the 2 twins. The second is a singleton pregnancy study on 11 women (GA: 27 – 40 weeks). For 10 of 11 subjects, scans were acquired in the supine and left lateral positions, yielding 20 different segmentations. MRI BOLD scans were acquired on a 3T Siemens Skyra scanner (GRE-EPI, 3 $\text{mm}^{3}$, TR=5.8–8s, TE=32–36ms, FA= 90∘). The placenta was manually segmented by a trained observer and input to the meshing software, which produced 6,500 tetrahedra and 2,800 surface triangles on average.

**Parameters: ** We used a grid search to determine the values of the hyper-parameters. We set shape distortion parameter $\lambda=1$ as it was in the optimal trade-off range between the template match and distortion (Fig. 1a). We set the template half-height $h$ to be half of the placenta thickness estimated from the histogram of the distance transform values inside the segmentation boundary. We set the spectral clustering parameter $\gamma=20$ and used a boundary geodesic distance of $5$ voxels as the width of the rim.

Evaluation: We use 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 [4]. We quantify metric distortion using the ratio of edge lengths $\log_{2}\left(x_{ij}/z_{ij}\right)$. We visually assess the quality of the transformation by mapping the BOLD MRI to the template coordinate system.

We compare with the prior parameterization approach in [6] where 2D surfaces spanning the placenta were independently parameterized to a disk in $\mathbb{R}^{2}$. The surfaces were derived by cutting Euclidean level sets to minimize local curvature changes. To emulate this method, we derive such 2D surfaces by intersecting the flattened placenta volume with planes spaced one voxel apart. We harmonically parameterize each surface to a disk [3], each with a corresponding point mapped to the north point of the disk. We scale the areas and edge lengths in the parameterized space to have mean [math] areal and metric log-distortion.

3.1 Results

For all cases, our algorithm achieves sub-voxel accuracy of matching the template (median of $0.09$  voxels, max. of $0.30$ voxels). The resulting transformations achieve close to minimum values of the symmetric Dirichlet energy (median $3.87\%$  higher, max. $9.04\%$ higher than the smallest possible value), where the energy is minimized by the identity transformation. Fig. 1b demonstrates the mapping achieves minimal local volumetric and metric distortions. We did not observe differences across twin and singleton pregnancies.

Fig. 2 illustrates the mapping results for four placentae, highlighting the variability in shape encountered in the dataset. We are able to map difficult structures robustly such as the fold in the first placenta and bowl-shape in the last. The estimated rim also effectively separates the maternal and fetal regions.

Fig. 3 demonstrates that the mapping enhances visualization of landmarks. The mapped BOLD intensity patterns clearly visualize local anatomy and function as is apparent in the honeycomb structure of the cotyledons, which are the circular structures that exchange oxygen and nutrients between the maternal blood and the fetal blood in the chorionic villi [5]. Cotyledons appear hyperintense in BOLD MRI. Similarly, contextual information is lost when visualizing a fabricated landmark (the letters “MIC”) in the curved volume. Since the inherent geometry of the placenta is flat, the letters could represent biomarkers. The curved geometry of the in vivo placenta demonstrates immediate difficulty in visualizing the details that are clearly seen in the flattened view. Several views in the original image space are required to identify anatomical landmarks.

Fig. 4 compares our method with the baseline 2D approach. Our method produces considerably lower distortion across all cases. Mapping to a fixed boundary results in higher distortion, and we observe image artifacts due to the lack of coupling across planes. These results confirm the need for a free boundary volumetric parameterization method. Finally, we note that the nested surfaces in the original space were derived from the volumetric parameterization to the flattened space. In practice, estimating these surfaces is challenging and requires a multi-step specialized pipeline as in [6].

4 Conclusion

We developed a volumetric mesh-based mapping of the placenta to a flattened template represented by two parallel planes, resembling the ex vivo flattened shape of the organ, thereby enabling visualization of local anatomy and function. An immediate next step is to assess the utility for clinical research by quantifying the improvement in identifying known key biomarkers and anatomical features. In future work, we will improve the template using anatomical data such as the umbilical cord insertion. We are collecting higher-resolution anatomical images to include landmarks that are not easily seen in the BOLD (functional) images. This work is the first step towards developing a common coordinate system to visualize, examine, and study the organ as well as to support statistical analysis across subjects and time. Such a framework promises to advance the state of the art in studies of the placenta and to provide MRI biomarkers of fetal health.

Acknowledgments

This work was supported in part by NIH NIBIB NAC P41EB015902, NIH NICHD U01HD087211, NSF IIS-1838071, Air Force FA9550-19-1-0319, Wistron, SIP, AWS, NSF Graduate Research Fellowship, and NSERC Post Graduate Scholarship.