Histopathology in 3D: From three-dimensional reconstruction to multi-stain and multi-modal analysis.

Light microscopy applied to the domain of histopathology has traditionally been a two-dimensional imaging modality. Several authors, including the authors of this work, have extended the use of digital microscopy to three dimensions by stacking digital images of serial sections using image-based registration. In this paper, we give an overview of our approach, and of extensions to the approach to register multi-modal data sets such as sets of interleaved histopathology sections with different stains, and sets of histopathology images to radiology volumes with very different appearance. Our approach involves transforming dissimilar images into a multi-channel representation derived from co-occurrence statistics between roughly aligned images.

INTRODUCTION

A number of authors have addressed the problem of reconstruction of volumetric data from serial histopathology sections. These approaches may be divided into those approaches that rely on the tissue only,[1] and those that use a form of three-dimensional (3D) imaging (radiology, or block face imaging) to aid the reconstruction.[23] Within both groups, image-based registration is usually based on either an iterative optimization of a similarity metric,[4] or feature detection and matching.[5] Either approach has drawbacks. Optimization-based approaches can find local optima of the similarity function. With volumetric reconstruction, it is necessary to either perform a large number of registrations (one per section),[5] or optimization in a very high dimensional space,[4] either increasing the likelihood of failure. With feature detection-based methods, the features must be appropriate to the data, and thus a truly generic method is not possible. We have previously presented an alternative approach based on combining multiple local rigid registrations into a single non-rigid transform using a robust statistical estimator.[6] We use a closed form method of rigid registration based on phase correlation[78] to perform local rigid registration. This was selected as it is computationally efficient (four Fourier transforms, and a fixed set of multiplies and adds), and is guaranteed to find the maxima of similarity as it is equivalent to an exhaustive search. The drawback of this method when applied to multi-modal data is that it is based on greyscale similarity. The conventional approach to multi-modal registration is to use mutual information (MI) as a similarity metric.[9] However, this implies an iterative method with associated computational complexity and risk of local optima. Our alternative approach is to transform pairs of multi-modal images into pairs of multi-channel “tissue class probability” images based on co-occurrence statistics of roughly aligned images. MI is used within the process of forming these emergent tissue classes from the image pair (a rather different use of MI than the conventional similarity metric in registration). Once the multi-channel tissue class probability images have been formed, registration proceeds as with the single stain registration, excepting for the fact that there are N sets of local registrations (one per emergent tissue class). These are combined as before within our multi-level robust statistics framework to form a single B-spline based registration. The remainder of this paper is as follows: Section 2 details the basic robust statistical framework used in both the single stain and multi-modal registrations; Section 3 details the formation of the multi-channel “tissue class probability” images from roughly aligned images; Section 4 presents some case studies of applications of the technology; and Section 5 presents a discussion and conclusion.

VOLUMETRIC RECONSTRUCTION FROM SERIAL SECTIONS USING ROBUST STATISTICS

The main idea behind our method is that a single non-rigid registration for a pair of large images may be performed as a set of rigid registrations on sub-images, which are subsequently combined. This has the dual advantage of computational efficiency (memory and processor usage) and robustness (a single registration failure is not catastrophic as there is redundancy). In order to implement this idea, images are padded to the same size and pre-aligned rigidly, in order to that (roughly) corresponding regions may be extracted by dividing the images into regular grids. The basic workflow of our robust statistical framework is:

For images n = 1:N-1:

Pad images to the same size.

Rigidly align images n and n + 1 using greyscale phase correlation (image n = static image, image n + 1=moving image).

Divide each image pair into equal sized 50% overlapping patches, and rigidly align corresponding patches using a form of phase correlation that recovers rotation and translation.[6]

For each local registration construct five transform vectors (one at each corner, and one in the middle) from each registration.

Approximate the set of vectors by a rigid transform using a least squares minimizing method, and subtract this transform from each vector.

Approximate the “residual transform vector” set using a B-spline using a robust least squares minimizing method.[10]

Use transformed image n + 1 as static image for image n + 2.

In practice, steps 3–6 are repeated at multiple scales (from coarse to fine) and increasing degrees of freedom of the B-spline. The first reference image is selected by hand as an image with minimal distortion (to avoid propagating distortions to subsequent images). This approach has been applied successfully to reconstruct several hundredvolumes of different tissue types and chemical stains. A number of examples are visualized in Figures 1 and 2.

Figure 1

Single stain reconstruction results (stack views, 1 line per image) (a) bowel cancer in human liver (50 μm slice spacing, paraffin-embedded, H & E stained), (b) rat glomerulus (0.5 μm slice spacing, plastic embedded, H&E stained)

Figure 2

Liver tissue quantification. (a) Left: Original data, right: “Stacks view” of reconstructed data (one row from each image). (b) Volume rendering of reconstructed liver tissue

MULTI-STAIN AND MULTI-MODAL REGISTRATION

In this section, the generation of mapping functions to map images of different appearance to multi-channel images of more similar appearance is described. The outline of the method is as follows:

Represent each pixel of each image by a feature vector derived from local intensity, color, and texture. These features in include the output of Gaussian filters on color and greyscale channels, and a novel derivative based texture feature.[11]

Quantize the set of features separately for each image such that each pixel is represented by a prototype label (L1x,y, L2x,y). Clustering is performed using a binary PCA-tree method.[11] This method was selected for its computational efficiency, and ability to work with variation of different scales.

Consider a pair of hypothesized mapping functions that map the prototypes to a finite set of common tissue classes: C1x,y = M1 (L1x,y) C2x,y = M2 (L2x,y).

For a given pair of mapping functions, it is possible to generate a tissue class co-occurrence matrix. The MI calculated from this matrix is a measure of the similarity between the two images under that pair of mapping functions (and that feature set). A greedy search of potential mapping functions is performed in order to maximize MI, and select the best mapping functions.

The method is repeated with different feature sub-sets in order to perform feature selection. The feature sub-set with highest MI is selected.

Once the mapping functions for each image have been determined, the construction of probability images for each tissue class for each image is simply a matter of considering the co-occurrence of prototype labels in one image with tissue classes in the other. Counting these co-occurrences and normalizing gives P(Tissue Class|Prototype), which is mapped to a pixel value by multiplying by 255. Figure 3 illustrates results of applying this process for both multi-stain histopathology pairs and histopathology: Magnetic resonance imaging (MRI) pairs. Once the images have been constructed registration is applied as described in Section 2, with 5 x Nc vectors per block (where Nc is the number of tissue classes). Initial rigid alignment is using the same greyscale phase correlation method as described previously, which works on such multi-modal data (at low resolution) because of the clear distinction between foreground and background at low resolution in histopathology images. Full details may be found in Song et al.[11]

Figure 3

Tissue class images: (a) Two histopathology images with different stains (left: Original images and sub image, right: 3 “tissue class probability images” corresponding to each image/tissue class, (b) histopathology image and magnetic resonance imaging image)

APPLICATION TO VOLUMETRIC RADIOLOGY DATA

For multi-stain histopathology data sets, the 3D correspondence (slice to slice) is explicit in the data set. For histopathology to radiology registration, a two-dimensional (2D) oblique slice needs to be determined in order to compute tissue class probability images and subsequently perform 2D:2D non-rigid registration. Initially, this is performed manually for a single histopathology image using an interactive tool (MIM Medical Image Manager, HeteroGenius Ltd, Leeds, UK http://www.heterogenius.co.uk). Once one section is aligned, its 3D location can be optimized locally by maximizing MI between prototype labels (MI [L1,L2]) over a 3D rigid transform using Levenberg Marquart (LM) optimization.[12] LM is an iterative gradient-based optimization method. Subsequent sections can be placed in 3D space with reference to their theoretical geometric relation to the initial slice(s) (i.e., parallel with known normal offset based on section thickness/separation). Again, their location can be optimized by local optimization. Once placed in 3D space registrations to slices above and below (as Section 3) can be performed, in addition to registration to the volumetric radiology data. Accuracy of the method is demonstrated in Figure 4a. Typically, registration is accurate to within 200 μ (evaluated by measuring the distance between corresponding landmarks, such as blood vessels, in 2D), although larger deformations (such as tissue folds, and severe deformation) cannot be corrected for.

Figure 4

Rat heart collagen quantification (a) histology to magnetic resonance imaging (MRI) registration (b) three-dimensional segmentation of MRI based on the AHA heart model

CASE STUDIES

Liver disease quantification

We used the original volumetric reconstruction algorithm (Section 2) to generate volumes from liver tissue with five different types of liver disease (alcoholic liver disease; hepatitis C virus; primary biliary cirrhosis; primary sclerosing cholangitis; and polycystic liver disease) plus a healthy control [Figure 2]. Two 1 cm3 tissue samples were taken for each disease and sectioned with a microtome to give approximately 100 sections per tissue sample (separation 100 μm). These were stained with picrosirius red and scanned using an aperio T2 or T3 scanner (Aperio Inc., San Diego) at ×20 objective. Figure 2 shows liver nodules (stained brown/yellow) are surrounded by patterns of fibrotic tissue (stained red). The size, shape, and connectivity of nodules were quantified by (i) interactively segmenting nodules from other pixels using our in house Volume Viewer software (Volume Viewer, University of Leeds, Leeds, UK) with inplane resolution 1/64 of native resolution, (ii) separating nearby nodules using a 3D sub-voxel anisotropic morphological opening procedure, and (iii) assigning statistics to connected 3D components (size, elongation, etc.) using C++ code based on the Insight Toolkit (Kitware Inc., NY). The number of connected components and size of connected components showed a statistically significant variation between diseases, which is an indication of the (loss of) liver function in different diseases. Full results will be presented elsewhere.

CARDIAC COLLAGEN QUANTIFICATION

The purpose of this study was to quantify the effect of sub-sampling sections on collagen quantification in rat hearts. Previous works[13] had used small number of sections 1–3 to quantify collagen density in different parts of the heart. The problem with using very sparse sections is twofold; (i) accurate identification of the cardiac regions and (ii) the collagen density is heterogeneous and as such, sampling by taking a single section could introduce an undersampling error. To quantify the degree of the undersampling error, we took 1000 5 μm serial sections from each of two rat hearts (a male normal Wistar rat and a male Wistar rat in right heart failure) and aligned them to high resolution MRI volumes of the same hearts pre-sectioning scanned using a FLASH (fast low angle shot) MRI sequence in a Bruker (Ettlingen, Germany) 9.4T spectroscope with spatial resolution of 50 μm × 50 μm × 50 μm[18]. Each MRI volume was manually segmented into regions as defined by a modified American Heart Association model [Figure 4b],[14] which enabled labeling of each co-registered histopathology image. Collagen quantification was carried out in 2D using a standard method,[15] using all 1000 sections, and also subsets of 100, 50, 10 sections. Results showed acceptable quantification down to 100 sections (100 μm spaced sections), but thereafter sub-sampling resulted in increased variance over different data sub-sets. The lesson to be learned from this is that quantification using a single (or small numbers of) 2D section is potentially subject to sampling noise, even if the whole 2D sample is analyzed. A more robust way of performing the quantification is to take a larger number of samples of the tissue and aggregate the results. Full results are presented by Gilbert et al.[16]

COMPUTATIONAL MODELING OF SPINAL DISCS FROM MULTIPLE DIFFERENT STAINS

Data driven computational models are a tool that can help understand the disease and the implication of clinical actions (e.g., surgery). In the domain of musculoskeletal medicine, physics-based models (e.g., Finite Element models) may be constructed based on MRI or MicroCT data.[17] Such data only provide one value per voxel (Density in the case of MicroCT). Chemical stains used in histopathology can provide a wealth of other functional information like collagen density (relating to elasticity), cell density, etc. However, such data are 2D only and as such, not suitable for use in building 3D models. We have run experiments to reconstruct 3D data sets from multiple interleaved sections stained with different chemical stains using the methods described in Section 4. Data used was from an ovine intervertebral disc. In all five stains were used (Alcian Blue, ERSR, FAST, Elastic Pico Sirius Red and Sirius Red) to build a multi-parametric 3D representation of the data. This volumetric data was also aligned with high resolution MRI (scanner as in the previous section) to provide further anatomical information. The aim is to build multi-scale physics-based models based on the anatomical and functional data provided by this rich data set. The modeling work is ongoing.

DISCUSSION

Reconstructing microscopic functional and anatomical datasets in 3D using multiple 2D digital images is a powerful tool in a number of research areas including disease quantification and computational modeling. In this paper, we have described how multiple sources of information (multiple chemical and immunohistochemical stains, radiology) may be combined in a similar manner to stacking single stain 2D datasets using an information theory based image pre-processing method. In order to facilitate such research, the process of volumetric reconstruction must be as robust as possible. We have tackled this challenge using a combination of fast local analysis, robust statistics, a multi-scale approach, and a minimal (but important) amount of manual intervention. The combined method has been demonstrated to outperform iterative optimization-based techniques both in terms of accuracy and run-time.[11] Our techniques have been applied in a number of different areas, and we continue to explore applications and collaborations in surgical planning, radiology sequence development validation, disease quantification, and a number of other areas.