1.1.

Application Domain

The type of imaging modality we are interested in is called whole slide imaging (WSI) or virtual microscopy. WSI is a recent innovation in the field of digital pathology in which digital slide scanners are used to create images of entire histologic sections. Traditionally, the use of the optical capabilities of a microscope to “focus” the lens on a small subsection of the slide (based on the field of view of the device) to review and evaluate the specimen is often carried out by a trained pathologist. This process may need to be repeated for other sections of the slide depending on the scientific/clinical question of interest, toward obtaining consistently well-focused digital slides. WSI essentially automates this procedure for the whole slide. The ability to do so, automatically for a large number of slides, ideally capturing as much information as the pathologist may have been manually able to glean from the histological specimen via a light microscope, offers an immensely powerful tool with many immediate applications in clinical practice, research, and education.^6–8 For instance, WSI makes it feasible to solicit subspecialty disease expertise regardless of the location of the pathologist, integration of patient medical records in their health portfolio, pooling data between research institutions, and reducing long-term storage costs of histological specimens. However, given the relatively recent advent of WSI technology, there are several barriers that still need to be overcome, before it is widely accepted in clinical practice. The chief among these are the fact that HR WSI scanners, which have been shown to match images obtained from light microscopy in terms of quality for diagnostic capability, are typically very expensive, even for LR usage. In addition, the size of the files produced also generates a bottleneck. Typically, a virtual slide acquired at HR is about 1600 to 2000 megapixels, which results in a file size of several gigabytes. Such files are typically much larger that image files used by other clinical imaging specialties such as radiology. If one has to transport, share, or upload multiple such files or 3D stack, it results in a consequential increase of storage capacity and network bandwidth. Notwithstanding these issues, WSI offers tremendous potential and numerous advantages for pathologists, which is why it is important to find a way to alleviate the existing issues with WSI, such that it can be widely applicable. One potential way to address these issues is to use images from low magnification slide scanners, which are widely available, easy to use, comparatively inexpensive, and can also quickly produce images with smaller storage requirements. However, such LR images can increase the chance of misdiagnosis and false treatment if used as the primary source by a pathologist. For example, cancer grading normally requires identifying tumor cells based on size and morphology assessments,⁹ which can be easily distorted in low magnification images. If such images were indeed to be used for research, clinical, or educational purposes, we need a way to convert such LR data and produce an output that closely resembles images acquired by a top-of-the-line WSI scanner, without substantial increase in storage and computational requirements.

1.2.

Solution Strategies and Associated Challenges

One way to address the above issues is to dynamically improve the quality and resolution of LR images to render them comparable in quality to those acquired from HR scanners, as and when needed by the end user. Such a proposed workflow would need a fraction of the time and can yield near real-time quantifiably accurate results at a fraction of the setup cost of a standard WSI system. The implementation of such a system would require an SR approach that works well for WSI images. But still, we find that there is no off-the-shelf deep network architecture that can be used for our application directly. The reasons for this are numerous. First, most existing methods have been trained on databases of natural images. However, the WSI images under consideration do not have the same characteristics as natural images, such as, textures, edges, and other high-level exemplars, which are often leveraged by the SR algorithms. Second, such deep learning models are often trained using large training datasets (usually consisting of synthetic/resized HR–LR pairs), in the order of millions. This does not directly translate to our application for two reasons: (a) large training datasets are typically harder to acquire since each image pertains to a unique acquisition that requires significant human involvement, and (b) in our case, the LR images are acquired from a different scanner and is not just a resized version of HR image. Third and perhaps the most important limitation is that existing deep SR methods typically only reconstruct successfully up to a resolution increase factor of 4, whereas in case of WSI, the resolution (from a low-cost scanner to an expensive HR scanner) can increase up to a factor of $10\times$, since there can be wide variance in resolution between an LR scanner ($4\times$) to HR scanners, which typically scan at $20\times$ or $40\times$. We discuss this issue in detail in the following paragraph.

1.3.

Our Contribution

Existing CNN-based methods have shown limited performance in scenarios when the resolution difference is high. The reason for this is that the complexity of the transformation that morphs the LR image to the HR one increases greatly in such situations, which in turn manifests in the CNN models taking longer to converge, or learning a function that generalizes poorly to unseen data. A typical way to address this issue is to make the CNN models deeper, by adding more layers ($>5$ layers) or increasing the number of examples required to learn a task to a given degree of accuracy while still keeping the network shallow. Both these directions pose challenges for our application: (a) a deeper network is associated with far more parameters, increasing the computational and memory footprint, to the extent that model may not be applicable in a real-time setup and (b) increasing the number of samples the extent needed would be impractical, due to associated time and monetary costs.

Our approach to solving this problem draws upon the inherent nature of the SR problem. While it is hard to acquire a large training dataset in this scenario, it is much more time and cost efficient to obtain the WSI images at different resolutions by varying the focus of the scanner. In this paper, we study how such multi-resolution data can be used effectively to make the transformation learning problem more tractable. Suppose $I_1$ and $I_h$ represent a particular LR and HR image pair. If $I_h$ is a significant resolution ratio higher than $I_1$, learning their direct transformation function $f(I_1)=I_h$ can be challenging leading to a overparameterized system. But if we had access to some intermediate resolutions say $I_2,\dots I_{h-1}$ (with a smaller resolution change between consecutive images), it makes intuitive sense that transformation that converts an image of a given resolution into the closest HR image would be roughly the same across all the resolutions considered, if we assume that resolution changes vary smoothly across the sequence. Having more image pairs $(I_{k-1},I_k)$ for $k=2\dots h$ to train, it may be computationally easier to learn a smooth function $\widehat{f}$, such that $\widehat{f}(I_{k-1})\approx I_k$ for all $k$. In this paper, we formalize this notion and develop a recurrent convolutional neural network (RCNN) [Note that the acronym RCNN is also used to refer to region-based CNNs,¹⁰ but in the context of this paper, we use it to refer to recurrent convolutional neural network.] to learn from multi-resolution slide scanner images. Our main contributions are as follows. (1) We propose a new version of SR problem motivated from this problem, multi-resolution SR (MSR), where the aim is to learn the transformation function, given a sequence of resolutions, rather than simply the LR and HR images. To the best of our knowledge, this is new problem scenario for SR that has not been studied before. (2) We propose an RCNN model to solve the MSR problem. (3) We demonstrate using experimental results on three WSI cell lines that the MSR idea can indeed reduce the need for large sample sizes and still learn the transformation that generalizes to unseen data.

2.1.

Deep Network Models for SR

Stacked collaborative local autoencoders are used¹¹ to construct the LR image layer by layer. Reference 12 suggested a method for SR based on an extension of the predictive convolutional sparse coding framework. A multiple layer CNN, similar to our model, inspired by sparse-coding methods, is proposed in Refs. 1, 2, and 13. Chen and Pock¹⁴ proposed to use multistage trainable nonlinear reaction diffusion as an alternative to CNN where the weights and the nonlinearity are trainable. Wang et al.¹⁵ trained a cascaded sparse coding network from end to end inspired by learning iterative shrinkage and thresholding algorithm¹⁶ to fully exploit the natural sparsity of images. Recently, Ref. 17 proposed a method for automated texture synthesis in reconstructed images by using a perceptual loss focusing on creating realistic textures. Several recent ideas have involved reducing the training complexity of the learning models using approaches, such as Laplacian pyramids,¹⁸ removing unnecessary components of CNN,¹⁹ and addressing the mutual dependencies of LR and HR images using deep back-projection networks.²⁰ In addition, generative adversarial networks (GAN) have also been used for the problem of single image SR, these include Refs. 2122.23.–24. Other deep network-based models for image SR problem include Refs. 2526.27.–28. We also briefly review SR approaches for sequence data such as videos. Most of the existing deep learning-based video SR methods using motion information inherent in video to generate a single HR output frame from multiple LR input frames. Kappeler et al.²⁹ warp video frames from the preceding and subsequent LR frames onto the current one using the optical flow method and pass them through a CNN that produces the output frame. Caballero et al.³⁰ followed the same approach but replaced the optical flow model with a trainable motion compensation network. Huang et al.³¹ used a bidirectional recurrent architecture for video SR with shallow networks but do not use any explicit motion compensation in their model. Other notable works include Refs. 32 and 33.

2.2.

Recurrent Neural Networks

A recurrent neural network (RNN) is a class of artificial neural network where connections between nodes form a directed graph along a sequence. This allows it to exhibit temporal dynamic behavior for a time sequence. Unlike feedforward neural networks such as CNNs, the input and outputs are not considered independent of each other, rather such models recompute the same/similar function for each element in the sequence, with the intermediate and final output of subsequent elements in the network being dependent on the previous computations on elements occurring earlier in the sequence. RNNs have most frequently been used in applications for language modeling,³⁴ speech recognition,³⁵ and machine translation.³⁶ But it can be applied to many learning tasks applied to sequence data, for more details, see survey paper by Ref. 37.

2.3.

CNN Architectures for Small Sample Size Training

In this regard, Erhan et al.³⁸ devised unsupervised pretraining of deep architecture and showed that such weights of the network generalize better than randomly initialized weights. Mishkin and Matas³⁹ have proposed layer-sequential unit-variance that utilizes the orthonormal matrices to initialize the weights of CNN layer. Andén and Mallat⁴⁰ proposed scattering transform network (ScatNet), which is a CNN-like architecture where predefined Morlet filter bank is used to extract features. Other notable architectures in this regard include PCANet,⁴¹ LDANet,⁴² kernel PCANet,⁴³ MLDANet,⁴⁴ DLANet,⁴⁵ to name a few.

2.4.

Deep Network Models in Microscopy

Since the application domain of this paper is in microscopy, we briefly review related papers that have used deep networks in this area. Most similar to our work is Rivenson et al.,⁴⁶ who showed how to enhance the spatial resolution of LR optical microscopy images over a large field of view and depth of field. But unlike our model, this framework is meant for single-image SR, where the model is trained on pairs of HR and LR images and provided a single LR image at test time. The design of this model includes a single CNN architecture, though they show that feeding the output of the CNN, back to the input, can improve results further. Wang et al.⁴⁷ proposed a GAN-based approach for super-resolving confocal microscopy images to match the resolution acquired with a stimulated emission depletion microscopes. Grant-Jacob et al.⁴⁸ designed a neural lens model based on a network of neural networks, which is used to transform optical microscope images into a resolution comparable to a scanning electron microscope image. Sinha et al.⁴⁹ used deep networks to recover phase objects given their propagated intensity diffraction patterns. Other related methods that have used deep learning-based reconstruction approaches for various applications in microscopy include Nehme et al.,⁵⁰ Wu et al.,⁵¹ and Nguyen et al.⁵² A more detailed survey of deep learning methods in microscopy can be found in Ref. 53.

3.1.

CNN Subnetwork

Here, we describe the basic structure of each CNN subnetwork and its constituent layers briefly. A more detailed description can be found in Sec. 7. Note that the components/layers of each CNN pipeline is kept the same, see Fig. 1. The first layer is a feature extraction-mapping layer that extracts features from LR images (denoted by $Y_1^j$ for the $j$’th pipeline), which are then served as input to the next layers. This is followed by three convolutional layers. We briefly elaborate on this, since they are useful to understand the RCNN terminology in the next section.

Fig. 1

Architecture of the proposed CNN for image super-resolution.

3.2.

Convolutional Layers

The feature extraction layer is followed by three additional convolutional layers. We also refer to these as hidden layers, denoted by $H_i^j$, which is the $i$’th hidden layer ($i\in \{\mathrm{2,3},4\}$) in the $j$’th CNN pipeline. The input to this layer is referred to as $Y_{i-1}^j$, and the output is denoted by $Y_i^j$. The filter functions in these intermediate layers can be written as

The last layer of the CNN architecture consists of a subpixel layer that upscales the LR feature maps to the size of the HR image.

3.3.

Recurrent Convolution Network

The recurrent convolution network is built by interconnecting the hidden units of the CNN subarchitectures, see Fig. 2. We index the CNN pipeline components with superscript $j\in \mathrm{0,1},2$ with the $j$’th pipeline being given image $I^j$ as input and reconstructing the image $I^{j+1}$. The basic premise of our RCNN model is that the hidden units processing the each of the images can be informed by the outputs of the hidden units in other CNN pipelines. We use one directional dependence (low to high) as it is more challenging to reconstruct HR images compared to lower resolution ones. We can introduce bidirectional dependencies as well, but in practice, this increases the number of parameters substantially and contributes to an increase in training time for the model.

Fig. 2

Connections between hidden units of the three CNN subarchitectures: feedforward connections (in black), recurrent connections (in red), and prelayer connections (in blue).

Besides the feedforward connections already discussed as a part of the CNN subarchitectures, we introduce two additional connections to encode the dependencies among the various hidden units, see Fig. 2. These are as follows.

3.4.

Training and Loss Function

The complete architecture of our network can be seen in Fig. 3. The output from Eq. (2) (for pipelines $j=1$ and $j=2$) is then passed on as an input to the subpixel layer [described in Eq. (4)], which outputs the desired prediction (let this be denoted by $R^j$). For the pipeline ($j=0$), the prediction is simply $R^0=Y_5^0$. This network is learned by minimizing a weighted function of mean square error (MSE) and structured similarity metric (SSIM) between the predicted HR and the ground truth at each pipeline

Fig. 3

Architecture of our proposed RCNN for image super-resolution.

The objective/loss function in Eq. (3) is optimized using stochastic gradient descent. During optimization, all the filter weights of recurrent and prelayer convolutions are initialized by randomly sampling from a Gaussian distribution with mean 0 and standard deviation 0.001, whereas the filter weights of feedforward convolution are set to 0. Note that one can also initialize these weights by pretraining the CNN pipeline on a small sized dataset. We experimentally find that using a smaller learning rate ($1e-5$) for the weights of the filters is crucial to obtain good performance.

4.1.

Datasets

We performed experiments to evaluate our RCNN approach on three previously published tissue microarray (TMA) cancer datasets, a breast TMA dataset consisting of 60 images,⁵⁴ and a kidney TMA dataset with 381 images,⁵⁵ and a pancreatic TMA dataset with 180 images.⁵⁶ The datasets are summarized in Table 1.

Table 1

Summary of datasets.

4.2.

Imaging Systems

In the datasets we analyze, highest resolution images were acquired and digitalized at $40\times$ using an Aperio CS2 digital pathology scanner (Leica Biosystems),⁵⁷ with $4\text{pixels}/\mu \mathrm{m}$, and lowest resolution images were acquired and digitized using PathScan Enabler IV⁵⁸ with $0.29\text{pixels}/\mu \mathrm{m}$ ($4\times$). Besides these, we have acquired images at two different resolutions ($10\times$ and $20\times$ from the Aperio CS2 digital pathology scanner, which served as our intermediate resolutions). These images are then resized (in the order of resolution) to $256\times 256$, $512\times 512$, $1024\times 1024$, and $2048\times 2048$, which provides us the dataset with four different resolutions.

4.3.

Evaluations

We evaluate various aspects of our RCNN model to determine the efficacy of our method. First, we evaluate the quality of reconstruction of our RCNN model compared to a single CNN pipeline (which is a baseline for our model) and other comparable methods for SR. We show both qualitative and quantitative results in this regard. Second, we study the advantage of having intermediate resolutions, by varying the number of resolutions available. We also analyze how useful our obtained reconstructions are toward end-user applications, such as segmentation and perform a user study, done by a pathologist to evaluate the utility of the reconstructed images for cancer grading. In addition, we study how the reconstruction accuracy is affected as a function of the training set size. Finally, we discuss the parameters used and the computational time requirements of our model. We discuss these issues next.

4.3.1.

Quality of reconstruction

Metrics

We evaluate the reconstruction quality of the obtained images by our approach by evaluating it relative to HR ground truth image and calculating seven different metrics: (1) root mean square error (RMSE), (2) signal-to-noise ratio (SNR), (3) SSIM, and (4) mutual information (MI), (5) multiscale structured similarity (MSSIM), (6) noise quality measure (NQM),⁵⁹ and (7) weighted peak signal-to-noise ratio (WSNR).⁶⁰ RMSE should be as low as possible, whereas SNR, SSIM (1 being the maximum), MSSIM (1 being the maximum), and the remaining metrics, should be high for good reconstruction.

Experimental setup

Note that our model was trained by providing three sets of images (of three different resolutions) as input. However, our testing experiments can be done in the following two settings: (a) in the first case, we provide the images of the same three resolutions $I^0$, $I^1$, and $I^2$ as input to the trained model, which then outputs the reconstructed highest resolution image $I^3$. We call this setup RCNN(full). (b) In the second case, we see how our model behaves given only the lowest resolution image $I^0$. In this case, we first generate the two intermediate resolutions as follows. First, we only activate pipeline $j=0$, which outputs $I^1$. Then, use this as input and activate both pipelines $j=0$ and $j=1$, which then reconstructs $I^2$ as well. Using $I^0$ and the reconstructed $I_1$ and $I_2$ as input, we activate all three pipelines and reconstruct $I_3$. We call this setup RCNN(1 input).

Comparable methods

To the best of our knowledge, we do not know of any other neural network framework to study the MSR problem presented in this paper. Therefore, to compare our method to other state-of-art methods, we choose other SR approaches that work in the two resolution setting (low and high). Our default baseline is the CNN architecture shown in Fig. 1. We refer to this method as CNN. In addition, we also compare with the following methods: (i) the CNN-based framework (FSRCNN) by Dong et al.,¹³ (ii) a CNN model that uses a subpixel layer (ESCNN), and (iii) a GAN-based approach for SR.²¹

Results

Results for the breast, kidney, and pancreatic datasets are shown in Tables 2Table 3–4, respectively. In each case, we see that the RCNN(full) setting outperforms all other methods, giving a significant improvement in all the metrics calculated. A qualitative analysis of the reconstructed images in Figs. 4Fig. 5–6 shows that the reconstructed images are indeed very similar to the HR images. The comparatively poorer performances of comparable methods, such as FSRCNN or SRGAN, are expected since these methods are not designed to learn a resolution ratio of 8 used in our experiments. Still we find that RCNN(1 input), which is trained on all three input resolutions but tested by providing the lowest resolution only, outperforms the other baselines, showing that the weights learned our model generalizes better than other neural network models for this difficult learning task. The qualitative results showing the performance of RCNN(1 input) is shown in Fig. 7.

Table 2

Quantitative results from reconstructed breast images.

Breast TMA
Metric	SRGAN	ESCNN	FSRCNN	CNN	RCNN(full)	RCNN(1 input)
RMSE	48.12	42.86	46.64	39.57	15.64	31.27
SNR	14.63	15.37	14.95	16.37	24.36	18.37
SSIM	0.40	0.35	0.42	0.34	0.98	0.51
MI	0.05	0.09	0.01	0.08	0.31	0.36
MSSIM	0.42	0.39	0.19	0.38	0.95	0.53
NQM	0.37	0.39	0.28	1.09	20.33	2.48
WSNR	13.83	14.61	14.78	15.77	26.59	18.04

Note: Best values are highlighted in bold.

Table 3

Quantitative results from reconstructed kidney images.

Kidney TMA
Metric	SRGAN	ESCNN	FSRCNN	CNN	RCNN(full)	RCNN(1 input)
RMSE	28.75	25.48	38.90	29.15	11.60	22.92
SNR	19.00	20.06	16.35	18.96	28.31	21.03
SSIM	0.82	0.72	0.39	0.76	0.98	0.85
MI	0.11	0.16	0.07	0.13	0.35	0.31
MSSIM	0.70	0.70	0.41	0.68	0.97	0.78
NQM	7.77	4.61	0.45	6.80	11.15	10.30
WSNR	20.55	20.34	15.71	19.58	19.75	21.89

Note: Best values are highlighted in bold.

Table 4

Quantitative results from reconstructed pancreatic images.

Pancreas TMA
Metric	SRGAN	ESCNN	FSRCNN	CNN	RCNN(full)	RCNN(1 input)
RMSE	84.59	37.30	39.56	35.39	20.32	33.50
SNR	10.0	16.78	16.26	17.26	22.07	17.78
SSIM	0.39	0.52	0.42	0.64	0.96	0.79
MI	0.07	0.16	0.16	0.16	0.29	0.33
MSSIM	0.39	0.50	0.42	0.56	0.93	0.69
NQM	0.14	7.10	5.95	10.24	16.94	10.97
WSNR	7.93	16.27	15.57	16.99	24.79	17.88

Note: Best values are highlighted in bold.

Fig. 4

Results of reconstruction of breast cancer TMA: columns 1 and 3 show HR and LR images and column 2 shows the reconstructed image. Rows 2 and 4 show a zoomed in region of interest from the corresponding images in row 1 and row 3, respectively.

Fig. 5

Results of reconstruction of kidney cancer TMA: columns 1 and 3 show HR and LR images and column 2 shows the reconstructed image. Rows 2 and 4 show a zoomed in region of interest from the corresponding images in row 1 and row 3, respectively.

Fig. 6

Results of reconstruction of pancreatic cancer TMA: columns 1 and 3 show HR and LR images and column 2 shows the reconstructed image. Rows 2 and 4 show a zoomed in region of interest from the corresponding images in row 1 and row 3, respectively.

Fig. 7

Results of reconstruction of all three cell types using RCNN(1 input): column 1 shows breast cells; column 2 shows kidney cells; and column 3 shows pancreatic cells. The HR images for the breast cells [images (a) and (d) in this figure] are shown in Figs. 4(a) and 4(d) and the corresponding LR images are in Figs. 4(c) and 4(f), respectively. Similarly for the kidney [images (b) and (c) of this figure], the HR images are in Figs. 5(a) and 5(d) and LR images are in 5(c) and 5(f), respectively. Finally, for the pancreatic cells [(c) and (f) in this figure], the HR images are in Figs. 6(a) and 6(d) and LRn images are in 6(c) and 6(f), respectively.

4.3.5.

Effect of training set size

The necessary size of the training set for a particular learning task is generally hard to estimate and depends mostly on the hardness of the learning problem as well as on the type of model being trained. Here, we provide empirical evidence of how our model behaves wrt to increasing the size of the training set. For this purpose, we use the Kidney dataset, since it is the largest dataset considered in this paper. We vary the training size from $\{\mathrm{50,100,150,200,300}\}$. The test set is set to 50 in each case. We analyze the quality of reconstruction by comparing the SNR in each case.

As seen in Fig. 8, the SNR improves only slightly when the training set is increased, indicating that a higher training size does not significantly improve the image reconstruction metrics.

Fig. 8

SNR as a function of training size.

7.1.

Feature Extraction-Mapping Layer

This layer is used as the first step to extract features from the LR input images ($I^j$) for the $j$’th CNN subarchitecture. The feature extraction process can be framed as a convolution operation and hence implemented as a single layer of the CNN. This can be expressed as

7.2.

Convolutional Layers

The feature extraction layer is followed by three additional convolutional layers. We also refer to these as hidden layers, denoted by $H_i^j$, which is the $i$’th hidden layer ($i\in \{\mathrm{2,3},4\}$) in the $j$’th CNN pipeline. The input to this layer is referred to as $Y_{i-1}^j$ and the output is denoted by $Y_i^j$. The filter functions in these intermediate layers can be written as

7.3.

Subpixel Layer

In our CNN architecture, we leave the final upscaling of the learned LR feature maps to match the size of the HR image, to be done at the last layer of the CNN. This is implemented as a subpixel layer similar to Ref. 74. The advantage of this is that is that layers prior to the last layer operate on the reduced LR image rather than HR size, which reduce the computational and memory complexity substantially.

The upscaling of the LR image to the size of the HR image is implemented as a convolution with a filter ${\theta }_{\text{sub}}^j$ whose stride is $\frac{1}{r}$ ($r$ is the resolution ratio between the HR and LR images). The size of the filter is denoted as $f_{\text{sub}}$. A convolution with stride of $\frac{1}{r}$ in the LR space with a filter ${\theta }_{\text{sub}}^j$ (weight spacing $\frac{1}{r}$) would activate different parts of ${\theta }_{\text{sub}}^j$ for the convolution. The patterns are activated at periodic intervals of $\mathrm{mod}(a,r)$ and $\mathrm{mod}(b,r)$ where $a$, $b$ are the pixel position in HR space. This can be implemented as a filter ${\theta }_5^j$, whose size is $n_4\times r^2\times f_5\times f_5$, given that $f_5=\frac{f_{\mathrm{sub}}}{r}$ and $\mathrm{mod}(f_{\mathrm{sub}},r)=0$. This can be written as