extfieldsofstudy:
Engineering
Computer Science
Physics
Medicine

text:
Deformable image registration with deep network priors: a study on longitudinal PET images

Longitudinal image registration is challenging and has not yet benefited from major performance improvements thanks to deep-learning. Inspired by Deep Image Prior, this paper introduces a different use of deep architectures as regularizers to tackle the image registration question. We propose a subject-specific deformable registration method called MIRRBA, relying on a deep pyramidal architecture to be the prior parametric model constraining the deformation field. Diverging from the supervised learning paradigm, MIRRBA does not require a learning database, but only the pair of images to be registered to optimize the network's parameters and provide a deformation field. We demonstrate the regularizing power of deep architectures and present new elements to understand the role of the architecture in deep learning methods for registration. Hence, to study the impact of the network parameters, we ran our method with different architectural configurations on a private dataset of 110 metastatic breast cancer full-body PET images with manual segmentations of the brain, bladder and metastatic lesions. We compared it against conventional iterative registration approaches and supervised deep learning-based models. Global and local registration accuracies were evaluated using the detection rate and the Dice score respectively, while registration realism was evaluated using the Jacobian's determinant. Moreover, we computed the ability of the different methods to shrink vanishing lesions with the disappearing rate. MIRRBA significantly improves the organ and lesion Dice scores of supervised models. Regarding the disappearing rate, MIRRBA more than doubles the best performing conventional approach SyNCC score. Our work therefore proposes an alternative way to bridge the performance gap between conventional and deep learning-based methods and demonstrates the regularizing power of deep architectures.

Introduction
Medical image registration is the precise overlaying of a fixed and moving image. It is notably used to create patient models (inter-patient monomodal registration), exploit information from different modalities acquired for a single patient (intra-patient multimodal registration), and to monitor tumor evolution (longitudinal registration).
However, the accuracy of image registration is a long-standing fundamental problem in medical image analysis ( [32,16]). In the past 30 years, registration methods have considered registration as an optimization problem: the goal is to minimize a dissimilarity term between moving and fixed images. It is usually measured with mean square error (MSE), mutual information (MI) or normalized cross-correlation (NCC). In addition, a term enforcing smooth and realistic transformations is often used as a regularizing term.
Most conventional methods manage registration by parameterizing the transformation between the fixed and moving image. They use for instance discrete cosine transforms ( [15]), 3D Fourier series ( [9]), cubic B-splines ( [26]) or optimized velocity fields ( [1]). The parameters of the model are typically optimized though efficient second-order minimization ( [44]) or stochastic gradient descent (SGD) ( [25]), within a single-or multi-level optimization strategy. The minimized cost function combines dissimilarity and regularization terms.
Recent deep learning (DL) supervised registration methods ( [16]) also rely on dissimilarity measurements and regularization terms that can be optimized with variants of SGD. However, there are two key differences compared to conventional approaches. DL methods require a training stage with a training database that gives the model prior knowledge regarding deformations. They also use a different type of transformation model, with convolutional neural networks (CNN) used as an over-parameterized but structured model of the deformation field.
Despite recent developments in DL-based registration ( [8]), conventional methods still perform better and obtain more accurate results in most applications ([21]). In addition, even though trained networks are faster than conventional methods, training requires large databases which often restricts their use to a specific therapeutic area.
With this work, we propose a new registration method formalized as a conventional registration approach, with a deformation field modeled by an untrained deep pyramidal network. We named our approach MIRRBA for Medical Image Registration method Regularized By Architecture.
We applied MIRRBA on longitudinal PET images acquired for the evaluation of treatment response in patients with metastatic breast cancer and compared the results to registrations obtained with standard approaches. We also studied the impact of different architecture configurations (14 configurations) on the deformation field.
The contributions of this paper are i) the proposition of a novel registration method regularized by architecture MIRRBA, ii) an extensive comparative study of the effects of different network components on the deformation field, and iii) a solution to register whole body PET images both at a global and local level without the need for prior registration to facilitate the simultaneous monitoring of multiple lesions.

Related work
Automatic longitudinal lesion monitoring Automatic approaches for lesion segmentation and/or image registration have been developed over the years to help with cancer monitoring and treatment response assessments. [36] used affine registration to align baseline and follow-up PET images, before subtracting them to identify tumor voxels showing significant changes between the two scans. However, since affine registration is a global method, it mainly performs well for single localized tumors. In a situation where there were multiple lesions to assess, [22] first segmented all lesions on liver-centered CT images, then used a longitudinal correspondence module to find matching pairs of lesions from the segmentation maps, and finally computed the lesions' evolution. Even if promising, multi-staged methods suffer from error propagation between stages. [7] removed the need for a prior segmentation step and relied on a conventional registration algorithm between longitudinal CT pairs to obtain the deformation fields and their Jacobian determinants. This information was then used to train a classifier network assessing systemic sclerosis interstitial lung disease.
Image registration can be a key step in lesion monitoring. As our long-term goal is to use registration results to develop tools for the monitoring of metastatic breast cancer, this work focuses on registration and on how to improve its accuracy. We present hereafter the most used registration approaches and the good practices which contribute to make them successful.
Established registration tools There is a large body of literature on conventional registration methods ( [27,41]). Some of today's most used registration tools include the Demons method ( [45]) using non-parametric diffeomorphic displacement fields, the Elastix toolbox ( [26]) based on cubic B-splines, and the Advanced Normalization Tools (ANTs) ( [1]) parametrizing the velocity field and relying on bi-directional diffeomorphisms. We compared our proposed approach to the last two methods, as they perform well with different datasets ( [27]) and their pyramidal coarse-to-fine optimization has inspired recent works on DL-based registration which is presented hereafter.
DL methods In recent years, different types of DL-based registration approaches have been proposed ( [8]). With the monitoring of metastatic breast cancer, deformation fields are not available and are difficult to obtain. We therefore focused on unsupervised registration methods. We categorize as unsupervised methods which do not require ground truth deformation fields, segmentations, or landmarks. Typically, they learn by minimizing a dissimilarity term between the fixed and the wrapped moving image. Among the unsupervised registration methods, various CNN architectures have been proposed. The widely used U-Net ( [38]) inspired the reference Voxelmorph network first trained on brain MRIs ( [4]). [42] added dilated convolutions to the encoder path of Voxelmorph to capture a wider field of view from lung MRI images. Later, [46] proposed a U-Net-based cascade network applied to cardiac cine MRI and chest CT data to perform affine and deformable registrations in stages, at the cost of an increase in complexity. On 3D PET images, [29] also proposed an iterative DL-based registration method to reconstruct motion compensated PET images. While still using a U-Net-shaped architecture, [14] proposed a coarse-to-fine training mimicking the best performing conventional approaches. With a similar training strategy, [33] won the Learn2Reg 2021 MICCAI challenge 1 using a pyramidal network named LapIRN. Contrary to the above cited methods, [33] did not only impose a smoothness constraint on the deformation field through its gradient, but also enforced diffeomorphic transformations using stationary vector fields under the Log-Euclidean framework, as in [12].
Despite recent efforts, when using DL-based registration techniques no major gain in registration performance was reported. Conventional iterative optimization methods still yield better accuracy in many tasks such as inter-patient alignment or intra-patient lung motion registration ( [21]). Large databases would be needed to learn the network parameters and produce accurate deformation fields. Moreover, we can argue that the generalization ability of the trained network is questionable when the deformation patterns do not repeat consistently across the dataset.
Deep Image Prior registration To alleviate the dependency to a dataset and the need for repeated deformation patterns, we adapted the Deep Image Prior (DIP) framework ( [43]) to propose a learning-free method for deformable medical image registration. Contrary to standard DL-based approaches, DIP does not learn from a database but relies on a single image. It uses a deep architecture not to summarize the information across samples but as a prior. Effectively, the architecture plays the role of a parametric model in an optimization problem restricting the solution space. While DIP was initially designed for denoising and inpainting tasks ( [43]), we adapted it here to medical image registration. [28] suggested image registration as a potential application of DIP, with preliminary results in the context of 2D brain MR images. DIP has also been used on medical data for the reconstruction of CT and PET images ( [17,3]). To perform these reconstructions, an untrained deep network was used as a denoiser and iteratively optimized by a conventional algorithm. Apart from the use of an initial reconstruction and classical regularization, images from other modalities were also given as input to condition the network output.
With this work, we establish a link between conventional, deep learning and DIP-based approaches. In particular, we focus on the role of the untrained network in parameterizing the displacement field, showing that each architecture induces an implicit regularization when used within a conventional optimization scheme. Similar observations have been made in the context of DL for inverse problems [31] and [13]. As a result, well-structured architectures (e.g. [33]) provide better over-parameterizations for the deformation fields, both in the supervised and in the untrained cases. We also investigate the role of the input (random, moving or moving and fixed images) and potential interactions of the untrained network with conventional supervised approaches.

Image registration
Given a pair of fixed and moving (F , M ) images, deformable image registration aims to estimate a dense deformation field Φ such that the wrapped image W = M (x + Φ(x)) is aligned with F for each voxel x. Given a metric S(., .) measuring the dissimilarity between two images, conventional registration methods choose a parametrization of transformation and optimize it through a cost function of the form: arg min However, the problem stated Eq. 1 is ill-posed: it belongs to the nonlinear transformation class ( [35]). By explicitly constraining the transformation, parametric approaches can become well-posed ( [40]), at the cost of limiting the range of admissible transformations. Constraining an over-parameterized or complex deformable transformation is effectively done by adding a regularization term R smooth (Φ) to the objective function Eq. 1. Regularization terms incorporate priors, which if correctly chosen, guide the optimization towards better estimates of the deformation fields. Weighted by λ smooth , the most common regularization term enforces smoothness onto the displacement field by penalizing the spatial derivatives of Φ.
DL-based registration approaches ( [16]) adapt the cost function in Eq. 2 by modeling Φ to be the output of a CNN trained on a dataset (See Eq. 3). The inputs of this network are the (F , M ) pair, and the warping operation to get W is performed using a spatial transformer layer ( [23]), handling the transformation, the sampling and the interpolation. Regarding the regularization, DLregistration based methods continue relying on a term enforcing the smoothness of the transformation R smooth , as done in conventional approaches. Smoothness is often completed by R dif f eo , which enforces the diffeomorphism of the transformation by penalizing the determinant of the Jacobian negative values. Regularization terms are weighted respectively by λ smooth and λ dif f eo . Moreover, we state that two additional regularizing priors are also implicitly added in the DL setup and influence the predicted deformation field Φ: the first one R dataset is induced by training on a domain-specific dataset, while the second one R archi is entailed by the network architecture choice (see Fig. 1.B). Due to the dataset dependency, one limitation of DL-based registration methods is the difficulty to generalize across organs or modalities.
arg min Next, we recall the deep image prior concept, which removes the dataset dependency in Section 3.2, and present our method for registration based on an untrained network in Section 3.3.

Deep Image Prior
The DIP method proposed in [43] uses a deep architecture to denoise images using a network, but without any prior learning. Supposing X 0 a distorted image and X the network output, the fitting process is characterized by Eq. 4, with R archi the implicit prior captured by the network architecture and S DIP a reconstruction function on a single image.
arg min DIP reconstructs a noisy image (e.g. with JPEG compression noise) from a white noise image by training a generator architecture to fit the noisy image. A denoised image is obtained by stopping the training before completely fitting the noise. We propose to adapt this idea to the registration of a pair of images, where we modify the moving image to match the fixed one. Contrary to the original DIP, we are interested in fitting all the way up to the finest deformations.

MIRRBA (Medical Image Registration method Regularized by Architecture)
In this paper, we argue that deep architectures, as parametric models with high capacity, are powerful representations for deformation fields. They can thus be exploited as implicit priors for image registration in iterative optimization schemes without a training stage ( [19,13,31]). Although DIP priors have been explored in the context of image reconstruction, there is no prior in depth study of such architecture priors in the context of image registration.
Since there is no learning step, the training dataset no longer influences the transformation. Thus, we transform Eq. 3 into Eq. 5 to directly optimize a patient-specific CNN for the pair of interest (F , M ), as traditionally done with iterative optimization methods (see Fig. 1 By optimizing Eq. 5, we find the best warped image allowed by the overparametrization of the architecture. Carefully designing the architecture to be used, it is possible to integrate several of the commonly used tricks in conventional approaches. In this paper, we rely on the LapIRN architecture ( [33]), which incorporates filtering through convolutions, pyramidal coarse-to-fine refinement, and interpolation steps with down-and transpose convolutions. Eq. 5 is optimized with an SGD for an input pair of images. Our idea applies to other architectures and optimization algorithms.

Dataset description
We ran our experiments on images from a private dataset, acquired in the context of the ongoing prospective multicentric EPICURE seinmeta study (NCT03958136) for metastatic breast cancer monitoring([10]) 2 . Patients underwent between two and three PET/CT acquisitions, corresponding to pre-, early-(after a month) and mid-treatment time points. A total number of 110 pairs of images were obtained, a pair being composed of a pre-treatment and either an early-or a mid-treatment image (58 and 52 images respectively). Images were acquired at two different centers. The 54 pairs of images from center A were obtained using a Philips Vereos or a GE Discovery PET/CT imaging systems, while the 56 pairs of images from center B were acquired on two different dual-slice Siemens Biograph PET/CT scanners.
Since we are interested in lesion monitoring in the context of metastatic breast cancer, we worked only on PET images, as shown useful in previous studies [6] and [2]. Expert physicians manually segmented all lesions. The brain and the bladder were also delineated, since they can be useful to mask irrelevant regions for patient response assessment. All PET images were normalized by the standardized uptake value (SUV) ( [24]) and resampled to obtain an isotropic resolution of 1 × 1 × 1 mm 3 with 200 pixels for each side. Besides, no prior registration of any kind was performed.

Architectural implementation details
Our MIRRBA method relies on the LapIRN network architecture proposed in [33]. LapIRN is a pyramidal network with N = 3 depth levels, each level being composed of a feature encoder, a set of residual blocks and a feature decoder, as shown in Fig. 1. For each level L i∈{1,2,3} , input images are downsampled by a factor 0.5 (N −1) using a trilinear interpolation. Hence, for the coarsest level L 1 the image resolution is divided by 4, while for the finest L 3 it remains identical. Moreover, a scaling and squaring module ( [12]) enforces diffeomorphic deformations.
Network levels were trained in a coarse-to-fine manner, meaning the coarsest level L 1 is first trained alone, and then higher levels are progressively trained to refine the registration. To avoid unstable starts when training levels L i>1 , lower levels weights were frozen for a fixed number of epochs. Regarding the optimization process, the learning rate was set to 10e−4, and the Adam optimizer was used for 1000 iterations on the two lower levels and 2000 iterations on the finest. As a dissimilarity metric, we used the NCC, regularized by smooth and diffeomorphic terms, for which λ smooth and λ dif f eo were determined with a grid search to 0.1 and 1.0 respectively to minimize the overall loss function over a test set (see Eq. 5), as described below for DL-based methods training. Indeed, since we aim at using registration-based features to monitor metastatic breast cancer in future work, we need to reach precise global and local registration. More details can be found in [33].
Starting from the architecture in Fig. 1, we performed an ablation study to measure the influence of each component. Hence, we i) changed the depth of the network (network with 1, 2, 3 or 4 depth levels), ii) computed the results after training each level during the coarse-to-fine registration i.e. the coarsest, the intermediate and the finest, iii) deleted the residual connections of residual blocks to transform them into simple convolutional blocks, iv) replaced the downand up-convolutions to respectively max-pooling and upsampling operations, v) used deformable convolutions ( [11]) in the finest level, and vi) used a Gaussian noise image N (0, 0.001) (as in [28]) or concatenated fixed and moving images as network inputs. Moreover, we also vii) set both λ smooth and λ dif f eo to 0 to remove the registration-specific regularization terms (see Fig. 1.C). All these implementations are detailed in Table 1.

Reference methods implementation details
Regarding the methods used as reference, we used a grid search on a test set (as described below for DL-based methods training) to find hyperparameters yielding best performance running the ANTs pipeline ([1]), i.e. with a three resolutions coarse-to-fine optimization, a gradient step of 0.2 and a symmetry transformation penalty. We used both MI and NCC as similarity metrics, leading to SyN and SyNCC methods respectively. We also ran the Elastix pipeline ( [26]) to perform successive rigid, affine and deformable image registration with four resolutions optimized with an adaptive SGD minimizing the NCC similarity term for 1000 iterations (Elastix). As for our method, we used a penalty on the bending energy as regularization to enforce smooth deformations.
To highlight the regularizing power of the architecture in MIRRBA, we ran the pipeline without the network, optimizing directly the deformation field (initialized as a Gaussian noise N (0, 0.001)) with the Adam optimizer, MIRRBA wo Archi (see Table 1 and Fig 1.D).
Regarding deep learning-based methods, DL LapIRN ( [33]) was ran with a similar regularization balance as MIRRBA, i.e. with λ smooth = 0.1, while all other recommended settings were used. For DL Voxelmorph, we used the latest diffeomorphic version to date ( [12]), with the NCC loss and recommended settings. For methods relying on a training stage, we split our dataset into five folds, paying attention to balance data from different acquisition centers among folds. For each fold, we refer to D train and D test as the train and test dataset, respectively. We trained the DL-based approaches on D train before testing them on D test for the five folds.
All architectures were implemented with PyTorch ( [37]) and trained from scratch on a Nvidia V100 32GB SXM2 GPU.

Evaluation metrics
The first criterion we used to evaluate the registration accuracy is the detection rate, defined for an individual as the percentage of lesions presenting an overlap greater than 50% between the wrapped and fixed segmentations ( [34]). A higher percentage refers to a better detection.
To evaluate the registration more locally, we computed Dice scores between fixed and wrapped i) brain and bladder and ii) lesions segmentation masks. A Dice score close to 1 means a precise local registration, while a Dice close to 0 is unsatisfactory. Since the detection rate represents the percentage of lesions correctly detected, it is positively correlated with the Dice score of the lesions. However, since some lesions (155) are cured over time and disappear on PET images, we removed them from the Dice score and detection rate computations to avoid erroneous null values. Instead, we evaluated the capacity of a method to effectively make lesions disappear by computing its disappearing rate, or percentage of volume reduction of a lesion induced by the registration, where a complete disappearance would mean a rate of 100%.
Registration smoothness was evaluated by measuring, for every deformation field, the standard deviation of its Jacobian determinant SDJDet. Null values indicate an identity transformation and high ones disorganized and incoherent displacements. Although an optimal value is difficult to define, we sought to obtain small positive values, characterizing smooth deformations ( [33]).
Finally, we evaluated the approximate running time of each approach, using a CPU for SyN, SyNCC and Elastix, and a GPU for all other methods.

Statistical analysis
To evaluate the statistical significance of our results, we first studied their distribution. According to the Shapiro-Wilk test (testing the null hypothesis that a sample comes from a normal distribution), we cannot reject the null hypothesis. Hence, we ran a paired sample t-test on our results and considered them of statistical significance if p < 0.05.

Regularization terms
In the first experiment, we looked at the influence of the different regularization terms from Eq. 3, i.e. R smooth , R dif f eo , R dataset and R archi . To this end, we compared the results to our MIRRBA approach optimized directly (without learning) on the D test pairs, using R smooth , R dif f eo and R archi terms (see Table 1, Fig. 1.A and Eq. 5). Dataset regularization We trained the Voxelmorph (DL Voxelmorph) ( [12]) and LapIRN (DL LapIRN) ( [33]) with the four regularization terms from Eq. 3 (see Table 1 and Fig. 1.B). Predictions were made on D test images. Even if both DLbased methods reached similar numerical results (see Table 2) and coherent visual ones (see Fig. 2 and supplementary material), MIRRBA performed better than both of them globally and locally. Indeed, from Table 2, we note that running a patient-specific optimization, i.e. removing R dataset and optimizing an untrained network MIRRBA compared to DL Voxelmorph and DL LapIRN, helps improve the results especially for the lesions. The organs' Dice score improved by 6% and 5% while the Dice of the lesions improved by 52% and 65% for DL Voxelmorph and DL LapIRN respectively. Moreover, MIRRBA presented lower SDJDet values, hence produced smoother deformations than training-based approaches. On the other hand, DL Voxelmorph and DL LapIRN had a higher disappearing rate.
Registration-specific regularization To understand the influence of the registrationspecific regularization terms R smooth and R dif f eo , we removed them from Eq. 5 in MIRRBA wo Regu (see Table 1) and Fig. 1.C). As quantitatively show in Table  3, MIRRBA presented Dice scores for the organs and lesions which are respectively 6% and 4% higher compared to MIRRBA wo Regu. Yet, the detection and disappearing rates were increased for MIRRBA wo Regu over MIRRBA.
Architecture regularization Finally, to also understand the impact of the regularization power of the network architecture on the registration, we looked at the results of MIRRBA wo Archi, which directly optimizes the deformation field without including any specific architecture, hence not including R archi (see Table 1 and Fig. 1.D). According to Table 3, MIRRBA wo Archi was outperformed by MIRRBA for both organ and lesion segmentations by 22% and 78% respectively. Although MIRRBA wo Archi disappearing rate was higher than MIRRBA's, the method without architectural regularization presented a low detection rate confirming the structural bias of convolutional generators.

Architectural choices
To study the regularization effect of a network architecture on registration, we compared the results of various architectural choices built from the pyramidal network presented in [33] (see Table 1). Quantitative results are presented in Table 3, while qualitative ones are visible Fig. 2

and in supplementary materials.
Depth of the pyramidal network First, we modified the architecture to optimize a simple U-Net-shaped network (MIRRBA Depth 1), and pyramidal ones with two (MIRRBA Depth 2), three (MIRRBA), and four (MIRRBA Depth 4) resolutions. While MIRRBA Depth 1 was optimized on full resolution images, all other architectures were trained using a coarse-to-fine strategy.
Results show that increasing the network depth improves the Dice results, as well as the detection and disappearing rates. Indeed, MIRRBA Depth 4 presented the second highest Dice scores for both the organs and the lesions among all MIRRBA setups. Regarding the SDJDet values, they were higher when using a depth 3 or 4, than when the network was trained on less resolutions.

Trained network level
To understand the amount of information brought by each network level during the coarse-to-fine training strategy, we computed the registration after optimizing only the lower level (MIRRBA Level 1) on coarse resolution images, both lower levels (MIRRBA Level 2) on coarse and medium resolution images, and the whole network (MIRRBA) with the complete coarseto-fine training strategy.
According to Table 3, training on all three levels of the network improved the registration accuracy over training only on low resolution images. Moreover, the disappearing rate was null when the training only occurred on low resolution images.
Max-pooling and upsampling operations For each network level, we replaced the down-convolution (convolution with stride 2) by a max-pooling operation in MIRRBA Max, and the transpose convolution by an upsampling in MIRRBA Up. As presented Table 3, using these learning-free operations improved the registration accuracy globally and locally. The disappearing rate was however slightly reduced.
Residual blocks We removed the residual connections of the residual blocks (MIRRBA wo RB) to understand their influence in our architecture. Without them, we obtained null Dice scores and detection rate, as well as very high SDJDet value. Besides, wrapped images did not look realistic (see Fig. 2 in the supplementary materials).
Deformable convolutions To adapt the receptive field of the convolutions to the local scale of objects to be registered, we replaced those of the highest resolution level (i.e. level 3) by deformable convolutions (MIRRBA DefConv). While global registration results were improved over MIRRBA, local ones were not. Moreover, the SDJDet was low.
Input images Regarding the network inputs, instead of conditioning the network with the moving image, we fed it with a Gaussian noise, as in [28]. In Table 3, MIRRBA NoiseImg shows worse results than MIRRBA, except for the SDJDet. We also provided more information to the network by concatenating the fixed to the moving image in MIRRBA FixImg, which significantly improved all MIRRBA results, especially the disappearing rate.
Combining best practices Finally, we combined the best architectural variations presented above to perform registration with four resolutions, both fixed and moving images as input, max-pooling and upsampling operations, as well as residual blocks: MIRRBA Best.

Comparison to conventional methods
Conventional registration According to Table 4, even after rigid and affine preregistration, locally precise deformable registration is challenging on whole-body images for the conventional Elastix ( [26]) pipeline (see Fig. 1

in supplementary material).
On the other hand, both SyN and SyNCC ([1]) statistically performed better than MIRRBA (which used the NCC metric) on organ segmentation, while the NCC similarity metric allowed SyNCC to also reach a better accuracy on lesion segmentation. Regarding the disappearing rate, MIRRBA performed better than both SyN-based methods, whereas their SDJDet was lower than ours.
Combination of DIP and conventional registration Finally, we pushed the analysis by combining our best MIRRBA-based method with the best performing conventional SyNCC approach. To do so, we optimized MIRRBA Best, using as input the deformation fields and already registered images obtained by SyNCC.
Our assumption was that the conventional pre-registration would be improved by our method.
Results of this Combined pipeline are shown in Table 4. In terms of global and local accuracy, the combined approach outperformed both SyNCC and MIRRBA Best. In addition, the disappearing rate was significantly improved for both methods.

Impact of training on a database
When comparing MIRRBA to the DL-based approaches (DL LapIRN and DL Voxelmorph), we show that not learning registration patterns from a dataset helps to obtain precise segmentations, especially at a lesion level. Indeed, while the size of organs and their locations are relatively consistent across patients, this is not the case with the lesions. Hence, there are fewer deformation patterns that can be learned from a database. Since performing locally precise registration with a DL-based method is very challenging in this situation, our method adapts to each subject. Moreover, according to its lower SDJDet value, MIRRBA produces smoother deformations than the DL-based approaches.

Architectural choices
Regularization terms As presented in Section 3.1, removing the registrationspecific regularization terms is equivalent to solving an ill-posed problem in conventional registration methods. In our case, since we use a network to perform registration, R archi is present when running MIRRBA wo Regu. For similar Dice scores and detection rate, MIRRBA wo Regu has higher SDJDet than MIRRBA, indicating R smooth and R dif f eo help smooth the deformation field, even though they have less impact than the regularization of the architecture.
On the opposite configuration MIRRBA wo Archi, where R archi is not used but R smooth and R dif f eo are, quantitative results show that the lack of architecture negatively impacts the registration smoothness, as well as its local and global precision, making questionable the convergence of the method without architecture. Hence, we confirm that the architecture has a regularization effect on the registration, which helps to find an admissible solution.
Coarse to fine strategies and image resolution Considering the overall network, [13] affirms that an architecture ran in a DIP pipeline acts as a low-pass filter in the beginning of the optimization, allowing higher frequencies to pass only after lower ones. We observe similar results over the coarse-to-fine training strategy (see Fig. 2), where low frequencies are registered first (global registration), followed later by higher frequencies (local registration). In the same way, looking at the intermediate results during the 4-depth pyramidal optimization of Elastix (see Fig. 1 in supplementary materials), we observe a tendency to register global features before local ones.
Modifying the number of coarse resolution levels, i.e. the depth of the pyramidal architecture, we studied the impact of global structural choices to optimize the whole network. The higher SDJDet values obtained when more resolutions are used could be explained by the generation of more local transformations, hence a globally less regular deformation field. Moreover, with four resolutions, the receptive field of the coarsest level of MIRRBA Depth 4 captures the whole image (i.e. 200 × 200 × 200), explaining the high Dice scores, as well as detection and disappearing rates. Indeed, successful conventional pipelines as Elastix ( [26]) or ANTs ( [1]) also uses this type of pyramidal strategy.
Network conditioning with input images Contrary to [43] and [28], who used Gaussian noise as input for their architectures, [17] and [3] respectively improved CT and PET DIP-based image reconstructions by providing acquisitions from other modalities to their networks. We made similar observations, as MIRRBA NoiseImg reached less accurate results than when we conditioned the model with PET images, as in MIRRBA. Moreover, feeding the network with more patient information, using the fixed image as additional input to the moving image in MIRRBA FixImg, improves the results by increasing the network conditioning on a single patient.
Other architectural choices Moving inside the architecture, residual blocks can be related to diffeomorphic registration according to [39]. Indeed, stacking residual blocks in ResNets ( [18]) aims to incrementally map the embedding space to a new unknown space, each block being defined as y = F (x) + x, with x and y the respective input and output of the residual blocks, and F the residual mapping to be learned. Similarly, diffeomorphic registration models ( [5,41]) address the registration issue by piling up incremental diffeomorphic mappings. Making the link between ResNets and registration, the function F can be seen as a parametrization of an elementary deformation flow, and training a series of residual blocks as learning continuous and integral diffeomorphic operator. With our deep architecture, MIRRBA wo RB results indicate that the registration without residual blocks fails to converge. Indeed, as explained above, residual blocks allow incremental diffeomorphic mappings, and removing them leads to gradient vanishing for our patient-specific method.
Regarding the convolutions, [20] suggested deformable convolutions ( [11]) to capture larger deformations. These convolutions add 2D or 3D offsets to the regular grid sampling of standard operations. If these offsets are set to zero, deformable convolutions become standard convolutions, otherwise they modify the receptive field. Since these offsets are learnable, deformable convolutions are trained to adapt their receptive fields, in order to focus on objects of interest in classification problems for instance. According to the high detection rate, this is performed globally by MIRRBA DefConv. Yet, the resulting low Dice scores show that MIRRBA DefConv did not achieve precise local registration, which can be explained by the low SDJDet value, showing that only smooth and regular deformations occurred instead of locally irregular ones. Indeed, deformable convolutions might need a more complex integration in the architecture ( [30]). Moreover, the additional learnable parameters and adaptive receptive field of deformable convolutions are to some extent redundant with the pyramidal structure, and make it harder to train.
Considering the amount of parameters to optimize in a network, we looked at the results of MIRRBA Max and MIRRBA Up compared to MIRRBA. Even if the learnable down-and up-convolutions are now common to respectively increase or decrease image dimensions within a network, max-pooling and upsampling operations were originally used. As in conventional registration methods, these operations are not learnable. Hence, their results do not depend on optimization parameters, and they help to control the overfitting in a traditional trainable setup. Therefore, the good results obtained using max-pooling and upsampling operations in MIRRBA Max, MIRRBA Up and MIRRBA Depth 4 Max Up may be explained by the fewer number of parameters to fit (see Fig. 2 in supplementary materials).

Lesion segmentation
Regarding lesion registration, Fig. 3 shows different lesion evolution scenarios. We can see that the moving masks adapt to the fixed ones when the lesion does not disappear. Although the disappearing lesion from the first column was not deleted, it was reduced by the registration algorithm. Indeed, lesion disappearance implies deformations which are not smooth nor diffeomorphic, hence the good disappearing rate of MIRRBA wo Regu and the probable need to adapt our registration strategy to the specific problem of lesion disappearance in future work.

Conclusion
In this paper, we propose an alternative method to perform image registration using a neural network without the typical learning stage on a database. We formalize the registration problem by following the conventional approaches relying on image-based similarities and regularization terms, but also explicitly consider the dataset and architecture bias. Indeed, our study was motivated by recent work on DIP, implying that neural networks create an inductive bias when learning from a database, but also an intrinsic structural bias induced by the architecture [19]. These biases are shown sufficient to solve certain image processing tasks. Our proposition also resonates with [13], who suggested that deep convolutional neural networks process low-frequency information first, to later focus on the finer deformations, both desirable properties for a registration algorithm. We integrate the LapIRN network from [33], who developed an effective pyramidal architecture, tested in the DL-based set-up. Here, we further demonstrate that beyond any prior coming from the dataset learning step, the architecture design has an important effect on the registration results, acting as an implicit regularizer. Our study also shows the impact of some of the architecture components, particularly the residual blocks, and we justify this behavior by making a link with findings from [39]. Moreover, we found that for our application, a pyramidal architecture capturing the whole image with a limited amount of parameters to optimize, as conventional registration methods, provides precise registration results.
The architectural prior seems to be a better option than learning from data in cases where there are no positional consistency, which is the case with metastatic breast cancer lesions, which arbitrarily vary in position, size and number. Indeed, finding a set of network parameters allowing precise registration for a whole dataset is a challenging task. Our approach enables to correctly register active organs such as the brain and the bladder, which could be used to automatically propagate annotations masking regions irrelevant for patient response evaluation. Although the Dice scores are relatively low for the lesions, we obtain good detection values and improve the disappearance rate.
Our approach makes a step in bridging conventional and DL-based methods for image registration, and provides several suitable approaches for the challenging 3D full-body longitudinal registration problem. We demonstrated the possibility to perform both global and local registration on whole body medical images using a network but without suffering from dataset bias. In future work, we would study the feasibility of extracting registration-based feature from our method to monitor lesion evolution without depending on manually performed segmentations.  Grayscale color indicates good overlapping. We can see that the higher the depth, the more precise the registration. A simple U-Net-shaped architecture as MIRRBA Depth 1 produces transformations of very low accuracy, while the more resolutions are used, the more precise the registration. The all green MIRRBA wo RB image is due to the non-convergence of the registration algorithm and a wrapped image not registered to the fixed one. All four other approaches produce realistic looking and coherent transformations, even if MIRRBA DefConv lacks a bit of precision around the bladder. MIRRBA NoiseImg does not reach local precise registration. Best viewed in color. Table 1. Loss terms and structural setups of the methods presented in the paper. "Pyr. net." stands for "Pyraminal network", "Sym. diffeo. transfo." for "Symmetric diffeomorphic transformation", "Max." for "Maxpooling", "Up." for "Upsampling", and "Def. conv." for "Deformable convolutions". An absence of information in the Other column indicates the use of setups described in Sections 4.2 and 4.3. Table 2. Comparison of MIRRBA to DL-based methods (Section 5.1) -Detection rate, Dice score of the organs and lesions (Dice organs and Dice lesions resp.), SDJDet, disappearing rate, and approximate computational time. Both training and inference computational times are indicated for DL-based methods. All pipelines were computed on GPU. Statistically significant improvement of our MIRRBA method over the others with p < 0.05 is indicated with *. Best results are marked in bold, and second best ones in bold-italic, except for SDJDet since no ideal value is defined.


Neural Representation Learning for Scribal Hands of Linear B

In this work, we present an investigation into the use of neural feature extraction in performing scribal hand analysis of the Linear B writing system. While prior work has demonstrated the usefulness of strategies such as phylogenetic systematics in tracing Linear B's history, these approaches have relied on manually extracted features which can be very time consuming to define by hand. Instead we propose learning features using a fully unsupervised neural network that does not require any human annotation. Specifically our model assigns each glyph written by the same scribal hand a shared vector embedding to represent that author's stylistic patterns, and each glyph representing the same syllabic sign a shared vector embedding to represent the identifying shape of that character. Thus the properties of each image in our dataset are represented as the combination of a scribe embedding and a sign embedding. We train this model using both a reconstructive loss governed by a decoder that seeks to reproduce glyphs from their corresponding embeddings, and a discriminative loss which measures the model's ability to predict whether or not an embedding corresponds to a given image. Among the key contributions of this work we (1) present a new dataset of Linear B glyphs, annotated by scribal hand and sign type, (2) propose a neural model for disentangling properties of scribal hands from glyph shape, and (3) quantitatively evaluate the learned embeddings on findplace prediction and similarity to manually extracted features, showing improvements over simpler baseline methods.

Introduction
Neural methods have seen much success in extracting high level information from visual representations of language. Many tasks from OCR, to handwriting recognition, to font manifold learning have benefited greatly from relying on architectures imported from computer vision to learn complex features in an arXiv:2108.04199v1 [cs.CV] 14 Jul 2021 automated, end-to-end trainable manner. However much of the emphasis of this field of work has been on well-supported modern languages, with significantly less attention paid towards low-resource, and in particular, historical scripts. While self-supervised learning has seen success in these large scale settings, in this work we demonstrate its utility in a low data regime for an applied analysis task. Specifically, we propose a novel neural framework for script analysis and present results on scribal hand representation learning for Linear B.
Linear B is a writing system that was used on Crete and the Greek mainland ca. 1400-1200 BCE. It was used to write Mycenaean, the earliest dialect of Greek, and was written on leaf or page-shaped clay tablets for accounting purposes. Linear B is a syllabary, with approximately 88 syllabic signs as well as a number of ideograms, which are used to indicate the commodity represented by adjacent numerals. The sites which have produced the most material to date are Knossos, Pylos, Thebes, and Mycenae (see Palmer [15] for more detail) but for the purposes of this paper, we focus specifically on Knossos and Pylos from which we collect a dataset of images of specific instances of glyphs written by 74 different scribal hands. Many signs exhibit slight variations depending on which scribe wrote them. Therefore, being able to uncover patterns in how different scribes may write the same character can inform us about the scribes' potential connections with one another, and to an extent even the ways in which the writing system evolved over time. This is a process that has previously been performed by hand [18], but to which we believe neural methods can provide new insights.
Rather than building representations of each scribe's writing style based on the presence of sign variations as determined by a human annotator, we propose a novel neural model that disentangles features of glyphs relating to sign shape and scribal idiosyncrasies directly from the raw images. Our model accomplishes this by modeling each image using two separate real-valued vector embeddings -we share one of these embeddings across glyphs written by the same scribe, and the other across glyphs depicting the same sign, thereby encouraging them to capture relevant patterns. These embeddings are learned via three networks. The first is a decoder that attempts to reconstruct images of glyphs from their corresponding embeddings using a series of transpose convolutional layers. We penalize the reconstructed output based on its mean squared error (MSE) vs. the original. In addition, we also use two discriminator networks, which attempt to predict whether a given image-embedding pair actually correspond to one another or not, both for scribe and sign embeddings respectively. These losses are backpropagated into those original vectors to encourage them to capture distinguishing properties that we care about. This model has many potential downstream applications, as paleography has long been an integral part of the study of the Linear B texts [14]. For a handful of examples from prior work, Bennett's study of the Pylos tablets [1] revealed the existence of a number of different scribal hands, all identified on the basis of handwriting, since scribes did not sign their work. This research was carried out even before the decipherment of Linear B. Later, Driessen [5] used the archaic sign forms on Linear B material from the Room of the Chariot Tablets at Knossos, in conjunction with archaeological evidence, to argue that these tablets pre-dated the remainder of the Linear B tablets from Knossos. Skelton [18] and later Skelton and Firth [20] carried out a phylogenetic analysis of Linear B based on discrete paleographical characteristics of the sign forms, and were able to use this analysis to help date the remainder of the Linear B tablets from Knossos, whose relative and absolute dates had already been in question. Our hope is that leveraging the large capacity of neural architectures can eventually lead to further strides in these areas.
Additionally, qualitative assessments of Linear B paleography have long been used to date tablets whose history or findplaces are uncertain. This has been vital in the case of the Knossos tablets, which were excavated before archaeological dating methods were well-developed [7], and the tablets from the Pylos megaron, which seem to date to an earlier time period than the rest of the palace [19]. While qualitative techniques have been fruitful, a key motivating question in our work is to what extent we can use quantitative methods to improve on these results.
Going back further, Linear B was adapted from the earlier writing system Linear A, which remains undeciphered. There is a hundred-year gap between the last attestation of Linear A and the first attestation of Linear B. With a more thorough understanding of the evolution of these two writing systems, can we reconstruct the sign forms of the earliest Linear B, and elucidate the circumstances behind its invention? With this in mind, our approach explicitly incorporates a decoder network capable of reconstructing glyphs from learned vector representations, that we hope can eventually be useful for such investigation.
In total, three writing systems related to Linear B remain undeciphered (Cretan Hieroglypic, Linear A, and Cypro-Minoan), and a fourth (the Cypriot Syllabary) was used to write an unknown language. One motivating question is whether a better understanding of the evolution of sign forms would help contribute additional context that could make a decipherment possible. In particular, scholars are not yet in agreement over the sign list for Cypro-Minoan, and how many distinct languages and writing systems it may represent [6]. Understanding variation in sign forms is foundational for this work [14], and it makes sense to test our methods on a related known writing system (i.e. Linear B) before attempting a writing system where so much has yet to be established.
Overall these are just some of the applications that we hope the steps presented here can be directed towards.

Prior Work
One of our primary points of comparison is the work of Skelton [18], who manually defined various types of systematic stylistic variations in the way signs were written by different scribal hands, and used these features to analyze similarity between scribes via phylogenetic systematics. In a sense this work can serve as our human baseline, giving insight into what we might expect reasonable While neural methods have not previously been used in analyzing this particular writing system, there is established literature on applying them to typography more broadly. For example, Srivatsan et al. [21] used a similar shared embedding structure to learn a manifold over font styles and perform font reconstruction. There has also been extensive work on writer identification from handwritten samples [2,17,25,4], as well as writer retrieval [3]. However these approaches rely on supervised discriminative training, in contrast to ours which is not explicitly trained on predicting authorship, but rather learning features.

Methods
In this section we describe the layout of our model and its training procedure. At a high level, our goal is to train it on our dataset of images of glyphs, and in doing so learn vector embeddings for each scribe that encode properties of their stylistic idiosyncrasies -for example which particular variations of each sign they preferred to use -while ignoring more surface level information such as the overall visual silhouettes of the signs we happen to observe written by them. To do this, we break the roles of modeling character shape and scribal stylistic tendencies into separate model parameters, as we now describe.
Our model, depicted in Figure 2, assigns each sign and each scribal hand a vector embedding to capture the structural and stylistic properties of glyphs in our data. These embeddings are shared by glyphs with the same sign or scribe, effectively factorizing our dataset into a set of vectors, one for each row and column. The loss is broken into two parts, a reconstruction loss and a discriminative loss. The reconstruction loss is computed by feeding a pair of sign and scribe embeddings to a decoder, which then outputs a reconstruction of the original glyph that those embeddings correspond to. A simple MSE is computed between this and the gold image. The discriminative cross-entropy loss is computed by the discriminator, a binary classifier which given an image of a glyph and a scribe embedding must predict whether or not the two match (i.e. did that scribe write  Fig. 2. Diagram of the model -each row is assigned a scribe embedding, and each column is assigned a sign embedding, which feed into a decoder and discriminator to collectively compute the loss function. that glyph). We train this on equally balanced positive and negative pairs. A similar discriminator is applied to the sign embeddings. In our experiments, we ablate these three losses and examine the effects on performance.
More formally, suppose we have a corpus of images of glyphs X, where X i ∈ [0, 1] 64×64 depicts one of J characters and is authored by one of K scribes. We define a set of vectors Y , where each Y j ∈ R d represents features of sign j. Similarly we define a set of vectors Z, where each Z k ∈ R d represents features of scribe k. Thus, a given X i depicting sign j and authored by scribe k is represented at a high level in terms of its corresponding Y j and Z k , each of which is shared by other images corresponding to the same character or scribal hand respectively.
For a given X i , the decoder network takes in the corresponding pair of vectors Y j and Z k , and then passes them through a series of transpose convolutional layers to produce a reconstructed imageX i . This reconstruction incurs a loss on the model based on its MSE compared to the original image.
Next, the scribal hand discriminator with parameters φ passes X i through several convolutional layers to obtain a feature embedding which is then concatenated with the corresponding Z k and passes through a series of fully connected layers, ultimately outputting a predicted p(Z k |X i ; φ) which scores the likelihood of the scribal hand represented by Z k being the true author of X i . Similarly, a sign discriminator with parameters ψ predicts p(Y j |X i ; ψ), i.e. the likelihood of X i depicting the sign represented by Y j . These discriminators are shown both one true match and one false match for an X i at every training step, giving us a basic cross-entropy loss.
These components all together yield the following loss function for our model: where j and k are randomly selected such that j = j and k = k .
The format of this model lends itself out of the box to various downstream tasks. For example, when a new clay tablet is discovered, it may be useful to determine if the author was a known scribal hand or someone completely new. The discriminator of our model can be used for exactly such as a purpose, as it can be directly queried for the likelihood of a new glyph matching each existing scribe. Even if this method is not exact enough to provide a definitive answer, we can at least use the predictions to get a rough sense of who the most similar known scribes are, despite never having seen that tablet during training.
Also, while our model treats the scribe and sign embeddings as parameters to be directly learned, one could easily imagine an extension in which embeddings are inferred from the images using an encoder network. Taking this further, we could even treat those embeddings as latent variables by adding a probabilistic prior, such as in Variational Autoencoders [9]. This would let us compute embeddings without having trained on that scribe's writing, taking the idea from the previous paragraph one step further, and letting us perform reconstructions of what we expect a full set of signs by that new scribe to look like.

Architecture
We now describe the architecture of our model in more detail. Let F i represent a fully connected layer of size i, T represent a 2x2 transpose convolutional layer, C be a 3x3 convolution, I be an instance normalization layer [23], R be a ReLU, M be a 2x2 blur pool [24]. The decoder network is then F 1024 * 4 * 4 → 4 × (T → 2 × (C → I → R)), where each convolutional layer C reduces the number of filters by half. The discriminator's encoder network is then C → 4 × (3 × (C → R) → M ) → F 16 , where each convolutional layer C increases the number of filters by double, except for the first which increases it to 64. Having encoded the input image, the discriminator concatenates the output of its encoder with the input scribe/sign embedding, and then passes the result through 7 fully connected layers to obtain the final output probability.

Collection
There has been previous work on transcribing the tablets of Knossos by Olivier [11] and Pylos by Palaima [13]. These authors present their transcriptions in form of a rough table for each hand containing every example of each sign written by them. While these are an important contribution and certainly parseable by a human reader, they are not of a format readily acceptable by computer vision methods. Since our model considers the atomic datapoint to be an individual glyph, we must first crop out every glyph from each table associated with a scribal hand. This is a nontrivial task as the existing tables are not strictly grid aligned, and many glyphs overlap with one another, which requires us to manually paint out the encroaching neighbors for many images. Ultimately we were able to convert the existing transcriptions into a digital corpus of images of individual glyphs, labeled both by scribal hand as well as sign type. This process was carried out by two graduate students under the guidance of a domain expert. Before being fed into our model, images were also binarized via Otsu's thresholding [12].
The digitization of this dataset into a standardized, well-formatted corpus is of much potential scientific value. While this process was laborious, it opens the doors not just to our own work, but also to any other future research hoping to apply machine learning methods to this data. Of course had this been a high-resource language such as English, we would have been able to make use of the multitude of available OCR toolkits to automate this preprocessing step. As such, one possible direction for future work is to train a system for glyph detection and segmentation, for which the corpus we have already created can serve as a source of supervision. This would be useful in more efficiently cropping similarly transcribed data, whether in Linear B or other related languages.

Augmentation
Due to the limited size of this corpus, there is a serious risk of models overfitting to spurious artifacts in the images rather than properties of the higher level structure of the glyphs. To mitigate this, we perform a data augmentation procedure to artificially inflate the number of glyphs by generating multiple variants of each original image that differ in ways we wish for our model to ignore. These are designed to mimic forms of variation that naturally occur in our data -for example, glyphs may not always be perfectly centered within the image due to human error during cropping, and strokes may appear thicker or thinner in some images based on the physical size of the sign being drawn relative to the stylus' width, or even as a result of inconsistencies in binarization during preprocessing. Following these observations, the specific types of augmentation we use include translation by 5 pixels in either direction horizontally and/or vertically, as well as eroding or dilating with a 2x2 pixel wide kernel (examples shown in Figure 4). All together we produce 27 variations of each image. Doing this theoretically encourages our model to learn features that are invariant to these types of visual differences, and therefore more likely to capture the coarser stylistic properties we wish to study.

Statistics
Our dataset represents 88 Linear B sign types, written by 74 scribal hands (of these 48 are from Knossos and 26 are from Pylos). For each scribal hand, we observe an average of 56 glyph images, some which may represent the same sign.
Conversely each sign appears an average of 47 times across all scribal hands. This gives us a total of 4, 134 images in our corpus, which after performing data augmentation increases to 111, 618.

Experiments
Having trained our model, we need to be able to quantitatively evaluate the extent to which it has learned features reflective of the sorts of stylistic patterns we care about, something which is nontrivial given the lack of any supervised ground truth success criteria. There are two metrics we use to accomplish this.
First, given that we expect scribes of the same findplace to share stylistic commonalities, we train a findplace classifier based on our scribal hand embeddings to probe whether our model's embeddings contain information that is informative for this downstream task. Secondly, since we have human generated scribal hand representations in the form of Skelton's manual features [18], we can measure the correlation between them and our model's learned manifold to see how well our embeddings align with human judgements of relevant features. We now describe these evaluation setups in more detail.

Findplace Classification
The first way we evaluate our model is by probing the learned embeddings to see how predictive they are of findplace, a metadata attribute not observed during training. Since our goal is to capture salient information about broad similarities between scribal hands, one coarse measure of our success at this is to inspect whether our model learns similar embeddings for scribes of the same findplace, as they would likely share similar stylistic properties. One reason for this underlying belief is that such scribes would have been contemporaries of each other, as can be infered from the fact that existing Linear B tablets were written "days, weeks, or at most months" before the destruction of the palaces [14]. For any scribes working together at the same place and time, we observe similar writing styles, likely as a result of the way in which they would have been trained. The number of scribes attested at any one palatial site is quite low, which would seem to rule out a sort of formal scribal school as attested in contemporary Near Eastern sources, and instead suggests that scribal training took place as more of a sort of apprenticeship, or as a profession passed from parents to children, as is attested for various sorts of craftsmen in the Linear B archives [14]. Thus, we expect that scribes found at the same findplace would exhibit similar writing styles, and therefore measuring how predictive our model's embeddings are of findplace serves as a reasonable proxy for how well they encode stylistic patterns.
To this end, we train a supervised logistic regression classifier (one per model) to predict the corresponding findplaces from the learned embeddings, and evaluate its performance on a held out set. In principal, the more faithfully a model's embeddings capture distinguishing properties of the scribal hands, the more accurately the classifier can predict a scribe's findplace based solely on those neural features. We perform k = 5 fold cross validation, and report average F 1 score over the 5 train/test splits. Our final numbers reflect the best score over 15 random initializations of the classifier parameters for each split.
Findplaces for the Knossos tablets were described in Melena and Firth [10], and for Pylos in Palaima [13]. However, there are some findplaces from which we only have one or two scribal hands in our corpus. As supervised classification is not practically possible in such cases, we exclude any findplaces corresponding to 2 or fewer scribal hands from our probing experiments. This leaves us with 63 scribal hands, labeled by 8 different findplaces to evaluate on.

QVEC
While findplace prediction is a useful downstream metric, we also wish to measure how well the neural features our model learns correlate with those extracted by expert humans. In order to do this, we use QVEC [22] -a metric originally developed for comparing neural word vectors with manual linguistic featuresto score the similarity between our learned scribal hand embeddings and the manual features of Skelton [18]. At a high level, QVEC works by finding the optimal one-to-many alignment between the dimensions of a neural embedding space and the dimensions of a manual feature representation, and then for each of the aligned neural and manual dimensions, summing their Pearson correlation coefficient with respect to the vocabulary (or in our case the scribal hands). Therefore, this metric provides an intuitive way for us to quantitatively evaluate the extent to which our learned neural features correlate with those extracted by human experts. As a side note, while the original QVEC algorithm assumes a one-to-many alignment between neural dimensions and manual features, we implement it as a many-to-one, owing to the fact that in our setting we have significantly more manual features than dimensions in our learned embeddings, implying that each neural dimension is potentially responsible for capturing information about multiple manual features rather than vice versa as is the case in the original setting of Tsvetkov et al. [22].

Baselines
Our primary baseline is a naive autoencoder, that assumes a single separate embedding for every image in the dataset, therefore not sharing embeddings across signs or scribal hands. Its encoder layout is the same as that of our discriminator, and its decoder is identical to that of our main model so as to reduce differences due to architecture capacity. In order to obtain an embedding for a scribal hand using this autoencoder, we simply average the embeddings of all glyphs produced by that scribe. We report several ablations of our model, indicating the presence of the reconstruction loss, scribal hand discriminator loss, and sign discriminator loss by +Recon, +Scribe, and +Sign respectively. Finally, we provide the F 1 score for a naive model that always predicts the most common class to indicate the lowest possible performance.

Training Details
We use sign and scribe embeddings of size 16 each for all experiments. Our model is trained with the Adam [8] optimization algorithm at a learning rate of 10 −4 and a batch size of 25. It is implemented in Pytorch [16] version 1.8.1 and trains on a single NVIDIA 2080ti in roughly one day. The classifier used for evaluation is trained via SGD at a learning rate of 10 −3 and a batch size of 15. Table 1. (Left) Best cross validated F1 achieved by logistic regression findplace classifier on embeddings produced by ablations of our model, compared to an autoencoder, and naive most common label baseline. (Right) QVEC [22] score between models' embeddings and manual features of Skelton [18].

Model
Findplace  Table 1 shows results of a findplace classifier trained on our system's output vs. that of an autoencoder, with the score for a most common baseline as a lower bound. We see that our full model with all three losses achieves the highest score, with reduced performance as they are successively removed. The worst score comes from the model trained on purely discriminative losses with no decoder, suggesting that the reconstruction loss provides important bias to prevent overfitting given the limited dataset size. This also makes sense given the strong performance by the vanilla autoencoder. For reference, Skelton [18] provides manual feature representations for 20 of the scribal hands in our corpus (17 of which are from the 6 findplaces shared by at least one other scribe in that set). If we perform a similar probing experiment with the hand crafted features for those 17 scribes, we get an F 1 score of 0.37. While this result is useful as an indicator of where the upper bound on this task may be, it is very important to note that this experiment is only possible on a small subset of the scribes we evaluate our model on and is therefore not directly comparable to the results in Table 1. A more direct comparison to the neural features of our model would require extracting manual features for the remaining 54 scribal hands, which is beyond the scope of this work. Table 1 also shows the QVEC similarity between each of the models' learned embeddings and the manual features of Skelton [18]. Once again, our system ranks highest, with the autoencoder also performing more strongly than our ablations. Interestingly, we see here that our discriminative losses appear more important than the reconstructive loss, suggesting that by encouraging the model to focus on features that are highly distinguishing between scribes or signs we may end up recovering more of the same information as humans would. Figure 5 shows example reconstructions from our decoder. These reconstructions are based solely on the scribe and sign representations for the corresponding row and column. We see that the model is able to realistically reconstruct the glyphs, showing that it is very much capable of learning shared representations for sign shape and scribal hand styles, which can then be recombined to render back the original images. This would imply that our embedding space is able to adequately fit the corpus, and perhaps requires some additional inductive bias in order to encourage those features to be more high level.
There are also additional downstream tasks we have not yet evaluated on that may be worth attempting with our system. The most obvious is that of phylogeny reconstruction. Now that we have trained scribal hand embeddings, it remains to be seen how similar a phylogenetic reconstruction based on these features would resemble that found by Skelton [18]. This may give us a better idea as to whether the embeddings we have computed are informative as to the broader hierarchical patterns in the evolution of Linear B as a writing system. Example reconstructions from our model for various scribal hands. The model manages to capture subtle stylistic differences, but we see a failure mode in the middle row right column where the output appears as a superposition of two images.

Conclusion
In this paper we introduced a novel approach to learning representations of scribal hands in Linear B using a new neural model which we described. To that end, we digitized a corpus of glyphs, creating a more standardized dataset more readily amenable to machine learning methods going forwards. We benchmarked the performance of this approach compared to simple baselines and ablations by quantitatively measuring how predictive the learned embeddings were of held out metadata regarding their findplaces. In addition, we directly evaluated the similarity of our scribal hand representations to those from existing manual features. We believe this work leaves many open avenues for future research, both with regard to the model design itself, as well as to further tasks it can be applied towards. For example, one of the key applications of manual features in the past has been phylogenetic analysis, something which our neural features could potentially augment for further improvements. This type of downstream work can further be used in service of dating Linear B tablets, understanding the patterns in the writing system's evolution, building OCR systems, or even analyzing other related but as of yet undeciphered scripts.


A Semi-Supervised Classification Method of Apicomplexan Parasites and Host Cell Using Contrastive Learning Strategy

A common shortfall of supervised learning for medical imaging is the greedy need for human annotations, which is often expensive and time-consuming to obtain. This paper proposes a semi-supervised classification method for three kinds of apicomplexan parasites and non-infected host cells microscopic images, which uses a small number of labeled data and a large number of unlabeled data for training. There are two challenges in microscopic image recognition. The first is that salient structures of the microscopic images are more fuzzy and intricate than natural images' on a real-world scale. The second is that insignificant textures, like background staining, lightness, and contrast level, vary a lot in samples from different clinical scenarios. To address these challenges, we aim to learn a distinguishable and appearance-invariant representation by contrastive learning strategy. On one hand, macroscopic images, which share similar shape characteristics in morphology, are introduced to contrast for structure enhancement. On the other hand, different appearance transformations, including color distortion and flittering, are utilized to contrast for texture elimination. In the case where only 1% of microscopic images are labeled, the proposed method reaches an accuracy of 94.90% in a generalized testing set.

I. INTRODUCTION
Apicomplexans comprise a group of intracellular protozoan parasites, including Plasmodium, Babesia, and Toxoplasma, which are amongst the most prevalent and morbiditycausing pathogens in humans and animals worldwide, and they live inside host cells Erythrocytes. Plasmodium is the causative agent of malaria, which impacts over 200 million individuals and kills over 300,000 children annually. Babesiosis caused by Babesia is a disease with many clinical features that are similar to those of malaria. Toxoplasma is estimated to infect one-third of the world's population, which is the leading cause of infectious retinitis in children and is life-threatening in pregnancy and to the immunocompromised [1].
The risk of illness and death could be significantly reduced with accurate and affordable diagnostic testing. Measuring parasite infection by direct microscopy observation remains relatively widespread as a point-of-care diagnostic in clinical and epidemiological settings [2].
Recently, deep learning techniques become a popular choice in both computer vision and the medical imaging community [3]. For parasites and Erythrocytes recognition, deep learning models have obtained impressive results. A *The authors are with the Graduate School of Information, Production and Systems, Waseda University, Kitakyushu-shi, 808-0135 Japan.
deep learning model using hand-craft features in the shape to classify healthy and abnormal Erythrocytes [4]. Wellknown Convolutional Neural Networks (CNNs), including Inception v3 [5], LeNet, AlexNet and GoogLeNet [6], are used to identify Plasmodium parasites. Furthermore, a solution for multiple apicomplexan parasites and Erythrocytes [7] is given based on the fact that Toxoplasma, Babesia, Plasmodium, and Erythrocyte are variant in morphology under microscopy. However, the outlined methods have been designed in the prerequisite that all the training data has accurate human annotations. In practical clinical scenarios, lack of labels and annotations is a common and fatal problem [8] due to the shortage of specialists in diagnostic imaging despite the increasing spread of equipment.
In addition, there are two challenges for microscopic image recognition. The first is that salient structures of microscopic images are more fuzzy and intricate than realworld objects', which leads to that microscopic images are not as distinguishable as real-world images. The second is that insignificant textures, like image background staining, lightness, or contrast level, are variant a lot in samples from different clinical scenarios, while the available training data [9] is of a similar pattern in each category, which may lead to poor generalization in the real-world application.
In this paper, we propose a semi-supervised learning (SSL) solution for the classification of apicomplexan parasites and host cell microscopic images, considering that a small number of reliable annotations are relatively affordable and it is easy to acquire a large amount of unlabeled data. The proposed method contains data pre-processing, a CNNbased feature extractor, and a multilayer perceptron (MLP) classifier.
To address these challenges, we learn a distinguishable and appearance-invariant representation by using the contrastive learning strategy. For structure enhancement, the macroscopic images which share similar shape characteristics in morphology are introduced for contrast. We encourage the similarity between the microscopic images and their corresponding macroscopic images, and the dissimilarity between the unrelated images. For texture elimination, different appearance transformations, including color distortion and flittering, are utilized for contrast. We maximize the similarity between two different appearance transformed views of the same image, while simultaneously minimizing the similarity between different images.

II. PROBLEM STATEMENT
Microscopic images from parasite infection samples of preserved slides stained with Giemsa [9] including 5758 Plasmodium, 5741 Toxoplasma, 5878 Babesia, and 6981 Erythrocytes, which are non-infected host cells, are used as Micro data. 1000 Micro data in each category are split as a testing set. To evaluate the generalization ability of a model, the testing set is transformed to a color dropped version, which will never be used as data augmentation in the training phase. For the training phase, we suppose that only 50 Micro data in each category are labeled, about 1% of the whole data set, in a practical clinical scenario. According to the parasitologist-level knowledge, Plasmodium is ring-shaped, Toxoplasma is generally banana-shaped, Babesia is typically double-pear-shaped and Erythrocytes resemble an apple. We photo 500 macroscopic images each of ring, banana, double-pear, and apple as Macro data, which best match the Micro objects of interest, instead of using images from the Internet [7], since there are usually irrelative items and confusing background, and collecting from the Internet takes much more effort than photoing. Macro data share the same label with their corresponding Micro data.
Therefore, we have 2000 Macro data which is fully labeled with y s ∈ {0, 1, 2, 3} as source data, 200 labeled Micro data X T = {(x t , y t )} T t=1 with y t ∈ {0, 1, 2, 3} as target data, and 20158 unlabeled Micro data X U = {x u } U u=1 in the training set. Fig. 1 shows some samples. As illustrated in Fig. 2, here are three modules in our proposed method. Firstly, Macro data is transformed to its adaptive version with the style of corresponding Micro data, as data pre-processing. Then a feature extractor is trained by all the training sets, where the contrastive learning strategy is employed to enhance the salient structure and eliminate the insignificant texture. Finally, with the learned feature extractor fixed, an MLP classifier is trained by the labeled data including adaptive Macro data and labeled Micro data.

A. Data Pre-processing
Real-world objects have clearer shapes and are easier to discriminate. However, if simply use full labeled Macro data and weakly labeled Micro data together in the training phase, the model may perform worse due to the different data distributions, known as "domain shift" [10].
Since the shape of Macro data is clearer, while Micro data accounting for a large proportion in quantity. We try to adapt Macro data to appear as if drawn from the Micro domain to realize the visual alignment, as shown in Fig. 2

(A).
In this subsection, we take rendering the ring image with Plasmodium style as an example. The low-level features of overall the labeled Plasmodium data should be separated as Plasmodium style, and the high-level content in terms of the object in the ring image should be extracted and recombined with Plasmodium style to produce an adaptive ring image.
Firstly a pre-trained CNN is used to extract and store the features. Content features are from the ring image (x s , 0), and style features are from all the labeled Plasmodium The features on CNN layer l could be stored in a matrix F l ∈ R N l ×H l ×W l , where N l is the number of distinct filters, H l and W l denotes the height and width of the feature map respectively. The style on CNN layer l could be represented by Gram matrix G l ∈ R N l ×N l , where G l,ij = F l,i F l,j is the inner product between the vectorized feature maps i and j in layer l [11].
Secondly, the overall loss function is minimized to update the image x a iteratively until it simultaneously matches the content features of x s and the style features of X t , where x a can be initiated as a white noise image and λ s is a coefficient. Since our goal is structure enhancement, which requires preserving the semantic content precisely, we empirically set a small weight of λ s = 10 −3 for this purpose. Fig. 3. Train the feature extractor e using a contrastive learning strategy. The network F , including the feature extractor e and a nonlinear projection head h, is updated by the back-propagation on a combination of supervised loss J S and unsupervised loss J U . An offline network Fm is used as a memory trick, which is an exponential moving average (EMA) of F , with a discrete dictionary Q on metric embeddings vm maintained as a queue. J S is defined by aligning the labeled Micro data with corresponding Macro data. J U is defined by aligning two views of different appearance transformation(at) of the same image, where at1 is random color distortion and at2 is Sobel flittering for example. For the learned representation z, J S contributes to the distinguishability and J U contributes to the appearance-invariance.
Based on the fact that deeper convolutional layers respond to higher semantics [12], J content and J style are defined as where L h are the higher CNN layers which capture the highlevel content in terms of objects, L o are all the other CNN layers, w l is the weight of layer l. Specifically, we useḠ l t = t G l t p as the average style features in the Micro domain instead of using the style of a single Micro data.
With the usage of labeled Micro data, all the Macro data X S = {(x s , y s ) S s=1 } is transferred to its adaptive version X A = {(x a , y s ) S s=1 } with corresponding Micro style.

B. Feature Extractor Learned by contrastive learning strategy
We aim to learn a distinguishable and appearance-invariant representation z for the downstream classification task using a small number of labeled Micro data X T = {(x t , y t )} T t=1 and a large amount of unlabeled Micro data X U = {x u } U u=1 , with the additional adaptive Macro data X A = {(x a , y s ) S s=1 }. We denote all the labeled data as X L = X T ∪ X A , and all the data except labels as X = X T ∪ X A ∪ X U . Fig.3 illustrates the proposed method to train the feature extractor e in Fig. 2(B). Our model denoted as F (x; Θ) : X → R |v| can be decomposed further into a CNN based feature extractor e(x; θ e ) : X → Z and a MLP projection head h(z; θ h ) : Z → V conceptually [13], where z is the representation for downstream tasks and v is the metric embedding for contrastive loss, with respect to an input x, and a(·) is a non-linear activation function using k-sparse strategy [14] a Γ = supp k {W T z} containing hidden units with top-k% activation values. In addition to the online network F , a offline momentum network F m (x; Θ m ) is used as a memory trick [15]. The update of Θ is by back-propagation on a combination of supervised loss J S and unsupervised loss J U while Θ m is an exponential moving average(EMA) of Θ where we follow the empirical experience and set α = 0.999. Considering the problem of similarity matching as a form of dictionary look-up, with similarity measured by the dot product of the metric embeddings, InfoNCE [16] is adopted as the form of contrastive loss function due to its efficiency and simplicity, with a well-grounded motivation from information theory. Given query v, the InfoNCE is the negative log-likelihood − log p v (6) where the likelihood is and σ is a temperature parameter that controls the concentration level of the distribution [17]. The value of InfoNCE is low when v is similar to its positive key v + , and dissimilar to all the other keys(considered as negative keys v − for v).
The way we define the positive keys and negative keys for all the training data realizes structure enhancement and texture elimination. 1) Structure Enhancement in J S : Instead of using a simple cross-entropy loss for labeled data X L , a contrastive loss is used to enhance the structure information by connecting Micro and Macro data directly. Every query embedding v t from Micro data, its positive keys are sampled from its corresponding adaptive Macro set, and its negative keys are the embeddings of all the other unrelated adaptive Macro data. Thus, Micro data is encouraged to be similar to its corresponding adaptive Macro data, while dissimilar to unrelated adaptive Macro data simultaneously. And it is the same way for adaptive Macro data. Therefore, the supervised loss for X L is and we define where t is the index of Micro data and Ω a are all the indexes of its corresponding adaptive Macro data, a is the index of adaptive Macro data and Ω t are all the indexes of its corresponding Micro data. Ω T are all the indexes of Micro data and Ω A are all the indexes of adaptive Macro data.
2) Texture Elimination in J U : To eliminate the insignificant texture, two kinds of appearance transformations are introduced, color distortion (which randomly changes the hue, lightness, and saturation of an image) and flittering (such as Sobel, Scharr, Laplacian). For every x u in X, it is randomly transformed to two different views x q , x m , and we define the positive pair p + (v q , v m ), naturally, all the other images are as the negative keys for x u .
Since InfoNCE benefits from more negative keys, we apply an offline momentum network F m (x; Θ m ) and maintain a discrete dictionary Q on metric embeddings v m as a queue [15]. The metric embeddings of the current minibatch are enqueued, and the oldest are dequeued. Introducing Q decouples the dictionary size from the mini-batch size, which can be much larger than the mini-batch size. Θ m is the exponential moving average copy from online network Θ which smooths the learning dynamics. The unsupervised loss for X is and we define where v q = F (x q ; Θ), v m = F m (x m ; Θ m ), and i is the index of the metric embeddings in the maintained queue under the current state.

C. MLP Classifier
All the labeled data X L is used to train a non-linear MLP classifier c in Fig. 2(C), which has a ReLU hidden layer, with the learned feature extractorê fixed.

A. Evaluation Metrics
We use accuracy (AC), F1-measure (F1) of sensitivity (SE) and specificity (SP), and Jaccard index (JA) to evaluate the performance AC = T P + T N T P + F P + T N + F N

B. Implementation
We adopt DenseNet-121 [18] architecture as the CNN backbone. In particular, we only include the five convolutional layers in the set, i.e., L = {conv1, dense1c, dense2d, dense3f, dense4c} for rendering the Macro data with Micro style, as the representations of these layers, in general, have the highest capability in each scale.
For a fair comparison, DenseNet-121 architecture is employed as the CNN backbone for all the methods, whose last fully-connected layer (after global average pooling) has a fixed-dimensional output (1024-D). The MLP projection head h leads to 128-D v with a 256-D hidden layer. The ReLU hidden layer of MLP classifier c is 256-D. In particular, we set temperature σ = 0.08 and k-sparse 20%. The batch size is 256, and the Queue size is 4096. The network was trained by SGD algorithm with a learning rate of 0.0001 and a momentum of 0.9. All the experiments are trained for 500 iterations and three repeat training are performed.

C. Result of Proposed Method
In this part, we report the performance of the supervised baseline and our method trained with only 50 randomly selected Macro data in each category, 200 labeled Micro data X T , and 20158 unlabeled Micro data X U in Table. I.  It is obvious that our semi-supervised method achieves higher performance on all the evaluation metrics, with 6.78%, 5.01%, and 11.76% improvements on AC, F1, and JA, respectively, compared with the supervised-only method. Our method, trained with 1% labeled Micro data and 10% Macro data, reaches an accuracy of 94.90% under a generalized testing set. It is comparable to the reported 95.7% (AC) in [7], which is trained in a full supervised scenario.
As illustrated in Fig. 4, a receiver operating characteristic (ROC) curve is created to visualize the classification performance of our proposed method for all four classes, and an overall area under the curve (AUC) value is computed to summarize the diagnostic performance. The proposed model achieves an overall AUC value of 0.94.

D. Effectiveness of Macro Data Introducing
We analyze how introducing Macro data helps the final performance, and the results are listed in Table. II.
The first three rows are the results of the supervised method, it helps a lot when introducing Macro data even in its original version. However, the results of the semi-supervision scenario in the last three rows tell a different story. There is a tiny improvement when introducing original Macro data, while an impressive improvement when the Macro data is transformed to its adaptive version. These results conform to our assumption that 'domain shift' is remarkable in a semi-supervision scenario since there are a large amount of unlabeled data being trained together with the labeled data.

E. Comparison With Other Semi-supervised Methods
As it is much easier to collect real-world objects naturally labeled than annotating microscopic parasites, we analyze the performances with the different numbers of Macro data, where two state-of-the-art semi-supervised learning methods are compared as counterparts.
Mean Teacher (MT) [19] used the same EMA strategy as our proposed method from a student network to a teacher network, but simply encourage consistent representations from the same input, without any comparison with others. With the same goal of consistency regularization, Virtual Adversarial Training (VAT) [20] directly approximated a tiny perturbation to add to input which would most significantly affect the output of the prediction function, instead of using data augmentation. Quantitative results are in Table. III and we draw the AC score of the results in Fig. 5. All these three semi-supervised methods perform better than the supervised baseline consistently, demonstrating that all the SSL methods effectively utilize the unlabeled data and bring performance gains. As expected, the performances of all the methods increase when more labeled training data are available, and the gap between supervised baseline and semi-supervised methods narrows as more labeled training data are available.
Compared with other methods, we achieve the greatest improvement over the supervised baseline, especially when there are only 200 labels. The comparison shows the effectiveness of our proposed method, compared with other semisupervised methods. And the success is naturally attributed to contrastive learning strategy.

V. CONCLUSION AND FUTURE WORK
In this work, we propose a semi-supervised method using a contrastive learning strategy to classify three apicomplexan parasites and non-infected host cell. The learned representation is more distinguishable by the alignment between Micro data and Macro data, while it is appearance-invariant by encouraging the consistent representation between two different appearance transformed views of the same image. The results verify the effectiveness of the proposed method. At present, the classifier is relatively simple, it worth further exploring how to take full advantage of the learned representation for classification.