Towards arbitrary-scale histopathology image super-resolution: An efficient dual-branch framework via implicit self-texture enhancement

Single-image super-resolution (SISR) refers to recovering an HR image from an LR image or an LR image sequence, which is a classical low-level computer vision task with a wide range of applications. Deep neural networks can achieve accurate mapping from LR images to HR images due to their powerful fitting ability. Thus, they have become the mainstream approach in current SR studies. Numerous methods based on convolutional neural networks (CNNs) have been proposed for natural image SR, including SRCNN [16], EDSR [9], RDN [17], and RCAN [18]. To further improve the performance of SR, some methods utilized residual modules [19, 20], densely connected modules [21, 22], and other blocks [23, 24] for the design of the CNNs. Subsequently, a series of attention-based SR methods have emerged, such as channel attention [25, 18], self-attention (IPT [26], SwinIR [27], HAT [28]), and non-local attention [29, 30]. However, these methods can only be trained and tested at a fixed integer scale, and the networks need to be retrained for new scaling factors.

In recent years, implicit neural representation (INR) has been proposed as a continuous data representation for various tasks in computer vision. INR uses a neural network (usually a coordinate-based MLP) to establish a mapping between coordinates and their signal values, which allows continuous and efficient modeling of 2D image signals. For example, Chen et al. [15] first used INR in the SR algorithm and proposed the local implicit image function (LIIF) for arbitrary-scale SR. Lee et al. [31] proposed the local texture estimator (LTE), which transforms coordinates into Fourier domain information to enhance the representation of the local implicit function. Although these methods can be directly applied to pathological images for continuous magnification, they fail to recover the special textures in pathological images effectively.

In recent years, deep learning-based SR algorithms have been widely used in pathological images to improve imaging resolution [32, 10, 33, 11, 34, 8, 35, 36]. Upadhyay et al. [32] developed a generative adversarial network that considered pathological image SR and surgical smoke removal tasks at the same time. Mukherjee et al. [10] implemented SR image generation using a CNN and up-sampling layer and augmented the outputs using the K-nearest neighbor algorithm. Chen et al. [11] accomplished the SR task through a spatial wavelet dual-stream network incorporating a refine context fusion module. Li et al. [8] utilized a generative adversarial network based on a multi-scale CNN for SR image generation and introduced a curriculum learning training strategy. Wu et al. [35] added a branch for magnification classification to the SR network and improved SR performance through multi-task learning. These studies demonstrate the promise of using SR to improve pathological image resolution in low-resource settings. However, they still have some limitations. For instance, they restrict training and testing to specific scaling factors, and the resultant SR outputs still exhibit scope for refinement. We attribute this primarily to a lack of adequate consideration for the unique textural characteristics of pathological images. In this paper, we introduce ISTE as a solution to overcome these challenges, aiming to achieve arbitrary-scale SR of pathological images with high quality.

Given a set of $N$ pairs of corresponding LR images and HR images $\left\{X_{LR}^{i},Y_{HR}^{i}\right\}_{i=1}^{N}$, the objective is to find the optimal parameters $\hat{\theta}$ of the SR model $F_{\theta}$:

where $X_{LR}^{i}$ is a LR image and $Y_{HR}^{i}$ is its corresponding ground truth (GT), and L is the L1 loss function to measure the difference between the ground-truth and the generated HR images. Fig.2 shows the overall framework of our proposed ISTE. We first utilize SwinIR [27] to perform feature pre-extraction on the input LR image $X_{LR}$ and then input the pre-extracted feature $F_{LR}$ into the upper feature aggregation branch and lower texture learning branch of ISTE, respectively. In the feature aggregation branch, we input the feature $F_{LR}$ into the local feature interactor (LFI) to enhance the interaction of features in the local region and obtain feature $F_{LFI}$. In the texture learning branch, we input the image feature $F_{LR}$ into the texture learner (TL) to enhance the learning of high-frequency information and extract the feature $F_{TL}$. Then we design a two-stage texture enhancement strategy for these two branches, where the first stage is feature-based texture enhancement, and the second stage is spatial domain-based texture enhancement. In the first stage, we designed the self-texture fusion (STF) module to leverage the interaction of similar regions of the pathological images in the feature space, thereby accomplishing feature-based texture enhancement to assist in reconstruction. In the second stage, we decode the $F_{STF}$ from the STF module to obtain the image $I_{LPD}$ through the local pixel decoder (LPD). Simultaneously, we decode the $F_{TL}$ from the TL module to obtain the image $I_{LTD}$ through the local texture decoder (LTD). Subsequently, we perform spatial summation of $I_{LTD}$ and $I_{LPD}$, obtaining the final reconstructed HR image $I_{Pred}$. The primary purpose of the second stage is to fully utilize the features $F_{TL}$ learned by the texture learner and decode them into the spatial domain for texture enhancement.

We propose the LFI to enhance the interaction of features within local regions, thereby capturing the correlation of features within local regions. As shown in Fig.3, the size of the feature map $F_{LR}$ is $h\times w\times 64$, and we denote each vector of $F_{LR}$ as $F_{LR}^{j}(j=1,2,\ldots,h\times w)$. The LFI first assigns a window of size $3\times 3$ to each vector of $F_{LR}$, and the eight neighboring vectors in the window around $F_{LR}^{j}$ form a set $F_{N}^{j}=\left\{F_{N_{i}}^{j}\mid i=3,4,\ldots,10\right\}$. The average pooling result of the vectors within a window is denoted as $F_{P}^{j}$. The feature map $F_{LFI}$ output by the LFI is calculated through self-attention so that each point on the feature map incorporates local features while paying more attention to itself. We denote each vector of $F_{LFI}$ as $F_{LFI}^{j}(j=1,2,\ldots,h\times w)$, and it is calculated through Eq.(2).

where $Q_{LR}^{j}$ is the query mapped linearly from $F_{LR}^{j}$, $K_{1}^{j}$ is the key mapped linearly from $F_{LR}^{j}$, $V_{1}^{j}$ is the value mapped linearly from $F_{LR}^{j}$, $K_{2}^{j}$ is the key mapped linearly from $F_{P}^{j}$, $V_{2}^{j}$ is the value mapped linearly from $F_{P}^{j}$, $\left\{K_{i}^{j}\mid i=3,4,\ldots,10\right\}$ is the key mapped linearly from $F_{N}^{j}$, $\left\{V_{i}^{j}\mid i=3,4,\ldots,10\right\}$ is the value mapped linearly from $F_{N}^{j}$, and $d$ is the dimension of these vectors. The parameters used by each window are shared in the self-attention calculation.

Inspired by LTE [31], we propose the TL for learning high-frequency texture information in pathological images. We employ sine activation to effectively enhance implicit neural representations for learning high-frequency details in images, thereby mitigating spectral bias issues stemming from the ReLU activation functions[12]. Specifically, we normalize the value of 2D pixel coordinate $\left(X^{\prime},Y^{\prime}\right)=\left\{\left(\mathrm{x}_{i}^{\prime},\mathrm{y}_{j}^{\prime}\right)\mid\mathrm{i}=1,2,\ldots,mw,\mathrm{j}=1,2,\ldots,mh\right\}$ in the continuous HR image domain and the value of 2D pixel coordinate $(X,Y)=\left\{\left(\mathrm{x}_{i},\mathrm{y}_{j}\right)\mid\mathrm{i}=1,2,\ldots,mw,\mathrm{j}=1,2,\ldots,mh\right\}$ nearest to $\left(X^{\prime},Y^{\prime}\right)$ in the continuous LR image domain between -1 and 1, and the Local Grid is defined as $\left(X^{\prime}-X,Y^{\prime}-Y\right)$. Since each pixel coordinate of the HR image has a corresponding coordinate in the LR image grid that is closest to it, the number of both the HR and LR image coordinates is equal to $mh\times mw$, where $m$ represents the scale factor. As shown in Fig.4(a), the TL module firstly outputs three feature maps $F_{Amp}\in h\times w\times 256$, $F_{FreqX}\in h\times w\times 256$ and $F_{FreqY}\in h\times w\times 256$ through three $3\times 3$ convolutional kernels respectively, and predicts the feature maps $Amp\in mh\times mw\times 256$, $FreqX\in mh\times mw\times 256$ and $FreqY\in mh\times mw\times 256$ corresponding to each pixel coordinate of the HR image through nearest-neighbor interpolation. Then we use linear projection based on an MLP and Sigmoid activation function to map $(2/\mathrm{mw},2/\mathrm{mh})$ to a 256-dimensional feature vector $Phase$ to simulate the effect of texture fragment offset when the image scaling factor changes. The output of the TL module is calculated by Eq.(3):

where $\otimes$ represents element-wise multiplication and $\odot$ represents inner product operation.

Inspired by SRNTT [37] and T2Net [38], we propose the cross-attention-based STF module, whose main idea is to globally retrieve texture features most similar to $F_{LFIC}$ in $F_{TL}$ and fuse the retrieved features to $F_{LFIC}$, thus completing the feature-based texture enhancement. As shown in Fig.4(b), we use the features sampled from $F_{LFIC}$ by nearest-neighborhood interpolation as the query (Q) and use $F_{TL}$ as the key ($K$) and value ($V$) of the cross-attention module. To retrieve the texture features that are most relevant to the pixel feature $F_{LFIC}$, we first compute the similarity matrix $R$ of $Q$ and $K$, where each element $r_{i,j}$ of $R$ is computed according to Eq.(4), where $q_{i}$ represents an element of $Q$, and $k_{j}$ represents an element of $K$. Then we obtain the coordinate index matrix $T$ with the highest similarity to $q_{i}$ in $K$. An element in $T$ is $t_{i}=\arg\max_{j}\left(r_{i,j}\right)$, and $t_{i}$ represents the position coordinates of the texture feature $k_{j}$ with the highest similarity to $q_{i}$ in $F_{TL}$. We pick the feature vector $a_{i}$ with the highest similarity to each element in $Q$ from $V$ according to the coordinate index matrix $T$ to obtain the retrieved texture feature $A$, which can be represented by $a_{i}=v_{t_{i}}$ where $a_{i}$ is an element in $A$ and $v_{t_{i}}$ represents the element at the $t_{i}$-th position in $V$. To fuse the retrieved texture feature $A$ with the feature $F_{LFIC}$, we first concatenate $F_{LFIC}$ with $A$ and obtain the aggregated feature $Z$ through the output of an MLP, that is $Z=MLP(Concat(F_{LFIC},A))$. Finally, we calculate the soft attention map $S$, where an element $s_{i}$ in $S$ represents the confidence of each element $a_{i}$ in the retrieved texture feature $A$, and $s_{i}=\max_{j}\left(r_{i,j}\right)$. $F_{STF}$ is calculated as below:

where $\langle\cdot\rangle$ represents inner product operation, $\|\cdot\|$ represents the square root operation, and $\oplus$ represents element-wise summation.

In the spatial domain-based texture enhancement, we decode the texture feature $F_{TL}$ directly into the spatial domain $I_{LTD}$ and add it to $I_{LPD}$, which is reconstructed from $F_{FLIC}$ by the LPD, to obtain the final output $I_{Pred}$. Firstly, we utilize the LPD to decode the feature $F_{STF}$ into the RGB value $I_{LPD}$. We parameterize the LPD as an MLP. As shown in Fig.4(c), ${u_{t}}$ denotes the coordinates of the ${F_{LR}}$ and ${x_{q}}$ denotes the coordinates of the ${F_{STF}}$ and $F_{TL}$. We use $u_{t}(t\in 00,01,10,11)$ to denote the upper-left, upper-right, lower-left, and lower-right coordinates of an arbitrary point $x_{q}$, respectively. The RGB value at the coordinate ${x_{q}}$ in the HR image decoded by the LPD can be represented by Eq.(6), where $c$ contains two elements, $2/mh$ and $2/mw$, representing the sizes of each pixel in the $I_{LPD}$. Similarly, we calculate the RGB values of the texture information $I_{LTD}$ at coordinate $x_{q}$ via Eq.(7), where the LTD is parameterized as an MLP $g_{\varphi}$. We use the LTD to decode the texture features into the spatial domain texture information $I_{LTD}$ and add it to the $I_{LPD}$ via Eq.(8) for spatial domain texture enhancement to obtain the prediction result $I_{Pred}$, where ${\varphi}$ is the network parameter of the MLP $g_{\varphi}$. $S_{t}(t\in 00,01,10,11)$ is the area of the rectangular region between $x_{q}$ and $u_{t}$, and the weights are normalized by $S=\sum_{t\in\{00,01,10,11\}}S_{t}$.

Following previous SR methods based on implicit neural representation [15, 31], we used the patches with the size of ${48\times 48}$ as the input for training. We first randomly sampled the scaling factor ${m}$ in a uniform distribution U(1, 4) and cropped patches with the size of $48m\times 48m$ from the raw HR images in a batch, where ${m}$ represents the scaling factor. Following [8, 36], we resized the patches to ${48\times 48}$ via bicubic downsampling and did a Gaussian blur to simulate degradation since it is difficult to acquire authentically downsampled images at arbitrary scales through scanners. The size of the Gaussian kernel was set to 1/2 of the scaling factor ${m}$. We sampled $48^{2}$ pixels from the corresponding cropped patches to form RGB-Coordinate pairs. We utilized the deep learning toolbox Pytorch to implement ISTE and Adam as the optimizer, setting the initial learning rate to 0.0001 and epochs to 1000. We employed structure similarity index measure (SSIM) and peak signal-to-noise ratio (PSNR) to evaluate the quality of reconstructed HR images.

We compared the performance of ISTE with state-of-the-art SR methods in both the pathological image domain: SWD-Net [11] and Li et al. [8], and the natural image domain: Bicubic, EDSR [9], SwinIR [27], LIIF [15] and LTE [31], where the latter two are methods based on implicit neural representation. For a fair comparison, the backbone used for LIIF [15] and LTE [31] is also SwinIR [27] without upsampling layers.

To validate the effectiveness of each module in our proposed method, including the LFI, TL, STF, and LTD, we designed several variant networks for ablation experiments at scaling factors of ×2, ×3, and ×4 on the TCGA dataset.

In this section, we experimentally demonstrate that the proposed SR method effectively enhances the performance of two downstream tasks: gland segmentation and malignancy classification. First, for the gland segmentation task, we trained and tested the state-of-the-art segmentation model U-Net [42] on the Glas dataset from the MICCAI 2015 Gland Segmentation Challenge [43]. The Glas dataset consists of a training set and two test sets, Test A and Test B. The training set contains 85 images and the corresponding labels, Test A contains 60 images and the corresponding labels, and Test B contains 20 images and the corresponding labels. We performed ×4 downsampling on HR images to generate LR images using bicubic interpolation. We compared segmentation results under the following settings: (1) Original High-Resolution: Train U-Net on the original HR GlaS dataset for segmentation of original high-resolution images; (2) SISR: Directly employing U-Net trained on the original HR GlaS dataset for segmentation of the reconstructed images generated by the SISR model. (3) HR U-Net: Train U-Net on the reconstructed images generated by the SISR model for segmentation of original HR images; (4) Bicubic: Train U-Net on LR images obtained after bicubic interpolation for segmentation of original HR images. Table 5 shows the quantitative test results, where larger values indicate better performance for the F1 score and Object Dice score, while smaller values indicate better performance for object Hausdorff distance. It can be seen that the U-Net model trained on the reconstructed images of the SISR model performs better than the UNet model trained on the LR image dataset after bicubic interpolation, showing higher F1 scores and object Dice scores, as well as lower object Hausdorff distances. In particular, when tested on the Test B dataset, our results for segmentation of reconstructed images using U-Net trained on the original HR GlaS training set are close to those for segmentation of the original HR image, both with an F1 score of 0.93. Fig.10 shows representative results for different experimental setups, and it can be observed that U-Net trained on LR images produced the worst results; not only did it fail to detect small glands, but also the segmentation results of large glands appeared to be crippled. In contrast, the U-Net trained on the reconstructed image could outline the boundaries of the macro glands and detect the tiny glands. Compared to using LR images for training, using the generated SR images for training can improve the segmentation accuracy during testing.

To further evaluate the contribution of the SR method to the malignancy classification task, we conducted tumor recognition on the PCam dataset [44]. The PCam dataset comprises 262,144 color images for training and 32,768 images for testing, with each image annotated with a binary label indicating the presence of metastatic tissue. We performed ×2 downsampling on HR images of the test set to generate LR images using bicubic interpolation. The ResNet-50 [45] was chosen as the classifier and trained on the original PCam dataset. We compared classification results under the following settings: (1) Original: Directly employing trained ResNet-50 model to test on the original HR images in the test set; (2) Low Resolution: Directly employing trained ResNet-50 model to test on the LR images of the test set; (3) Bicubic: Directly employing trained ResNet-50 model to test on the bicubic images of the test set; (4) LIIF: Directly employing trained ResNet-50 model to test on the images generated by LIIF from the LR test set images; (5) ISTE: Directly employing trained ResNet-50 model to test on the images generated by our method ISTE from the LR test set images; Table 6 illustrates the enhancement in diagnostic performance by the SR method. By introducing additional prior knowledge, our ISTE leads to a performance improvement, with an accuracy increase of 4.06% compared to Bicubic. These results indicate that ISTE can improve classification performance by recovering more distinctive details.