Unraveling The Past: Applying GNN Techniques to Ancient Scroll Reconstruction

By Yasmina Abukhadra, Frank D’Agostino, and Aakash Mishra

Background

After the volcanic eruption of Vesuvius in 79 AD, hundreds of scrolls of antiquity text were buried and preserved in Pompeii and Herculaneum. Unable to be unrolled without destroying the contents, non-invasive approaches have been developed to read what is inside. This has involved CT scanning the scrolls to generate a dataset composed of 2D CT slices [2]. The effort to decode the scrolls from this data is led by the Vesuvius Challenge, which is employing machine learning methods to detect the ink in the CT scans and virtually unwrap the scrolls.

Motivation

Current methods and models of virtually unwrapping an entire scroll would cost 1 to 5 million dollars given the manual verification still needed. With 300+ scrolls yet to be read, step function improvements in model iteration time and accuracy are needed, with an emphasis on unwrapping at scale [1]. Reconstructing the remaining scrolls would unveil lost works of ancient literature, benefiting both historians and classicists.

A portion of the virtual unwrapping process involves using a graph representation to stitch together scroll patches segmented from the CT data. Current approaches use random walks or optimization-based approaches to solve this problem, but are slow and difficult to work with. We aimed to apply GNN techniques to this graph-based patch-stitching problem, allowing for a more efficient, streamlined, and generalizable approach. Additionally, inductive GNN models could more easily be applied to quickly stitch scroll-patch graphs from unseen scrolls.

Figure 1: To the left is PHerc. 332 (Scroll #3) from the Biblioteca Nazionale di Napoli, which we analyze and use for training/testing in this project. To the right is a slice of the CT scan of PHerc. 332 [1].

Data

The EduceLab-Scroll dataset contains volumetric X-ray micro-CT scans with resolutions of 2.24μm to 27.17μm of several scrolls and scroll fragments. Data used in the preparation of this project were obtained from the EduceLab-Scrolls dataset [2]. The dataset is made publicly available (conditional to terms agreement) by the Vesuvius Project: https://scrollprize.org/data.

To convert this data to a graph-representation for the patch-stitching problem, a detailed process must be followed (as will be outlined below). For this reason, we additionally provide a precomputed graph-representation of a portion of one of the scrolls, scroll 3, and pipeline for converting it for a format compatible with torch_geometric. In our pipeline, we convert from the existing ScrollGraph object to a NetworkX graph and from NetworkX graph to PyG Data object. We ensure that we include relevant node features and labels from the ScrollGraph object, and split the nodes of the dataset according to training and test proportions.

def convert_scrollgraph_to_networkx(scroll_graph):
    graph = nx.DiGraph(name="scroll graph")
    nodes = scroll_graph.nodes
    nodes_to_add = []

    print("length of node items", len(nodes.items()))
    for key, val in tqdm(nodes.items()):
        # include a subset of the features
        centroid = val['centroid']
        assigned_k = val['assigned_k']
        winding_angle = val['winding_angle']
        nodes_to_add.append((key, {'winding_angle': winding_angle, 'centroid': centroid, 'assigned_k': assigned_k}))
    graph.add_nodes_from(nodes_to_add)

    edges = scroll_graph.edges
    # convert the edges into connections between nodes
    edge_list = []
    for edge in tqdm(scroll_graph.edges):
      #ignore the k offset being the key
      for number, edge_info in scroll_graph.edges[edge].items():
        # only add edges between non-deleted nodes
        if (not nodes.get(edge[0])) or (not nodes.get(edge[1])):
          continue
        # add edge but do not append edge features
        edge_list.append((edge[0], edge[1]))
    graph.add_edges_from(edge_list)

    return graph

As will be described in the Graph Problem section, we have an initial graph, and a solved graph. The initial graph contains noisy nodes that are removed in the solved graph, and the solved graph contains ground truth labels for the nodes. In our pipeline, we combine the information from the initial and solved graphs into Data objects for training and testing.

Preprocessing

To get to the graph, there are a variety of steps from going from the CT scan images to the point cloud data which we will briefly describe.

First, we hosted all our work on a AWS EC2 instance to store the scroll data and models. Using the Vesuvius Data Download repo, we used rclone to grab all the 2D X-ray CT scan data that together represent the volumes depicting the scroll. With the data, we then want to generate our graph object, where we follow the steps and scripts outlined in the ThaumatoAnakalyptor repository. After starting a docker container, we first downsample the resolutions of the images to 8μm to make segmentation easier. With the downsampled grids, we then detect high gradient regions to detect surfaces, and output 3D point cloud data on the detected papyrus surface. Finally, with the point cloud data, we group patches of point cloud data together using the Mask3D segmentation model [3].

Figure 2: PHerc_332 (scroll 3) graph being traversed by the ThaumatoAnakalyptor solver to determine the winding numbers for each node.

In previous iterations (Chesler [4], Schilliger [5, 7, 8], Mou and Ahmed [6]), they implemented a patch stitching algorithm which takes segmented surface patches and treats them as nodes, with edges based on a similarity score based on the degree of overlap of the patches with other patches. Past approaches utilized a random walk over the subgraph in order to filter out noisy nodes and assign a winding number (scroll layer) to each node. More recent efforts use a MAP estimation solver over the graph to assign winding angles and numbers that minimize an objective related to the physical constraints of the graph [5, 7, 8]. Once the patches were stitched together using the assigned winding number, the pipeline would use Poisson surface reconstruction to create a non-manifold mesh for each half-winding. The resulting mesh is then flattened using Symmetric L1 energy Minimization (SLIM) for texture parameterization and then rendered using TIFF surface reconstruction [8].

Graph Problem

For each segmented patch, we have the 3D coordinate of its centroid (the “balancing point” for all of the points in its point cloud representation). From its location relative to the “umbilicus” or center of the scroll, we also have an angle in the range of -180 degrees to 180 degrees representing its angular position relative to the center. We will call this the initial winding angle.

Each patch is considered a node, and is connected by edges to patches that are considered close/similar to it, as determined by scores computed in pre-processing.

Figure 3: Animation we made showing a synthetic 3D scroll with 3 different sheets spiraled up in the scroll. We can see that our aim is to garner the winding number so that we can project the sheets into a 2D side-by-side view. Link to notebook.

The goal of the graph problem is to both filter out noisy patches and determine an updated winding angle, which indicates which layer of the scroll the patch is in. The initial winding angle mentioned above provides no information about what layer of the scroll the patch is in. The updated angles can smaller than -180, and are roughly a — 360*k, where k represents the layer of the patch and a represents the initial winding angle. [5][7][8]

We split this layer portion of the problem into two approaches:

1. Determining the layer of the patch: This is a node classification problem with the number of classes equal to the number of layers of the scroll (or slices taken by the CT scan).
2. Determining the updated winding angle directly: This is a node regression problem, with the goal of predicting an angular value similar to the ground truth “updated winding angle”.

There is also the noisy patch problem, where many patches may be artifact point cloud data that is not physically consistent with the rest of the patches. The existing algorithm will delete noisy nodes based on inconsistencies with their neighbors. This is simply a binary node classification problem. We classify a patch as 0 if it is a true patch and 1 if it is a noisy patch.

Overall, we predict 3 tasks with a variety of GNNs: layer multi-class classification, winding angle regression, and noise classification.

Existing Approach

Past approaches used random walks to assign winding numbers and angles, however the most recent approach uses a solver that optimizes objectives based on physical constraints. It uses a spring-based energy model over many iterations, progressively reducing the spring constants to converge in a low-energy state. Basically, it treats every edge connecting patches in the scroll as a spring, and tries a variety of winding numbers for each node such that it reduces the energy in the system where it is stable. Conceptually, nodes that have a high confidence edge between them should not differ substantially in winding number, meaning the solver finds a decent heuristic for the winding number assignments based on these physical constraints. To assign winding numbers, they use a ring solver that enforces monotonic growth of winding numbers along subsequent neighbors of the starting seed node. During this stage, it also filters edges and deletes nodes that are noisy and disrupt the system beyond some threshold by detecting the largest connected components. [8, 9]

Limitations of this approach are its high barrier-to-entry C++ implementation, long computation time, and heuristic based learning. The implementation is also incomplete in that it often encounters memory management issues. We were able to modify the implementation to get a few cycles of this approach so that we could get ground truth winding angle, winding number, and noisy node data (our changes can be found in our Github). We aimed to train GNNs using this data to create models that are more generalizable and efficient than the current graph solver.

Applying GNNs to This Problem

Feature Augmentation

The graph begins with the centroid (x,y,z) coordinates and starting winding angle (in [-180, 180]) as the only node attributes. We also have the coordinates of the umbilicus at each level, which is the theoretical line going through the center of the scroll. For feature engineering, we computed node features that better capture the structure of the point cloud graph and better represent each node’s geometry with respect to the umbilicus.

Node features.

We also computed a variety of edge features to augment the graph with structural information.

Edge features.

We know that the eigenvectors associated with larger eigenvalues are associated with more local structure. In our ablation studies, we compare the performance of models with just the original node features, versus with combinations of our additional edge and node features.

Figure 4: Higher eigenvectors depict more local structures of the graph. We can see a heatmap of the first 4 eigenvectors of the normalized Laplacian matrix D^(-1/2)LD^(-1/2) and its entries of the sampled graph. We use eigenvectors 1, 2, and 3 for the node features.

Models

We used PyG’s built-in GNN models to quickly compare the efficacy of different model types for these tasks. We compared Graph Convolutional Network (GCN) [13], GraphSAGE [14], Graph Isomorphism Network (GIN) [15], and Graph Attention Network (GAT) [16] models.

As GIN, GAT, and GraphSage are inductive models, and therefore generalize well to unseen graphs, we hoped they would provide good solutions for future applications to unseen scrolls. Because each model architecture uses slightly different message passing, aggregation, and update methods, we aimed to use experimentation to determine which was best suited for each task.

In addition to these built-in models, we implemented a DCGNN (Wang et al., 2019) [10], a GNN specialized for learning on point cloud data. Since our graph is derived from point cloud data, we suspect this model could outperform other more general GNN models. The main innovation of DGCNN is the EdgeConv layer. The paper applies convolution-like operations on edges connecting pairs of points. They also allow the graph to be dynamic where the k-nearest neighbors of a node changes layer to layer. This operation is permutation invariant and partially translation invariant, “balancing local information of the patches while keeping global shape information” [10]. Formally, they define an asymmetric edge function with shared MLP weights:

Which is mathematically equivalent to learning a linear layer over the concatenation as shown in the PyG documentation [9]:

We can see our code implementation of the DGCNN using the PyG supported PyG EdgeConv layer. We employ 5 layers with ReLU activation and dropout regularization, with cross entropy loss for the multi-classification task (predicting which layer each node is in).

class DGCNN(nn.Module):
  def __init__(self, in_channels, hidden_channels, out_channels, num_layers=5):
    """
    DGCNN with EdgeConv layers

    Args:
      in_channels (int)
      hidden_channels (int)
      out_channels (int)
      num_layers (int)
    """

    super().__init__()

    self.convs = nn.ModuleList()

    for i in range(num_layers):
      in_c = in_channels if i == 0 else hidden_channels
      self.convs.append(
          EdgeConv(
              nn=nn.Sequential(
                  nn.Linear(2 * in_c, hidden_channels), 
                  nn.ReLU(), 
                  nn.Linear(hidden_channels, hidden_channels)
              )
          )
      )

    self.mlp = nn.Sequential(
        nn.Linear(hidden_channels, hidden_channels), 
        nn.ReLU(), 
        nn.Linear(hidden_channels, out_channels)
    )

  def forward(self, x, edge_index, batch=None, return_embed=False):
    embeddings = []

    for conv in self.convs:
      x = conv(x, edge_index)
      embeddings.append(x)

    if batch is not None:
      x = global_max_pool(x, batch)

    out = self.mlp(x)
    if return_embed:
        return out, embeddings
    else:
        return out

Smoothing Loss Term

The regression task of predicting the winding angle was proven to be difficult for our models and subject to outlier predictions because in our dataset, it can range up to ~2500 degrees (since every spiral is 360 degrees and there are 7 layers in the scroll we worked on). Inspired by the solver’s constraint of monotonicity, we worked on adding a potential loss term to our models capturing this physical constraint. We call it the smoothness loss, as it penalizes the model for giving very different winding angles than its neighbors, enforcing a smooth change in winding angles as you follow the patches’ neighbors through the scroll. We define it as:

Formulation for the smooth loss and how we combine it with the MSE loss.

We can see that we penalize the square of the difference of the predicted winding angles with the pairs of nodes and add it to the MSE loss, weighted by some term, which empirically we set to 0.30. Below is how we incorporated it in Python:

def smoothness_loss(pred, edge_index):
  i, j = edge_index
  diff = pred[i] - pred[j]
  return (diff**2).mean()

...

if use_smoothness_loss:
      loss += lambda_smooth * smoothness_loss(out, data.edge_index)

Jumping Knowledge LSTM

We compared GraphSAGE with and without enabling the built-in PyG Jumping Knowledge mode (using an LSTM). Jumping knowledge enables intermediate representations to inform the final node embedding. The LSTM-enabled version uses attention to determine which neighborhood ranges are the most relevant for the final representation. We ablated against this modification, as Jumping Knowledge has previously been shown to provide performance improvements for a variety of GNN architectures. [17]

Performance

Layer Classification

These results for the multi-class layer classification task are with no edge/node features and have the same number of layers and hyperparams for each model type.

• GraphSAGE Accuracy: 0.84
• GAT Accuracy: 0.71
• GIN Accuracy: 0.82
• GCN Accuracy: 0.83
• DGCNN Accuracy: 0.81

Figure 5: We can see the training accuracy over time for the multi-class layer classification task. We can see that GraphSAGE performs best.

Winding Number Regression

These results are with the smoothness loss but no edge/node features. We can see that GraphSAGE and DGCNN perform best for the regression task.

• GraphSAGE Mean absolute error: 135.62
• GAT Mean absolute error: 192.19
• GIN Mean absolute error: 408.33
• GCN Mean absolute error: 257.20
• DGCNN Mean absolute error: 93.58

Figure 6: Training loss over time compared across the models. We can see that DGCNN and GraphSAGE perform the best for the winding angle regression task.

Noisy Patch Classification

For binary classification of noisy patches we see that GIN outperforms all other models on this task.

• GraphSAGE Accuracy: 0.69
• GAT Accuracy: 0.66
• GIN Accuracy: 0.75
• GCN Accuracy: 0.51
• DGCNN Accuracy: 0.64

Figure 7: ROC AUC plot of our binary classification GIN model to determine which nodes are noise and should be deleted or not.

Ablation Analysis for Regression Task

We wanted to understand how different components of the model development contributed to the regression task. We incorporated node features (including eigenvector features), edge features, a smoothing loss term, and multiple model types. To conduct the ablation analysis, we created an ablation config to easily turn each of these components on or off and could call them as follows:

@dataclass
class AblationConfig:
    use_node_features: bool = True
    use_edge_features: bool = True
    use_smoothness_loss: bool = False
    lambda_smooth: float = 0.3
    model_name: str = "DGCNN"

example_ablation_settings = [
    AblationConfig(model_name="GraphSAGE",
                   use_node_features=False,
                   use_edge_features=False,
                   use_smoothness_loss=False),

    AblationConfig(model_name="DGCNN",
                   use_node_features=True,
                   use_edge_features=False,
                   use_smoothness_loss=False)
]

Figure 8: Training loss of a variety of models in our ablation analysis. We can see that GraphSAGE with the LSTM jumping knowledge seemed to perform the best, and whether it had edge features or the smoothing loss did not impact the model significantly.

We can see that node features are critical for the model to perform well on the regression task, which makes sense since the eigenvector features and other geometric computations provide crucial structural information for the model to better predict the angles. We also see that the smoothing loss and LSTM jumping knowledge improve model performance significantly.

Latency

We can see that the training and inference time of the DGCNN model is almost 10x faster than the physics-based solver. Both models were GPU-enabled for the purposes of comparison.

Table comparing the computation time to run the Thaumato graph solver for Scroll 3 compared to training a DGCNN model on Scroll 3 pointcloud data. We can see ~10x reduction in latency when both runs were GPU-enabled.

Embedding Visualization

We employed TSNE [12] to represent and analyze a subsample of the high dimensional node embeddings after a single EdgeConv layer in the DGCNN model. The class labels represent each layer. We can see the visualization below and the ability of the EdgeConv layer to represent the nodes in a meaningful way for the task at hand (classification or regression).

Figure 9: TSNE visualization of the node embeddings after an EdgeConv layer in the DGCNN model. We can see that the sampled points mainly include layers 3, 4, and 5, and we can see an interesting separation/pattern across the different nodes.

Takeaways

GNNs provide a promising alternative to solver-based approaches to the patch stitching problem in ancient scroll reconstruction.

For determining which layer a sheet is in, all GNNs performed similarly and were able to assign the correct layer often.

For regressing on the winding angle, the models had more trouble and our model developments helped improve it substantially. For model types, GraphSAGE and our implemented DGCNN model performed the best. From the ablation analysis, we found that the node features were crucial in allowing the model to make reasonable predictions. From there, our smoothing loss function added further physical constraints to improve the quality of predictions, and the LSTM jumping knowledge framework improved the model performance even more.

Not only were our models able to achieve high performance on these tasks, but it is substantially faster than the solver (~10x faster), making it a promising candidate for scaling. Right now, single digit scrolls have been unfolded, and there are still 300+ scrolls yet to be read. Using GNNs to generalize and bring the quality insights experts and the existing solvers have to many new scrolls will make this feat more attainable.

For detecting noisy sheets, our model was also able to get an AUC of 0.88, showing the GNNs were able to use the point cloud relationships to be able to filter out noisy nodes for other downstream tasks.

We chose to use supervised learning methods by training on ground truth results derived from a computationally intensive physics-based solver. For future directions, we also hypothesize that unsupervised or self-supervised approaches could be beneficial for the winding number and noise classification tasks. Architectures such as the adversarially regularized graph autoencoder for graph embedding (ARGVA) use graph-based variational autoencoders to generate node embeddings that can be used for k-means clustering [11]. This self-supervised learning technique could leverage structural properties of the graph to cluster nodes based on winding angle or noisiness. Because all data we can obtain about the scrolls and ground truth values are derived computationally, self-supervised methods that use inherent structural properties could be a more direct way to train these predictive models.

In sum, patch-stitching for scroll reconstruction is a challenging and interesting task, for which GNNs are a promising approach. Accurate inductive GNN models could speed up the process of scroll reconstruction for unseen scrolls, and increase the feasibility of automation. We hope that future work can build off these approaches, helping to reconstruct the remaining 300+ scrolls and reveal lost history.

Code

Our Colab with all models and experiments is publicly available here: https://colab.research.google.com/drive/1gtvPtIeVhciZiyPaHtLf_7lz-b_lAQcf#scrollTo=WeikYCah1mH0

To run, you will need to use the data (precomputed graph files) from this private folder: https://drive.google.com/drive/folders/1F3kXssGQz_D0mvZUpJYVBnlOrFOyxTJC

For access to the repositories we used for processing the data, please contact yasabukh@stanford.edu, frankdag@stanford.edu, or aamishra@stanford.edu.

They are available in private repositories here:

• Data download repository: Frankdag20/vesuvius-224w-data-download
• Modified preprocessing repository: https://github.com/yasabukh/scroll-project
• Colab with models and figures: https://colab.research.google.com/drive/1gtvPtIeVhciZiyPaHtLf_7lz-b_lAQcf?usp=sharing

References

[1] Vesuvius Challenge. https://scrollprize.org

[2] Parsons, S., Parker, C. S., Chapman, C., Hayashida, M., & Seales, W. B. (2023). EduceLab-Scrolls: Verifiable Recovery of Text from Herculaneum Papyri using X-ray CT. ArXiv [Cs.CV]. https://doi.org/10.48550/arXiv.2304.02084

[3] Schult, J., Engelmann, F., Hermans, A., Litany, O., Tang, S., & Leibe, B. (2022). Mask3d: Mask transformer for 3d semantic instance segmentation. arXiv preprint arXiv:2210.03105.

[4] A. L. A. T. Ryan Chesler, Ted Kyi. Solution for kaggle vesuvius ink detection challenge, 2023. URL https://github.com/ainatersol/Vesuvius-InkDetection.

[5] J. Schilliger. Thaumato anakalyptor, 2024. URL https://github.com/schillij95/ThaumatoAnakalyptor.

[6] F. S. Mou and T. Ahmed. Ink detection from carbonized herculaneum papyri using deep learning. In 2023 26th International Conference on Computer and Information Technology (ICCIT), pages 1–6. IEEE, 2023.

[7] J. Schilliger. Sheet Stitching Problem Definition, 2024. https://github.com/schillij95/ThaumatoAnakalyptor/blob/main/documentation/Sheet_Stitching_Problem_Definition.pdf

[8] J. Schilliger. Thaumato Anakalyptor Technical Report and Roadmap, 2024. https://github.com/schillij95/ThaumatoAnakalyptor/blob/main/documentation/ThaumatoAnakalyptor___Technical_Report_and_Roadmap.pdf

[9] PyG Team. (2024). EdgeConv layer. PyTorch Geometric Documentation. https://pytorch-geometric.readthedocs.io/en/latest/generated/torch_geometric.nn.conv.EdgeConv.html

[10] Wang, Yue, et al. “Dynamic graph cnn for learning on point clouds.” ACM Transactions on Graphics (tog) 38.5 (2019): 1–12.

[11] Pan, Shirui, et al. “Adversarially regularized graph autoencoder for graph embedding.” arXiv preprint arXiv:1802.04407 (2018).

[12] Maaten, L. V. D., & Hinton, G. (2008). Visualizing data using t-SNE. Journal of machine learning research, 9(Nov), 2579–2605.

[13] Kipf, T. N. “Semi-supervised classification with graph convolutional networks.” arXiv preprint arXiv:1609.02907 (2016).

[14] Hamilton, Will, Zhitao Ying, and Jure Leskovec. “Inductive representation learning on large graphs.” Advances in neural information processing systems 30 (2017).

[15] Veličković, Petar, et al. “Graph attention networks.” arXiv preprint arXiv:1710.10903 (2017).

[16] Xu, Keyulu, et al. “How powerful are graph neural networks?.” arXiv preprint arXiv:1810.00826 (2018).

[17] Xu, Keyulu, et al. “Representation learning on graphs with jumping knowledge networks.” International conference on machine learning. pmlr, 2018.