Photoacoustic-MR Image Registration Based on a Co-Sparse Analysis Model to Compensate for Brain Shift

Brain shift is an important obstacle to the application of image guidance during neurosurgical interventions. There has been a growing interest in intra-operative imaging to update the image-guided surgery systems. However, due to the innate limitations of the current imaging modalities, accurate brain shift compensation continues to be a challenging task. In this study, the application of intra-operative photoacoustic imaging and registration of the intra-operative photoacoustic with pre-operative MR images are proposed to compensate for brain deformation. Finding a satisfactory registration method is challenging due to the unpredictable nature of brain deformation. In this study, the co-sparse analysis model is proposed for photoacoustic-MR image registration, which can capture the interdependency of the two modalities. The proposed algorithm works based on the minimization of mapping transform via a pair of analysis operators that are learned by the alternating direction method of multipliers. The method was evaluated using an experimental phantom and ex vivo data obtained from a mouse brain. The results of the phantom data show about 63% improvement in target registration error in comparison with the commonly used normalized mutual information method. The results proved that intra-operative photoacoustic images could become a promising tool when the brain shift invalidates pre-operative MRI.


Introduction
Maximal safe resection of brain tumors in eloquent regions is optimally performed under image-guided surgery systems [1,2]. The accuracy of the image-guided neurosurgery system is drastically affected by intra-operative tissue deformation, called brain shift. Brain shift is a dynamic and complex spatiotemporal phenomenon that happens after performing a craniotomy and invalidates the pre-operative images of patients [3,4]. The brain shift, which is known as brain deformation, is a combination of a wide variety of biological, physical, and surgical causes and occurs in both cortical and deep brain structures [2,[5][6][7]. Brain shift calculation and compensation methods are based on updating the pre-operative images with regard to the intraoperative tissue deformation. These methods fall into two main categories: biomechanical models and intra-operative imaging approaches. Biomechanical model-based approaches are time and computation-consuming methods; however, they can be highly accurate [8][9][10]. The main drawback of model-based techniques ies benefits from simplicity and fast implementation. However, it is not customized for different types of data. In the second approach, an over-complete dictionary is constructed via learning methods. Among learning methods, the K-singular value decomposition (K-SVD) method has been widely used for image registration [58]. There are some studies which used synthesis sparse models for multimodal image registration [59]. However, a learned dictionary includes a large number of atoms. This leads to the increased computational complexity of multi-modal image registration, which is not suitable for the real-time compensation of brain shift. The analysis sparse model, named the co-sparse analysis model, represents a powerful alternative to the synthesis sparse representation approach in order to reduce the computational time [60]. Co-sparse analysis models can yield richer feature representations and better results for image registration in real-time processes. As a result of richer feature representation using co-sparse analysis models, better results for image registration can be obtained in real-time processes [61,62]. There are a few studies for multi-modal image registration via a co-sparse analysis model, but none of them were in the medical field. Kiechle et al., proposed an analysis model in a joint co-sparsity setup for different modalities of depth and intensity images [63]. Chang Han et al., utilized the analysis sparse model for remote sensing images [64] and Gao et al., used it to register multi-focus noisy images with higher quality images [65]. In our previous work, we applied an analysis sparse model for US-MR image registration to compensate for the brain shift [66].
To date, a few research studies have investigated PA and MR image registration. Ren et al., proposed a PA-MR image registration method based on mutual information to yield more insights into physiology and pathophysiology [67]. Gehrung et al., proposed the co-registration of PA and MR images of murine tumor models for the assessment of tumor physiology [68]. However, these studies were dedicated to solving the rigid registration problems and did not focus on the intra-operative application of PA imaging, and therefore did not face any complicated brain deformation.
To the best of our knowledge, in this study, for the first time, PA and MR image registration was used for the purpose of compensating complicated brain shift phenomena. The co-sparse analysis model is proposed for PA-MR image registration, which is able to capture the interdependency of two modalities. The algorithm works based on the minimization of mapping transform by using a pair of analysis operators which are learned by the alternating direction method of multipliers (ADMM).

Brain-Mimicking Phantom Data
To assess the performance of the multi-modal image registration algorithm to compensate for brain shift, a phantom that mimics brain tissue was prepared. The phantom was made of polyvinyl alcohol cryogel (PVA-C) which has been successfully used for mimicking brain tissue in previous studies [19]. The PVA-C material also has been applied in the fabrication of phantoms for ultrasound, MRI, and, recently, PA imaging [69]. A 10% by weight PVA in water solution was used to form PVA-C, which is solidified through a freeze-thaw process. The dimensions of the phantom were approximately 150 × 40 mm, with a curved top surface mimicking the shape of a head as shown in Figure 1a. Two plastic tubes with 1.2 and 1.4 mm inside diameters were inserted randomly into the mold before the freeze-thaw cycle to simulate blood vessels. Figure 1b shows the 3D model of the phantom including random vessels. Two types of chromophores, copper sulfate pentahydrate (CuSO 4 (H 2 O) 5 ) and human blood (1:100 dilution), were used to fill the embedded vessels before PA imaging (Figure 1c). To acquire MR images of the phantom before any deformations, the phantom was scanned using a Siemens 1.5 Tesla scanner using a standard T1-and T2-weighted protocol. Pulse-sequence parameters were set to TR = 600 ms, TE = 10 ms, and Ec = 1/1 27.8 kHz for T1-weighted and TR = 8.6, TE = 3.2, TI = 450, and Ec = 1/1 31.3 kHz for T2-weighted considering a 1 mm slice thickness with full brain phantom coverage and a 1 mm isotropic resolution.
PA images were achieved by using an ultrasound scanner (Vantage 128, Verasonics Inc., Kirkland, WA, USA) with a 128-element linear array US transducer (L11-4v, Verasonics, Inc., Kirkland, WA, USA) operating at a frequency range between 4 and 9 MHz. A pulsed tunable laser (PhocusCore, Optotek, California, USA) and Nd:YAG/OPO nanosecond pulsed laser (Phocus core system, OPOTEK Inc., Carlsbad, CA, USA), with a pulse repetition rate of 10 Hz at wavelengths of 700, 800, and 900 nm, were used to illuminate the phantom. The scan resolution was 1 mm, and the laser fluence was ~1 mJ/cm 2 ( Figure  2). It is notable that we used frame averaging for de-noising and spectral un-mixing as an image reconstruction algorithm to obtain high quality PA images.  To acquire MR images of the phantom before any deformations, the phantom was scanned using a Siemens 1.5 Tesla scanner using a standard T1-and T2-weighted protocol. Pulse-sequence parameters were set to TR = 600 ms, TE = 10 ms, and Ec = 1/1 27.8 kHz for T1-weighted and TR = 8.6, TE = 3.2, TI = 450, and Ec = 1/1 31.3 kHz for T2-weighted considering a 1 mm slice thickness with full brain phantom coverage and a 1 mm isotropic resolution.
PA images were achieved by using an ultrasound scanner (Vantage 128, Verasonics Inc., Kirkland, WA, USA) with a 128-element linear array US transducer (L11-4v, Verasonics, Inc., Kirkland, WA, USA) operating at a frequency range between 4 and 9 MHz. A pulsed tunable laser (PhocusCore, Optotek, CA, USA) and Nd:YAG/OPO nanosecond pulsed laser (Phocus core system, OPOTEK Inc., Carlsbad, CA, USA), with a pulse repetition rate of 10 Hz at wavelengths of 700, 800, and 900 nm, were used to illuminate the phantom. The scan resolution was 1 mm, and the laser fluence was~1 mJ/cm 2 ( Figure 2). It is notable that we used frame averaging for de-noising and spectral un-mixing as an image reconstruction algorithm to obtain high quality PA images. To acquire MR images of the phantom before any deformations, the phantom was scanned using a Siemens 1.5 Tesla scanner using a standard T1-and T2-weighted protocol. Pulse-sequence parameters were set to TR = 600 ms, TE = 10 ms, and Ec = 1/1 27.8 kHz for T1-weighted and TR = 8.6, TE = 3.2, TI = 450, and Ec = 1/1 31.3 kHz for T2-weighted considering a 1 mm slice thickness with full brain phantom coverage and a 1 mm isotropic resolution.
PA images were achieved by using an ultrasound scanner (Vantage 128, Verasonics Inc., Kirkland, WA, USA) with a 128-element linear array US transducer (L11-4v, Verasonics, Inc., Kirkland, WA, USA) operating at a frequency range between 4 and 9 MHz. A pulsed tunable laser (PhocusCore, Optotek, California, USA) and Nd:YAG/OPO nanosecond pulsed laser (Phocus core system, OPOTEK Inc., Carlsbad, CA, USA), with a pulse repetition rate of 10 Hz at wavelengths of 700, 800, and 900 nm, were used to illuminate the phantom. The scan resolution was 1 mm, and the laser fluence was ~1 mJ/cm 2 ( Figure  2). It is notable that we used frame averaging for de-noising and spectral un-mixing as an image reconstruction algorithm to obtain high quality PA images. Figure 2. Schematic of the PA imaging setup, which includes a tunable pulsed laser and a programmable ultrasound data acquisition system. Figure 2. Schematic of the PA imaging setup, which includes a tunable pulsed laser and a programmable ultrasound data acquisition system.

Murine Brain Data
For further evaluation of the proposed image registration method, we used ex vivo mouse brain data which was provided by Ren et al., in a previous study [67]. After removal of the mouse brain skull, the whole brain of the mouse was embedded in 3% agar in phosphate-buffered saline and was then imaged ex vivo. To acquire T2-weighted MR images of the mouse brain, a 2D spin-echo sequence with imaging parameters of TR = 2627.7 ms, TE = 36 ms, a slice thickness of 0.7 mm, a field of view of 20 × 20 mm, and a scanning time of 12.36 min were used. For PA imaging, the laser excitation pulses of 9 ns were delivered at five wavelengths (680, 715, 730, 760, 800, and 850 nm) in coronal orientation with a field of view of 20 × 20 mm, step sizes of 0.3 mm moving along the horizontal direction, and a scan time of 20 min. To validate these data, five natural anatomical landmarks were manually selected as registration targets ( Figure 3).

Murine Brain Data
For further evaluation of the proposed image registration method, we used ex vivo mouse brain data which was provided by Ren et al., in a previous study [67]. After removal of the mouse brain skull, the whole brain of the mouse was embedded in 3% agar in phosphate-buffered saline and was then imaged ex vivo. To acquire T2-weighted MR images of the mouse brain, a 2D spin-echo sequence with imaging parameters of TR = 2627.7 ms, TE = 36 ms, a slice thickness of 0.7 mm, a field of view of 20 × 20 mm, and a scanning time of 12.36 min were used. For PA imaging, the laser excitation pulses of 9 ns were delivered at five wavelengths (680, 715, 730, 760, 800, and 850 nm) in coronal orientation with a field of view of 20 × 20 mm, step sizes of 0.3 mm moving along the horizontal direction, and a scan time of 20 min. To validate these data, five natural anatomical landmarks were manually selected as registration targets ( Figure 3).

Inducing Brain Deformation
The proposed algorithm was designed to compensate for brain deformation during neurosurgery. Since the brain deformation is a complicated non-linear transformation, it is a challenging task to implement it physically on the phantom or mouse brain data. To evaluate our proposed registration algorithm, we performed brain deformation numerically by applying pre-defined pixel shifts to images. For this purpose, we used pre-operative and intra-operative MR images of brain tissue. The intra-operative MR image was considered as a gold standard. The deformation matrix was obtained by mono-modal registration of these images using the residual complexity algorithm [70] (Figure 4). Then, the obtained brain deformation matrix was applied on PA images of the brain phantom and mouse brain data.

Inducing Brain Deformation
The proposed algorithm was designed to compensate for brain deformation during neurosurgery. Since the brain deformation is a complicated non-linear transformation, it is a challenging task to implement it physically on the phantom or mouse brain data. To evaluate our proposed registration algorithm, we performed brain deformation numerically by applying pre-defined pixel shifts to images. For this purpose, we used pre-operative and intra-operative MR images of brain tissue. The intra-operative MR image was considered as a gold standard. The deformation matrix was obtained by mono-modal registration of these images using the residual complexity algorithm [70] (Figure 4). Then, the obtained brain deformation matrix was applied on PA images of the brain phantom and mouse brain data.

PA-MR Image Registration Framework
The workflow for automatic multi-modal image registration to compensate for the brain deformation was shown in Figure 5. After preparing two data sets, including brainmimicking phantom data and murine brain data, pre-deformation MR images were set as reference images, and pre-deformation PA images were set as float images. Then, a real

PA-MR Image Registration Framework
The workflow for automatic multi-modal image registration to compensate for the brain deformation was shown in Figure 5. After preparing two data sets, including brainmimicking phantom data and murine brain data, pre-deformation MR images were set as reference images, and pre-deformation PA images were set as float images. Then, a real brain deformation matrix which was achieved by the registration of intra-operative and pre-operative patient MR images using the residual complexity method was applied on PA images to generate deformed PA images. Then, by using the proposed registration method based on joint co-sparse analysis, registration of the MR image and deformed PA image was done. Finally, the image registration results were evaluated and visualized for brain shift calculation. To evaluate the registration algorithm, root mean square error (RMSE) was calculated for the phantom and mouse image registration. Additionally, target registration error (TRE) was calculated for defined targets in the phantom and mouse brain data. Furthermore, we used the Hausdorff distance (HD) between the PA and MR images. The HD between two point sets is defined as where d(., .) is the Euclidean distance between the locations and a smaller value of HD indicates a better alignment of the boundaries. To avoid the effect of outliers [71], we used 95% HD instead of maximum HD. Figure 4. (a) Pre-operative MR image, (b) intra-operative MR image, and (c) brain deformation field was achieved by registration of intra-operative and pre-operative MR images using residual complexity method.

PA-MR Image Registration Framework
The workflow for automatic multi-modal image registration to compensate for the brain deformation was shown in Figure 5. After preparing two data sets, including brainmimicking phantom data and murine brain data, pre-deformation MR images were set as reference images, and pre-deformation PA images were set as float images. Then, a real brain deformation matrix which was achieved by the registration of intra-operative and pre-operative patient MR images using the residual complexity method was applied on PA images to generate deformed PA images. Then, by using the proposed registration method based on joint co-sparse analysis, registration of the MR image and deformed PA image was done. Finally, the image registration results were evaluated and visualized for brain shift calculation. To evaluate the registration algorithm, root mean square error (RMSE) was calculated for the phantom and mouse image registration. Additionally, target registration error (TRE) was calculated for defined targets in the phantom and mouse brain data. Furthermore, we used the Hausdorff distance (HD) between the PA and MR images. The HD between two point sets is defined as where (. , . ) is the Euclidean distance between the locations and a smaller value of HD indicates a better alignment of the boundaries. To avoid the effect of outliers [71], we used 95% HD instead of maximum HD. The workflow for automatic multi-modal image registration to compensate for brain deformation. MR and PA images including pre-defined targets were set as a reference and float images, respectively. After applying brain deformation on PA images, registration of MR and deformed PA was conducted and evaluated.

Co-Sparse Analysis Model
Image (I) can be approximated via the sparse representation x ∈ R n , which is a linear combination of a few non-zero elements (named atoms) in an over-complete dictionary matrix D ∈ R n×k (n << k).
x ≈ Dα where α ∈ R k is a sparse vector with the fewest k non-zero elements. The sparse coefficients describe the salient features of the images. Therefore, the sparse representation problem could be solved as the following optimization problem:

of 15
Here, α 0 is the zero norm of α that represents the number of non-zero values in a vector (α). The sparse representation of an image considers that a synthesis dictionary represents the redundant signals.
There is also another representation of an image based on the co-sparse analysis model [60]. This alternative assumes that for a signal of interest (x), there exists an analysis operator Ω ∈ R k×n such that Ωx ≈ α as an analyzed vector is sparse for all x ∈ R n . The rows of Ω represent filters that provide sparse responses and indices of the filters with zero response determine the subspace to which the signal belongs. This subspace is the intersection of all hyperplanes to which these filters are normal vectors, and therefore, the information of signals is encoded in its zero responses. The index set of the zero entries of Ωx is called the co-support of x as below: As the key property of analysis sparse models, these models put an emphasis on the zeros in the analysis representation rather than the non-zeros in the sparse representation of the signal. These zeros in the analysis representation model inscribe the low-dimensional subspace which the signal belongs to. Consequently, analysis operator learning procedures find the suitable operator Ω for signal x as below: where Ω * is the optimized operator Ω. In order to relax the co-sparsity assumption, the logsquare function as a proper approximation of zero norm is used for large values of ν as below: where ν is the positive weight. Therefore, Equation (4) could be converted to One should consider that there are three main constraints on the Ω * to avoid trivial solutions as below [72]: 1.
The rows of Ω * have the unit Euclidean norm; Ω * ∈ obliquemani f old.

2.
The operator Ω * has full rank, i.e., it has the maximal number of linear independent rows. 3.
The rows of the operator Ω * are not trivially linearly dependent.

Multi-Modal Image Registration Algorithm
In this study, we formulated the multimodal image registration problem in terms of a co-sparse analysis model. There are different co-sparse models that can be used in multimodal image registration approaches [73]. In our approach, a joint analysis co-sparse model (JACSM) was proposed for the registration of PA and MR images. JACSM indicates that different signals from different sensors of the same scene form an ensemble. The signals in an ensemble include a common sparse component, shared between all of them, and an innovation component which represents individual differences [74].
Consider two images, I PA and I MR , which are provided through PA and MR imaging, respectively, from a brain simulated phantom as the input data. The interdependency of the two image modalities was modeled via JACSM and common sparse components were considered in this study. This image pair has a co-sparse representation with an appropriate pair of analysis operators, (Ω PA , Ω MR ) ∈ R k×n PA × R k×n MR . By considering the structures of images encoded in their co-supports based on Equation (3), there is a pair of analysis operators so that the intersection of the co-supports of Ω PA I PA and Ω MR I MR is large. In particular, we attempted to learn the pair of co-sparse analysis operators (Ω PA , Ω MR ) for two different image modalities.
On the other hand, the PA and MR images should be matched with a transformation T such that where x determines homogeneous pixel coordinates in PA images. The goal of multimodal image registration problem in this approach is to optimize T by using the pair of analysis operators (Ω PA , Ω MR ). We consider that, for an optimized transformation, there is a coupled sparsity measure to be minimized. Thus, by considering Equation (6) and constraints based on Equation (7), we are searching for T * such that To tackle the problem of Equation (9), we proposed the ADMM. In other words, the analysis operators were learned by optimizing a JACSM via an ADMM. The ADMM is a candidate solver for convex problems, breaking our main problem into smaller subproblems as below: where x ∈ R n , y ∈ R m , A ∈ R p×n , and B ∈ R p×m . The Lagrangian for augmentation Equation (10) can be written as where the term ρ is a penalty term that is considered positive and λ is the Lagrangian multiplier. Equation (11) is solved over three steps-x-minimization and y-minimization are split into N separate problems and followed by an updating step for the multiplier λ as follows: x k+1 := arg min x L p x, y k , λ k , y k+1 := arg min y L p x k+1 , y, λ k ,

Results and Discussion
To implement the proposed image registration algorithm, 20,000 pairs of square sample patches of size 7 pixels from the total images in the training set were randomly selected. It is notable that in our experiments, the patch sizes of 3, 5, 7, 9, and 11 pixels were applied. Based on our experience, a small patch size would cause an over-smoothing effect, and a larger patch size would lead to more computation. Therefore, based on our results, the patch size of 7 × 7 was selected to balance the two effects.
The performance of the JACSM-based registration method was evaluated using a phantom with simulated vessels and using ex vivo mouse brain data with anatomical landmarks. In Figure 6, the performance of the proposed registration method for PA-MR, US-MR, and MR-MR images on the phantom data is shown and compared. Also, for further evaluation the results of our proposed method were compared to the commonly used normalized mutual information (NMI) registration method. In the first row, the MR image and its corresponding US and PA images are shown. Dashed yellow circles show the same fields of view in three different modalities (MRI, US, and PA). Corresponding structures which were used to calculate target registration error are labeled with numbers 1 to 3 in the three imaging modalities. The brain deformation field was applied to the images in the first row, and the second row represents deformed MR, US, and PA images. As shown in Figure 6d-f, labeled targets have been displaced due to inducing deformation. Finally, the images in the third and last rows show the image registration results of MR, US, and PA after deformation (second row) with the original MRI before deformation (Figure 6a) using two different algorithms, NMI and JACSM, respectively. The results of the registration between the original MR and deformed MR, deformed US, and deformed PA using NMI as a common multimodal registration method are shown in the third row of Figure 6g-i, respectively. Also, the results of the registration between the original MR and deformed MR, deformed US, and deformed PA using our proposed method are shown in the last row of Figure 6j-l, respectively. The result of the registration between the original MR image and deformed MR image (Figure 6j) was used as a gold standard to evaluate the proposed algorithm. Also, the registration result of the deformed PA image (Figure 6l) was compared to the registration result of the deformed ultrasound image (Figure 6k) as a commonly used intra-operative imaging modality for brain shift compensation. As we have shown in the last row, images registered more accurately in the MR-MR image registration compared to the PA-MR image registration. Also, images registered more accurately in the PA-MR image registration compared to the US-MR image registration. As we have shown with the blue arrow in the last row of images, the surface of the phantom was matched accurately in the result of the MR-MR image registration. The registration of US-MR had the worst performance in matching the surface of the phantom in two modalities, and the registration of PA-MR had an acceptable performance in matching the surface of the phantom in two modalities, PA and MRI. Comparing the blue arrows in the third and last row images, it is clear that our proposed algorithm was more accurate than NMI. Also, white arrows in Figure 6i,l show that the PA-MR registration results for vessels were located in the depth of the phantom.
To quantitatively evaluate the proposed registration method, the RMSE, TRE, and HD for the PA-MR, US-MR, and MR-MR image registration were calculated and shown in Table 1. In total, we used 23 phantom data. The registration accuracy of MR and MR images was considered as a gold standard. The algorithms were implemented in MATLAB and tested on an Intel Core i7 3.2 GHz CPU with 8GB RAM.
The results of the phantom study showed that the PA-MR image registration had a better RMSE, TRE, and HD by about 60%, 65%, and 59%, respectively, compared to the US-MR image registration as a common imaging modality for brain shift compensation. On the other hand, the proposed method reached an RMSE of about 0.73 mm, which is acceptable in comparison with the MR-MR image registration as a gold standard, with an RMSE of about 0.62 mm. The proposed method improved the results of RMSE and TRE by about 60% and 63% (on average) compared to NMI.
For further evaluation of the proposed method, ex vivo mouse brain data were used. In Figure 7, the performance of the JACSM-based registration method for the PA-MR image registration for mouse brain data is shown and compared with the MR-MR image registration. Figure 7a,b represents MR and PA images of the mouse brain before any deformation, respectively. The PA image after applying non-linear deformation is shown in Figure 7c, and the registration result of the deformed PA and original MR of the mouse brain images is shown in Figure 7. The registration of MRI images before and after deformation was shown in Figure 7e as a gold standard. Also, in Figure 7f, the mean of the RMSE, TRE, and HD of the PA-MR image registration for all data of the mouse brain was calculated and compared to the result of the MR-MR image registration. To quantitatively evaluate the proposed registration method, the RMSE, TRE, a HD for the PA-MR, US-MR, and MR-MR image registration were calculated and sho in Table 1. In total, we used 23 phantom data. The registration accuracy of MR and M images was considered as a gold standard. The algorithms were implemented MATLAB and tested on an Intel Core i7 3.2 GHz CPU with 8GB RAM.  registration. Figure 7a,b represents MR and PA images of the mouse brain before any deformation, respectively. The PA image after applying non-linear deformation is shown in Figure 7c, and the registration result of the deformed PA and original MR of the mouse brain images is shown in Figure 7. The registration of MRI images before and after deformation was shown in Figure 7e as a gold standard. Also, in Figure 7f, the mean of the RMSE, TRE, and HD of the PA-MR image registration for all data of the mouse brain was calculated and compared to the result of the MR-MR image registration. The results acquired from the ex vivo mouse brain also proved the ability of the proposed registration method to recover non-linear deformation, with a calculated mean of the RMSE, TRE, and HD of 1.13, 0.98, and 0.85 mm, respectively. The results are acceptable The results acquired from the ex vivo mouse brain also proved the ability of the proposed registration method to recover non-linear deformation, with a calculated mean of the RMSE, TRE, and HD of 1.13, 0.98, and 0.85 mm, respectively. The results are acceptable when compared to the results of the MRI-MRI registration as a gold standard, with an RMSE, TRE, and HD of about 0.98, 0.85, and 0.77 mm, respectively. In fact, intra-operative PA imaging as a real-time imaging with about a 15% RMSE increase could be a good alternative to intra-operative MR imaging. Additionally, with a 60% improvement in registration accuracy, PA imaging could be an alternative for intra-operative ultrasound imaging. Generally, it cannot be concluded that PA imaging is an alternative to US imaging due to its insufficiency in providing structural and anatomical information. However, for brain shift calculation in neurosurgery, the blood vessel meshwork of brain tissue is considered as the main landmark during surgery, which is better visualized using PA imaging.
Having a closer look at the comparison between the synthesis and analysis models, the synthesis model contains very few low-dimensional subspaces and an increasingly large number of subspaces of higher dimension. In contrast, the analysis model includes a combinatorial number of low-dimensional subspaces with fewer high-dimensional subspaces. The co-sparse analysis models can yield richer feature representations, and joint co-sparse analysis models consider the common sparse components of different signals from different sensors. Therefore, the JACSM-based registration method was found to be more suitable for multi-modal image registration. Despite the promising results that were obtained for multimodal image registration based on the joint co-sparse analysis model, there are certainly limitations to adopting our proposed approach for other multimodal medical image registrations. The joint co-sparse analysis model is based on local assumptions and thus fails where large areas of one modality are not available (such as an existing gap in one of the modalities). It seems we could overcome this limitation by developing a co-sparse analysis model for each modality separately and proposing an optimized cost function in co-sparse space. This is something we will do in the future.
It is noteworthy that the quality of PA images also affects the registration accuracy. In our previous works, we also focused on improving the quality of PA images using advanced methods in image de-noising [75] and image reconstruction [21]. On the other hand, recently, there has been a growing interest in using low-fluence based photoacoustic imaging systems such as LED-based systems for guiding real-time interventions [76][77][78]. In fact, the development of methods to improve the quality of LED-based PA images [79] as well as advantages such as high frame rates and low-cost imaging have made PA imaging promising to achieve a higher SNR compared to the system used here for the purpose of brain shift compensation.

Conclusions
There has been a growing interest in intra-operative imaging approaches to update the pre-operative images with real-time data when tissue deformation occurs during surgery. In particular, accurate and real-time brain shift compensation remains a challenging problem during neurosurgery. For the first time in this study, we proposed the application of PA imaging as an interventional solution during neurosurgery in combination with pre-operative modalities such as MRI to track brain deformation. However, the accurate combination of PA and MR images requires the development of a real-time and robust image registration algorithm. Accurate registration of intra-operative PA images with preoperative MR images of brain tissue could calculate and compensate for brain deformation. In this study, the JACSM-based registration, which can capture the interdependency of two modalities, was proposed for the PA-MR image registration. The proposed algorithm works based on the minimization of mapping transform by using a pair of analysis operators in PA and MR images which are learned by the ADMM. The algorithm was tested on two data sets of phantom and mouse brain data and the results showed a more accurate performance for PA imaging versus US imaging for brain shift calculation. Furthermore, the proposed method showed about a 60% improvement in TRE in comparison with the common NMI registration method. The co-sparse analysis models can yield richer feature representations and better accuracy for medical image registration in real-time processes, which is crucial for surgeons during neurosurgery to compensate for brain shift. Finally, by using this JACSM-based registration, the intra-operative PA images could become a promising tool when the brain shift invalidates pre-operative MRI.