Topology Estimation of Simulated 4D Image Data by Combining Downscaling and Convolutional Neural Networks

Data design experiments were focused on choosing radius parameters (namely, $a$, $r$, and $R$) that would generate samples with non-trivial topologies in a resolution that would ensure that holes were represented clearly (in a homological sense). Since the toxels of an image were attributed integral coordinates, persistent homology software would, theoretically, be capable of correctly detecting holes of any dimension, provided that the diameter of a hole was greater than the distance between diagonal points of a 4D unit-cube ($\sqrt{4}$ units). Otherwise, calculations would suffer as a result of there being insufficient resolution to describe a hole with such a small radius. Conversely, if parameters were too large, then the cavities would be so big that we would be limited to samples with fewer holes and less interesting topologies.

Table 1 provides several formulas that may be used to depict a variety of 4D objects in an $(x,y,z,w)$-system, and were used to design the cavities for the dataset. These objects vary in their geometry (relative to each other), and this can be demonstrated by experimenting with the parameters in each formula. For example, $B^{4}$ does not have a tunnel, whereas $S^{1}\times S^{1}\times B^{2}$ does have a tunnel and can vary from being quite ‘flat and expansive’ to being more ‘round’, depending on the choice of the parameter $\alpha$, which sets the orientation of the $S^{1}\times B^{2}$ factor and ranges from $0$ to $\pi/2$. Figure 1 shows several 2D visualisations that offer some intuition of how varying $\alpha$ can impact the resulting embedding of $S^{1}\times S^{1}\times S^{1}$ (and, equivalently, $S^{1}\times S^{1}\times B^{2}$) in $\mathbb{R}^{4}$. The idea is to begin with a (dark-grey) torus $S^{1}\times S^{1}$ that is oriented according to $\alpha$, and positioned $R$ units from the origin along the $x$-axis in the 3D $xzw$-hyperplane. The torus is then rotated around the origin, through the $xy$-plane, in order to introduce the third $S^{1}$ factor. Although the figures may suggest otherwise, the implicit formula for $S^{1}\times S^{1}\times B^{2}$ that is used to generate our data guarantees that overlaps or self-intersections do not occur. This is because we are working in $\mathbb{R}^{4}$, rather than $\mathbb{R}^{3}$, as the figures may also suggest; the extra dimension cannot be shown easily in 2D. Figure 1a demonstrates a construction in which $\alpha=0$, and Figures 1b and 1c demonstrate constructions in which $\alpha=\pi/2$; notice that Figure 1c is, in fact, a $\pi/2$ radians $zw$-rotation of Figure 1b. Because of the symmetry of the torus that we begin with, any $zw$-rotation of the construction in Figure 1a is inconsequential, as is a $\pi$ radians $zw$-rotation of the remaining examples. In practice, we only need to consider when $\alpha$ ranges from $0$ to $\pi/2$; the remaining options arise freely from the random rotations that are applied during data generation.

In order to maintain some homogeneity in the range of sizes of these objects, we selected parameters for each object that would produce cavities with a similar range of hypervolume (4D-volume). Formulas for the hypervolumes of these objects, along with some simplifications that arise by setting $R=2a$ (for $S^{1}\times B^{3}$ and $S^{2}\times B^{2}$) are also provided in Table 1. For $S^{1}\times S^{1}\times B^{2}$, we assume that $r=2a$, and that the value of $R$ depends on $\alpha$ and ranges from $2a$ to $4a$. Therefore, for $a\in[a_{\mathrm{min}},a_{\mathrm{max}}]$, the hypervolume of $S^{1}\times S^{1}\times B^{2}$ ranged from $16\pi^{3}a_{\mathrm{min}}^{4}$ to $32\pi^{3}a_{\mathrm{max}}^{4}$.

The hypervolumes of the remaining objects were scaled into the same range by finding suitable values for $a_{object}$. For example, if the hypervolume of $B^{4}$ were to also fall within this range, it was necessary to choose $a_{B^{4}}$ such that $\frac{\pi^{2}}{2}a_{B^{4}}^{4}\in[16\pi^{3}a_{\mathrm{min}}^{4},32\pi^{3}a_{\mathrm{max}}^{4}]$. Rearranging this expression leads to Equation 1, and the remaining ranges were deduced in the same way.

The parameters that we selected in order to generate a 4D dataset with which to investigate our approach are summarised in Table 2. Our study focused on recognising the global topology of compact manifolds. Hence, we restricted our experiments to single-component samples ($\beta_{0}=1$), and ensured that both a cube’s boundaries were not disturbed and that cavities did not intersect each other. These rules were enforced by implementing a 1-toxel boundary and minimum 6.5-unit spacing between cavities. The self-intersection of a cavity was prevented via the simplifications that were explained in Section 3.2. The choices for $a_{\mathrm{min}}$ and $a_{\mathrm{max}}$, as listed in the $B^{2}$ row of the $S^{1}\times S^{1}\times B^{2}$ column of Table 2, were sufficient to produce cavities with a nice resolution.

Each sample comprised of a 4D $128^{4}$ cube, with the combined number of 1D, 2D, and 3D holes ranging from 1 to 48. A $128^{4}$ cube was used because it was large enough to contain a non-trivial range of cavities, which afforded the analysis of samples with interesting topologies. This data was also slightly bigger than what would have been feasible to be directly analysed using persistent homology methods or CNNs with our hardware (NVIDIA DGX Station, with an Intel Xeon E5-2698 v4 CPU, 256GB RAM, and four V100-32GB GPUs). Hence, downscaling became a requirement in order to process this data. Note also that the cavity dimensions were small enough that the application of downscaling could potentially disturb the homology of a sample by closing up holes (see Section 4.1).

The data was generated in parallel on a High Performance Computing (HPC) Grid in 100-sample batches, over 320 nodes. An average of 8.60 hours was required to complete each batch, and the entire process utilised approximately 2753.21 HPC hours.

A visualisation of the cavities within a $128^{4}$ sample is shown in Figure 2. This is achieved by inverting the toxel values (setting 0 to 1, and vice versa) so that the cube itself is stripped away in order to reveal the objects that were used to produce its cavities, and then taking 3D slices along some axis. In this case, 18 equally-spaced slices have been taken along the $w$-axis. Taking finer slices allows one to see a more continuous-looking evolution of the cavities (see Appendix A).

Each label was produced on-the-fly during the generation of a sample. The homology of each cavity that was introduced into a sample was algebraically derived by, firstly, observing that the 4D cube $I^{4}$ is homeomorphic to the 4D ball, and therefore shares the same homology. Secondly, the homology of both the object being considered for removal $M$ and its boundary $\partial M$ were computed by using the Künneth theorem. Thirdly, the Mayer-Vietoris Sequence was applied to find the homology of the cube with the cavity $I^{4}-M$; an introduction to the theorems that were used in this derivation can be found in Hatcher (2002). The results of these calculations are provided in Table 3. Since the Betti numbers that each cavity contributed to a sample were known, they could simply be summed over their dimensions in order to produce a label for the sample. Furthermore, the parameters outlined in Table 2 were large enough to always produce a homologically correct representation in the $128^{4}$ cube setting. Consequentially, there was no requirement to analyse the samples using persistent homology software in order to produce a label. A more detailed description of the mathematics involved in producing the labels, can be obtained from Hannouch and Chalup (2023).

The data were downsampled from $128^{4}$ to $32^{4}$ by taking the voxels at every 4^th coordinate along each axis; the resulting images were 256 times smaller. Figure 3 gives an indication of how course-looking a 3D slice of a 4D sample becomes after the application of this degree of downsampling. Labels were not modified, however, the impact of downsampling was investigated using GUDHI. This was performed in parallel on a HPC Grid in 50-sample batches, over 640 nodes. An average of 0.66 hours was required to analyse each batch, and the entire process utilised approximately 420.14 HPC hours. As summarised in Table 4, all of the samples remained as single-component images ($\beta_{0}=1$), however, the number of samples with a structure that was consistent their label for $\beta_{1}$, $\beta_{2}$, and $\beta_{3}$ was markedly affected, with only 16065 (50.2%), 14154 (44.23%), and 4121 (12.88%) samples matching their label for each respective Betti number. Only 2175 (6.8%) of the samples still possessed a structure that was completely consistent with their label; this result is recorded in the ‘complete match’ column.

The downsampled data was then used to train a 4D CNN, which was implemented with PyTorch, and had a similar architecture to the 2D and 3D CNNs that were used by Paul and Chalup (2019); these consisted of several convolution layers, followed by a max-pooling layer, and finishing with a sequence of fully connected linear layers. Our models began with four iterations of a module consisting of: a 4D convolution layer, followed by a ReLU activation function, and then a 4D max-pooling layer. The convolutional layers output 8, 16, 32, and 64 channels, respectively, with 2 units of padding and a $5^{4}$ kernel. The pooling kernel was $2^{4}$ units. After the fourth convolution module, the result was flattened, and then passed through two fully connected layers that were separated by a ReLU operation. Finally, the result was output to a sparsely coded vector by reserving one output neuron for each possible value of $\beta_{n}$. Based on the design of our dataset, we accommodated for the values 0 and 1 for $\beta_{0}$, and accommodated for the values 0 to 16 for $\beta_{1}$, $\beta_{2}$, and $\beta_{3}$.

Full-scale deep learning experiments were performed on the above-mentioned NVIDIA DGX Station using a multi-GPU arrangement. A high-bandwidth connection between the GPUs, and between the GPUs and CPU, via NVLink, made it possible to consider these devices as a single, larger, computing element, which was able to accommodate the CNN and allow for a 192-sample batch size.

The dataset was randomly divided into 90% training, 5% validation, and 5% test sets at the beginning of each experiment. For each epoch, the samples were rotated in random multiples of 90 degrees through a randomly selected coordinate plane as they were fed into the Pytorch dataloader; this offered twelve possible variations on each sample. The Cross Entropy loss function was appropriately set up to handle the four separate outputs for $\beta_{0},\ \beta_{1},\ \beta_{2}$, and $\beta_{3}$. The Adam optimiser was initialised with a learning rate of 0.001, and a scheduler was employed to reduce the learning rate by a factor of 10 at epochs 160 and 190 over a 200 epoch training schedule.

Table 5 presents the test set accuracy average $\mu$ and standard deviation $\sigma$ that were achieved in five repeats of the experiment. Each experiment required just over 4 days (approximately 98 hours) to complete, and utilised its own randomly selected training set, validation set, and test set in the proportions set out above. An average combined test set accuracy of 82.41% was achieved.

Similar experiments were performed with another downscaled dataset that was produced by reducing the same $128^{4}$ samples to $32^{4}$ via 4D average-pooling; the software that was used to perform this version of downscaling was implemented with Pytorch. The labels were, again, left unaltered. The resulting samples were a smaller, blurry, grey-scale image of the original cube. Figure 4 depicts a sequence of 2D slices, each taken from the $32^{4}$ cubes that were produced by downsampling and average-pooling the same 4D sample that was the subject of Figure 3. For comparison, Figure 4 shows equally-spaced 2D slices (taken from bottom to top) of the same downsampled $32^{3}$ 3D slice that is seen in the left of Figure 3, and compares these slices with their corresponding average-pooled slice. The average-pooled slices appear to be less coarse and richer in features, versus their corresponding downsampled slice. For example, the first downsampled slice is empty, however, the average-pooled slice contains evidence of a cavity; the fifth, sixth, and eighth slices demonstrate something similar. The seventh and tenth downsampled slices contain pairs of features that are actually part of the same cavity, which we deduce by inspecting the corresponding average-pooled slice. This occurs because, while downsampling only collects the value (0 or 1) of the voxels at every fourth coordinate, the $4^{4}$ average-pooling kernel (with a 4-unit stride) takes an average of the $256$ voxel values that it sees.

The dataset division percentages, CNN architecture, and scheduler that were detailed in Section 4.1 were also used for these experiments, however in this case, the Adam optimiser was initialised with a learning rate of 0.0001 in response to inconsistent results that were observed during several abbreviated preliminary test runs with a learning rate of 0.001. Despite the smaller learning rate, the scheduler performed the same 200 epochs. The final results of these experiments are presented in Table 6, along with some statistics. An average combined accuracy of approximately 78.52% was observed in these experiments.

In order to understand how well the CNN was able to perform, several 1000-sample datasets of $32^{4}$ cubes were again generated from $128^{4}$ cubes via the downsampling approach described in Section 4.1; these data were therefore new to the trained models. As expected, an analysis of these datasets, using both GUDHI and a trained CNN, produced similar results to those presented in Tables 4 and 5. Despite the significant inconsistency between the structure of these downsampled images and their labels (only 7.5% of the new samples exhibited a complete label match after downsampling, as demonstrated by GUDHI), the CNN was still able to achieve a complete match for over 65% of the samples, and more accurately estimated the respective Betti numbers of each sample versus the results produced via persistent homology (see Table 7).

Benchmarking the two proposed downscaling options was not the aim of this work, however, the CNN architecture, and many of the training hyperparameters, were left unchanged across the training experiments. We observed that average-pooling produced a slightly better $\beta_{3}$ estimation, which, possibly, could have been a consequence of the richer-looking data (as demonstrated in Figure 4), although, this CNN also performed worse when estimating $\beta_{1}$. We also note that the randomness that was applied during dataset generation and subsetting, or changing the CNN architecture, could have produced different results. Overall, the test set accuracy that was achieved in the downsampling experiment was comparably similar to the results that were achieved in the average-pooling experiment (82.41% versus 78.52%), therefore, both options could be equally as useful.

The results and insights that have been presented in this work extend the approach that was pioneered by Paul and Chalup (2019) in the 2D and 3D setting. Collectively, these results have, primarily, been restricted to synthetic, single-component samples. More complicated features, such as multiple-components, links, and connected sums, or the use of topology-preserving deformations to vary the geometric appearance of cavities, are yet to be considered in this line of research. Addressing this would be essential in order to apply the CNN approach more generally to the task of estimating homology; it is likely that real data may possess these features. For example, the random deformation of the canonical embeddings that are used in our experiments could be accomplished by developing a 4D version of the repulsive tangent-point energy algorithm that was proposed by Yu et al. (2021b, a) (although, we concede that this method itself would be computationally expensive, and would not protect against the homological consequences of producing deformed cavities with self-intersections).

Although we have demonstrated that it is possible to analyse large 4D samples with CNNs where using existing options may not be feasible, it is apparent from our work that the process of generating synthetic data with which to train a CNN comes at a significant resource and time cost that may not be accessible to everyone. Choosing how to generate data that models real data may also be a necessary preliminary step, which brings with it a new set of challenges. In our case, we faced the complication of finding a balance between the size of the data that we considered, and the size of the CNN that we used. The decision to use a relatively simple CNN architecture in our experiments was made in order to demonstrate how readily CNNs could be applied to our task, but was also a result of being confined to the VRAM capacity of the available hardware; introducing additional layers, connections, or training parameters would potentially increase the hardware requirements. In the near term, performance gains could be achieved through improved software engineering. In the long term, improvements may come in the form of cost reduction and hardware advances, which would certainly make exploring deeper, wider, or more sophisticated 4D CNN architectures, similar to the many well-established options that are available in lower-dimensions (Simonyan and Zisserman, 2014; Szegedy et al., 2015; He et al., 2016), more tractable; the hope would be to determine which architectures cope best with topological applications, such as estimating Betti numbers.

As previously alluded to, there are several immediate directions into which the results that are presented in this work could be extended in future work. The downscaling results could be expanded by implementing, and then exploring the use of, the 4D equivalents of other image downscaling algorithms, such as those that were mentioned in Section 2. The downscaling approach that we took in this work could also be compared and combined with multi-view methods, which employ lower-dimensional representations of data that are produced by gathering lower-dimensional images from different angles (perspectives), or by projecting 4D data onto a lower-dimensional space; this could potentially offer some speed or memory advantage. This would contrast volumetric approaches that process 4D data using 4D operations, such as those applied through a 4D CNN (see Section I.B of Cao et al. (2020) for a discussion about multi-view and volumetric approaches). In the 3D setting, Su et al. (2015) demonstrate how CAD and voxel data may be projected onto two dimensions by using an ‘outside’ perspective to capture several images of the data from different angles in a similar way to taking X-ray scans of a subject along three orthogonal axes. Alternatively, Shi et al. (2015) project 3D samples outwards from their centre onto an annulus that wraps around the object. Similar ideas have been employed by Qi et al. (2016) and Kanezaki et al. (2018). Nevertheless, Wang et al. (2019) suggest that standard volumetric approaches may be more capable of gathering information when compared to multi-view models, and argue that this is because multi-view strategies often fail to encode information from different views. It may, however, be easier to implement larger models in lower dimensions, or exploit pre-trained models by fine-tuning (Russakovsky et al., 2015; Brock et al., 2016), which could potentially be useful when handling large 4D data.

More generally, several other matters also require deeper investigation in order to fully assess the capabilities of CNNs in estimating topology. For one, the task of modelling a real dataset for the purpose of training a CNN would need to be addressed. Secondly, generating this data would need to be carefully controlled, possibly via some form of topology-preserving algorithm; this itself raises the question of whether generating a diverse training dataset from the outset, that is, a dataset that comprises of more complicated features, such as multiple components, links, and connected sums, could be enough to train a CNN that is capable of general homology estimation, without the need to understand any properties of the real dataset under consideration. These aspects could, potentially, benefit from a (Bayesian) statistical approach, such as those considered in zero-shot, one-shot, or $N$-shot learning models (Fe-Fei et al., 2003; Palatucci et al., 2009), which are capable of generalising knowledge to unfamiliar cases after seeing little, or no, training examples without requiring extensive retraining, or where complete training may not be possible due to dataset limitations, or because real data may be infinitely-variable (as may be the case with homeomorphic deformations).