Privacy Preserving Image Registration

Image registration is a key task in medical imaging applications, allowing to represent medical images in a common spatial reference frame. Current approaches to image registration are generally based on the assumption that the content of the images is usually accessible in clear form, from which the spatial transformation is subsequently estimated. This common assumption may not be met in practical applications, since the sensitive nature of medical images may ultimately require their analysis under privacy constraints, preventing to openly share the image content.In this work, we formulate the problem of image registration under a privacy preserving regime, where images are assumed to be confidential and cannot be disclosed in clear. We derive our privacy preserving image registration framework by extending classical registration paradigms to account for advanced cryptographic tools, such as secure multi-party computation and homomorphic encryption, that enable the execution of operations without leaking the underlying data. To overcome the problem of performance and scalability of cryptographic tools in high dimensions, we propose several techniques to optimize the image registration operations by using gradient approximations, and by revisiting the use of homomorphic encryption trough packing, to allow the efficient encryption and multiplication of large matrices. We demonstrate our privacy preserving framework in linear and non-linear registration problems, evaluating its accuracy and scalability with respect to standard, non-private counterparts. Our results show that privacy preserving image registration is feasible and can be adopted in sensitive medical imaging applications.


Introduction
Image Registration is a crucial task in medical imaging applications, allowing to spatially align imaging features between two or multiple scans.Registration methods are today a central component of state-of-the-art methods for atlas-based segmentation (Shattuck et al., 2009;Cardoso et al., 2013), morphological and functional analysis (Dale et al., 1999;Ashburner and medical imaging software and applications, including public web-based services for automated segmentation and labelling of medical images.Using these services generally requires uploading and exchanging medical images over the Internet, to subsequently perform image registration with respect to one or multiple (potentially proprietary) atlases.Besides these classical medical imaging use-cases, emerging paradigms for collaborative data analysis, such as Federated Learning (FL) (McMahan et al., 2017), have been proposed to enable analysis of medical images in multicentric scenarios for performing group analysis (Gazula et al., 2021) and distributed machine learning (Kaissis et al., 2021;Zerka et al., 2020).However, in these settings, typical medical imaging tasks such as spatial alignment and downstream operations are generally not possible without disclosing the image information.
Due to the evolving juridical landscape on data protection, medical image analysis tools need to be adapted to guarantee compliance with regulations currently existing in many countries, such as the European General Data Protection Regulation (GDPR)1 , or the US Health Insurance Portability and Accountability Act (HIPAA)2 .Medical imaging information falls within the realm of personal health data (Lotan et al., 2020) and its sensitive nature should ultimately require the analysis under privacy preserving constraints, for instance by preventing to share the image content in clear form.
Advanced cryptographic tools hold great potential in sensitive data analysis problems (e.g., (Lauter, 2021)).Examples of such approaches are Secure-Multi-Party-Computation (MPC) (Yao, 1982) and Homomorphic Encryption (HE) (Rivest et al., 1978).While MPC allows multiple parties to jointly compute a common function over their private inputs and discover no more than the output of this function, HE enables computation on encrypted data without disclosing either the input data or the result of the computation.
This work presents privacy-preserving image registration (PPIR), a new methodological framework allowing image registration under privacy constraints.To this end, we reformulate the image registration problem to integrate cryptographic tools, namely MPC or FHE, thus preserving the privacy of the image data.Due to the well-known scalability issues of such cryptographic techniques, we investigate strategies for the practical use of PPIR.In our experiments, we evaluate the effectiveness of PPIR on a variety of registration tasks and medical imaging modalities.Our results demonstrate the feasibility of PPIR and pave the way for the application of secured image registration in sensitive medical imaging applications.

Background
Given images I, J : R d → R, image registration (IR) aims at estimating the parameters θ of a spatial transformation W θ ∈ R d → R d , either linear or non-linear, maximizing the spatial until ∥∆θ∥ ≤ ϵ 10: return I(W θ ), e 11: end function overlap between J and the transformed image I(W θ ), by minimizing a registration loss function f : (1) The loss f can be any similarity measure, e.g., the Sum of Squared Differences (SSD), the negative Mutual Information (MI), or normalized cross correlation (CC).Equation (1) can be typically optimized through gradient-based methods, where the parameters θ are iteratively updated until convergence.In particular, when using a Gauss-Newton optimization scheme (Algorithm 1), the update of the spatial transformation can be computed through Equation ( 2): where G = ∂ f ∂θ is the Jacobian and H = ∂ 2 f ∂θ 2 the Hessian of f .Besides the Gauss-Newton schemes proposed in the field of IR (Pennec et al., 1999;Modersitzki, 2009), gradient-based techniques are classically adopted to solve the IR task, for example in diffeomorphic image registration problems (Ashburner, 2007;Avants et al., 2011).
In all these cases we consider a scenario with two parties, party 1 , and party 2 , whereby party 1 owns image I and party 2 owns image J.The parties wish to collaboratively optimize the image registration problem without disclosing their respective images to each other.We assume that only party 1 has access to the transformation parameters θ and that it is also responsible for computing the update at each optimization step.In what follows, we introduce the basic notation to develop PPIR based on different registration frameworks.We focus on registration methods of increasing complexity, including (i) rigid, (ii) affine, and (iii) non-linear registration based on cubic splines or diffeomorphisms parametrized by time-varying velocity fields (large deformation diffeomorphic metric mapping, LDDMM) (Beg et al., 2005).In all these settings, we demonstrate how the registration problem can be naturally adapted for accounting to privacy-preserving operations, and illustrate the effectiveness of PPIR on a variety of registration tasks.

Analysis of classical IR loss functions under a privacypreserving perspective
In this section we review typical loss functions used in image registration, and analyze the related requirements for privacypreserving optimization.

Optimization of SSD loss
A typical loss function to be optimized during the registration process is the sum of squared intensity differences (SSD) evaluated on the set of image coordinates: with Jacobian: where the quantity quantifies image and transformation gradients, and is the second order term obtained from Equation (3) through linearization (Pennec et al., 1999;Baker and Matthews, 2004).
The solution to this problem requires the joint availability of both images I and J, as well as of the gradients of I and of W θ .In a privacy-preserving setting, this information cannot be disclosed, and the computation of Equation ( 2) is therefore impossible.We note that to calculate the update of the registration ∆θ of Equation ( 2), the only operation that requires the joint availability of information from both parties is the term R = x S (x) • J(x), which can be computed a matrix-vector multiplication of vectorized quantities R = S T • J.

Mutual Information
Mutual Information quantifies the joint information content between the intensity distributions of the two images.This is calculated from the joint probability distribution function (PDF): where, given N r and N t the maximum intensity for respectively I and J, we define r as the range of discretized intensity values of I and J, respectively.A Parzen window (Parzen, 1961) is used to generate continuous estimates of the underlying intensity distributions, thereby reducing the effects of quantization from interpolation, and discretization from binning the data.Let ψ 3 I : R → [0, 1] be a cubic spline Parzen window, and let ψ 0 J : R → [0, 1] be a zeroorder spline Parzen window.The smoothed joint histogram of I and J (Viola and Wells III, 1997;Mattes et al., 2003) is given by: In the above formula, the intensity values of I and J are normalized by their respective minimum (denoted by I(W θ ) • and J • ), and by the bin size (respectively ∆b r and ∆b t ), to fit into the specified number of bins (b r or b t ) of the intensity distribution.The final value for p(r, t; θ) is computed by normalizing by N x , the number of sampled voxels.Marginal probabilities are simply obtained by summing along one axis of the PDF, that is, p(r) = t p(r, t; θ) and p(t) = r p(r, t; θ).Let the matrices A 3 I ∈ R N x ×N r and B 0 J ∈ R N x ×N t be defined as: the discretized joint PDF can be rewritten in a matrix form via the multiplication: The first derivative of the joint PDF is calculated as follows (Dowson and Bowden, 2006): is the input of the cubic spline.We also introduce the tensor C 3 I ∈ R N x ×N r ×|θ| defined as ∂θ , to write the discretized first derivative as: where P ′ ∈ R N t ×N r ×|θ| .The Jacobian of the MI is obtained from the chain rule and takes the form: while the linearized Hessian (Dowson and Bowden, 2007) can be written as: where P I = t p(r, t, θ) is a vector which defines the discretized marginal PDF of the moving image.The derivatives can be easily calculated from the properties of B-splines since we have . In a privacy-preserving scenario, to calculate the update of the registration ∆θ of Equation ( 7), two operations require the joint availability of information from both parties, which are the matrix P of Equation ( 8) and the matrix P ′ of Equation (9).

Cross Correlation with Advanced Normalization Tools
In the Advanced Normalization Tools (ANTs) introduced by (Avants et al., 2008), the normalized cross-correlation (CC) loss was specified in the context of diffeomorphic image registration.Let ϕ define a diffeomorphism over the domain Ω = R d , parameterized by a time-varying velocity field v(x, t).In the ANTs setting, inverse consistency is obtained by optimizing the CC loss with respect to both forward and backward (inverse) transformations, here denoted by ϕ 1 (x, t) and ϕ 2 (x, t), and parameterized by velocity fields v 1 (x, t) and v 2 (x, t) respectively.In particular, both images I and J are simultaneously warped towards a "half-way" space, to obtain I 1 = I(ϕ 1 (x, 0.5)) and J 2 = J(ϕ 2 (x, 0.5)).The CC loss is thus defined as: quantify the images appearance at location x i , with respect to the average intensity µ I 1 and µ J 2 measuread in a local window of size M. Coherently with the LDDMM formulation, the variational optimization problem is defined as: where L is a linear operator prescribing a norm on the velocity fields acting as a regularizer.The equation for the derivative of the forward update is given by: while the derivative of the backward update is analogously given by: In privacy-preserving scenario, the sensitive terms carrying private image information are 2D EF , ( J(x) , which must therefore be computed in a privacypreserving regime.

Building blocks for Secure Computation
After introducing in Section 1 the IR optimization problem and the related functionals addressed in this work, in this section, we review the standard privacy-preserving techniques that will be employed to develop PPIR.

Secure Multi-Party Computation
Introduced by (Yao, 1982), MPC is a cryptographic tool that allows multiple parties to jointly compute a common function over their private inputs (secrets) and discover no more than the output of this function.Among existing MPC protocols, additive secret sharing consists of first splitting every secret s into additive shares ⟨s⟩ i , such that n i=1 ⟨s⟩ i = s, where n is the number of collaborating parties.Each party i receives one share ⟨s⟩ i and executes an arithmetic circuit in order to obtain the final output of the function.In this paper, we adopt the twoparty computation protocol defined in SPDZ (Damgard et al., 2011), whereby the actual function is mapped into an arithmetic circuit and all computations are performed within a finite ring with modulus Q. Additions consist of locally adding shares of secrets, while multiplications require interaction between parties.Following (Damgard et al., 2011), SPDZ defines: MPC-Mul to compute element-wise multiplication, MPCDot to compute matrix-vector multiplication and MPCMatMul to compute matrix-matrix multiplication.These operations are performed within an honest but curious protocol.

Homomorphic Encryption
Initially introduced by Rivest et al. in (Rivest et al., 1978), Homomorphic Encryption (HE) enables the execution of operations over encrypted data without disclosing either the input data or the result of the computation.Hence, party 1 encrypts the input with its public key and sends the encrypted input to party 2 .In turn, party 2 evaluates a circuit over this encrypted input and sends the result, which still remains encrypted, back to party 1 which can finally decrypt them.Among various HE schemes, CKKS (Cheon et al., 2017) supports the execution of all operations on encrypted real values and is considered a levelled homomorphic encryption (LHE) scheme.The supported operations are: Sum (+), Element-wise multiplication ( * ) and DotProduct.With CKKS, an input vector is mapped to a polynomial and further encrypted with a public key in order to obtain a pair of polynomials c = (c 0 , c 1 ).The original function is further mapped into a set of operations that are supported by CKKS, which are executed over c.The performance and security of CKKS depend on multiple parameters including the degree of the polynomial N, which is usually sufficiently large (e.g.N = 4096, or N = 8192).

Methods: from IR to PPIR
In this section, we describe the privacy preserving variants of the three IR methods described in Section 3. We propose two versions of PPIR for SSD according to the underlying cryptographic tool, namely PPIR(MPC), integrating MPC, and PPIR(FHE), integrating FHE.In the case of MI and NCC, we focus on the design of the MPC-variant, due to the nonnegligible computational overhead of FHE in these applications.Finally, we also study rigid point cloud registration and describe its privacy-preserving variant in Appendix A 1.

PPIR based on SSD
As mentioned in Section 3.1, when optimizing the SSD cost, the only sensitive operation that must be jointly executed by the parties is the matrix-vector multiplication: R = S T • J, where S T is only known to party 1 and J to party 2 .Figure 1 illustrates how cryptographic tools are employed to ensure privacy during registration.
With MPC (Figure 1a), party 1 secretly shares the matrix S T to obtain (⟨S ⟩ 1 , ⟨S ⟩ 2 ), while party 2 secretly shares the image J to obtain (⟨J⟩ 1 , ⟨J⟩ 2 ).Each party also receives its corresponding share, so that party 1 holds (⟨S ⟩ 1 , ⟨J⟩ 1 ) and party 2 holds (⟨S ⟩ 2 , ⟨J⟩ 2 ).The parties execute a circuit with MPCMul operations to calculate the 2-party dot product between S T and J.The parties further synchronize to allow party 1 to obtain the product and finally to calculate ∆θ (Equations ( 4) and ( 6)).
When using FHE (Figure 1b), party 2 first uses a FHE key k to encrypt J and obtains ⟦J⟧ ← Enc(k, J).This encrypted image is sent to party 1 , who computes the encrypted matrix-vector multiplication ⟦R⟧.In this framework, only vector J is encrypted, and therefore party 1 executes scalar multiplications and additions in the encrypted domain only (which are less costly than multiplications over two encrypted inputs).The encrypted result ⟦R⟧ is sent back to party 2 , which can obtain the result by decryption: R = Dec(k, ⟦R⟧).Finally, party 1 receives R in clear form and can therefore compute ∆θ.
Thanks to the privacy and security guarantees of these cryptographic tools, during the entire registration procedure, the content of the image data S and J is never disclosed to the opposite party.

PPIR based on MI
In the MPC variant to calculate the joint PDF P in a privacypreserving manner, the quantity A 3 I is only known to party 1 , while the quantity B 0 J to party 2 (Section 3.2).With reference to Supplementary Figure A1a, party 1 secretly shares matrix (A 3 I ) T = (⟨A 3 I ⟩ 1 , ⟨A 3 I ⟩ 2 ), while party 2 does the same with . Each party also receives its corresponding share: party 1 now holds (⟨A 3 I ⟩ 1 , ⟨B 0 J ⟩ 1 ), and party 2 holds (⟨B 0 J ⟩ 2 , ⟨A 3 I ⟩ 2 ).The parties execute a circuit with MPCMatMul operation to calculate the 2-party matrix multiplication between (A 3 I ) T and B 0 J .The next operation carried out in the privacypreserving setting is the computation of the first derivative of Equation ( 9).Supplementary Figure A1b illustrates the MPC variant, where party 1 only knows C 3 I and party 2 only knows B 0 J .Initially, party 1 secretly shares the matrix . Each party also receives its corresponding share, namely: party 1 holds (⟨C 3 I ⟩ 1 , ⟨B 0 J ⟩ 1 ) and party 2 holds (⟨C 3 I ⟩ 2 , ⟨B 0 J ⟩ 2 ).The parties also execute a circuit with the MPCMatMul between (B 0 J ) T and C 3 I .

PPIR based on ANTS CC loss
According to Section 3.3, party 1 has access to Ī and E, whereas party 2 has access to J and F. The computation begins with the optimization of the CC term, specifically the quantity  A2a, party 1 and party 2 secretly share Ī = (⟨ Ī⟩ 1 , ⟨ Ī⟩ 2 ), and J = (⟨ J⟩ 1 , ⟨ J⟩ 2 ), respectively.The subsequent multiplication of these shared values results in the computation of the shares of D = (⟨D⟩ 1 , ⟨D⟩ 2 ), which is never reconstructed.In the next step of the protocol, party 1 secretly shares 1 E , and party 2 secretly shares 1 F .Through a multiplication of the shares of D, 1 E , and 1 F , both parties collectively obtain the final value of 2D EF .The other two terms that need to be jointly computed are the third term of Equation ( 10) and ( 11), namely ( J − D E Ī) and ( Ī − D F J), reported in Supplementary Figure A2b.In the case of ( J − D E Ī), party 1 secretly shares Ī and 1 E , while party 2 secretly shares J. Since both parties already have access to the share of D from the previous protocol, they proceed to multiply the secret shares of D, 1 E , and Ī.Subsequently, they subtract the result from J to obtain the final value of ( Ī − D F J).The computation of ( J − D E Ī) is analogous to the one of ( Ī − D F J), where party 2 secretly shares 1 E , and both parties participate to the multiplication of D, 1 E , and J.The resulting value is finally subtracted from Ī, yielding the final result ( J − D E Ī).

Protocols enhancement for SSD loss
Effectively optimizing Equation (1) with MPC or FHE is particularly challenging, due to the computational bottleneck of these techniques when applied to large-dimensional objects (Haralampieva et al., 2020;Benaissa et al., 2021), notably affecting the computation time and the occupation of communication bandwidth between parties.Because cryptographic tools introduce a non-negligible overhead in terms of performance and scalability, in this section we introduce specific techniques to optimize the underlying image registration operations.

Gradient sampling
Since the registration gradient is generally driven mainly by a fraction of the image content, such as the image boundaries in the case of SSD cost, a reasonable approximation of Equations (4) and ( 6) can be obtained by evaluating the cost only on relevant image locations.This idea has been introduced in medical image registration (Viola and Wells III, 1997;Mattes et al., 2003;Sabuncu and Ramadge, 2004), and here is adopted to optimize Equation (3) by reducing the dimensionality of the arrays on which encryption is performed.We test two different techniques: (i): Uniformly Random Selection (URS), proposed by (Viola and Wells III, 1997;Mattes et al., 2003), in which a random subset of dimension l ≤ d of spatial coordinates is sampled at every iteration with uniform probabilities, p(x) = 1 d ; and (ii): Gradient Magnitude Sampling (GMS) (Sabuncu and Ramadge, 2004), which consists of sampling a subset of coordinates with probability proportional to the norm of the image gradient, p(x) ∼ ∥∇I(x)∥.We note that gradient sampling is not necessary for computing the MI since in Equation ( 8) the computation is already performed on a subsample of the image voxels.

Matrix partitioning in FHE
We now describe additional improvements dedicated to PPIR(FHE) and propose two versions of this solution: (i) PPIR(FHE)-v1 implements an optimization of the matrixvector multiplication by partitioning the vector image into a vector of submatrices, whereas (ii) PPIR(FHE)-v2 enhances the workload of the parties by fairly distributing the computation among them.
PPIR(FHE)-v1.We introduce here a novel optimization dedicated to PPIR with FHE, in particular when the CKKS algorithm is adopted.CKKS allows multiple inputs to be packed into a single ciphertext to decrease the number of homomorphic operations.To optimize matrix-vector multiplication, we propose to partition the image vector J into k sub-arrays of dimension l, and the matrix S T into k sub-matrices of dimension |θ| × l.Once all sub-arrays J i are encrypted, we propose to iteratively apply DotProduct as proposed by Benaissa et al. (2021), between each sub-matrix and corresponding sub-array; these intermediate results are then summed to obtain the final result, namely: PPIR(FHE)-v2.In addition to packing multiple inputs into a single ciphertext, the application of FHE to PPIR can be optimized by more equally distributing the workload among the two parties.We note that in PPIR(FHE)-v1, party 1 is in charge of computing the matrix-vector multiplication entirely while party 2 only encrypts the input and decrypts the result.Following Supplementary Figure A4, PPIR(FHE)-v2 starts by splitting the matrix S of party 1 into two sub-matrices S 1 and S 2 using the operation split h,K .This operation partitions the matrix into K equally-sized sub-matrices.Next, the operation f latten is applied to S 1 and S 2 , obtaining vectors S ′ 1 and S ′ 2 respectively.Then, party 1 encrypts S ′ 2 and sends it to party 2 , which subsequently applies to its vector J the operation split v,K , obtaining J 1 and J 2 .This operation splits J into K equally-sized partitions, and it executes the operation replicate d to J 1 and J 2 , obtaining J ′ 1 and J ′ 2 respectively.party 2 encrypts J ′ 1 and sends it to party 1 .Both parties then iteratively perform the element-wise multiplication * and sum up the results of the different partitions using the primitive S um i,k .The protocol doesn't rely anymore on DotProduct and leads to a significant gain in computational load.

Experiments & Results
We demonstrate and assess the different versions of PPIR illustrated in Section 3 on a variety of image registration problem, namely: (i) SSD for rigid transformation of point cloud data, (ii) SSD with linear and non-linear alignment of whole body positron emission tomography (PET) data; (iii) SSD and MI for mono-and multimodal linear alignment of MRI and PET brain scans; (iv) diffeomorphic non-linear registration with CC of multimodal abdomen data from CT and MRI scans.Experiments are carried out on 2D (mainly for the SSD case) and 3D imaging data.In Table 1 is reported an overview of the datasets used specifying their dimensions, modality, registration type, loss functions, and the Privacy Enhancing Technologies (PETs) employed.

Experimental data
Point Cloud Data.We showcase rigid registration on 2D point cloud data representing the corpus callosum, as presented in Vachet et al. (2012), with a set size n = 193.The registration loss here considered is SSD between point coordinates (additional details are provided in Appendix A 1).
Whole body PET data.The dataset considered for linear and non-linear registration with SSD consists of 18-Fluoro-Deoxy-Glucose ( 18 FDG) whole body PET scans.The images are a frontal view of the maximum intensity projection reconstruction, obtained by 2D projection of the voxels with the highest intensity across views (1260 × 1090 pixels).
Brain MRI and PET data.This dataset regroups brain MRI and PET images obtained from the Alzheimer's Disease Neuroimaging Initiative (Mueller et al., 2005).MRI data were processed via a standard processing pipeline to estimate gray matter density maps (Ashburner and Friston, 2000).Non-linear registration was carried out on the extracted mid-coronal slice, of dimension 121 × 121 pixels.For 3D multimodal linear registration with MI, we use both MRI images and PET images, with respective dimension of 180 × 256 × 256 and 160 × 160 × 96 voxels.
Abdomen MR and CT data.The multimodal dataset Abdomen-MR-CT (Hering et al., 2022) was used for experiments with ANTs registration based on CC.The data was compiled from public studies of the cancer imaging archive (TCIA) Clark et al. (2013) that contains 8 paired scans of MRI and CT from the same patients.The data have an isotropic resolution of 2mm and a voxel dimension of 192 × 160 × 192.They also provide 3D segmentation masks for the liver, spleen, and left and right kidney.All scans were pre-aligned by groupwise affine registration.

Experimental Details
In order to avoid local minima and to decrease computation time, we use a hierarchical multiresolution optimization scheme.The scheme involves M resolution steps, denoted as r 1 . . .r M .At each resolution step r m , the input data is downsampled by a scaling factor m, where m ∈ [1 . . .M].The quality of PPIR is assessed by comparing the registration results with those obtained with standard registration on clear images (Clear).The metrics considered are the difference in image intensity at optimum (for SSD), the total number of iterations required to converge, and the displacement root mean square difference (RMSE) between Clear and PPIR.We also evaluate the performance of PPIR in terms of average computation (running time) and communication (bandwidth) across iterations.In the multimodal Abdomen data, the quality of the registration result was assessed by the overlap across the labeled anatomical regions, quantified by the DICE score.
Point Cloud Data.The registration protocol here adopted is detailed in Appendix A 1. For MPC we set as the prime modulus Q = 2 32 .For PPIR(FHE), we define the polynomial degree modulus as N = 4096.
Whole body PET data.Whole-body PET image alignment was first performed by optimizing the transformation W θ in Equation ( 1) with respect to affine registration parameters.The multiresolution steps used are r 1 , r 5 , r 10 , r 20 .A second wholebody PET image alignment experiment was performed by nonlinear registration, without gradient approximation based on a cubic spline model (one control point every four pixels along both dimensions), with multiresolution steps r 1 ,r 2 ,r 5 ,r 10 ,r 20 and r 30 .Concerning the PPIR framework, transformations were optimized for both MPC and FHE by using gradient approximation (Section 5) using the same sampling seed for each test.For MPC we set as the prime modulus Q = 2 32 .For PPIR(FHE), we define the polynomial degree modulus as N = 4096, and set the resizing parameter D to optimize the trade-off between run-time and bandwidth.Since D needs to be a divisor of the image size image data we set D = 128.
Brain MRI and PET data.The registration of brain gray matter density images was performed by non-linear registration based on SSD, without gradient approximation, based on a cubic spline model (one control point every five pixels along both dimensions), with multiresolution steps r 1 and r 2 .For PPIR(MPC) we use the same configuration defined in the previous section, while for PPIR(FHE) we use the same N and we set D = 121.
We tested PPIR with MI for multi-modal 3D affine image registration between PET and MRI brain scans where, in addition to varying the multi-resolution steps, a Gaussian blurring filter is applied to the images with a kernel that narrows as multi-resolution proceeds.The kernel size at different resolutions, denoted with σ 1 . . .σ M , is used to control the amount of blurring applied to the image at each step of the multi-resolution process.The multiresolution steps applied are r 5 and r 10 , with 10% of the image's pixels utilized as the number of subsample pixels (N x ).Gaussian image blurring is applied with a degree of σ 5 = 1 and σ 10 = 3.For MPC we set as the prime modulus Q = 2 64 .
Abdomen MRI and CT data.We tested PPIR with CC for multi-modal 3D ANTs image registration between MRI and CT abdomen scans.The multiresolution steps applied are r 3 , r 2 and r 1 and the CC window size M = 5 × 5 × 5. and TenSeal (Benaissa et al., 2021), which implements the CKKS protocol3 .PPIR based on MI and CC is implemented by extending the Dipy framework of Garyfallidis et al. (2014) 4 .Finally, PPIR for point cloud data is released in a separated repository 5 .
All the experiments are executed on a machine with an Intel(R) Core(TM) i7-7800X CPU @ (3.50GHz x 12) using 132GB of RAM.For each registration configuration, the optimization is repeated 10 times to account for the random generation of MPC shares and FHE encryption keys.

Results
Point Cloud Data.In Supplementary Table A1 we present the registration metrics for PPIR(MPC) and PPIR(FHE)-v1.The registration shows that PPIR(MPC) achieves the best results compared to PPIR(FHE), which exhibits not only a longer computation time but also requires higher bandwidth, thanks to its non-iterative algorithm.However, to carry out MatMul, a sufficiently large N (4096) is required, and in this scenario, it leads to a significant loss of chipertext slots compared to the dimension of the point set n = 193.Finally, the qualitative results reported in Figure A7 show negligible differences between point cloud transformed with Clear, PPIR(MPC) and PPIR(FHE)-v1.
Whole body PET data: affine registration (SSD).Table 2 compares Clear, PPIR(MPC), PPIR(FHE)-v1 and v2, showcasing metrics resulting from the affine transformation of whole-body PET images.Notably, registration through PPIR(MPC) yields negligible differences compared to Clear in terms of the number of iterations, intensity, and displacement.In contrast, registering with PPIR(FHE) is not feasible when considering entire images due to computational complexity.Nevertheless, Supplementary Figure A5 shows that neither MPC nor FHE decreases the overall quality of the affine registered images.A comprehensive assessment of the registration results is available in the Appendix.Table 2 (Efficiency metrics) shows that PPIR(MPC) performed on full images requires higher computation time and required communication bandwidth compared to PPIR(MPC)+URS/GMS.These numbers improve sensibly when using URS or GMS (by factors 10× and 20× for time and bandwidth, respectively).Concerning PPIR(FHE)-v1, we note the uneven computation time and bandwidth usage between clients, due to the asymmetry of the encryption operations and communication protocol (Figure 1).PPIR(FHE)-v2, which shares the computational workload between the two parties and avoids DotProduct, allows obtaining an important speed-up over PPIR(FHE)-v1 (100× faster).Notably, this gain is obtained without affecting the quality of the registration metrics, and improves the execution time of PPIR(MPC).On the other side, although PPIR(FHE)-v1 is able to improve communication with respect to PPIR(MPC), it still suffers from the highest communication among the three proposed solutions.This is due to the fact PPIR(FHE) protocols can find different applications depending on the requirements in terms of computational power or bandwidth.
Brain MRI data and whole body PET data: nonlinear registration (SSD).Fig. 2: Qualitative results for affine registration with MI over 3D medical images using ADNI dataset (Mueller et al., 2005).The images are presented in a 3 × 4 grid, with the first row representing the axial axis, the second row the coronal axis, and the third row the sagittal axis.In the first column of each row, the moving image obtained using PET modality is shown, while in the second column, the fixed image obtained using MRI modality is displayed.The third column shows the checkerboard alignment result using Clear, while the fourth column shows the result using PPIR(MPC).The different protocols are highlighted by red and green frames, respectively.
matter density images without the application of gradient approximation.Additionally, the table includes results for the registration between whole-body PET images when the gradient approximation is applied.
Brain MRI data without gradient approximation.Regarding the registration accuracy, we draw conclusions similar to those of the affine case, where PPIR(MPC) leads to minimum differences with respect to Clear, while PPIR(FHE)-v1 seems slightly superior.PPIR(MPC) is associated with a lower execution time and a higher computational bandwidth, due to the larger number of parameters of the cubic splines, which affects the size of the matrix S .Although PPIR(FHE)-v1 has a slower execution time, the demanded bandwidth is inferior to the one of PPIR(MPC), since the encrypted image is transmitted only once.PPIR(FHE)-v2, as in the affine case, outperforms PPIR(FHE)-v1 (still about 100× faster) leading to comparable values between registration metrics and is still inferior to PPIR(MPC).Here, the limitations of PPIR(FHE)-v1 on the bandwidth size are even more evident than in the affine case, since the bandwidth increases according to the number of parameters.This result gives a non-negligible burden to the party 1 , due to the multiple sending of the flattened and encrypted submatrices of updated parameters.Furthermore, in this case, PPIR(FHE)-v1 performs slightly worse than PPIR(MPC) in terms of execution time.
Whole Body PET Data With Gradient Approximation.Incorporating gradient approximation for handling whole-body PET data leads to similar conclusions as for the experiments on brain data.Qualitative results, reported in Supplementary Figure A6, show negligible differences between images transformed with Clear+GMS, PPIR(MPC)+GMS, and PPIR(FHE)-v1+GMS.
Brain MRI and PET data: affine registration (MI).Table 4 provides information on the computation cost of the protocol and the registration metrics for both the joint PDF and the First Derivative of the Joint PDF, using both Clear and PPIR(MPC).The results demonstrate a reasonable execution time (completed in under 5 minutes for the entire process) and noteworthy data transfer, totaling less than 4GB.Qualitative results for the image registration are shown in Figure 2, indicating that there is no difference between the moving transformed using Clear and PPIR(MPC).
Abdomen MRI and CT data: diffeomorphic non-linear registration (CC).Table 5 presents the metrics for ANTs registration using both Clear and PPIR(MPC) methods.The initial DICE scores for both forward and backward transformations are consistent between Clear and PPIR(MPC), with slight variations in the final DICE scores.The number of iterations and displacement RMSE values also exhibit similar trends.In terms of communication and computation times, PPIR(MPC) demonstrates comparable performance with Clear, showcasing its feasibility for secure registration.The computation times and communication bandwith between the two parties are well within reasonable limits.These qualitative result in Figure 3 indicates that the proposed PPIR(MPC) solution maintains the quality of registration metrics while ensuring secure and private communication between parties.

Discussion
This current work builds upon Taiello et al. (2022) by extending its application to include MI with linear registration (3D images), CC using the ANTs framework (3D images), enhancing the FHE for the SSD loss function, namely PPIR(FHE)-v2, and also integrating rigid point cloud registration.We recognize specific limitations associated with the use of FHE in our PPIR framework, which limited the effective use of this techniques besides the optimization of the SSD loss.FHE's computational cost during homomorphic operations poses challenges, limiting the scalability and real-time applicability of PPIR(FHE).This is particularly evident when dealing with large datasets, such as 3D images, or when employing advanced image registration cost functions that demand significant computational resources.To address this gap, researchers are actively exploring the optimization of FHE through the integration of hardware accelerators (Boemer et al., 2021).
For the SSD loss function, we provide comparison experiments with both URS and GMS (Viola and Wells III, 1997;Mattes et al., 2003;Sabuncu and Ramadge, 2004) for sake of completeness and compatibility with subsampling approaches in IR.We recognized that URS doesn't bring substantial improvements with resepct to GRS, and this latter method should be preferred in the considered application or testing scenario.
While PPIR focuses on the privacy-preserving formulation of classical image registration methods based on gradient-based optimization, throughout the past years the research community has been steering the attention towards deep learning (DL)based image registration (Simonovsky et al., 2016;Krebs et al., 2017;Yang et al., 2017;Balakrishnan et al., 2019).Among the medical imaging application of privacy-preserving methodologies, Kaissis et al. (2021) discussed privacy-preserving FL with Secure Aggregation (Bonawitz et al., 2017) and Differential Privacy (Abadi et al., 2016) for 2D medical image classification tasks.However, as highlighted by (Kaissis et al., 2021), deploying DL models for privacy-preserving inference nowadays is predominantly achievable through Multi-Party Computation (MPC).This process necessitates multiple servers and incurs significant overhead, primarily attributed to the size of the DL model, especially when handling 3D image registration tasks within a DL-based framework (Balakrishnan et al., 2019).To the best of our knowledge both DL and non-DL registration methods available in the literature do not satisfy the PPIR requirements investigated in our work, as they always require the disclosure of the target and moving images in clear.

Conclusion and future works
This study introduces the novel paradigm of Privacy Preserving Image Registration, designed for allowing image registration in privacy-preserving scenarios where images are confidential and cannot be shared in clear.Leveraging both secure multi-party computation (MPC) and Fully Homomorphic Encryption (FHE), we propose in PPIR effective strategies integrating cryptographic techniques into a variety of state-of-theart registration frameworks, encompassing different parameterization and loss functions.We evaluate the framework's per-formance across various registration benchmarks, conducting quantitative and qualitative assessment for all the considered image registration problems.Our future direction involve extending PPIR to encompass additional cost functions commonly used in image registration, aiming to enhance the framework's versatility and applicability.Fig. 3: Qualitative results for diffeomorphic registration with CC between 3D medical images from the AbdomenMRCT dataset (Hering et al., 2022).The images are presented in a 3 × 4 grid, with the first row representing the axial axis, the second row the coronal axis, and the third row the sagittal axis.First and second column show respectively MRI and CT images.The third column shows the MRI transformed using Clear, while the fourth column shows the MRI transformed using PPIR(MPC).The transformed images are highlighted by red and green frames, respectively.
Fig. A2: Optimization of ANTS NCC loss: proposed framework to calculate 2D EF and ( J − D E Ī) based on PPIR(MPC).
Fig. A3: Optimization of rigid point cloud: proposed framework to compute matrix multiplication L = Q • Q ′T based on PPIR(MPC) and PPIR(FHE).

Fig. A5 :
Fig. A5: Qualitative results for affine registration with SSD between 2D medical images.The red frame is the transformed moving image using Clear+URS registration.Green and Yellow frames are the transformed images using respectively PPIR(MPC)+URS and PPIR(FHE)v1+URS.

Fig. A6 :
Fig. A6: Qualitative results for Cubic splines registration with SSD between 2D medical images.The red frame is the transformed moving image using Clear+GMS registration.Green and Yellow frames are the transformed images using respectively PPIR(MPC)+GMS and PPIR(FHE)v1+GMS.

Fig. A7 :
Fig. A7: Qualitative results for rigid point cloud registration between 2D corpus callosum point sets.The red frame is the transformed moving image using Clear registration.Green and Yellow frames are the transformed images using respectively PPIR(MPC) and PPIR(FHE)v1.

Table 1 :
Overview of the datasets used in the study.PETs: Privacy Enhancing Technologies.
2DEF .In the initial phase, as illustrated in Supplementary Figure

Table 2 :
Implementation Details.The PPIR framework for the SSD is implemented using two state-of-the-art libraries: PySyft (Ryffel et al., 2018), which provides SPDZ two-party computation, Affine SSD registration test, comparison between Clear, PPIR(MPC), PPIR(FHE)-v1 and PPIR(FHE)-v2.Registration metrics are reported as mean and standard deviation.Efficiency metrics in terms of average across iterations.RMSE: root mean square error.

Table 3 :
Non-Linear SSD registration test comparison between Clear, PPIR(MPC), PPIR(FHE)-v1 and PPIR(FHE)-v2.The registration metrics are reported as mean and standard deviation.Efficiency metrics in terms of average across iterations.RMSE: root mean square error.

Table 4 :
Table 3, comparing Clear and PPIR(MPC), PPIR(FHE)-v1 and v2, showcases the metrics resulting from spline-based non-linear registration between grey MI Affine registration metrics Solution Intensity Error (MI) Num.Iteration Displacement RMSE (mm) Time party 1 (s) Time party 2 (s) Comm.party 1 (MB) Comm.party 2 (MB) Affine MI registration test, comparison between Clear and PPIR(MPC).Registration metrics are reported as mean and standard deviation.Efficiency metrics in terms of average across iterations.RMSE: root mean square error.CC ANTs registration registration metrics Solution Initial DICE score Final DICE score Num.Iteration Displacement RMSE (mm) Time party 1 (s) Time party 2 (s) Comm.party 1 (MB) Comm party 2 (MB)

Table 5 :
ANTs registration with CC, comparison between Clear and PPIR(MPC).Registration metrics are reported as mean and standard deviation.Efficiency metrics in terms of average across iterations.RMSE: root mean square error.