The IQT problem can be mathematically formulated as follows: Given an original vectorised high-quality volume ${\mathbf{X}}\in\mathbb{R}^{M}$, its corresponding low-quality version is denoted as ${\mathbf{Y}}\in\mathbb{R}^{P}$, where the relation between the two volumes can be modeled as

where ${\mathbf{H}}$ is the matrix representing a linear blurring operator, ${\mathbf{L}}$ is the downsampling operator, and ${\mathbf{W}}$ stands for additive noise, modelling observation noise and model mismatch and is assumed to be a white Gaussian noise sequence. This equation states that ${\mathbf{Y}}$ is a blurred and down-sampled version of the original high-quality volume ${\mathbf{X}}$.

In IQT, the goal is to recover a high-quality volume $\hat{{\mathbf{X}}}$ given its blurred and down-sampled version ${\mathbf{Y}}$, such that $\hat{{\mathbf{X}}}\approx{\mathbf{X}}$. The problem of estimating ${\mathbf{X}}$ from ${\mathbf{Y}}$ in Eq. (1) is an ill-posed linear inverse problem (LIP), i.e., the matrix ${\mathbf{L}}{\mathbf{H}}$ is singular and/or very ill-conditioned, since for a given low-quality input, infinitely many high-quality volumes satisfy the above equation. Consequently, this problem requires additional regularisation (or prior information from Bayesian perspective) in order to reduce uncertainties and improve estimation performance.

Figure 1 shows a schematic diagram to the IQT problem using a sparse representation model and dictionary learning. The proposed model consists of two separate stages. First, the coupled low-quality and high-quality dictionaries, ${\mathbf{D}}_{\ell}$ and ${\mathbf{D}}_{h}$ respectively, are constructed from training data set. Then, a reconstruction algorithm is applied to upscale a test low-quality volume to recover its high-quality version. This algorithm considers the patch-based sparse prior model to recover an estimate to the high-quality volume in a patch-by-patch basis. The following sections provide more details about the two stages mentioned above.

Constructing the high-quality and low-quality dictionaries requires a set of matched high- and low-quality volume patches. The training set is composed by a set of high-quality and the corresponding low-quality volumes. As proposed by Zeyde et al. (2010), the high-quality volumes are processed to obtain only the high-frequency information, whereas the intensity maps are used for the low-quality volumes. Each of the high- and low-quality volumes are then split into a set of 3D patches which are vectorised and training pairs are generated. Patches containing $>80\%$ background voxels are excluded from the patch library. The coupled-dictionary training algorithm proposed by Zeyde et al. (2010) is then used in order to obtain the low- and high-quality dictionaries ${\mathbf{D}}_{\ell}$ and ${\mathbf{D}}_{h}$ respectively. For this local model, the two dictionaries ${\mathbf{D}}_{h}$ and ${\mathbf{D}}_{\ell}$ are trained such that they share the same sparse representations for each high- and low-quality volume patch pair. Finally, the dimensionality of ${{\mathbf{D}}_{\ell}}$ may be reduced to speed up the subsequent computations, given the intrinsic redundancy of the multi-scale edge analysis. For doing so, a Principal Component Analysis (PCA) is applied to this matrix, searching for a set of projection coefficients that represents at least $90\%$ of the original variance. All patches are collected together to form the reduced low-quality dictionary ${\mathbf{D}}_{\ell}$, whereby the number of atoms in the dictionary has not changed.

The low-quality volume ${\mathbf{Y}}$ can be split into a set of overlapping 3D patches ${\mathbf{y}}$, each of size $\sqrt[3]{m}\times\sqrt[3]{m}\times\sqrt[3]{m}$. With the sparse generative model, each patch ${\mathbf{y}}$ can be represented by a linear combination of a few atoms drawn from a dictionary ${\mathbf{D}}_{\ell}$, which characterises the low-quality patches. This can be written as

where ${\boldsymbol{\alpha}}\in\mathbb{R}^{K}$ is a sparse vector and $\lVert{\boldsymbol{\alpha}}\rVert_{0}\ll K$. The corresponding high-quality patch ${\mathbf{x}}$, with size $\sqrt[3]{p}\times\sqrt[3]{p}\times\sqrt[3]{p}$, can be computed by again applying the following sparse generative model

From Eq. (2) and (3), it can be assumed that the sparse representation of a low-quality patch in terms of ${\mathbf{D}}_{\ell}$ can be directly used to recover the corresponding high-quality patch from ${\mathbf{D}}_{h}$, namely, that ${\boldsymbol{\alpha}}_{\ell}={\boldsymbol{\alpha}}_{h}$. Therefore, the reconstructed high-quality image $\hat{{\mathbf{X}}}$ can be built by applying the sparse representation to each patch ${\mathbf{y}}$ in ${\mathbf{Y}}$ and then using the estimated ${\boldsymbol{\alpha}}$ with ${\mathbf{D}}_{h}$ to obtain each ${\mathbf{x}}$, which together form the image $\hat{{\mathbf{X}}}$.

The aim is to estimate a high-quality version $\tilde{{\mathbf{X}}}$ from a given low-quality volume ${\mathbf{Y}}$. Given a test low-quality volume, for each input low-quality patch ${\mathbf{y}}$, we find a sparse representation with respect to ${\mathbf{D}}_{\ell}$. The corresponding high-quality patch bases ${\mathbf{D}}_{h}$ will be combined according to these coefficients to generate the output high-quality patch ${\mathbf{x}}$. The problem of finding the sparsest representation of ${\mathbf{y}}$ can be formulated as

where $\lambda$ balances sparsity of the solution and fidelity of the approximation to ${\mathbf{y}}$, and ${\mathbf{F}}$ is a linear feature extraction operator as in Zeyde et al. (2010). Given the optimal solution ${\boldsymbol{\alpha}}^{*}$ of Eq.(4), the high-quality patch ${\mathbf{x}}$ can be reconstructed as ${\mathbf{x}}={\mathbf{D}}_{h}{\boldsymbol{\alpha}}^{*}$. This optimisation problem can be solved using the Basis Pursuit algorithm Chen and Donoho (1994).

The complete IQT process is summarised in Algorithm (1). In this algorithm, the input low-quality volume ${\mathbf{Y}}$ is up-sampled using bicubic interpolation to provide a preliminary high-resolution volume. For each cubic patch of size $\sqrt[3]{m}\times\sqrt[3]{m}\times\sqrt[3]{m}$ from the up-sampled volume, starting from the top left corner with an overlap $p$, the mean intensity $\mu$ of the patch is computed to ensure that the dictionary represents image textures rather than absolute intensities. The sparse representation ${\boldsymbol{\alpha}}^{*}$ of the patch is then obtained by solving an optimisation problem that minimises the difference between the transformed low-quality patch and its sparse representation in the low-quality dictionary ${\mathbf{D}}_{\ell}$, subject to a sparsity constraint controlled by the regularisation parameter $\lambda$. Using the high-quality dictionary ${\mathbf{D}}_{h}$ and the sparse coefficients ${\boldsymbol{\alpha}}^{*}$, a high-quality patch ${\mathbf{x}}$ is generated. This high-quality patch, with the mean intensity restored, is placed in the corresponding location in the initial high-quality volume ${\mathbf{X}}_{0}$. After processing all patches, the final high-quality volume ${\mathbf{X}}$ is obtained.

The high-quality dictionary ${\mathbf{D}}_{H}^{(N\times s_{i}\times s_{j}\times s_{k})}$ is constructed from random noise using a standard Gaussian distribution, where $s_{i},s_{j}$ and $s_{k}$ are up-scaling factors in $i$, $j$ and $k$ directions, and $N$ is the number of dictionary atoms. The high-quality dictionary ${\mathbf{D}}_{H}$ is then encoded by $s_{i}\times s_{j}\times s_{k}$ convolution with groups $N$, followed by ReLU Nair and Hinton (2010) and $1\times 1\times 1$ convolution. Each $N$ element of the resultant code ${\mathbf{C}}_{H}^{N\times 1\times 1\times 1}$ represents each $s_{i}\times s_{j}\times s_{k}$ atom as a scalar value. Note that low-quality dictionaries can be naturally replaced by convolutional operations, and therefore only ${\mathbf{D}}_{H}$ is constructed. The ${\mathbf{D}}_{H}$ generator has a tree-like structure, where the nodes consist of two $1\times 1\times 1$ convolutional layers with ReLU activation. The final layer has a Tanh activation followed by a pixel shuffling layer. To produce $N$ atoms, depth $d$ of the generator is determined as $d=\log_{2}(N)$.

We use the UNet++ Zhou et al. (2018) as a deep feature extractor in Fig. 2, with depth of three, and a long skip connection is added. For an input image $I\in\mathbb{R}^{i\times j\times k}$, the deep feature extractor generates a tensor of size $f\times i\times j\times k$. The per-pixel predictor then takes as input a concatenation of the extracted feature and the expanded code of ${\mathbf{D}}_{H}$, such that $C_{H}^{N\times i\times j\times k}=R_{1\times i\times j\times k}(C_{H}^{N\times 1\times 1\times 1})$, where $R_{a\times b\times c(\cdot)}$ denotes the $a\times b\times c$ repeat operations. The per-pixel predictor is composed of ten bottleneck residual blocks followed by a softmax function that computes the $N$ coefficients of ${\mathbf{D}}_{H}$ for each input pixel. Both the deep feature extractor and per-pixel predictor contain batch normalisation layers Loffe and Normalization (2014) before the ReLU activation. The resultant prediction map $M^{N\times i\times j\times k}$ is further convolved with a $2\times 2\times 2$ convolution layer to produce a complementary prediction map $M^{\prime N\times(i-1)\times(j-1)\times(k-1)}$, that compensates the patch boundaries when reconstructing the final output. The detail of the compensation mechanism is described in the next subsection.

The prediction map $M^{N\times i\times j\times k}$ is upscaled to $N\times s_{i}i\times s_{j}j\times s_{k}k$ by nearest-neighbor interpolation, and the element-wise multiplication of that upscaled prediction map $U_{s_{i}s_{j}s_{k}}(M^{N\times i\times j\times k})$ with the expanded dictionary $R_{1\times i\times j\times k}$ $(D_{H}^{N\times s_{i}\times s_{j}\times s_{k}})$ produces $N\times s_{i}i\times s_{j}j\times s_{k}k$ tensor $T$ consists of weighted atoms. The $U_{s_{i},s_{j},s_{k}}(\cdot)$ denotes $s_{i}\times s_{j}\times s_{k}$ nearest neighbor upsampling. Finally, the tensor $T$ is summed over the first dimension, producing the output $x$ as

The same sequence of operations is applied to the complementary prediction map to obtain the output $x^{\prime}$. The final high-field prediction is obtained by centering $x$ and $x^{\prime}$ on top of each other and concatenating the overlapping parts of the centered $x$ and $x^{\prime}$, and applying a $5\times 5\times 5$ convolution. For non-overlapping parts, $x$ is simply used as the final output.

The number of atoms and patch-sizes in dictionaries ${\mathbf{D}}_{\ell}$ and ${\mathbf{D}}_{h}$ has impact on two important aspects of the proposed IQT-SRep model; that are the reconstruction accuracy and reconstruction time. Larger dictionaries include more image patterns, and therefore more accurate super-resolved volumes. However, the drawbacks are the computational complexity of solving the optimisation problem and the longer time required for patch extraction. Following this, from an initial set of $100,000$ 3D-vectorised patches, we learned compact dictionaries of different atom numbers, including 150, 256, 512, 1024 and patch-sizes of $3\times 3\times 3$, $5\times 5\times 5$ and $7\times 7\times 7$. We first present those of $1024$ atoms using $7\times 7\times 7$ patch-size, which provide best construction quality, and the effect of different atom number is presented afterwards.

In this work, we adopt a model using different atoms numbers of 64 and 128 atoms. The number of filters of the models is adjusted according to the number of atoms. The scaling factors $s_{i},s_{j},\text{ and }s_{k}$ are set to $s_{i}=1,s_{j}=1,\text{ and }s_{k}=4$. The network is trained using low-quality patch size of $32\times 32\times(32/s_{k})$ with a mini-batch size of 32. Random flipping and rotation augmentation is applied to each training sample. An Adam optimiser Kingma and Ba (2014) with $\beta_{1}=0.9,\beta_{2}=0.999,$ and $\epsilon=10^{-8}$ is used. The learning rate of the network except for the ${\mathbf{D}}_{H}$ generator is initialised as $2e^{-4}$ and halved at $[200k,300k,350k,375k]$. The total training iterations is $400k$. The learning rate of the ${\mathbf{D}}_{H}$ generator is initialised as $5e^{-3}$ and halved at $[50k,100k,200k,300k,350k]$. In addition, to stabilise training of the ${\mathbf{D}}_{H}$ generator, we randomly shuffle the order of output atoms for the first $1k$ iterations. The results of the 128 atoms dictionary are first presented, followed by a comparison with those of 64 atoms dictionary.

As in Lin et al. (2022), we use a default patch size of $32\times 32\times(32/s_{k})$ and $32\times 32\times 32$, respectively for low-field and high-field volumes, and a step size of $8$, $16$, and $16/s_{k}$ along $x$-, $y$-, and $z$-directions, which provide best construction quality. Training model is then constructed using the training procedure explained in Section 2.3.

Table 1 provides NRMSE and SSIM results of InD1, InD2 and OOD using the three methods IQT-SRep, IQT-DDL and IQT-DL. We can observe that the supervised deep learning approach IQT-DL provides better results (lowest NRMSE and highest SSIM) using the in-distribution datasets (InD1 and InD2), compared to the unsupervised learning IQT-SRep algorithm, revealing that supervised learning is more robust for super-resolving images that follow the same distribution of the training data set compared to unsupervised learning. However, when testing using out-of distribution data that is different from the distribution of the training samples, the unsupervised learning approach IQT-SRep provides lower NRMSE and higher SSIM compared to the supervised deep learning model IQT-DL. This highlights the importance of unsupervised learning models since the full data distribution cannot be represented in the training data set. On the other hand, we can observe that the supervised deep learning model IQT-DL performs better (lower NRMSE and higher SSIM) than the blended supervised and unsupervised learning IQT-DDL approach using InD1, whereas the IQT-DDL provides better results using both InD2 and OOD datasets. This reveals the robustness of the blended learning IQT-DDL approach in super-resolving datasets differ from that the model was trained on, in addition to data that slightly deviates from InD1 but still part of the training samples.

Figure 4 shows examples of coronal T1 weighted images from the HCP data set, corresponding to synthesised low-field images using InD1, InD2 and OOD, and results of IQT-SRep, IQT-DDL and IQT-DL. Figure 5 shows corresponding absolute error maps between high-quality ground truth images and corresponding low-quality images, and results of IQT-SRep, IQT-DDL and IQT-DL. Moreover, the binary maps of regions (in red label) where the IQT-SRep and IQT-DDL models provide closer estimates to ground truth high-quality images compared to IQT-DL are also presented. The qualitative results in general follow the same behaviour of the quantitative results described earlier: although IQT-DL provides better visual results of brain structure compared to IQT-SRep using InD1, and InD2, the IQT-SRep model shows better visual results using OOD data compared to IQT-DL. This is clearer in the absolute error maps in Fig. 5, between high-quality ground truth images and results of both IQT-DL and IQT-SRep, and in the binary maps where there are more image regions where IQT-SRep and IQT-DDL performs better than IQT-DL. This implies that the IQT-SRep approach is more robust for image enhancement using out-of-distribution data, which are created using a different distribution to that of the training samples mimicking real-world examples. Moreover, the IQT-SRep approach provides smoother outputs compared to that of the IQT-DL approach where artifacts arising from patch construction are very obvious. On the other hand, the blended learning IQT-DDL approach provides better visual results than the supervised deep learning approach IQT-DL using InD2 and OOD datasets. This is also clear in the absolute error maps and in the binary maps where there are more image regions where IQT-DDL performs better than IQT-DL in Fig. 5. On the other hand, while in this work we process data volumes by splitting them into overlapping patches, the proposed approaches ensure that information from the borders of each patch is preserved and integrated into the subsequent patches, thereby there is no information loss and the continuity of image features across the entire volume is maintained. Moreover, we synthesise the low-quality volumes from the high-quality ones, which ensures that there are no pixel alignment problems, as both low-quality and high-quality volume pairs are inherently aligned during the synthesis.

To summarise, the blended learning IQT-DDL approach provides best visual results compared to the supervised deep learning IQT-DL and the unsupervised leaning IQT-SRep approaches using both InD2 and OOD datasets, whereas the unsupervised learning approach IQT-SRep provides better visual results than the supervised deep learning IQT-DL approach using OOD which is generated using a different distribution to that of the training samples. There are widespread regions where the errors are lower for IQT-SRep and IQT-DDL compared to IQT-DL highlighting bias in the regression model estimates, which both IQT-SRep and IQT-DDL can avoid.

Now, we evaluate the effect of atom number and patch size on both approaches. From the sampled 100,000 image patch pairs, and for the IQT-SRep approach, we train four dictionaries of size 150, 256, 512, 1024, and use each to estimate the high-field image from the low-field counterpart. Moreover, for the IQT-DDL approach, in order to assess the performance of the algorithm using different atoms numbers, we construct dictionaries using atom number of 64, in addition to that of 128 whose results are presented in the previous section. Table 2 shows NRMSE of the IQT-SRep and IQT-DDL approaches using different atom number using the three testing datasets InD1, InD2 and OOD. We can observe that in general, as atom number increases, construction quality improves (NRMSE decreases and SSIM increases), but saturates for atom number higher than 512 for the IQT-SRep approach. For the IQT-SRep approach, all tested atom numbers still provide better construction results using OOD data as compared to the IQT-DL approach. On the other hand, we tested several patch sizes for both the IQT-SRep and IQT-DDL models. Specifically, for the IQT-SRep model, we tested patch sizes of P3: $3\times 3\times 3$, P5: $5\times 5\times 5$, and P7: $7\times 7\times 7$ using a dictionary size of 1024 (which provides the best results). For the IQT-DDL model, we tested patch sizes of P16: $16\times 16\times 16$, P32: $32\times 32\times 32$, and P48: $48\times 48\times 48$ using a dictionary size of 128 (which provides the best results). Table 3 provides the NRMSE and SSIM for two in-distribution datasets (InD1 and InD2) and one out-of-distribution (OOD) dataset. We can observe that for the IQT-SRep model, the reconstruction results improves (NRMSE decreases and SSIM increases) as the patch size increases. Conversely, for the IQT-DDL model, a patch size of P32: $32\times 32\times 32$ outperforms both P16: $16\times 16\times 16$ and P48: $48\times 48\times 48$, indicating that it is a good operating point, balancing structural information content with the ability to learn and generalise from a finite training set. Fig. 6 shows an example of super resolved images using the OOD data set using dictionaries of different sizes at patch sizes providing best results (P7 for IQT-SRep and P32 for IQT-DDL). While there are no substantial visual differences, we can observe in the binary error maps that the number red pixels (improvement over interpolated low-field image) gradually increase with larger dictionaries until saturation for dictionary size of 1024. In terms of computation time, the IQT-SRep algorithm is implemented in MATLAB and the experiments are carried out on a laptop with a 2.8 GHz processor CPU, with 16 GB of RAM, under Microsoft Windows 10. Dictionary construction times ranges from $\sim 25$ min for a 150-size dictionary to $\sim$ 80 min for a 1024-size dictionary. During the testing, in terms of test image reconstruction time, the computation is approximately linear to the size of the dictionary, that larger dictionaries will result in heavier computation. For example, smaller dictionaries, such as those with 150 atoms, yield reconstructions in an average time of $\sim 7$ min, while larger dictionaries, such as those with 1024 atoms, yielded image reconstructions in an average time of $50$ min. On the other hand, for the IQT-DDL algorithm, as shown in Table 2, the NRMSE is slightly lower using dictionary with atom number of 128 compared to that of 64 atoms for all testing datasets InD1, InD2 and OOD, as it retains more image patterns. Fig. 6 shows visual results of the IQT-DDL approach using an OOD example. Similar to the IQT-SRep approach, while there is no substantial visual difference, we indeed observe the increase in more super-resolved pixels (in red) in the binary error map images compared to the interpolated low-field image for atom number of 128 compared to that of 64. In terms of computation time, the IQT-DDL algorithm is implemented in PyTorch, and the testing construction time ranges from $\sim 4$ to 7 min for atom numbers of 64 to 128, respectively.