## 2 Materials and Methods
### Deformable Registration of Point Clouds  
Let us consider two point clouds: the fixed point cloud \(\mathcal{P}=\mathbf{p}_{1},\mathbf{p}_{2},\ldots,\mathbf{p}_{n}\) and the moving point cloud \(\mathcal{Q}=\mathbf{q}_{1},\mathbf{q}_{2},\ldots,\mathbf{q}_{m}\), where \(\mathbf{p}_{i},\mathbf{q}_{j}\in\mathbb{R}^{3}\). The goal is to find a transformation \(\mathcal{T}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}\) that aligns \(\mathcal{Q}\) to \(\mathcal{P}\) by minimizing the alignment error. In the case of nonrigid registration, a global transformation is insufficient to recover point correspondences accurately due to the deformable nature of the shapes. Therefore, we approximate the transformation function locally by grouping points into \(c\) clusters, where neighboring points are likely to undergo similar transformations. The deformation of the moving point cloud is modeled as follows:  
\[\mathbf{q}i=\sum_{j=1}^{c}w_{ij}\mathcal{T}_{j}(\mathbf{p}_{i}\ \mathbf{1}),\quad \forall\mathbf{p}_{i}\in\mathcal{P}, \tag{1}\]  
where \(\mathcal{T}_{j}\in\mathbb{R}^{3\times 4}\) represents the affine transformation matrix for cluster \(j\), and \(w_{ij}\) is a weight indicating the influence of cluster \(j\) on point \(\mathbf{p}_{i}\). The homogeneous coordinate (factor 1) is used to facilitate affine transformations. The choice of the number of clusters \(c\) is crucial; if \(c=n\), where \(n\) is the number of points, the problem becomes mathematically and computationally intractable. Conversely, a very low \(c\) may not capture the complexity of the deformation. We typically select a \(c\) such that \(c<n\) is sufficient to model the deformation accurately. To determine the clusters, we use an initial set of correspondences, possibly obtained through a matching strategy or feature descriptors, to group each point \(\mathbf{p}_{i}\) with the key points \(\mathbf{n}_{1},\ldots,\mathbf{n}_{k}\). We define the surface representation of \(\mathcal{P}\) as \(\mathcal{D}\) and of \(\mathcal{Q}\) as \(\mathcal{S}\). We then partition the surfaces into the patches \(\mathcal{D}_{1},\ldots,\mathcal{D}_{k}\subset\mathcal{D}\) and \(\mathcal{S}_{1},\ldots,\mathcal{S}_{k}\subset\mathcal{S}\) based on geodesic distances from the key points. Each patch \(\mathcal{D}_{i}\) and \(\mathcal{S}_{i}\) is further divided into \(r\) neighboring rings, \(\mathcal{D}_{i}^{1},\ldots,\mathcal{D}_{i}^{r}\) and \(\mathcal{S}_{i}^{1},\ldots,\mathcal{S}_{i}^{r}\), respectively, based on geodesic distances. The entire surface can thus be represented as \(\mathcal{D}=\bigcup i=1^{k}\mathcal{D}_{i}\), and each patch as \(\mathcal{D}i=\bigcup j=1^{r}\mathcal{D}_{i}^{j}\). The same applies to \(\mathcal{S}\). We estimate the deformation hierarchically using different geometric levels: local, intermediate, and global. At the local level, we estimate the transformation \(\mathcal{T}_{ij}\) for each neighboring ring \(\mathcal{D}_{i}^{j}\). At the intermediate level, we estimate the transformation \(\mathcal{T}_{i}\) for each patch \(\mathcal{D}_{i}\) based on the deformed neighboring rings. Finally, at the global level, we estimate the overall transformation \(\mathcal{T}\) based on the deformed patches. The estimation of the transformation parameters is performed iteratively, updating one level at a time while keeping the others fixed, as summarized in Algorithm 1.