Image recovery from unknown network mechanisms for DNA sequencing-based microscopy

Imaging-by-sequencing methods are an emerging alternative to conventional optical micro- or nanoscale imaging. In these methods, molecular networks form through proximity-dependent association between DNA molecules carrying random sequence identifiers. DNA strands record pairwise associations such that network structure may be recovered by sequencing which, in turn, reveals the underlying spatial relationships between molecules comprising the network. Determining the computational reconstruction strategy that makes the best use of the information (in terms of spatial localization accuracy, robustness to noise, and scalability) in these networks is an open problem. We present a graph-based technique for reconstructing a diversity of molecular network classes in 2 and 3 dimensions without prior knowledge of their fundamental generation mechanisms. The model achieves robustness by obtaining an unsupervised sampling of local and global network structure using random walks, making use of minimal prior assumptions. Images are recovered from networks in two stages of dimensionality reduction first with a structural discovery step followed by a manifold learning step. By breaking the process into stages, computational complexity could be reduced leading to fast and accurate performance. Our method represents a means by which diverse molecular network generation scenarios can be unified with a common reconstruction framework.


1.A.2 Binary unweighted networks versus networks with weighted edges
Spatial graphs are built from a set of points located in a space with a particular metric, most commonly Euclidean. In this work, we studied several ways to build such graphs by establishing proximity rules. This resulted in unweighted graphs because it captured a binary type of information: while edges linking nodes indicate neighborhood, absence of edges indicates absence of neighborhood ( Fig. S1 A). As explored in this work, there are several ways to define neighborhood. For example, the Heaviside step function f (dist) = H(threshold) can be used to return 1 for nodes that are within a certain threshold distance and 0 otherwise (edge absence). However, edges can take values different from 0 and 1, therefore providing an additional layer of information. For instance, Euclidean graphs ( Fig. S1 B) consist of spatial graphs where edges are weighted with the pairwise distance between nodes. Similarly, one can construct graphs where edges are weighted with a monotonically decreasing function with distance such as the inverse distance (Fig. S1  C). This case can model experiments where higher numbers are associated with closer distances. An example is the count of interactions between two molecules: the closer they are, the more interactions.

Unweighted
Weighted cases are of special interest to the Imaging-by-sequencing field if experiments are designed to quantize proximity beyond a binary manner. For example, DNA Microscopy by amplicon diffusion 5 leverages the fact that there is a relation between expected diffusion encounters and distance between molecules. Therefore, it is possible to weigh edges by counting the number of encounters. In this case, similarly to Fig. S1C, the closer a pair of molecules are, the greater the weight.

1.B Dimensionality reduction
The imaging-by-sequencing computational problem is one in which a high dimensionality dataset, e.g. an adjacency matrix of connections or a matrix of pairwise node-to-node distances, is used as input to recover the original spatial locations of the points that gave rise to proximity-dependent network topology. This progression from a large matrix to a small Cartesian coordinate vector thus falls under the umbrella of dimensionality reduction tasks, where the goal is to take a high-dimensional data set and map it onto a lower-dimensional space while preserving properties of interest in the source data. Dimensionality reduction techniques can be roughly broken down into linear and nonlinear categories:

1.B.1 Linear projection methods
Linear projection techniques, including principal component analysis (PCA), linear discriminant analysis (LDA), and random projections, aim to find a linear combination of input features that maximizes variance in the data. Projection methods assume that data points in a data set are correlated, and that the most important features are ones that explain the variance in the data. If the data points are sampled from a low dimensional manifold and embedded in a high dimensional space, as in the case of a pairwise inverse distance matrix or matrix of shortest path graph distances 1 , then the directions with the greatest variance in the data points may align with those that preserve pairwise distances or the local neighborhood structure of the data.

1.B.2 Nonlinear methods
Nonlinear techniques can capture non-linear relationships in the input matrix and can be more effective at preserving important features of the original data. These include manifold learning techniques like Isomap 8 , locally linear embedding (LLE) 9 , Laplacian Eigenmaps 10 , and node embedding methods. These techniques are based on the idea that high dimensional data lies on or close to a lower-dimensional manifold and use non-linear functions to discover a representation of the data that preserves the geometry of the data by maintaining pairwise distances between points.

Methods
Of particular interest to us are the node embedding methods, including DeepWalk 11 , Large-scale Information Network Embedding (LINE) 12 , PyMDE 13 , and Node2vec 14 . These are methods that learn a low-dimensional representation of the nodes in a graph. By pairing this technique with a manifold learning technique, we are able to set up an autoencoder format whereby a high-dimensional data set, such as the graph adjacency matrix, is first reduced to an intermediate medium-dimensionality encoding using node embedding and subsequently decoded by using manifold learning or other dimensionality reduction to achieve the final low-dimensionality embedding comparable to the original points.

2.A Node2Vec
Node2Vec is a node embedding technique that learns graph structural features by performing biased random walks and embedding nodes into low-medium-dimensional vectors. It effectively compresses the information of the adjacency matrix into a N × D matrix. Storing information in such a way improves the overall computational complexity for the subsequent manifold learning and the robustness of our approach. The node embedding is optimized by maximizing the log-likelihood of observing a node given a particular random walk context. In short, several random walks are performed starting at each node while storing every visited node. The sequence of nodes visited is used to train a skip-gram model. This results in feature vectors that are low dimensional representations of each node. The embedding algorithm optimizes for an objective function that maximizes neighborhood proximity likelihood, and thus feature vectors populated with similar values are likely to reflect nodes that are close in spatial proximity. Negative sampling is used to reduce the optimization's computational complexity.

2.B UMAP
Uniform Manifold Approximation and Projection (UMAP) 15 is a nonlinear manifold learning technique that can be used to reduce the number of dimensions in a data set while preserving important features. UMAP is designed to be fast and efficient and furthermore is focused on preserving both global and local relationships between points by minimizing the distance between the original highdimensionality data and its low-dimensionality embedding. Because of this, it is a more desirable method when one wishes to preserve spatial proximity information accurately. It is worth noting that UMAP can directly recover images without requiring a previous node embedding step, as done in Fig.3 C. A dense version of the adjacency matrix can be obtained by computing the shortest-path distance between every node, where d i j is the shortest-path distance between node i and node j. The result is a N × N distance matrix, assigning N-dimensional vectors to each node. Here, the linear relation between Euclidean distance and graph shortest-path distance is exploited to obtain high-dimensional vectors that preserve Euclidean geometric information.

2.C Hyperparameters
We use the values shown in Table S1 as default hyperparameters when using Node2Vec and UMAP, having found them to be a good compromise between reconstruction accuracy and low computational complexity, as will be discussed below from studying Tables S2-S6.

2.C.1 Grid Search study
We evaluated the accuracy and computational complexity (time, memory) for different hyperparameter values to establish default values (Table S1). Evaluations were done for a system size of N = 5000 and can be found in Tables S2-S6.
As can be seen in Table S2, a low embedded dimension results in poor accuracy. However, after a dimensionality threshold of 16, accuracy values remain similar while computational complexity increases for this system size.
As can be seen in Table S3, an increase in the estimated number of neighbors has, in general, a positive effect on the accuracy scores. However, accuracy can decrease if such a value is too large (from 512 to 2048). We conclude that a value of the order of magnitude ∼ 10 produces accurate enough results while its space complexity is low for this system size.
Both 1/q and 1/p variations in hyperparameter value are not particularly impactful on accuracy (Tables S4, S5), although higher values contribute to a better local quality metric. Furthermore, choosing a value of the order of magnitude below or above 1 contributes to a higher time complexity. We suggest keeping such hyperparameters at value 1. This is equivalent to pure random walks with no bias towards Breadth-First Search or Depth-First Search strategies.
Finally, the random walk length has a positive effect on the local quality metric KNN score but also leads to an increase regarding time and space complexity (Table S6). We consider the trade-off to be adequate for a random walk length of the order of 10, although high system sizes might require greater values. This is further explored in section 2.C.2.

2.C.2 Accuracy impact of the Embedded Dimension and Random Walk Length
We observe that accuracy results are particularly sensitive to the length of the random walks and the choice of embedded dimension (Tables S2, S6), and elaborate further in this subsection. Regarding the embedded dimension, an educated assumption is to consider that the best accuracy is achieved when the embedded dimension is equal to the number of nodes in the network. However, an optimal dimension can be found by comparing different embedding dimensions 16 . Because of the low accuracy resulting from using an embedding dimension of 4 (Table S2), we find empirically that the optimal embedding dimension has to be greater than 2 or 3, the dimension of the original positions. Further empirical results supporting this can be found in Tables S8-S13, where Node2Vec was implemented without a subsequent dimensionality reduction, i.e., directly embedding in 2 or 3 dimensions respectively. The accuracy results are low in comparison to other methods, suggesting that choosing a low embedding dimension does not provide good performance regarding image reconstruction accuracy.
A possible interpretation of why this happens is that the information encoded in the feature vectors is complex enough that it has to be represented through more than 2 or 3 coordinates. This is not unreasonable as the information encoded in the spatial graph is related to the neighborhood of nodes rather than its exact Euclidean distance. This means that the underlying dimension of the graph is likely to be higher than the underlying dimension of the original points. For instance, if 4 points are originally close to each other they will form a K 4 graph (unweighted complete graph with size 4). However, it can be proven that a K 4 graph cannot be perfectly embedded in 2D while exactly preserving graph distances. Therefore, the network structure is in general more accurately depicted through higher embedding dimensions. After the features are encoded in D-dimensional vectors, it is possible to use a manifold learning procedure such as UMAP under the assumption that similar feature vectors correspond to points that are close in the original space.
Moreover, Fig.2 B shows a notable decline in performance with a system size increase. Consequently, properly choosing hyperparameter values with regard to the system size might be a key aspect of preserving accurate results. We investigate the accuracy impact of independently changing the embedding dimension and the random walk length with different system sizes (Fig.S3). In general, increasing hyperparameter values led to an increase in accuracy. For the local KNN metric we set a standard accuracy threshold of 0.8, dividing the grid points that reached this standard from the ones that did not (Fig.S3 A). A similar procedure is done in (Fig.S3 B) with an accuracy standard of 0.98. This reveals that the embedding dimension and the random walk length have to be raised in order to keep the same standards, with a few exceptions for the global metric and the random walk length.
We further study the behavior of changing both the embedding dimension and the random walk length in Fig.S4 by setting the same accuracy thresholds. For the local KNN quality metric, an increase in size requires raising both the value of the embedding dimension and of the random walk length to preserve the same accuracy standards. Conversely, a lower random walk length and a higher embedding dimension are required to preserve the same global CPD quality metric. This suggests, on the one hand, that higher system sizes require raising the embedding dimension to be as accurate as lower system sizes. On the other hand, the random walk length acts as a trade-off between the local and global structure, as raising it can improve the local quality metric but also worsen the global quality metric when the embedding dimension is high enough.

2.C.3 Modeling Accuracy with a non-linear squares fit
In this section, we study STRND's performance to understand better how can we can maintain a given accuracy standard by changing hyperparameters. In particular, we want to find the accuracy function A(x, y, z) that depends on three variables: embedded dimension, random walk length and system size. For the sake of simplicity, they will be referred to as x, y, z respectively. We chose these parameters because they have a considerable impact on accuracy scores (in particular, the KNN local metric) as can be observed from Tables S2, S6 and Fig.2 B. While an increase in embedded dimension and random walk length have a positive impact on accuracy, an increase in system size has been observed to have a negative impact.
An educated guess for the accuracy function is to assume that each variable has a logarithmic impact and are independent from one another: where A is the accuracy score, a, b, c, d are parameters and x, y, z are the embedded dimension, the random walk length and the system size respectively. In order to find the parameters a, b, c, d and obtain the corresponding fit, we use data containing several combinations of x, y, z and their associated accuracy. In particular, the embedded dimension values range from 4 to 32, the random walk length also from 4 to 32 and the system size from 1000 to 18000. The fit we find will be bounded to such values, meaning that predictions below and over them might not be valid. Furthermore, the assumptions that x, y, z are not coupled and that they have a logarithmic impact are not necessarily true. Nonetheless, we find a fit using the non-linear least squares method and test it. As seen in Fig.S5 A, the values predicted by the fit and the ground truth do not have a linear correspondence, possibly to the nature of the assumptions made to find the fit. However, they have a logistic correspondence, both for the KNN and for the CPD quality metric. We decide to perform a second fit, this time a logistic fit with its corresponding parameters: L, k, x 0 . As seen in Fig.S5 B, the predicted and ground truth accuracy scores now have a linear relationship with noise. The accuracy equation becomes: This equation can predict accuracy in an approximate manner for a bounded range of its variables x, y, and z. For example, it is possible to set an accuracy preservation condition, A 1 (x 1 , y 1 , z 1 ) = A 2 (x 2 , y 2 , z 2 ) and obtain: This is particularly useful for the task of choosing the embedded dimension x 2 and the random walk length y 2 that preserves accuracy when changing to a system size z 2 , given a previous accuracy observation with variables x 1 , y 1 , z 2 . The results obtained with this approach are an approximation. However, it can provide educated guesses that may lead to faster results when compared to solving the task of accuracy preservation manually.  (KNN and CPD) for a logarithmic fit. The data can be modeled with a further logistic fit with parameters L, k and x 0 . B. Ground truth and predicted accuracy (KNN and CPD) after applying the logistic fit. The result is a (noisy) linear correspondence between ground truth and prediction values. C. Predicted random walk dependence with size in order to preserve predicted accuracy.

2.C.4 Computational complexity of constant accuracy reconstructions
In the previous section, we built an accuracy predictor that predicts performance given certain reconstruction hyperparameters. In this section, we investigate the computational impact of using the values suggested by such a predictor. In particular, we focus on the problem of finding configurations that preserve accuracy given size variations for the 2D case. We use the local KNN quality metric as an accuracy measure and the random walk length as the variable that can raise computational complexity. The reason we choose random walk length and not embedded dimension is that changes in random walk length will lead to different empirical complexities for low values (ranging from 4 to 32), which is a behavior that is not observed for the embedded dimension (Tables S2, S6).
Fig.S5 C shows the random walk length value that is suggested by the predictor given a certain system size in order to achieve an accuracy standard of 0.7. We find that the predictor suggests that for a size increase Z = z 2 z 1 , the random walk length increase should be Y = Z 0.56 . We reconstructed images with 7 different sizes ranging from 500 to 15000 following this random walk length behavior, which resulted in random walk lengths ranging from 3 to 19. Furthermore, we performed 50 reconstructions for each size in order to increase the statistical validity of the results. We annotated the KNN accuracy, time and space complexity and computed their average value for every size. Fig.S6 A shows that the empirical time complexity that preserves accuracy is of order O(N 1.08 ), which is greater than the complexity we found for a constant random walk length in Fig.2 C (O(N 1.02 )). This supports the intuition that in order to preserve accuracy the time complexity must increase. Here we present a quantitative measure for such an increase. However, Fig.S6 shows that an increased random walk length presents a space complexity is O (N 1.04 ), while a constant random walk length presented a complexity of O(N 1.05 ). Therefore, peak memory values are affected by the system size rather than random walk length, at least for values ranging from 3 to 19 and sizes from 500 to 15000. Fig.S6 C shows the mean accuracy values for each size. The reconstructed images had mean accuracy values ranging from 0.68 to 0.81, which are close to the accuracy standard 0.7 used to choose the random walk length hyperparameter.

2.D.1 The local and global structure
Local structure is preserved when points that are close to each other in the original image are also close in the reconstructed image. In this work, local structure is investigated through the KNN quality metric, which evaluates the fraction of neighbors preserved in the reconstruction. Conversely, global structure is preserved when points that are separated far apart in the original image are equally far in the reconstructed image. In this case, the quantity to be preserved is separation distance, and it can be investigated by computing the correlation between pairwise distances.

2.D.2 KNN
The KNN quality metric can be formally defined as (S1) Here, N is the total number of elements and k is the number of considered neighbors. We used k = 15 for quality metrics presented in this manuscript. Effectively, KNN computes the fraction between original and reconstructed neighbors (inner summation) for every node (outer summation) and averages it over the total number of elements.

2.D.3 CPD
The CPD quality metric is the Pearson correlation between original and reconstructed pairwise distances. As opposed to KNN, CPD measures concern not only neighboring points but also distant points. The motivation to measure correlation is to evaluate distance preservation. If pairwise distances in the original image are linearly correlated with pairwise distances in the reconstructed image, then the reconstructed image is equivalent to the original image up to an affine transformation. Dealing with scalability often require workarounds: for N points, there are N(N−1) 2 ∝ N 2 pairwise distances, which can be computationally prohibitive for large N. The adopted solution is to randomly sample N = 1000 points and compute their 499500 pairwise distances for the original and reconstructed points, to later obtain a representative Pearson correlation.

2.D.4 Distortion
Point cloud alignment is necessary to compute distortion, which is defined as the distance between every original and reconstructed point pair. Reconstructed points need to be aligned with the original points because scale, rotation, translation, or chirality are not necessarily preserved. Consequently, point cloud alignment consists in finding the affine transformation such that distortion between original and reconstructed points is minimized. To find such a transformation, we apply the coherent point drift algorithm 17 .
Nonetheless, coherent point drift is usually not enough to successfully align reconstructed points because for two reasons. First, it does not consider reflection transformations which account for chirality corrections. Second, if the original and reconstructed point clouds are initially far apart, the algorithm might not converge in the optimal transformation. To solve this, we first perform a set of transformations over the reconstructed points involving scaling obtained by computing the correlation between pairwise distances, point-to-point translation and rotation, and the set of coordinate reflections that minimize distortion.

2.D.5 Quality metric variation with proximity graph
To examine the robustness of STRND we studied the quantitative quality metric variation to proximity graph change. The results are presented in Table S7.

3.A Other approaches
As there are multiple possible image recovery routes, different dimensionality reduction techniques, and combination pipelines, we compared a selection of approaches with STRND in terms of quality metrics, computational complexity and robustness to network patterns.

3.A.1 Shortest-path distance matrix
Shortest-path distance matrices are obtained by applying a Breadth-First-Search algorithm to a sparse adjacency matrix. The result is a densely-populated N × N matrix where the element N i j contains the shortest-path distance between node i and node j. Such distance matrix can be reduced to a N × 2 or N × 3 matrix using e.g. PCA or UMAP (parameters in Table S1), as shown in section 3.B.

3.A.2 Spring relaxation
Force-directed approaches are used as a means to layout graphs in aesthetically pleasing ways. They are inspired by physical systems: nodes in a graph are considered particles that attract each other if they are close and repel otherwise. In particular, we implement the Fruchterman-Reingold algorithm 18 , where edges between nodes can be considered physical springs storing elastic potential energy. The goal is to find the set of node positions that reaches equilibrium, i.e., to find the global minimum of the energy function. In particular, we use the Networkx 19 implementation with 300 iterations to ensure convergence.

3.A.3 Landmark Isomap
Landmark Isomap was implemented by randomly sampling D landmark nodes, computing the squared shortest-path distance between every landmark and every node in the graph, and finally triangulating nodes positions through Multidimensional Scaling (MDS) 20 . The squared shortest-path computation results in a N × D matrix, as there are N nodes and D landmarks. Such computation is done by iterating the single-source Breadth-First-Search algorithm for every landmark. This approach is considerably better than computing the all-pairs shortest-path distance matrix in terms of empirical computational complexity, as it avoids computing a N × N matrix. The next step is to store the D × D squared distance matrix between all landmarks, as it is used to map the landmark nodes into the reconstructed space by applying MDS. Finally, the rest of the nodes are triangulated by solving an eigenvalue problem on the Moore-Penrose pseudo-inverse of the D × D landmark matrix 3 . As the default embedding dimension for STRND is 32, we also use D = 32 landmarks. Nonetheless, we study the accuracy impact of the number of landmarks in the section 3.B and in particular in Fig.S8.

3.A.4 DNA Microscopy by Amplicon Diffusion
Amplicon diffusion imaging is an image recovery approach tailored for a particular experimental setting that uses a weighted scheme 5 . It relies on molecular diffusion to gather proximity information. Proximity information can be extracted because there is a physical link between the expected diffusion encounters between molecules and their distance. It distinguishes from unweighted approaches which only gather binary proximity information regarding neighborhood, i.e., establishing if molecules are neighbors or not.
As a benchmark for the quality metrics used in this work, amplicon diffusion was simulated and reconstructed according to the recovery method from Weinstein and colleagues' code. The simulation was done for a system size of N = 10000 randomly positioned beacon and target molecules and 100 Unique Event Identifiers (UEI) per molecule 5 . A large enough number of UEI per molecule was chosen for a more accurate reconstruction. The total number of molecules is not the number of target (imaged) molecules, which is N ≈ 5000 in this simulation. The original and recovered image through the amplicon diffusion algorithm (which produces its own weighted proximity graph) can be seen in Fig.S7 A and B, along with the respective quality metrics. The global CPD quality metric is of similar magnitude to STRND (Tables S8, S9, S10). Close inspection reveals the presence of outlier points that influence the KNN local quality metric. These outliers can also be seen in Fig.S7 C, the distortion resulting from aligning original and recovered points, and

3.B Benchmark comparison
The robustness of each approach is first studied by generating unbiased network pattern rules in Fig.3 D. Such rules were obtained from Gaussian mixture distributions with randomly selected parameters (10 individual distributions each). For every reconstruction method, 100 images were recovered with system size N = 100. To compute reconstruction accuracy we use the local KNN quality metric and the global CPD quality metric, abstaining from using point cloud alignment (distortion) as it is not guaranteed to find the optimal affine transformation when reconstruction accuracy is low. Although a low system size enables the reconstruction of many images in a limited time, it may be subject to boundary effects. Therefore, we further study a system size of N = 5000 in Tables S8-S13. However, we restrict the network patterns to three representative scenarios for the 2D and 3D case: the KNN graph, the Delaunay graph and a graph with a (stochastically) decaying probability with distance. To compute the empirical computational complexity we evaluate the total time and peak memory used by the image recovery algorithm. Corresponding reconstructions from the benchmark are shown in Fig.S9.
Tables S8, S9, S10 contain the accuracy and computational complexity for three different proximity graphs for the 2D case, whereas Tables S11, S12, S13 are the 3D case equivalent.
Shortest Paths + PCA obtains satisfactory results for both the local and global quality metric in specific cases, e.g. for the KNN proximity graph (Tables S8, S11). PCA is reported to perform better in terms of global structure preservation than local 21 as was observed in most proximity graphs we tested. Second, Shortest Paths + UMAP slightly outperforms PCA in every case both for the local and global structure, except for the 3D Delaunay graph (Table S12) where the global quality metric is marginally lower.
Spring Relaxation stands out for its high empirical running time. Its worst-case computational complexity is quadratic with the number of nodes, O(N 2 + E), where N is the number of nodes and E is the number of edges. Interestingly, spring relaxation performs better for the 3D case (Tables S11, S12, S13) than for the 2D case (Tables S8, S9, S10). Similarly to the Shortest Paths + PCA approach, Landmark Isomap displays a high level of accuracy for the KNN proximity graph using D = 32 landmarks. Additionally, its running time and peak memory consumption are the lowest of all approaches, which makes it a suitable candidate for fast and scalable reconstructions. Nonetheless, it is not as robust as STRND regarding performance in different proximity graphs. For a stochastic decay graph (Tables S10, S13) its local quality metric is slightly worse and for the Delaunay graph (Tables S9, S12) both local and global quality metric are substantially lower.
We further investigate the effect of the number of landmarks on the accuracy scores in order to see if a poor parameter choice led to non-representative results. Although the purpose of landmark nodes is to decrease overall computational complexity, increasing the number of landmarks can have an impact on how well the images can be reconstructed. Indeed, we find that a greater number of landmarks positively affect the accuracy scores (Fig.S8). This effect is clear for Delaunay proximity graphs and in particular for the global CPD quality metric, where an increase in landmarks benefits the accuracy score. However, such scores are still considerably low when compared to other approaches. Conversely, high accuracy scores are shown for the KNN and Stochastic Decay graphs, but the impact of a higher number of landmarks seems more subtle. Table S14 shows the effect of the number of landmarks on empirical computational complexity. Although increasing landmark nodes might have a slight positive effect on accuracy scores, at the same time it can substantially increase computational complexity.
All approaches involving shortest-path computations show a high space complexity as a N × N shortest-path distance matrix has to be stored and processed by a manifold learning algorithm. The system size quadratic dependence becomes problematic as N grows, leading to memory allocation issues and increasing running times. An exception to this behavior is Landmark Isomap, as it only requires a markedly less computationally demanding N × D matrix, which can be computed through single-source shortest path algorithms. Possible workarounds to this problem might involve distance matrix pruning and sampling, e.g., establishing a maximum shortest-path threshold or randomly sampling matrix values.
STRND shows the best robustness to different proximity graph types, as it has a similar accuracy and complexity score for different proximity graphs both for 2 and 3 dimensions (Tables S8-S13). This behavior is not as prevalent for other approaches as they show lower robustness concerning accuracy scores and computational complexity.

Code Availability
The full code including STRND and other approaches can be found at https://github.com/DavidFernandezBonet/ImageRecovery, including a demo Jupyternotebook that demonstrates use cases and acts as a tutorial for new users. It is also possible to install it as a Python package using pip install ImageRecovery.
Fig. S10 Distortion plots for every system sized displayed in Fig.3 B for 2 and 3 dimensions.