Deep Hierarchical Super Resolution
for Scientific Data

Super resolution techniques are used to increase the resolution of an image[9, 28], video [29], volume [8, 2], or time-varying volumetric data [4]. The related works below do not include image SR unless otherwise noted, and we refer readers instead to a survey by Wang et al. [30] for a comprehensive review on image SR. In general, there are two categories of SR: spatial super resolution (SSR) and temporal super resolution (TSR). TSR increases the temporal resolution of input data, such as TSR-TVD by Han et al. [4], which uses a convolutional LSTM to learn the recurrence between timesteps and recovers interpolated timesteps with higher accuracy than linear interpolation. TSR is also useful for video frame interpolation [29].

Our approach is an SSR technique, which increases the spatial resolution of input data. Convolutional neural networks (CNNs) were adopted for image SR, with specific techniques introduced to improve quality such as residual learning [29, 31], adversarial training [32, 9], and removing bias from interpolation techniques used in the networks [28, 33]. Similar to our approach, LapSRN [28] estimates multi-level upscaled output to allow $2\times,4\times,8\times,...$ scale factor SR. However, their network trains end-to-end using a Charbonnier loss function and assumes the input is at the lowest resolution, whereas our hierarchy trains each model individually which allows SR from any starting resolution below the training data resolution. Weiss et al. use an SR network to increase the resolution of isosurfaces by upscaling the depth and normal maps [34]. Weiss et al. also use an SR approach for adaptive sampling SR for isosurface and volume rendering images [35].

Image-focused models were adapted to accommodate SSR of 3D scientific data as well, but no methods have yet been proposed for hierarchical data SSR. Zhou et al. [1] use a 3D convolutional neural network (CNN) to perform SSR on volumes with better feature reconstruction than trilinear interpolation or cubic-spline interpolation. Guo et al. [3] use three NNs to perform SSR on 3D vector fields by training a single network for each variable. Xie et al. create tempoGAN [36], which upscales fluid flows for temporally consistent HR output. Fukami et al. compare two ML-based SSR methods for 2D fluid flow [6]. Höhlein et al. [37] use CNNs to perform wind field downscaling, which learns the mapping from coarse input data to fine predicted output data. SSRTVD by Han and Wang [2] uses a generative adversarial network (GAN) for upscaling volumes such that the upscaled output frame is temporally coherent with adjacent timesteps. With the recent advances in SR, Jakob et al. have created a public 2D flow dataset with varying Reynolds number, citing that ML methods are a powerful method for interpolating flow maps and data reduction [38].

There are also works that use ML for spatiotemporal SR, wherein both space and time dimensions are upscaled [7, 5]. These networks use a combination of convolutional layers for spatial upscaling, and may use convolutional and/or recurrent layers for temporal upscaling.

Hierarchical data structures have been used in applications including rendering, scientific simulations, data reduction, and ML, to reduce computation and storage requirements by more efficiently allocating resources available. A survey of octree-based rendering methods is available by Knoll [18]. Rendering pipelines have utilized octrees [39] for efficient and scalable volume rendering by dynamically loading data from different octree depths depending on the current camera position and angle. Gobbetti et al. [40] organize large scalar fields into an octree structure and only feed the GPU information relevant to the current viewpoint and transfer function at each frame. Knoll et al. [16] create a technique for isosurface ray tracing for large octree volumes. Wilhelms and Van Gelder [17] use octrees for isosurface generation when the regions of interest are unknown. Hadwiger et al. [41] enable visualization of petascale volumes of microscopy data using an out-of-core multiresolution approach. Fogal et al. [42] use an out-of-core adaptive multi-grid rendering algorithm that samples densely only where needed to enable rendering of large scientific data on consumer graphics cards. A key problem for rendering hierarchical data is the seam artifact that can occur between data blocks. Methods for reducing the effect of this issue are presented by Santos et al. [24], Wald et al. [25], and Wang et al. [26].

Aside from rendering, hierarchical data formats have also been used during simulation to improve performance. Losasso et al. [43] create a method for simulating water and smoke on an unrestricted octree to capture small scale visual detail and allow for efficient solving of the Poisson equation. Popinet [44] uses quadtrees and octrees for a flexible and efficient approach to solving time-dependent incompressible Euler equations. Adaptive mesh refinement (AMR) schemes, a similar hierarchical data representation, have been utilized in numerous scientific software packages to improve performance by using a coarse mesh where there is little detail and fine meshes where precision is more important. Berger and Colella first design the AMR scheme for use with hyperbolic partial differential equations [45] and two dimensional shock dynamics [46]. Different implementations of AMR are implemented in modern simulation packages such as Enzo[15], Chombo[13], LAVA[14], Nyx [11], and AMReX[12].

Hierarchical representations can also be used for data reduction. Knoll et al. [16] use an octree representation that only uses a fraction of the original data’s memory footprint, similar to the approach of Velasco and Torres [47], to accelerate volume rendering. Bhatia et al. [48] introduce AMM, an adaptive multilinear mesh framework that reduces the memory footprint of raw data using tensor products of linear B-spline wavelets and allow a tradeoff between numerical precision and spatial resolution. Hoang et al. [23] use a precision-resolution tree, and an encoding of the tree, to perform data reduction on scalar field data. Ainsworth et al. [49] create a multilevel technique to compress univariate data that allows users to select a level that minimizes memory use while meeting some tolerance level. Another multilevel technique is presented in MGARD [50], which uses adaptive error based coefficient quantization to enable different tolerances at different levels of detail, and proposes to use the multilevel data decomposition as a preconditioner to terminate at an appropriate level.

Hierarchical data representations have also been used in ML. However, our approach differs as it is a SR approach for hierarchical data representing scalar fields. Riegler et al. [51] create octree convolution, pooling, and upscaling layers for a NN to classify occupancy octree data representing geometric objects. This method is limited to only occupancy octrees, and not octrees for scalar fields, which we use. Tatarchenko et al. [52] create a deep convolutional decoder that learns to decode dense input to octree representations of 3D objects, which faces the same limitation as occupancy octrees. Martel et al. [22] and Takikawa et al. [53] use hierarchical data representations to assist with implicit neural representation.

In this section, we detail the construction of an SR NN hierarchy. Since the hierarchical data may be at varying levels of detail, we require a design that can upscale various scales of data. For example, some regions may need $2\times$ SR, and others may need $16\times$. Hierarchical NN approaches, such as LapSRN [28] or EDSR [54] may offer multiple scale factors, but expect specific input/output levels of detail, which is incompatible with our SR method described in Section 3.4.

We design a hierarchy compatible with our hierarchical SR approach that is used for uniform scale factor SR at varying scale factors such as $2\times,4\times,...$, scaling with the number of NNs in the hierarchy. Our design works with our hierarchical SR algorithm outlined in Section 3.4 by allowing arbitrary input and output level of detail. Our design creates a hierarchy of SR networks that allows a volume of any level of detail as input, and produces an upscaled volume of any level of detail between the input and the original HR as output. Instead of level of detail, we define “downscaling level”, such that downscaling level $i$ corresponds to raw data that has been downscaled by a factor of $2^{i}$, or simply data that is a factor of $2^{i}$ more coarse than the finest data in the volume. As the downscaling level increases, the data become coarser. We create a hierarchy $G$ of $N$ SR networks $G=\{G_{0},G_{1},...,G_{N-1}\}$, where each network $G_{i}$ is trained to perform $2\times$ SR from downscaling level $i+1$ to $i$. By using multiple networks in tandem, larger scale factors like $4\times,8\times,...$ up to $2^{N}\times$ are possible. An example of an SR hierarchy is shown in Figure 1, which shows how coarse input is super-resolved as we move up the hierarchy to the finer scales. Any SR architecture capable of $2\times$ SR can be used in our defined hierarchy for hierarchical SR, and every network within $G$ is identical in terms of architecture, but no weights are shared between models. In this paper, we evaluate our hierarchy with ESRGAN [9], SSRTVD [2], and STNet [5].

To model the hierarchical data, we define an octree data structure called an SR-octree. Our SR-octree is designed to minimize the number of nodes in the tree and keep data with the same downscaling level contiguous until spatially divided. Though everything will be defined in three spatial dimensions (for an octree), the same methods apply to 2D (for a quadtree). For brevity, we will use the terms “octree” and “SR-octree” regardless of the dimensionality of the data.

Like most octrees, children of a parent node in an SR-octree represent the 8 octants that compose the parent. The octree data structure only defines the spatial partitioning and not spatial averages, where the root node refers to the entire domain. Unique to the SR-octree is a dwnscl_lvl attribute on each node that represents the factor that data have been downscaled in each spatial dimension within that node. This attribute allows us to represent contiguous large regions with a uniform downscaling ratio with one node. Specifically, each leaf node in the SR-octree holds a volume of data defined by the node’s spatial partition (position and size) and downscaling level. This is necessary to upscale the data properly in our hierarchical SR algorithm.

Figure 2 shows the concept of an SR-octree using a quadtree example in 2-D. Non-leaf nodes hold no data and only point to their child octants. Nodes are colored according to their downscaling level shown by the key on the left. We also show an equivalent “typical” octree (as defined by Meagher [39]) in the bottom right that has more dense nodes in the finer regions. Notice that the mesh density corresponds to the color (downscaling level) in the equivalent SR-octree on the left compared with the typical octree. Our SR-octree construction groups together octree nodes at the same downscaling level into a single node.

In some use cases and for testing purposes, it is necessary to downscale a uniform grid HR volume to an SR-octree, such as for our uses cases for data reduction (Section 4.4.1) and hierarchical SR for FTLE data (Section 4.4.2). To create an SR-octree from a uniform full resolution volume, we hierarchically reduce a full-resolution volume according to a user set error bound $\epsilon$, as follows. Starting with the original volume, downscale the volume by $2\times$ in each spatial dimension, and see if the $L^{\infty}$ norm of the downscaled version compared with the original volume is below the threshold $\epsilon$. If so, keep the volume downscaled by $2\times$ and increment the downscaling level of the volume’s octree node by 1. Otherwise, return the volume to the higher resolution version, split the volume into its octants, and repeat the above process for each suboctant if each dimension of the suboctant is larger than some minimum size min_chunk. This process will maximally hierarchically downscale the volume according to the error $\epsilon$. We add other downscaling parameters max_dwnscl_lvl, which controls the maximum downscaling factor, and min_dwnscl_lvl, which will automatically downscale the original volume to that scale factor before performing the hierarchical downscaling. After reduction, we join any adjacent nodes that have the same downscaling level so that the hierarchical data is represented with the fewest number of nodes as possible.

The resulting SR-octree has volume data in leaf nodes, each having an associated dwnscl_lvl attribute for the volume it represents. This allows us to easily determine what SR scale factor is needed to restore an octree node to its original full resolution. We note that in visualizations of an SR-octree, the coarseness/fineness of a region is not determined by the density of the grid in a region (such as with typical octrees), but instead is represented with the dwnscl_lvl attribute, which we visualize with a discrete color scale in Figure 2.

Since no methods currently exist for using SR NNs for upscaling hierarchical data, we contrast our algorithm with another baseline approach. The baseline solution to upscale an SR-octree upscales each octree node individually based on the dwnscl_lvl attribute each octree node is saved with (i.e. upscales by $2^{\texttt{dwnscl\_lvl}}\times$), and then stitches together the resulting patches, illustrated in the top portion of Figure 3.

There are two drawbacks to this method that lead to poor performance with trained neural networks and the introduction of seam artifacts. The first problem is that convolutional neural networks use padding within convolution layers, which introduces incorrect neighborhood information to the trained neural network, causing error in the upscaled result. In convolutional neural networks, padding is a technique that will surround an input volume with a layer of values around the boundaries before a convolution operation. In practice, padding is used to ensure that the output resolution of the convolutional layer within the NN is the same as the input resolution. Without padding, a convolution layer with a stride of 1 would reduce the spatial resolution by $2\times\lfloor\frac{\texttt{kernel\_size}}{2}\rfloor$ voxels in each direction. The network’s output is then incorrectly influenced by the padded values that are not actually there. When composing multiple convolutional layers in sequence, the effect is increased. Due to the error from padding, some evaluations for proposed SR NNs choose to ignore some number of border pixels when calculating the recovered image’s PSNR, SSIM, or perceptual quality [55, 56, 57, 58, 59].

Closely tied to the padding problem is the second problem, which is the lack of shared neighbor information. Deep convolutional neural networks will often have large receptive fields, which indicate the window size that influences each output voxel. The larger the receptive field, the more global information is used within the neural network to determine the estimated output voxel value. The receptive field is determined by the neural network architecture (number of convolutions and the stride/kernel size/padding, etc.). Therefore, a hierarchical upscaling algorithm that shares or copies a static number of neighboring voxels/neighbor octree nodes, such as work by Ljung et al. [24] or Wang et al. [26], would only work for networks that have a receptive field smaller than the number of shared voxels.

In our setting specifically, padding and limited shared neighborhood information will have a pronounced impact for SR-octree nodes that have few voxels. For example, a $2^{3}$ sized leaf node in an SR-octree would be dwarfed by the total number of padded voxels from up to 30 padding operations in our deep convolutional neural networks. Additionally, in this baseline approach, each node will use its own padding operations, so the upscaled results between each node will not align well, introducing seams between data chunks which can be distracting and/or misleading in visualizations.

In this section, we present our approach for hierarchical SR. We design an approach that uses the entire domain at each downscaling level during SR so that padding does not create seams within the volume, and so that the receptive field of the neural network is not limited when upscaling, both of which are limitations outlined in Section 3.3. Though data are distributed across many different downscaling levels, the only level that can uniformly represent the hierarchical data is the coarsest level. Converting the hierarchical data to this uniform representation is the first part of our algorithm and is called the downscaling process, outlined in Section 3.4.1. Next, our upscaling process (detailed in Section 3.4.2) will uniformly super resolve the uniform LR data without introducing seam artifacts and overwrite voxels that have more accurate information in the SR-octree, one level at a time, beginning with the coarsest level. In the end, we are left with an approximation of the HR uniform volume.

In our proposed NN hierarchy, any spatial SR model capable of $2\times$ upscaling is compatible. In this paper, we use ESRGAN [9], SSRTVD [2] and STNet [5], which are all state-of-the-art neural networks for SR. Small changes were made to the models for compatibility and training stability. For all three models, we find they train more stably and with better reconstruction without the discriminator network(s) or losses, so each of the three architectures are only regressors that upscale input. Additionally, we adjust the models to only perform $2\times$ upscaling instead of $4\times$ upscaling. This is done by using only a single pixel/voxel shuffle [33] at the end of the model to upscale the features by a factor of 2 in each dimension instead of what the original architectures suggests. Regardless of the architecture chosen for the NN hierarchy, each generator network $G_{i}$ in our constructed neural network hierarchy is responsible for learning $2\times$ upscaling on input data of downscaling level $i+1$, generating downscaling level $i$. Our loss function for any network during training is the L1 loss $\mathcal{L}_{1}(V,V^{GT})$, where $V$ is the upscaled result from the model, and $V^{GT}$ is the ground truth data.

We experiment with seven scalar field datasets. Three are from Johns Hopkins Turbulence Databases (JHTDB) [60] that are time-varying results from a direct numerical simulation (DNS) for turbulent fluid flow. Information about each dataset is shown in Table I. Isotropic3D is a 3D velocity magnitude dataset from a DNS of isotropic turbulent fluid flow at Reynolds number $R_{\lambda}\sim 433$, hosted in the “isotropic1024coarse” dataset in the JHTDB. We take the middle z-axis slice (z=512) of the Isotropic3D magnitude dataset to create an isotropic2D dataset for testing our approach with 2D data. The mixing3D dataset is the velocity magnitude field from the “mixing” dataset in JHTDB. The solar plume dataset is a velocity magnitude field from a simulation that examines the effect of the solar plume on the heat, momentum, and magnetic field of the sun. The vorts dataset is a $128^{3}$ vorticity magnitude scalar field dataset from a pseudo-spectral simulation of vortex structures. Nyx is a cosmological simulation which is run with parameter settings uniformly sampled within suggested ranges by domain scientists: total matter density $OmM\in[0.12,0.155]$, total density of baryons $OmB\in[0.0215,0.0235]$, and Hubble constant $h\in[0.55,0.85]$ [11, 61]. Heated flow is a 2D fluid simulation from the Gerris flow solver [62] created by Günther et al. [63] that simulates flow around a heated cylinder with Boussinesq approximation. All datasets are linearly scaled to [0.0,1.0] except for Nyx which is log scaled before being scaled to [0.0, 1.0]. All datasets are represented at single precision floating point accuracy to accommodate the architectures used.

Three NN hierarchies are trained per dataset using the ESRGAN, STNet, and SSRTVD architectures. To train $G_{i}$ in a hierarchy $G$, a ground truth volume $V$ is downscaled to LR $V^{\textbf{LR}}$ and HR $V^{\textbf{HR}}$ (if necessary, i.e. $i$ is not 0) by scale factors $2^{i+1}$ and $2^{i}$, respectively. Then, the loss function is performed on $G_{i}(V^{\textbf{LR}})$ and $V^{\textbf{HR}}$.

The number of models trained in each hierarchy is shown in Table I, column $|G|$, and the total number of iterations of training (with batch size 1) is listed in column “iters”. We use Python 3.9 with PyTorch to train on an NVidia A100 40GB Tensor Core graphics card. The full resolution 3D data cannot fit in the GPU’s memory during training, so we crop training volumes to $96^{3}$ and use a random starting position for the cropping, augmented with random flipping in any subset of the spatial dimensions to further increase training diversity.

The SSRTVD models are constructed as implemented by the paper with no changes to the number of kernels per convolution. For our ESRGAN and STNet models, all layers (aside from output) use 96 kernels for a good mix between accuracy, GPU memory use, and training time. All convolution operations in all networks use reflection padding during training and inference.

The training time per iteration and the storage of each model are listed in Table II. Since each architecture is the same size across all datasets, and the inputs are cropped to a fixed resolution, the models train at a similar time per iteration regardless of which level of detail they are training on. STNet trains the quickest, as it is the most lightweight architecture. ESRGAN trains quicker than SSRTVD, yet is a larger model with more parameters, likely due to SSRTVD’s decision to be a wider network with more kernels per convolution instead of deeper network, which uses more convolution layers like ESRGAN does.

We compare our proposed approach in two steps. In Section 4.3.1, we isolate the neural network hierarchies by comparing the uniform grid SR performance of the hierarchies against bilinear/trilinear interpolation. In Section 4.3.2, we isolate our hierarchical SR algorithm by comparing the results of hierarchical SR with our method against the baseline blockwise approach.

We provide examples of three use cases for our hierarchical SR approach. In Section 4.4.1, we show how our hierarchical SR method can be used to improve STNet’s data reduction capabilities by decomposing the spatial domain based on the presence of features instead of uniform downscaling. Section 4.4.2 shows how computation resources can be saved on the expensive FTLE field computation by increasing tolerated error in hierarchical FTLE calculation algorithms and using our hierarchical SR algorithm to approximate the HR uniform FTLE. Lastly, Section 4.4.3, we demonstrate that a trained hierarchy can be used to upscale LR simulation output for visualization, which can be used to save computation resources before running the higher resolution simulation.