3D Reconstruction and Segmentation of Dissection Photographs for MRI-free Neuropathology

Neuroimaging to neuropathology correlation (NTNC) promises to enable the transfer of microscopic signatures of pathology to in vivo imaging with MRI, ultimately enhancing clinical care. NTNC traditionally requires a volumetric MRI scan, acquired either ex vivo or a short time prior to death. Unfortunately, ex vivo MRI is difficult and costly, and recent premortem scans of sufficient quality are seldom available. To bridge this gap, we present methodology to 3D reconstruct and segment full brain image volumes from brain dissection photographs, which are routinely acquired at many brain banks and neuropathology departments. The 3D reconstruction is achieved via a joint registration framework, which uses a reference volume other than MRI. This volume may represent either the sample at hand (e.g., a surface 3D scan) or the general population (a probabilistic atlas). In addition, we present a Bayesian method to segment the 3D reconstructed photographic volumes into 36 neuroanatomical structures, which is robust to nonuniform brightness within and across photographs. We evaluate our methods on a dataset with 24 brains, using Dice scores and volume correlations. The results show that dissection photography is a valid replacement for ex vivo MRI in many volumetric analyses, opening an avenue for MRI-free NTNC, including retrospective data. The code is available at https://github.com/htregidgo/DissectionPhotoVolumes.

definitive diagnoses can only be obtained via neuropathology. To overcome this, neuroimaging to neuropathology correlational (NTNC) science seeks to establish imaging phenotypes that correlate with gold standard pathological diagnoses, in order to port these signatures to in vivo imaging as biomarkers. One candidate for NTNC is in vivo or ex vivo MRI. Reliable matching of histology and in vivo MRI requires a premortem scan acquired a short time before death. Unfortunately, these are difficult to obtain for precisely the most interesting individuals -asymptomatic, early-stage cases. This problem can be overcome with ex vivo MRI, which has been successfully used in NTNC (e.g., [7,13]), but is also challenging to perform: it requires scanning and sample preparation expertise that is not present at many research centres, cannot easily be done on the frozen tissue required in many genetics analyses, and is expensive.
Meanwhile, a wealth of information exists in brain banks that is hidden in existing images from routine dissection photography. Here we present algorithms to 3D reconstruct and segment imaging volumes from this underutilised modality, enabling morphometric NTNC studies without MRI at almost no cost.
Related work: Building 3D images from dissection photographs requires alignment of a stack of 2D photographs into a 3D consistent volume via image registration [15,27]. Registration of image pairs is a well-studied problem but, to the best of our knowledge, literature on joint registration of dissection photographs for 3D reconstruction is nonexistent. The closest related work is a method for volumetric reconstruction from printed films of MRI [9], which is not suitable for our task, as it requires a reference MRI volume (which we wish to avoid).
One step further removed is 3D histology reconstruction [20]. Despite the peculiarities of histological data in terms of contrast, resolution, and sectioning distortions, many of the challenges we face in this work are similar. Without an additional reference, recovering the 3D shape of a 2D stack of images is a heavily underconstrained problem. A common approach is to iteratively align each 2D image to its neighbours, possibly with an outlier rejection strategy [31]. This approach yields smooth volumes but can result in straightening of curved structures (often known as "banana effect" [30]), and accumulation of errors along the stack (z-shift, [19]). These can both be overcome with a reference MRI scan [2,1] -a requirement which, again, we are trying to avoid. For brain segmentation, the neuroimaging literature has long been dominated by multi-atlas segmentation (MAS), Bayesian segmentation and, more recently, deep convolutional neural networks (CNNs). MAS [23,11] nonlinearly registers several labelled atlases to a target scan, deforms the corresponding segmentations, and merges these warped label maps into a robust estimate of the segmentation with a label fusion algorithm. Segmentation CNN architectures [16,12,14], best represented by the ubiquitous U-Net [24], yield state-of-the-art accuracy and runtimes (seconds). Being supervised methods, MAS and CNNs share the disadvantage that performance quickly decreases when the training and test domains do not match. Despite progress in data augmentation [33] and transfer learning [26], manual labels are often needed for every new segmentation task.

3D Reference volume
Input photos Masked photos 3D reconstruction Segmentation Fig. 1. Diagram of the proposed processing pathway for dissection photography. Photographs are coarsely corrected for pixel size and perspective; arranged into a volume; registered to each other and a reference surface mask; and automatically segmented.
Bayesian segmentation with probabilistic atlases uses a generative model combining a supervised prior model of anatomy (the atlas) and a model of image formation (likelihood). Segmentation is then posed as a Bayesian inference problem, estimating the most probable hidden segmentation that generated the observed image, given the atlas. A subset of Bayesian segmentation methods [29,32,4,21] use an unsupervised likelihood model, usually a Gaussian mixture model (GMM) whose parameters are estimated specifically for each volume to segment, making the segmentation adaptive to different contrasts. Contribution: The contribution of this paper is twofold. First, we propose a joint registration algorithm for reconstruction of 3D imaging volumes from stacks of dissection photographs, requiring only a reference mask. This mask can be "hard", e.g., measured directly with a 3D surface scanner [10], or "soft", e.g., a probabilistic atlas of the whole brain; neither of these require MRI acquisition. Second, we present a Bayesian algorithm to segment 3D reconstructed stacks of photographs into 36 brain structures. We have designed an unsupervised likelihood term that models photography-specific artifacts, readily adapting to photographic hardware and brain fixation differences, thus making our publicly available code immediately usable with data from any institution.

General workflow
The workflow of our framework is outlined in Fig. 1. The inputs of the algorithm are dissection photographs of brain slices, and a reference volume R describing the exterior shape of the brain. In addition, the user provides two inputs. First, three landmarks on the photographs (e.g., on rulers which are commonly used in dissection photographs), which we use for coarse pixel size calibration and perspective correction. Second, segmentations for the different slices, which isolate tissue from background and encode the order of the slices. These can often be obtained with minimal interaction -we use a simple GMM [34] requiring one click per slice. The resulting slices and masks are ordered manually and arranged into sets of stacked slices S = {S n }, and corresponding masks M = {M n }.

3D reconstruction from dissection photographs
Let {Φ n } be a set of 2D affine geometric transforms for the brain slices, which brings them into alignment. These transforms correct for slice displacement and rotation, as well as perspective distortion. Then, S n [x; Φ n ] and M n [x; Φ n ] denote a resampling of slice n and its corresponding mask, to a discretised grid given by x, where the in-plane coordinates have been transformed according to the parameters in Φ n . Similarly, we define a 3D rigid transform for the reference volume R that brings it into alignment with the stack of slices, such that R[x; Ψ ] is a resampling of R parameterised by the 3D transform in Ψ . If the reference R is hard and directly represents the target shape, we also include an additional scaling in the direction of slicing (typically the anterior-posterior axis, for coronal slices) to account for deviations from the nominal slice thickness; with a soft atlas as reference, this is not possible. We jointly register the reference volume and the slices by maximising the following objective function F: where N s is the number of slices, D is the Dice score, C is the normalised cross correlation, f is a regulariser, and {α, β, γ, ν} are relative weights for each term. Equation 1 corrects for overall shape by encouraging a high Dice similarity coefficient between M [x; {Φ n }] and R, i.e., the 3D reconstructed mask and the reference volume. Smoothness within the reconstructed photography volume S[x; {Φ n }] is encouraged with two terms: the normalised cross correlation between successive slices in S, and the Dice coefficient between the corresponding masks. The final term, f (Φ n ), is a regulariser used to constrain the 2D spatial transforms not to be excessively scaled or sheared, which is particularly useful for the first and last slices in the stack, as they often contain little tissue.
The registration is solved in a hierarchical fashion using two levels with increasing complexity, combined with a multi-scale approach, to help avoid local maxima and increase convergence speed. At the first level, we limit the registration to correcting for slice displacement and rotation only, which is achieved by constraining each Φ n to be rigid rather than affine. At this level, no regularisation is needed. At the second level, we use the full model (i.e., affine {Φ n }). Undue scaling or shearing is avoided at this level by penalising transformations that excessively modify the area of a pixel, with f (Φ n ) = log |Φ n | in Equation 1.
Optimisation is performed at three levels of resolution (1/4, 1/2 and 1) with the L-BFGS algorithm [5,18]. The 2D transforms are initialised by aligning the centres of gravity (COGs) of the masks, and the 3D transform by matching the COG of the 3D mask with the COG of the initialised stack. Model parameters were set via visual inspection on a separate dataset. If the reference R is soft, we set α = 10, β = 1, γ = 2, ν = 0.1. If R is hard (i.e., measured directly), we give it a higher weight in the reconstruction (α = 50) and regularise less (ν = 0.05).

Segmentation
The ultimate purpose of the photographic volumes is for morphometric analyses, most of which require image segmentation. Since our goal is to make our code available to other researchers, supervised CNNs or MAS may not be appropriate as they may not generalise well to photographs of brains that have been fixed with potentially very different protocols. Instead, we propose a Bayesian algorithm with an unsupervised likelihood that includes a model of artefacts specific to photography, and thus adapts to cases fixed and imaged with any protocol.
Specifically, we maximise the probability of a 3D label map L given the image data D using Bayes' rule p(L|D) ∝ p(D|L)p(L). Both the prior and the likelihood have an associated set of parameters, θ L and θ D , respectively, with prior distributions p(θ L ) and p(θ D ). The prior p(L|θ L ) is a publicly available probabilistic atlas of anatomy [21] with K = 36 neuroanatomical classes, encoded as a tetrahedral mesh endowed with a deformation model [28]. Each voxel of the segmentation is assumed to be an independent sample of the discrete distribution defined by the deformed atlas at the location of the voxel.
The likelihood p(D|L, θ D ) combines a GMM with a model for brightness variations. Specifically, each of the K classes has an associated set of GMM parameters (weights, means, covariances), such that the intensity of a voxel is assumed to be a sample of the GMM associated with its label. These intensities are further corrupted by a slice-specific, smooth, multiplicative field (henceforth "brightness field"), which we assume to be a linear combination of smooth basis functions allowing bilinear variation in plane, independently for each slice.
It is typical in Bayesian segmentation to first compute point estimates of the model parameters (θ L ,θ D ) by maximising p(θ L , θ D |D), and then to estimate the segmentation as the maximum of p(D|L,θ D )p(L|θ L ). Let Γ = {Γ k } K k=1 be the GMM parameters of the different classes, C the matrix of brightness field coefficients (with 3 rows and as many columns as basis functions), and x the atlas mesh position, such that θ D = (Γ , C) and θ L = x. Then, taking the logarithm of p(θ L , θ D |D) yields the following objective function for the model parameters: x,Γ ,Ĉ = arg max (2) where d i is the vector with the log-transformed RGB intensities of voxel i and N is the number of voxels; a flat prior is assumed for Γ . The likelihood term is: where G k is the number of components of the GMM of class k; w k,g is the weight of component g of class k; N is the Gaussian distribution; φ i is a vector with the values of the brightness field basis functions at voxel i; and µ k,g , Σ k,g are the mean vector and covariance matrix associated with component g of class k.
Segmentation is achieved by maximising Equation 2 with coordinate ascent. We numerically optimise θ L with L-BFGS, initialised with an affine transform computed by registering the atlas to R with a robust approach [22]. We optimise θ D with the Generalised EM algorithm [8]. GEM involves iteratively: (i) constructing a lower bound of the objective function that touches it at the current estimate of the parameters, which amounts to a soft classification of each voxel (E step); and (ii) improving this bound to update θ D (generalised M step). Upon convergence, the probabilistic segmentation p(D|L,θ D )p(L|θ L ) is given by the soft classification of the final E step. We implement this optimisations by adapting routines from the public SAMSEG repository [21]. All parameters (number of Gaussians G k , mesh stiffness, etc.) are set to default SAMSEG values.

Datasets
We used a dataset consisting of dissection photography and matched ex vivo MRI for 24 cases, including only the cerebrum (i.e., no cerebellum or brainstem). Photographs were acquired of slices cut in the coronal plane with 4 mm thickness, using a 35 Megapixel camera. After initial perspective and pixel size adjustments, the calibrated photos were resampled to a pixel size of 0.1 mm. The ex vivo MRI were acquired using a FLAIR sequence with 0.8 mm isotropic resolution.
To evaluate the algorithms, we used two sets of reference segmentations. The first set consists of sparse manual delineations made on one slice photograph per volume. The slices were chosen to be close to the mid-coronal plane, while maximising visibility of seven representative subcortical structures: lateral ventricle, thalamus, caudate, putamen, pallidum, hippocampus and amygdala. The second set consists of dense segmentations of 36 brain structures, estimated from the FLAIR scans using SAMSEG [21]. Leaving aside two cases in which SAMSEG failed, we used these segmentations as a silver standard to evaluate the methods using every available voxel. Moreover, the cerebral tissue labels from these segmentations were also used to simulate a 3D surface scan (which in the future we plan to achieve with an inexpensive device) for hard reconstruction. We also tested a version with a soft reference, using the LPBA40 atlas [25].

Experimental Setup
Two volumes were reconstructed, using hard and soft references respectively, for each case in the dataset. Reconstructions were computed at 0.5 mm in-plane resolution, and segmented using the Bayesian method in Section 2.3. Since direct evaluation of registration error is very difficult, we used two measures of segmentation quality as surrogates. First, Dice scores were calculated against manual delineation of cerebral cortex, white matter, and the seven subcortical structures listed in Section 3.1. The evaluation was then extended to the whole brain by computing the correlation (and associated p value) between the volumes of the nine structures, derived from the ex vivo MRI and from the 3D reconstructed photographic volumes with our method. Figure 2 shows representative images from each stage in the registration and segmentation process, using a hard reference. The proposed procedure successfully aligns the photographs to the reference surface and estimates a brightness field that clearly increases the homogeneity of the image intensities, enabling accurate segmentation. Further qualitative results are shown in Figures 3(a-c), which compare the manual delineations with our automated segmentations, using both the hard and soft references. The corresponding quantitative results (Dice scores) are shown in Figures 3(d-e). The segmentations are quite accurate for most structures, except for the hippocampus and amygdala, whose interface is difficult to separate in this particular coronal plane. The method also commits minor mistakes that are common in Bayesian segmentation, e.g., including the  claustrum in the putamen. But overall, the Dice scores are competitive (above 0.8 for many structures), which is very encouraging given that they are computed from photographs. Particularly high scores are achieved for the cerebral cortex, since there is no extracerebral tissue in the images.

Results
While Dice scores on a single slice are informative, in order to show that measured trends from photo volumes are transferable to clinical imaging modalities it is crucial to test whether the volumes computed with our method on thick slices correlate well with the volumes derived from the isotropic MRI. For this reason we compare the silver standard volumes derived from the ex vivo MRI, with the volumes given by our proposed method. Table 1 shows the correlation coefficients and associated p values for the nine representative structures of interest. The results are consistent with the Dice scores on the sparse slices, showing strong correlations and significance for all structures except for the pallidum and the hippocampus. The pallidum is notoriously difficult to segment, even in MRI, due its low contrast. The hippocampus seems to be particularly affected by the large slice thickness in our reconstructed photography volumes.

Discussion and conclusion
We have presented the first algorithm for the construction of registered dissection photography volumes using only an external boundary shape as reference. For this proof of concept study we tested our methods on 24 cases, and assessed accuracy and sensitivity using Dice scores and correlations to silver standard volumes, respectively. The results are promising, both with the hard and the soft reference, and pave the way for inexpensive, large-scale NTNC studieseven retrospectively, using the soft version.
Future work will follow several directions. First, we will more thoroughly validate the methods, using additional cases, metrics, and manually traced images. We also plan to extend our method to slices with uneven thickness, and explore imputation algorithms (e.g., [6]) to increase the resolution of the 3D reconstructed scans. Access to super-resolved isotropic volumes is expected to enhance the quality of the segmentations (e.g., for the hippocampus), and also has the potential to enable other volumetric analyses that underperform with insufficient resolution (e.g., registration, cortical thickness). These additional analyses will likely benefit from a hard external reference: while drifts in the 3D reconstruction towards an average shape due to the probabilistic atlas do not seem to penalise segmentation, we hypothesise that using a 3D surface scan (an increasingly inexpensive technology) like the one in Fig. 1 will increase the precision of cortical measurements and registration, and enable the discovery of new imaging markers to study neurodegenerative diseases.