HNSF Log-Demons: Diffeomorphic demons registration using hierarchical neighbourhood spectral features

Many biomedical applications require accurate non-rigid image registration that can cope with complex deformations. However, popular diffeomorphic Demons registration algorithms suffer from difﬁculties for complex and serious distortions since they only use image greyscale and gradient information. To address these difﬁculties, a new diffeomorphic Demons registration algorithm is proposed using hierarchical neighbourhood spectral features namely HNSF Log-Demons in this paper. In view of three important properties of hierarchical neighbourhood spectral features based on line graph such as rotation invariance, invariance of linear changes of brightness, and robustness to noise, the hierarchical neighbourhood spectral features of a reference image and a moving image is ﬁrst extracted and these novel spectral features are incorporated into the energy function of the diffeomorphic registration framework to improve the capability of capturing complex distortions. Secondly, the Nystr ̈ om approximation based on random singular value decomposition is employed to effectively enhance the computational efﬁciency of HNSF Log-Demons. Finally, the hybrid multi-resolution strategy based on wavelet decomposition in the registration process is utilised to further improve the registration accuracy and efﬁciency. Experimental results show that the proposed HNSF Log-Demons not only effectively ensures the generation of smooth and reversible deformation ﬁeld, but also achieves better performance than state-of-the-art algorithms.


INTRODUCTION
Medical image registration refers to the process of finding a plausible spatial transformation for moving images, so that it can reach the same spatial relationship with the corresponding anatomical points or at least diagnostic points on reference images [1,2]. According to different spatial transformations, medical image registration is generally categorised into two groups, i.e. rigid registration and non-rigid registration. The rigid registration only describes the motion that is limited to global rotations and translations while non-rigid registration usually includes very complex local and global elastic deformations. Practically, in the process of image acquisition, same or different patients are often influenced by other factors, such as lung movement and bladder filling, resulting in complex non-rigid deformation. The non-rigid registration often plays an impor-the incorporation of new features. Hierarchical neighbourhood spectral features based on line graph can describe the unique and intrinsic geometry shape of images due to three useful properties: rotation invariance, invariance of linear changes of brightness and robustness to noise [8]. Therefore, it is suitable to guide the non-rigid registration procedure to capture complex distortions. Aiming to address the above issues, we proposed a diffeomorphic Demons registration algorithm using hierarchical neighbourhood spectral features in this paper, called HNSF Log-Demons for short. The proposed HNSF Log-Demons makes the following contributions: • A novel spectral matching similarity based on the hierarchical neighbourhood spectral features is presented and incorporated into the diffeomorphic registration framework, which not only solves complex distortions, but also ensures the diffeomorphic transformation and the registration accuracy. • The Nyström approximation based on random SVD is employed to reduce the spectral decomposition time of the high-dimensional weighted adjacency matrix and improve the registration efficiency. • By introducing the hybrid multi-resolution strategy based on wavelet decomposition, the registration procedure is executed from coarse to fine to further improve the registration accuracy and efficiency of the HNSF Log-Demons.
The rest of this paper is organised as follows. Section 2 illustrates the related work. Section 3 proposes a spectral matching similarity based on the hierarchical neighbourhood spectral features and describes the HNSF Log-Demons in detail. Section 4 demonstrates the registration performance of the HNSF Log-Demons in comparison to other well-established non-rigid registration algorithms. Section 5 discusses the practical properties of the HNSF Log-Demons. Finally, Section 6 makes the conclusion for this paper.

RELATED WORK
In recent years, a large number of non-rigid registration models for medical images has been reported, such as the free form deformation (FFD) models based on B-spline [9][10][11], finite elements model (FEM) [12][13][14][15], viscous fluid model [16,17], and Demons [4,5,[18][19][20][21][22] etc. Since FFD based on B-spline employs the cubic B-spline to model the elastic deformation and each B-spline curve is only related to four adjacent control points, it can provide a high degree of flexibility for estimating the local motions of human tissues and organs. Although FFD based on B-spline does not require assumptions about the elastic properties of the tissues and organs compared to physics-based deformation models, it is unable to effectively describe the global and large-scale deformation of human tissues and organs [9][10][11]. Image registration based on FEM transforms complex deformation into discrete and simple deformation elements, and then utilises these elements as the basic unit combined with the overall topological structure to fit the whole complex elastic deformation. Though FEM shows good adaptability to irregular geometry, it needs to solve many unknown parameters resulting in high computational cost [12][13][14][15]. Different with FFD and FEM, the viscous fluid model simulates the image registration as the fluid flow process, it can therefore better describe large-scale deformation. However, it is sensitive to greyscale changes of images, and the registration process suffers from serious time-consuming [16,17]. In these models, the Demons is a popular non-rigid registration algorithm based on the theory of optical flow field, and it has been widely used in medical image registration due to its complete mathematical theory [4,5,18]. However, the Demons algorithm suffers from three shortcomings. First, since only the gradient information of reference images is used to drive the deformation, the Demons algorithm has a very slow convergence speed. Secondly, when the value of gradients is close to zero in the textureless areas of reference images, the algorithm easily generates the wrong results of registration. Thirdly, Demons algorithm is based on the assumption that there is small deformation between images, therefore, it is difficult to accurately estimate the large deformation only using grey information of images. To speed up the convergence of registration procedure, Wang et al. [23] proposed the active Demons algorithm that introduces the gradient information of moving images into diffusion equation and utilises homogenisation coefficient to adjust the intensity of driving forces. However, the active Demons algorithm only adjusts the intensity of driving forces through the homogenisation coefficient, which cannot achieve a balance between the large and the small deformation. Xue et al. [24] presented an improved active Demons algorithm that introduces a new parameter called balance coefficient into the active Demons to adjust the driving force between the large and small deformation. To keep the topology of the deformation field and avoid the physical unreasonable deformation, Vercaiteren et al. [4,5] proposed the diffeomorphic Demons registration algorithm based on the Lie group theory, the algorithm regards the registration procedure as an optimisation process of energy function. Although the diffeomorphic Demons is a well-established intensity-based registration framework, it cannot capture the complex and large-scale distortions [6,7]. To tackle the problem and enhance the registration accuracy, many researchers presented improved algorithms that can be categorised coarsely into two different strategies for the diffeomorphic Demons. The two strategies for improving the diffeomorphic Demons are shown in Figure 1.

Employing a new similarity measure
Based on Vercauteren's work, many improved algorithms that employ new similarity measures are reported in [25][26][27][28]. For example, Avants et al. [25] incorporated the cross-correlation into the symmetric diffeomorphic registration framework, which improves the performance of medical image non-rigid registration. Lorenzi et al. [26] added stationary velocity fields (SVFs) constraints and used a new similarity measure named local correlation coefficient (LCC) in the Log-Demons registration framework for the deformation field. Authors indicated The recent advances of the application of diffeomorphic Demons for non-rigid image registration that the LCC-Demons was a generic, flexible, efficient and robust algorithm for non-rigid registration of images. Furthermore, Mehdi et al. [27] implemented a masking of the similarity term on the non-rigid registration, which enhances the robustness of algorithms with respect to intensity artefacts of the brain boundaries. Han et al. [28] extended the mono-modality, diffeomorphic form of the Demons algorithm to multi-modality registration using pointwise mutual information (pMI) as a similarity metric. The proposed pMI-Demons achieves registration accuracy comparable to MI-FFD and MI-SyN, and maintains diffeomorphic transformation similar to MI-SyN.

Incorporating new features into the driving force
Popular features mainly include similarity of grey gradient field [29], multi-support-region-order-based gradient histogram (MROGH) [30], geometric constraint [31,32], the well-known robust feature descriptors such as SIFT, SURF or ASIFT [33], unsupervised deep features using stack auto-encoder (SAE) extraction [34], modal independent neighbourhood descriptor (MIND) [35,36], spectral features [37] and Gabor features [38]. Reaungamornrat et al. [35] put forward a deformable image registration algorithm using the MIND and the Huber metric, the algorithm is robust against realistic anisotropic resolution characteristics in MR and yields high registration accuracy. However, the algorithm does not satisfy rotation and scale invariance properties, and requires a high computational cost. The spectral Demons is a new framework for capturing large and complex deformations in image registration, which achieves clear improvements in accuracy and robustness to large-scale deformations compared to conventional Demons algorithms. However, the final transformation obtained by spectral Demons is suitable for the registration of images with similar topological structure but unsuitable for the registration of images with occlusions and holes, which is bad for the monitoring and evaluation of postoperative rehabilitation. Additionally, the spectral Demons is time-consuming caused by spectral decomposition of high dimensional matrix, which limits its practicability in many clinical applications [37]. Wen et al. [38] extracted Gabor features of the reference and moving images to construct similarity metric since Gabor filters can effectively extract image texture information. Furthermore, they proposed an inertial constraint strategy to provide guided information for updating the current velocity field. As a result, both Gabor features and inertial constraint strategy improve the performance of Demons algorithm [38].
However, many improved non-rigid registration algorithms utilising the above two strategies are usually time-consuming due to the introduction of new similarities or new features. In practice, there are two popular ways to accelerate these improved algorithms: the multi-resolution strategy and the parallel computing technology. Researchers applied multiresolution strategy in the registration process to speed up the registration convergence [39,40]. In the coarse stage of multiresolution registration, the similarity measure is calculated for images with low resolution, which is not only fast but also can suppress local extrema. And then, registration parameters from the coarse stage are set as initial parameters of the fine stage, which can shorten the time of optimisation process and effectively improve the registration efficiency. In addition, a lot of non-rigid image registration algorithms are accelerated based on parallel computing techniques such as GPU and other multicore hardwares, which become an important way to improve the registration performance [41,42]. For instance, Gu et al. [43] implemented the classical Demons registration algorithm and other five improved Demons algorithms based on CUDA on GPU, and evaluated their performance. Shackleford et al. [44] used GPU's flow processing model to accelerate the Demons algorithm, where Nvidia 680 GTX GPU is selected leading to a 150 times speedup ratio for large medical volumes (e.g. 250 3 voxels).

METHODOLOGY
In this section, we propose the HNSF Log-Demons that is a hierarchical neighbourhood spectral features-based Log-Demons algorithm for image registration with complex deformation. First, we introduce briefly the proposed registration framework. Secondly, we design the extraction of hierarchical neighbourhood spectral features. Thirdly, we propose a spectral matching algorithm to calculate the similarity measure based on hierarchical neighbourhood spectral features and integrate it into the energy function of the proposed registration framework. At last, we summarise the HNSF Log-Demons.

The proposed registration framework
To maintain the reversibility of displacement field and the topological structure of images, the HNSF Log-Demons is based on diffeomorphic registration framework. Inspired by the diffeomorphic Log-Demons, we design the registration framework as shown in Figure 2.
In this framework, we first adopt the hybrid multi-resolution strategy based on wavelet decomposition, and the registration process is performed from coarse to fine, so as to improve the execution speed of registration process and prevent the registration process falling into local optimum solutions. More importantly, in view of the advantages of the hierarchical neighbourhood spectral features in geometry structure representation, we extract hierarchical neighbourhood spectral features of a reference image and a moving image, respectively. Furthermore, we combine the features with other important information of the reference and moving images, such as grey scale, gradient and space coordinate, to obtain robust and effective features used for non-rigid registration with complex distortions. Finally, we incorporate these combined features into the energy function of the proposed diffeomorphic registration framework to obtain better registration models.
Although HNSF Log-Demons is inspired by the diffeomorphic Demons, there are two difference between them. On the one hand, the energy function of diffeomorphic Demons only utilises the grey scale, gradient, and spatial information of the reference and moving images to guide the registration process, but HNSF Log-Demons employs much more information, i.e. hierarchical neighbourhood spectral features of the reference and moving images, which can describe geometry structures of the reference and moving images. Therefore, HNSF Log-Demons possesses a more powerful driving force than the diffeomorphic Demons and can capture more complex distortions. On the other hand, HNSF Log-Demons incorporates the hybrid multi-resolution strategy into the registration process, it can thus improve the accuracy and efficiency of registration better than diffeomorphic Demons.

Hierarchical neighbourhood spectral features
For the image I of size m × n, the feature of any non-boundary pixel p can be determined by the relationship between its neighbourhood pixels, so the feature of pixel p is represented based on the graph [8]. Supposing the eight neighbourhood pixels of p as q i (i = 1, 2, … , 8), we can construct the star graph S . In this graph, all these neighbourhood pixels are connected by the edge pq i and the weight r i of an edge is: where I (p) and X (p) are the greyscale value and the space coordinate of the pixel p, respectively. To further reflect the characteristic of the pixel p, the line graph L(S ) and the weighted adjacency matrix A corresponding to L(S ) are constructed. The element of matrix A is: where 1 ≤ i, j ≤ 8, i, j ∈ R + . Since the weighted adjacency matrix A is a real symmetric matrix, eight real eigenvalues can be obtained after spectral decomposition, and 1 (p) ≥ 2 (p) ≥ ⋯ ≥ 8 (p). Using i (p) to replace grey value of the pixel p and regularising it to [0, 255], we can get the eigenvalue image If we only use the pixel p and its eight neighbourhood pixels to construct a graph and obtain its spectral features, the features of the pixel p cannot be fully described. By enlarging the neighbourhood range of the pixel p, we can get (2k + 1) 2 − 1 neighbourhood pixels, where k is the number of layers of neighbourhood pixels. However, if all the selected neighbourhood pixels are used to construct line graph, the construction of the adjacency matrix and the spectral decomposition will require a high computational cost when the value of k is large. To reduce computational cost, we divide the neighbourhood pixels into k layers according to the distance to the pixel p, and utilise the hierarchical neighbourhood method to construct spectral features as shown in Figure 3.
In Figure 3, the feature construction consists of four steps. The first step is to create the hierarchical neighbourhood layer by layer and construct the corresponding line graphs. The second step is to compute spectral features of the pixel p from the t -layer(1 ≤ t ≤ k): The third step is to select main features ′ t from t to reduce the computational cost of feature matching. The last step is to combine the main spectral features of each layer to obtain the final spectral feature vector of the pixel p: where ‖.‖ is the L 2 norm. In the process of extracting hierarchical neighbourhood spectral features, it is necessary to construct weighted adjacency matrix and implement the spectral decomposition of the matrix. With the increase of the value of k(i.e. the number of layers), the size of the weighted adjacency matrix increases rapidly. Therefore, it will take much time to implement the spectral decomposition of the matrix. To reduce this expensive cost, we employ the Nyström approximation based on random SVD [45] to realise the low rank approximation of the weighted adjacency matrix. The Nyström approximation based on random SVD finally outputs t that is the approximate eigenvalue of A. Therefore, it does not require to directly decompose the high dimensional weighted adjacency matrix. We can significantly reduce the computational cost of spectral decomposition and improve the final registration efficiency.
There are three important properties of the hierarchical neighbourhood spectral features based on the line graph [8]: rotation invariance, invariance of linear changes of brightness, and robustness to noise.

Similarity measure based on HNSF
The registration process can be degraded to the problem of minimising the following energy function: where v is the regular diffeomorphic non-rigid transformation, and c is a non-parametric spatial transformation and an exact realisation of the spatial transformation v, which allows some error at each image pixel. The hidden variable c is introduced into the energy function, so that the diffeomorphic Log-Demons algorithm is regarded as an optimisation problem for a suitable criterion. The dist(c, v) = ‖c − v‖ represents the distance error before and after regularisation, and the Reg(v) is a regular term. The x and T are the traditional Demons parameters controlling the step size and the regularisation, respectively.
To deal with the complex deformation, we consider to combine the greyscale information, the spatial information and spectral features into the feature space of pixels to implement image registration. Meanwhile, R = ( i I R , s X R , g R ) and M = ( i I M , s X M , g M ) are obtained by weighting grey information, spatial information and hierarchical neighbourhood spectral features, so the similarity measure of the HNSF Log-Demons is given as follows: where integrates the norms of the displacements between corresponding points and acts as a form of spatial regularisation, and the ‖ R − M•exp(c ) ‖ 2 is a hierarchical neighbourhood spectral matching similarity between the reference and moving images, and parameters i , s and g are used to control the consistency in intensity, space and geometry, respectively. In this work, we propose the hierarchical neighbourhood spectral matching algorithm to calculate the spectral feature similarity between the reference and moving images, which is shown in Algorithm 1. First, the extraction of spectral features for each pixel includes four steps: constructing line graphs, calculating the spectral features for each layer, selecting main features and concatenating all feature vectors of all layers. Secondly, we need to reorder spectral features of the moving image with respect to spectral features of the reference image. Finally, we embed hierarchical neighbourhood spectral features of the reference and moving images into their feature spaces, respectively, and find a map from the feature space of the reference image to the feature space of the moving image. Owing to hierarchical neighbourhood spectral matching, the newly updated field not only takes the grey-scale and gradient information into account, but also utilises local structural features of images.

The HNSF Log-Demons
To keep the reversibility of displacement field and the topological structure of images, we select the diffeomorphic registration framework in this work. According to the theory of Lie group, a diffeomorphic transformation , which resides on a Lie group structure, is related to the exponential map of a velocity field v, and is associated with Lie algebra: = exp(v). However, in the registration process, every iteration needs to calculate the exponential map of the velocity field, so a simple and efficient approximate method is needed to calculate the exponential map.
In the case of stationary velocity fields, a practical and simple approximation is possible with the scaling-and-squaring method [5]. In the proposed HNSF Log-Demons algorithm, we employ the scaling-and-squaring to implement the exponential map.
To improve the accuracy and efficiency of Demons registration for medical images with complex distortions, we put forward a diffeomorphic Demons registration algorithm based on hierarchical neighbourhood spectral features as shown in Algorithm 2. In the Line 1 of Algorithm 2, the hybrid multiresolution strategy based on wavelet decomposition is utilised to generate the image pyramids of the reference and moving images. The velocity field applicable to the first layer of image pyramids is initialised in the Line 2. From the Line 3-16, for the image pyramids, the HNSF Log-Demons executes the registration process from the coarse layer to fine layer. From the Line 4-6, the HNSF Log-Demons initialises the velocity field of the current layer according to the velocity field of the previous layer. From the Line 7-14, first, the HNSF Log-Demons computes the updating fields u R→M and u M→ R , respectively using the Algorithm 1. And then, it computes the average updating field and smooths the updating field using Gaussian filter. At last, the velocity field is updated by adding the updating field and smoothed using a diffusion kernel. In the Line 17, the HNSF Log-Demons outputs the optimal transformation * using the exponent map. The HNSF Log-Demons extracts hierarchical neighbourhood spectral features of the reference and moving images and introduces them into the energy function of the proposed registration framework, ensuring that the generated deformation field has smoothness and reversibility, and provides better registration accuracy. Meanwhile, the Nyström approximation based on random SVD is utilised to reduce the computational burden of spectral decomposition and improve the registration efficiency of the HNSF Log-Demons. In addition, in the registration process, the hybrid multi-resolution strategy based on wavelet decomposition is combined to further improve the registration accuracy and efficiency. In the HNSF Log-Demons, the corresponding updates u can be computed using hierarchical neighbourhood spectral matching at the coarser levels of resolution, while local details and deformations can be computed more efficiently with the conventional gradient-based updates at the finer levels of resolution. Hierarchical neighbourhood spectral features are computed only at the coarser levels of resolu-tion, and the coarser level has a smaller image size. Therefore, it can complete image registration quickly.

EXPERIMENTS
We conducted experiments on two types of images: medical images that are randomly selected from the database of retrospective image registration evaluation (RIRE) [46] at the Vanderbilt university, as well as real images. Especially, to demonstrate the generalisation of the HNSF Log-Demons on medical image modality and body part, the test dataset of medical images includes different image modalities, such as CT, MRI-T1, MRI-T2 and MRI-T2-rectified, and includes different parts of human body, such as head, chest and abdomen. All algorithms and experimental evaluations are performed on a PC with an Intel(R) Core (TM)CPU, i7-9700, 4.20 GHz×8, and 32 GB RAM.

Parameter settings
To evaluate the effectiveness and efficiency of the HNSF Log-Demons, six popular non-rigid registration algorithms are considered for comparisons. They are the FFD based on B-spline [9], Demons [18], active Demons [23], improved active Demons [24], Log-Demons [5] and spectral Demons [37]. For the sake of fairness, the FFD algorithm based on B-spline changes parameters until the best performance. The active Demons, Log-Demons, spectral Demons and HNSF Log-Demons choose same parameters. Moreover, all algorithms adopt the three-level multi-resolution strategy, and the maximal number of iterations from high-resolution to low-resolution are 100, 50 and 20, respectively. To balance efficiency and accuracy of the HNSF Log-Demons, hierarchical neighbourhood spectral features are computed only at coarse levels of resolution, and the number of layer of hierarchical neighbourhood is set to 10 to capture complex large-scale distortions. Additionally, we use r = 2 eigenmodes with the HNSF Log-Demons to reduce the complexity of similarity calculation. All Demons algorithms employ subsequently the same parameters, fluid = 1, diffuse = 1, and x = 1. The experiment in [37] demonstrates that a larger value of weight on spatial consistency actually degrades the registration accuracy, since s penalises large deformations. Therefore, the weight value on spatial consistency should be low but non-zero, and the best empirical value of parameters should meet the condition of 1∕10 < s ∕ g < 1∕3. In practice, i = 0.8, s = 0.05 and g = 0.15 are often are used for spectral Demons to favour consistency in geometry and intensity, and to prevent unreasonable displacements during each iteration [37]. Compared to ref. [37], we use hierarchical neighbourhood spectral features instead of global spectral features in the HNSF Log-Demons. However, the similarity between the HNSF Log-Demons and ref. [37] is both algorithms make use of grey, spatial information and spectral features. Among them, the weight of spatial information should not be too large, otherwise it will lead to the decline of registration accuracy. Therefore, in the HNSF Log-Demons, the values of three weight parameters are consistent with ref. [37], especially, the weight value on spatial consistency should be low but non-zero, and the best value of parameters should meet the condition of 1∕10 < s ∕ g < 1∕3.
Meanwhile, for the sake of fairness, the HNSF Log-Demons choose the same weighting parameters i = 0.8, s = 0.05 and g = 0.15. In this work, registration results are evaluated by subjective visualisation and objective quantitative evaluation. There are four metrics for objective quantitative evaluation [47]: sum of squared differences (SSD), sum of absolute differences (SAD), normalisation mutual information (NMI) and normalisation cross correlation (NCC). For these evaluation metrics, fine registration results correspond to large values of NMI and NCC but small values of SSD and SAD.

Results on medical images
In this subsection, the non-rigid registration experiment is conducted on medical images using the above seven algorithms. Before this experiment, medical images are first registered rigidly, and then, the normalisation processing is carried out to effectively avoid the influence of registration results due to data inconsistency caused by different equipment positioning. Registration experiments are conducted by using the FFD based on B-spline, Demons, active Demons, improved active Demons, Log-Demons, spectral Demons and HNSF Log-Demons. Figure 4 shows registration results on six pairs of medical images with different modalities and body parts using the above seven algorithms. In Figure 4, owning to the complex large-scale deformation in moving images, the FFD based on B-spline, Demons, active Demons, improved active Demons and Log-Demons easily fall into local optimum in the registration process, which limits deformation correction of registration results. Therefore, these algorithms only provide registration results with low registration accuracy. Compared with the FFD based on B-spline, Demons, active Demons, improved active Demons and Log-Demons, the spectral Demons provides better registration results. However, these results still include obvious mismatched areas. Compared with other algorithms, there is no obvious registration error in difference images by using the HNSF Log-Demons, which achieves the best smoothness and the smallest registration error.
In addition, we quantitatively evaluate the registration accuracy of the HNSF Log-Demons. The above seven algorithms are used to match six pairs of medical images, and the values of four objective evaluation metrics for registration results are shown in Table 1. It can be seen that the HNSF Log-Demons achieves the best registration result in term of four metrics, which is also consistent with the observation result of the previous subjective visualisation. Therefore, this experiment concludes that the proposed HNSF Log-Demons can achieve better registration accuracy than other six algorithms for different modality medical images and different parts of human body.

Results on real images
To further evaluate the registration accuracy and robustness of the HNSF Log-Demons, we conducted experiments on four pairs of real images from public datasets. Experimental results are shown in Figure 5. From Figure 5, we can see that the Demons, active Demons, improved active Demons, Log-Demons and FFD based on B-spline achieve poorer registration results than the spectral Demons and HNSF Log-Demons. Therefore, it is difficult to obtain fine transformations for registration results. An important reason is that these algorithms only employ image grey and gradient information to provide driving force, and then use a Gaussian filter to update the displacement deformation field smoothly. Consequently, some obvious mismatched areas appear in Figure 5(c-g) due to insufficiently driving force for complex deformation. Spectral Demons can capture the global geometric similarity of the moving image by using spectral features, and can deal with global large-scale non-rigid deformation, but the improvement effect of local small deformation is limited. The HNSF Log-Demons can deal with the complex deformation better in the moving image. Meanwhile, the HNSF Log-Demons can also effectively deal with the local non-rigid deformation and improve the registration accuracy, which benefits from the employment of hierarchical neighbourhood spectral features. Furthermore, the quantitative result in Table 2 is consistent with subjective visualisation evaluation of registration results shown in Figure 5. From Table 2, it can be seen that the HNSF Log-Demons that integrates the hierarchical neighbourhood spectral features into the energy function of proposed registration framework, is superior to other six algorithms on real images in term of the four evaluation metrics.

DISCUSSION AND ANALYSIS
In this section, we discuss HNSF Log-Demons from four aspects, which mainly includes convergence, efficiency, computational complexity and limitations of HNSF Log-Demons.

Convergence analysis
In this subsection, we compared the HNSF Log-Demons with the Log-Demons and spectral Demons to further analyse and verify the convergence of the HNSF Log-Demons. We performed non-rigid registration for a pair of medical images and a pair of real images, respectively. For the sake of showing the trends of energy function with iteration number, we normalised the energy function value of each iteration in this experiment. Again, since we analyse the convergence of the HNSF Log-Demons, the same three-layer multi-resolution strategy is applied in these algorithms. In the third resolution image, we  purposely calculate the change of energy function with the number of iterations. The curves of the energy function value are shown in Figure 6. As can be seen from Fig. 6, the HNSF Log-Demons shows a fast convergence rate for two pairs of images, and it closes to the optimal registration result after about 12-13 iterations.
We propose a new spectral matching similarity based on hierarchical neighbourhood spectral features and incorporate it into the energy function of diffeomorphic registration framework, which provides the deformation field with a prior information that is invariant to rotation and brightness, and robust to noise. This prior information is used to guide the whole registration process, which overcomes the problem of insufficient driving force, accelerates the convergence speed of the energy function and improves the registration performance.

Efficiency analysis
Here we demonstrate the registration efficiency of the HNSF Log-Demons. For the above ten pairs of images, seven registration algorithms are implemented on the same experimental environment. The parameter settings of all algorithms are the same as above. To eliminate the system error, all algorithms are executed ten times for each pair of images, and then we compare the average running time. Average running time of seven algorithms is shown in Figure 7. Since the grey and gradient information are only considered in the energy functions of Demons, active Demons and improved active Demons, these three algorithms take relatively short running time. The Log-Demons requires much time than the Demons to obtain the reversible deformation field. The FFD based on B-spline involves many parameters to deal with the complex deformation and takes longer running time than the Demons, active Demons and improved active Demons. The spectral Demons and HNSF Log-Demons introduce spectral features into the registration process, and the spectral matching similarity is computed in each iteration. Therefore, both the spectral Demons and HNSF Log-Demons take longer running time than other comparative algorithms in each iteration. However, the spectral Demons and HNSF Log-Demons introduce spectral information into their energy functions, which makes the registration process have stronger driving force and achieve faster convergence. According to convergence analysis, the HNSF Log-Demons can approximate the best registration result after about 12-13 iterations, so we can get the same results with Log-Demons by reducing the number of iterations, and the HNSF Log-Demons provides a better registration accuracy than the Log-Demons. Compared with the spectral Demons, the HNSF Log-Demons employs the Nyström approximation based on random SVD to decompose high dimensional adjacency matrix, which effectively reduces the computational complexity of spectral decomposition. Therefore, the HNSF Log-Demons achieves better registration performance than the spectral Demons.

Computational complexity
For the HNSF Log-Demons, although extracting hierarchical neighbourhood spectral features is very important for improving the registration accuracy, it needs to pay the price of time. In the next, we briefly analyse the computational complexity of extracting hierarchical neighbourhood spectral features. The procedure of extracting hierarchical neighbourhood spectral features mainly includes two important components: the construction of adjacency matrix and the spectral decomposition. There are 8i points at the ith layer, therefore, considering the symmetry, the calculation cost of constructing the adjacency matrix of the ith layer is 8i(8i − 1)∕2 = 4i(8i − 1), the calculation cost of spectral decomposition is O((8i ) 3 ) = O(512i 3 ), and the calculation cost of generating total features is 4 ). According to above analysis, it can be indicated that the spectral decomposition is relatively timeconsuming than the construction of line graphs, especially when the value of k is large. For this problem, in the HNSF Log-Demons, we adopt the Nyström approximation based on random SVD to reduce the time of spectral decomposition when the value of k is larger than five. Here, the computational complexity of the Nyström approximation based on random SVD is analysed and compared with that of spectral decomposition of weighted adjacency matrix. For a weighted adjacency matrix A with the size n × n, the computational complexity of spectral decomposition is O(n 3 ). However, in the Nyström approximation based on random SVD, the computational complexity of generating the approximation matrix of A is O(m 2 k), the computational complexity of QR decomposition is O(mk), the computational complexity of SVD is O(k 3 ). Therefore, the overall computational complexity is O(m 2 k + k 3 ) due to n >> m ≥ k. To summarise, the computational complexity of the Nyström approximation based on random SVD is less than that of the spectral decomposition of weighted adjacency matrix. Consequently, the Nyström approximation based on random SVD can reduce the execution time of the HNSF Log-Demons and thus improve the efficiency of registration process.

Limitations of the HNSF Log-Demons
Although the proposed HNSF Log-Demons shows that hierarchical neighbourhood spectral features can improve the registration accuracy for moving images with complex deformation, it still suffers from two limitations as follows. From the perspective of practice, since hierarchical neighbourhood spectral features are usually extracted from the neighbourhood information of central pixels, the HNSF Log-Demons depends on an assumption that the reference and moving images have a similar topology without missing tissues or organs. This assumption is generally true in many practical applications, especially, we often perform image registration for same or similar objects. However, the assumption could be a limitation in some specific clinical applications. For example, registering some medical images before and after tumour resection surgery, where tissues or organs are clearly missed.
From the perspective of clinical application, the HNSF Log-Demons takes too much execution time to meet the high real-time requirement of several clinical applications, such as the online adaptive radiotherapy and the image-guided surgery. Since the limitation of the proposed HNSF Log-Demons currently resides in the computational cost of extracting hierarchical neighbourhood spectral features for each pixel of reference and moving images, especially, the reference and moving images have a high resolution. Therefore, there is still much room for improving the registration efficiency of the HNSF Log-Demons.

CONCLUSION
In this work, we have studied non-rigid image registration algorithms using spectral features. We proposed a novel diffeomorphic Demons registration algorithm for images with complex distortions by incorporating hierarchical neighbourhood spectral features. The proposed HNSF Log-Demons was used to match medical and real images, respectively, demonstrating that it can capture complex distortions of moving images due to the incorporation of the hierarchical neighbourhood spectral features. Meanwhile, the proposed algorithm reduces the computational cost of hierarchical neighbourhood spectral features in registration process by employ the Nyström approximation based on random SVD and the hybrid multi-resolution strategy. Moreover, the HNSF Log-Demons can effectively achieve the smooth and reversible displacement field for images with complex distortions and provide better performance of non-rigid registration than conventional Demons algorithms.
In future work, further reducing greatly the computational burden of the HNSF Log-Demons is very urgent and significant for several real-time clinical applications. Fortunately, owing to the intrinsic pixel-wise nature of extracting hierarchical neighbourhood spectral features in the HNSF Log-Demons, parallel computing techniques such as GPU or many-core architecture will be utilised to significantly improve the computational efficiency. Simultaneously, the future work may focus on deformation invariant feature descriptors and analysis of sensitivity to noise, as well as incorporation of some prior knowledge (e.g. edge information of key objects) to further constrain the solution space of the registration process and raise the registration accuracy. In addition, it is also the subject of future investigation to extends hierarchical neighbourhood spectral features for 3D reference and moving images and apply the HNSF Log-Demons to 3D medical images, such as intra-operative CBCT images.