Galindo-Cardiel

Fernández-Jiménez

Badiola Jj

Novel group A rotavirus G8 P[1] as primary cause of an ovine diarrheic syndrome outbreak in weaned lambs

Rotavirus is a worldwide major cause of diarrhea outbreaks in neonatal ruminants. An outbreak of ovine diarrheic syndrome (ODS) in 50-75 days-old lambs (weaned lambs) is described. Fecal immunochromatography and intestinal immunohistochemistry for rotavirus group A were performed. In addition, semi-nested multiplex RT-PCR for G

A c c e p t e d M a n u s c r i p t 2 and P rotavirus genotyping in combination with sequencing were performed, to support the diagnosis and identify the viral strain. A novel ovine rotavirus group A G8 P [1] strain was determined as the main cause of the ODS observed, whereas other pathogens were ruled out.

Keywords: Ovine diarrheic syndrome; rotavirus group A; weaned lambs

The economical management of ovine and other ruminant farms is seriously compromised by postnatal (less than 10 days-old) mortality associated to several infectious agents such as rotavirus [START_REF] Holland | Some infectious causes of diarrhea in young farm animals[END_REF][START_REF] Muñóz | Prevalence of neutralizing antibodies to 9 rotavirus strains representing 7 G-serotypes in sheep sera[END_REF], Dahma et al, 2009). The rotavirus genus is divided into seven antigenically distinct groups (A to G) [START_REF] Desselberger | Rotaviruses and rotavirus vaccines[END_REF]. They are classified by sequencing of the major structural protein VP6, which have a role in viral transcription [START_REF] Dhama | Rotavirus diarrhea in bovines and other domestic animals[END_REF]. The protective immunity obtained by neutralizing antibodies is elicited by two outer capsid proteins, VP4 and VP7, which act as major virulence factors. Both proteins are associated with P and G serotype specificity, respectively [START_REF] Desselberger | Rotaviruses and rotavirus vaccines[END_REF][START_REF] Dhama | Rotavirus diarrhea in bovines and other domestic animals[END_REF]. VP4 helps the virus to attach to host cells and VP7 provides the protective outer capsid layer [START_REF] Desselberger | Rotaviruses and rotavirus vaccines[END_REF]). Sheep are one of the few species in which group B rotaviruses appear to be more prevalent in neonatal lambs than group A rotaviruses (Fitzgerald et al, 1995).

Rotavirus infections have a very short incubation period and affected animals can manifest profuse diarrhea and severe dehydration. Moreover, concurrent infection with secondary pathogens may increase clinical severity [START_REF] Dhama | Rotavirus diarrhea in bovines and other domestic animals[END_REF].

Experimental infections with rotavirus in lambs have been performed [START_REF] Van Pinxteren | Characterisation of the primary local and systemic immune response in gnotobiotic lambs against rotavirus infection[END_REF], but descriptions of natural rotavirus infections are scarce and results are not A c c e p t e d M a n u s c r i p t 3 always conclusive [START_REF] Tzipori | The clinical manifestation and pathogenesis of enteritis associated with rotavirus and enterotoxigenic Escherichia coli infections in domestic animals[END_REF][START_REF] Moon | Comparative histopathology of intestinal infections[END_REF][START_REF] Van Metre | Diagnosis of enteric disease in small ruminants[END_REF]. In this report, we describe a new ovine group A rotavirus G8 P[1] strain (OVgAR) associated to an outbreak of ovine diarrheic syndrome (ODS) in 50 to 90 days old Rasa Aragonesa lambs. Immunochromatography (IC), immunohistochemistry (IHC), and a multiplex reverse transcriptase-polymerase chain reaction (RT-PCR) were performed to support the diagnosis and to characterize the new viral strain.

An outbreak of acute ODS was reported in November 2008 in an high sanitary status farm of Rasa Aragonesa sheep in Zaragoza (North-eastern Spain) under a semiintensive 3 lambing/2 years production system. All affected lambs belonged to the September 2008 lambing period and they were housed in pens up to 1 square meter per animal. No ODS outbreaks were noticed in the farm or in neighborhood flocks during the 5 years previous to the present outbreak. No management or feeding practices changes were recorded in the year previous to the disease outbreak.

A total of 192 out of 203, 50 to 60 day-old lambs (95%) were mild to severely ODS-affected in an outbreak which lasted 10 days in total. Mortality associated with the ODS was 17% (35 out of 203 lambs). A total of 25 ODS-affected and 5 apparentlyhealthy lambs were studied as follows. Clinical severity was scored by using the following parameters: (1) attitude/ lamb activity (0 = apparently healthy, 1 = slightly decreased, 2 = moderately decreased, 3 = severely decreased), (2) rectal temperature (0 = below 40 C, 1 = above 40 C) and (3) diarrhea (0 = not present, 1 = present). These partial scores were summed up, yielding a total cumulative clinical index that reflected apparent health (score = 0), the presence of mild (score = 1-2), moderate (score = 3-4) or severe (score = 5) clinical severity. Fecal swabs were cultured (aerobes, anaerobes and fungi) and a zinc sulphate flotation test for larvae and oocysts were performed.

Seven protozoa-positive fecal samples were speciated (sporulation test) for Eimeria spp.

A c c e p t e d M a n u s c r i p t S.A.U, Girona, Spain) was performed in feces following manufacturer´s instructions.

Clinical severity and pathogen detection results in these 30 lambs are shown in Table 1.

Five lambs severely affected by ODS were necropsied and samples of tissues were formalin fixed and paraffin-embedded. Consecutive sections of each paraffinembedded block were processed for haematoxylin&eosin and for IHC. Main pathological diagnosis consisted in segmental (proximal) moderate to severe catharral enteritis (5/5 lambs) with diffuse severe epithelial stunting and blunting (5/5), diffuse mild to moderate epithelial apical vacuolation (5/5), multifocal mild epithelial cryptal hyperplasia (3/5) and diffuse mild lacteal dilatation (3/5). No enteric bacteria adhered to microvilli were observed in any of these lambs. Moreover, diffuse severe mesenteric lymphoid hypertrophy was observed (5/5). Other secondary lesions recorded were multifocal moderate subacute necrotizing hepatitis (2/5), and multifocal mild subacute catharral-purulent bronchopneumonia (2/5). Feces from the 5 necropsied lambs were OVrGA-positive by IC. In addition, Cryptosporidium parvum was detected in two lambs and no oocysts of Coccidia spp. were cultured.

The small intestines from all necropsied lambs were processed for VP6 rotaviral antigen by IHC. Briefly, OVrGA antigen retrieval was performed using EDTA Buffer (1mM EDTA, 0.05% Tween 20, pH 8.0). Polyclonal rabbit anti-rotavirus (Dako, ref.

A0082) served as the primary antibody (1:250), peroxidase-conjugated avidine complex (ABC peroxidase staining kit, Thermo Scientific, ref 32020, 1:100) was utilized as the visualization system and 3,3′-diaminobenzidine (DAB) was used as the chromogen.

Negative control procedures included omission of primary antibody. PCR-and IHCpositive and negative intestinal tissues from lamb and pig were used as positive and negative controls respectively. OVrGA-positive diagnosis was assessed based on the Isolated Escherichia coli strains from the 5 ODS-affected necropsied lambs were serotyped and genotyped for enterotoxigenic (LT, STa, STb, K88, P987, K99 and F41), verotoxigenic (VT1 and VT2), atypical enteropathogenic (eae) and adhesine F17 virulence factors (ETEC, VTEC, aEPEC and F17EC) as described elsewhere [START_REF] Rey | Serotypes, phage types and virulence genes of shiga-producing Escherichia coli isolated from sheep in Spain[END_REF]. Serotyping and genotyping results of the detected Escherichia coli strains are presented in Table 2.

Two frozen OVrGA-positive fecal samples (by IHC and IC) from necropsied lambs were used for genotyping the rotavirus strain. RNA was extracted by using RNAwiz TM reagent (Ambion), following the manufacturer's instructions. Reverse transcription with random primers and PCR with rotavirus VP6 gene-specific primers were performed [START_REF] Gouvea | Polymerase chain reaction amplification and typing of rotavirus nucleic acid from stool specimens[END_REF]. Both fecal samples were VP6-positive PCR (type I2). G-and P-genotyping was carried out by semi-nested multiplex PCR with G-and Ptypes oligonucleotide primers previously described [START_REF] Gouvea | Polymerase chain reaction amplification and typing of rotavirus nucleic acid from stool specimens[END_REF][START_REF] Gentsch | Identification of group A rotavirus gene 4 types by polymerase chain reaction[END_REF][START_REF] Iturriza-Gómara | Diversity within the VP4 gene of rotavirus P[8] strains: implications for reverse transcription-PCR genotyping[END_REF] and the OVrGA isolate was identified as a G8 P[1] strain. One PCR amplicon per each isolate (3 in total) were purified with the QIAquick Gel Extraction Kit (Qiagen) and sequencing was performed by the dideoxynucleotide chain termination method with the ABI Prism BigDye terminator cycle sequencing reaction kit (Applied Biosystems) on an ABI Prism 3100 automated sequencer (Applied Biosystems). The sequencing was performed with the forward or reverse primers used for the RT-PCR. BLAST analysis was carried out with the VP6 (GenBank Accession Sequence alignment was performed with Clustal X (version 2.0) and the phylogenetic relationships of the strain described in the present study with other reference rotavirus strains were assessed employing the software MEGA version 4.0 using the neighbor-joining method [START_REF] Tamura | MEGA4: Molecular Evolutionary Genetics Analysis (MEGA) software version 4.0[END_REF]. The robustness of the taxonomic relationships was evaluated by bootstrap analysis (1,000 replications).

Reference strains for different VP7 G-, VP4 P-and VP6 I-genotypes were selected from GenBank (Figure 1). SAS 9.1.3 software was used (SAS Institute Inc., Cary, NC, USA) to perform the ANOVA test. Clinical severity and pathogen detection were considered as variables.

Significant association between group A rotavirus-infection and presence of ODS was the only association detected (p = .026, significant at p < .05; Table 1).

This ODS outbreak significantly correlates with the detection of the newlydescribed rotavirus strain whereas other pathogens isolated did not. Based on this correlation and the quick spread of the clinical process, high morbidity, limited mortality, compatible pathology and pathogen detection (IC, IHC, RT-PCR), the rotaviral etiology is strongly supported as the primary cause of this outbreak, following accepted evidence-based veterinary medicine approach [START_REF] Gay | Determining cause and effect in herds[END_REF].

The phylogenetic analyses confirmed the VP6, VP4 and VP7 genotypes identification previously obtained by RT-PCR, sequencing and BLAST analysis.

Interestingly, 96% identity was found between the OVrGA VP6 gene sequence, a rotavirus group A sable antelope G6P[14] (Accession No. FJ495131.1) and a sheep rotavirus G8P[14] strain previously isolated in Zaragoza (Spain) (Accession No.

EF554152.1) [START_REF] Ciarlet | Genomic characterization of a novel group A lamb rotavirus isolated in Zaragoza, Spain[END_REF][START_REF] Matthijnssens | Are human P[14] rotavirus strains the result of interspecies transmissions from sheep or other ungulates that belong to the mammalian order Artiodactyla?[END_REF]. A 90% of identity A c c e p t e d M a n u s c r i p t 7 between OVrGA VP4 gene sequence obtained and the VP4 sequenced genes selected for the analysis (Figure 1) was found. Moreover, OVrGA VP7 gene sequence obtained and the sheep rotavirus G8P[14] showed the higher inter-sequences identity founded in this study, a 98%. VP4 and VP7 are the main antigenic determinants for rotaviruses, but only one of the OVrGA-positive isolate was positive for VP4 gene. The different sensitivity of the primers used in the semi-nested multiplex PCR performed was considered to be the cause of such findings [START_REF] Gouvea | Polymerase chain reaction amplification and typing of rotavirus nucleic acid from stool specimens[END_REF][START_REF] Gentsch | Identification of group A rotavirus gene 4 types by polymerase chain reaction[END_REF].

However, genomic diversity within the VP4 gene of other rotavirus strains has been related with failures of VP4 amplification using RT-PCR [START_REF] Iturriza-Gómara | Diversity within the VP4 gene of rotavirus P[8] strains: implications for reverse transcription-PCR genotyping[END_REF].

The main difference was that these described strains were isolated from animals which were into the described range of age for rotavirus infection (i.e., less than 30 days-old) [START_REF] Holland | Some infectious causes of diarrhea in young farm animals[END_REF][START_REF] Dhama | Rotavirus diarrhea in bovines and other domestic animals[END_REF]. On the contrary, the strain described in this paper mainly affects 50-75 days-old lambs (weaned lambs). However, few ovine strains have been isolated and characterized to date [START_REF] Ciarlet | Genomic characterization of a novel group A lamb rotavirus isolated in Zaragoza, Spain[END_REF][START_REF] Chen | Sequence and evolutionary analysis of VP1 gene of ovine rotavirus NT[END_REF] possibly due to the lack of routine surveillance for rotaviral infections in lambs.

Further epidemiological studies should be carried out to clearly demonstrate the pattern of variation between rotaviral strains.

Escherichia coli strains isolated in this study were not considered as the most likely ODS-related primary etiology due to the lack of significant correlation with the clinical process and to the high variability among the isolates. This phenomenon has already been pointed out by others authors (Nagy and Fekete, 2005). However, Escherichia coli strains forming attaching-effacing lesions (eae-positive) can cause diarrhea due to their intimate association with intestinal epithelial cells [START_REF] Moon | Comparative histopathology of intestinal infections[END_REF] but these lesions have not been observed in our work. However, its role as possible complicating pathogen should be studied in further studies. Other pathogens detected in A c c e p t e d M a n u s c r i p t 8 this study (Cryptosporidium parvum, Coccidia spp., Campylobacter spp.) did not demonstrate to be the most likely ODS-outbreak related primary pathogen, although it is assumed that they likely cause severe clinical complications [START_REF] Holland | Some infectious causes of diarrhea in young farm animals[END_REF][START_REF] Van Metre | Diagnosis of enteric disease in small ruminants[END_REF].

In conclusion, epidemiological, pathological and laboratory evidence supports that a novel group A rotavirus G8 P[1] strain described in this paper was the main cause of an ODS outbreak in 50-75 days-old lambs (weaned lambs). Further experimental work using rotavirus and E. coli in lambs is being carried out by our group at the moment.

  A c c e p t e d M a n u s c r i p t 5 amount and distribution of antigen in positive-stained apical enteric epithelium and macrophages of the lamina propia. All samples were OVrGA-positive IHC. The number of OVrGA-positive cells found by IHC was directly related with the severity of the individual enteritis.

  ), VP4 (Accession No. FN675939) and VP7 (Accession No. FN675937) gene sequences obtained.

FigureFigure 1 .

 1 Figure

Figure

  Figure

Table 2 .

 2 Serotyping and genotyping of Escherichia coli strains isolated in five severe ODS-affected necropsied lambs.

	Lamb Strain	Serotyping*	eae #	F17 ‡	Reference number §
	1	S1	ONT	+	-	eae-12073
	2	S2	ONT	+	-	eae-12074
	3	S3	-	-	+	F17-CNL
		S4	ONT	+	-	eae-12075
	4	S5	O23	-	+	F17-12076
		S6	O4 (1/1280), O15 (1/640)	+	-	eae-12077
	5	S7	-	-	+	F17-CNL
		S8	O4 (1/640), O15 (1/640)	+	-	eae-12078
	* Serotyping: Results given as serotype detected (ONT, O4, O15 and O23) with serotype dilution added
	in brackets when it is relevant.

# eae: intimin gene responsible for attaching effacing lesions ‡ F17: adhesine fimbriae 17 gen § Reference number given by the Reference Laboratory of Escherichia coli at the University of Santiago de Compostela in Lugo (Spain)

Table A

 A

	FN675939
	M a n u s c r i p t
	c c e p t e d

Acknowledgements

The authors thank M. Pérez and E. Salazar for their support and observations in the pathological analysis of the data.


Sean Jy-Shyang Chen

Ingerid Reinertsen

Pierrick Coup

Charles X B Yan

Laurence Mercier

Rolando Del Maestro

Louis Collins

Sean Jy

Shyang Chen
email: <sjschen@bic.mni.mcgill.ca

Validation of a hybrid Doppler ultrasound vessel-based registration algorithm for neurosurgery

Keywords: Adult, Aged, Algorithms, Brain, blood supply, surgery, Female, Humans, Intraoperative Period, Magnetic Resonance Angiography, Male, Middle Aged, Neuronavigation, methods, Neurosurgical Procedures, Phantoms, Imaging, Preoperative Period, Ultrasonography, Doppler Volumetric registration, Brain-shift, Doppler ultrasound, Digital Phantom, Validation

Purpose

We describe and validate a novel hybrid non-linear vessel registration algorithm for intraoperative updating of preoperative magnetic resonance (MR) images using Doppler ultrasound (US) images acquired on the dura for the correction of brain-shift and registration inaccuracies. We also introduce an US vessel appearance simulator that generates vessel images similar in appearance to that acquired with US from MR angiography data.

Methods

Our registration uses the minimum amount of preprocessing to extract vessels from the raw volumetric images. This prevents the removal of important registration information and minimizes the introduction of artifacts that may affect robustness, while reducing the amount of extraneous information in the image to be processed, thus improving the convergence speed of the algorithm. We then completed 3 rounds of validation for our vessel registration method for robustness and accuracy using (i)a large number of synthetic trials generated with our US vessel simulator, (ii)US images acquired from a real physical phantom made from polyvinyl alcohol cryogel (PVAc), and (iii)real clinical data gathered intraoperatively from 3 patients.

Results

Resulting target registration errors (TRE) of less than 2.5mm are achieved in more than 90 of the synthetic trials when the initial % TREs are less than 20mm. TREs of less than 2mm were achieved when the technique was applied to the physical phantom, and TREs of less than 3mm were achieved on clinical data.

Conclusions

These test trials show that the proposed algorithm is not only accurate but also highly robust to noise and missing vessel segments when working with US images acquired in a wide range of real-world conditions.

Introduction

Modern neurosurgery relies heavily on computer assistance and image guidance to provide topological and locational information to surgeons, enabling them to accurately navigate within the patient s brain. Known as image guided neurosurgery (IGNS), these techniques ' typically use the preoperative image of the patient s head acquired through magnetic resonance imaging (MRI) for surgical planning and ' navigation. Significant errors may be introduced into surgical navigation from image misregistration or through intraoperative deformation of brain tissues in neurosurgery, also known as brain-shift. Displacements from brain-shift can vary greatly between 5 50mm on the cortical surface 3,8,15,18, 22 24,35,34,36 . However, even prior to dural opening, displacements caused by brain-shift can nevertheless be [ -] quite significant. [START_REF] Maurer | Investigation of intraoperative brain deformation using a 1.5-t interventional mr system: preliminary results[END_REF] 15 have shown that even with the dura intact, the magnitude of brain-shifts observed was as large as

[ ]
3.4mm pushing outwards and 8.1 mm collapsing inwards.

These tissue displacements, both and after dural opening, are quite significant and enough to render guidance using the before preoperative magnetic resonance (MR) image useless for surgical use, even after accurate initial patient to image registration. In order to compensate for such displacements in the soft tissue, the preoperative MR image can be corrected by registering it to an intraoperative image, which could be acquired using either intraoperative computed tomography (CT) 13 , intraoperative MR 14,[START_REF] Nimsky | Intraoperative compensation for brain shift[END_REF][START_REF] Roche | Rigid registration of 3-D ultrasound with mr images: a new approach combining intensity and gradient information[END_REF] or

[ ] [ ]
intraoperative ultrasound (US) imaging [START_REF] Lange | Vessel-based non-rigid registration of MR/CT and 3D ultrasound for navigation in liver surgery[END_REF][START_REF] Reinertsen | Clinical validation of vessel-based registration for correction of brain-shift[END_REF] .

[ ] US is a good candidate for intraoperative imaging. However, images acquired using US may be difficult to interpret by untrained users. This is due in part to its relatively small imaging field-of-view, its non-standard anatomical views, speckling in its acquired images, and the need for a craniotomy in order to do intraoperative brain imaging. Nevertheless, US imaging has many advantages that overcome these shortcomings. For instance, US imagers have the benefit of being significantly less expensive, capable of real-time imaging, and more portable than either MR and CT imaging devices. As well, US imaging gives good resolution and contrast in soft tissues without the need to compromise surgical ergonomics or exposing the patient to ionizing radiation when compared to intraoperative MR and CT, respectively. Intraoperative US images can be used effectively to correct and update the preoperative MR images to the intraoperative state of the imaged tissues.

The challenge of registering images between US and MR modalities is that they have different image contrasts, noise patterns, and artifacts. This is due to the fact that the two imaging modalities are acquired using very different physical principles, with MR images depending on proton spin interaction and US images depending on the material s acoustic properties. Two types of strategies are used in ' the literature to overcome the challenge of registering the images:

Salient anatomical features that can be imaged on both MR and US such as sulci 5 , tumour outlines 21 , or blood Feature-based: This strategy typically employs image processing methods to relate the intensities of MR to US. For instance, Intensity-based: gradient information in the MR image 37 may be used directly for registration, or the preoperative MR image may be segmented and [ ] processed to generate a set of pseudo-ultrasound images, which allows for better image correspondence and registration 1 . Mutual

" " [ ]
information could also be used to relate the intensities found in MR and US 16 .

[ ] Blood vessels are important and useful features for the registration of deforming soft tissues due to several of their properties. First, vessels are easily identifiable and also usually well distributed in most surgical regions-of-interest (ROI). Second, many image processing methods exist to extract vessels from the acquired images. For instance, the moving blood in a blood vessel is easily imaged using Doppler US and can be implicitly segmented from the image background based on colour saturation alone. As well, vasculature from different types of MR angiography can be extracted preoperatively using segmentation algorithms such as that of Descoteaux et al. [START_REF] Descoteaux | A multi-scale geometric flow for segmenting vasculature in mri : Theory and validation[END_REF]

and Frangi

[ ] [ 11 . Finally, the vessel s unique branching and stable relative paths and positions to surrounding tissues allow for high specificity in ] '

registration. Therefore as vessels move and change shape with deforming tissues undergoing brain-shift, the correction of these vessel movements enable the detection and correction of the corresponding tissue deformations.

To date, blood vessels have been successfully used for linear and non-linear registration of deforming soft tissues such as the cerebral cortex 25,33 and the liver 19 . From previous work in our lab, brain shift compensation using vascular images from Doppler US and MR

[ ] [ ]
was characterized and shown to be highly accurate 33 .

[ ] Purely feature-based registration algorithms are fast and accurate techniques that typically rely on image preprocessing to parameterize and reduce the complexity in the raw data. For vessel registration, this preprocessing is typically done through skeletonization, which reduces and parameterizes the vessel image into a set of vessel centerpoints. This process dramatically reduces computation time for registration since the amount of data that needs to be processed is often greatly reduced.

However, such methods can remove important information for registration, such as vessel lumen diameter or voxel intensity differences in the vessels images themselves. As well, thin, low contrast vessels and important vascular features such as branching and continuity can be eliminated altogether. Removal of this information introduces ambiguities that allows skeletonized branches of thick vessels to be incorrectly matched with thin vessels. Likewise, centerpoint fragments from different vessels can be incorrectly matched to a single vessel. As well, the extracted center-points of densely situated vessels can be reduced to featureless point clouds, making purely iterative closest point (ICP) algorithms 2 fail and rendering the registration process impossible. For these reasons, purely feature-based

[ ]
registration can produce poor results on clinical data.

Indeed our previous experience using ICP-based techniques yielded good results using vessel traces that were relatively sparse 34 , [ ] with distances between the vascular branches being greater than the magnitude of misregistration. However our recent experience with more complex anatomies containing many more blood vessels resulted in failed registrations.

Intensity-based registration methods are very accurate and quite robust to noise due to the amount of information available for use in alignment. As such, they are widely used for atlas, inter-subject, and various inter-modality registration tasks. Although more sophisticated intensity-based methods may be developed to improve processing time, intensity-based registration methods capable of correcting free-formed non-linear deformations required to properly register cerebral vasculature are quite computationally intensive and can take several minutes for images of modest size. For this reason, their use is limited to offline registration tasks since they are not suitable for intraoperative registration, which should ideally complete in at most 2 3 minutes during the time needed to prepare for tissue dissection.

-These issues inherent to feature or intensity-based registration methods, necessitate the investigation and creation of a robust vessel registration method capable of overcoming the reliance on extensive preprocessing to allow for the registration of MR and US vessels with a high degree of accuracy in clinically acceptable time.

Motivation for Hybrid Registration

In this paper, we describe an algorithm for a hybrid feature and intensity-based non-linear registration of MR to US images through cerebral blood vessels using ANIMAL 4 for registration. In contrast, our registration method attempts to eliminate the background by extracting only the vessel information while at the same time minimizing other forms of preprocessing to retain as much of the original vessel information as possible. To this extent, we use the whole volumetric vessels from the angiographic intensity images instead of relying on discrete points or lines produced through skeletonization of vessel data. In the following sections, we present our hybrid non-linear volumetric registration technique developed for the registration of US vessels acquired on the dura with MR vessels. Our goal for hybrid registration was to be able to automatically and robustly correct tissue displacements up to 20mm; a threshold that is more than twice the maximum recorded brain-shift tissue displacement of 8.1mm prior to dura opening 15 .

[ ]

Validation data types

Our registration method was then validated using three types of data with increasing realism:

Multiple synthetic US data sets were generated numerically from a clinical MR Synthetic data generated from a digital phantom: image. This enables us to evaluate the behavior under controlled conditions where the ground truth is known. However, the simulated data are not completely realistic and do not account for all possible imaging situations.

Real MRI and US data were acquired from a PVA phantom. This set of experiments Real data acquired from a physical phantom: permits evaluation of the algorithm with an object having a known gold standard, which produces images that are more realistic. However, the spatial distribution of the image signal is quite simplistic and does not accurately replicate data gathered in clinical scenarios.

Real preoperative MRI and intraoperative US data were acquired from 3 patients. This Real clinical data acquired from 3 patients: is the most realistic set of data and represents real surgery and clinical imaging. However the image registration results are difficult to evaluate since the ground truth is not known exactly.

Overview

We will describe in detail our hybrid registration algorithm in Section 2. We will then describe the tool we created and images we acquired to validate the robustness and accuracy of our registration algorithm, including the digital phantom (Section 3.1), the physical PVA-based phantom (Section 3.2), and clinical data from 3 patients (Section 3.3). The methodology for validating our registration algorithm is detailed in Section 3.5 followed by the result of the validation in Section 4. We then conclude with a discussion of our test results and future work in Section 5.

Hybrid Registration

Our hybrid registration method takes the raw MR and US clinical images and preprocesses them to extract vascular features while retaining the intensities of the vessels in volumetric form. This combination of intensity and feature-based cues are used for hybrid registration. The source image for registration consists of a 3D contrast enhanced MR vascular imaging processed with Frangi s vesselness ' filtering 11 , while the target image consists of an intraoperatively acquired and reconstructed 3D Doppler US volume with the volumetric [ ] vessels extracted from the background tissues by using a colour saturation threshold. Detailed descriptions of the preprocessing done to the target and source images can be found in Sections 2.1.

Using the minimally processed vessel data, registration proceeds in two phases, with multiscale linear registration prior to non-linear registration. We noticed that although a large amount of registration error can be corrected with rigid-body transformations, the 6 degrees of freedom of the transform do not allow for an adequate correction of all the non-rigid displacements observed between the preoperative and intraoperative images. It is for this reason that we register the images using full affine transformations after an initial round of rigid-body registration for a first alignment, and then subsequently using non-linear transformations to correct the remaining brain-shift not captured by the affine linear transformation.

Preprocessing

The first step in our registration method is to preprocess our source and target images to extract their vasculature. Our source images were T1 weighted Gadolinium MR images that were processed using Frangi s vesselness method 11 , which effectively enhances tube and

' [ ]
vessel-like structures in the image. The parameters used for the filtering are , and at 0.5, 0.5, and 0.5 ( ) and five α β γ ×max HessianNorm filter scales with values 0.4, 0.74, 1.1, 1.64, and 2.5mm.

Our target images were 2D power Doppler US images, which were fused with the raw B-mode US images. The latter were filtered by removing all pixels under a colour saturation threshold. The filtered 2D images were then reconstructed in a 3D volume using the tracking information contained in each image. This was done by nearest-neighbour regridding of each pixel in the ultrasound image slices to a 3D volume with the same resolution as the MR image in a similar fashion as that of 40 . Although other interpolation methods could have [ ] been used rather than the basic nearest-neighbour regridding method, we found that the other methods increased the 3D reconstruction time and did not perceptibly improve registration accuracy. A thorough review of 2D ultrasound image reconstruction techniques and their considerations can be found in 41 .

[ ]

Linear Registration Phase

The linear registration algorithm takes as input the aforementioned source and target intensity volumes which have been preprocessed for vessel features, then blurs them with a Gaussian kernel prior to each round of linear registration, using a smaller kernel for each subsequent round. We used the reconstructed US vessel volume as the registration target image and the MR vessel volume as the source image.

The starting estimates for registration depend on the method used for data acquisition. For our clinical cases (Section 3.3), the starting estimates are based on the transforms acquired from the facial skin landmark registration performed by a surgeon. In the case of the physical phantom data (Section 3.2), the initial registration is done by point-based registration on the container of the physical phantom. In any case, the transforms puts the acquired US images into the same coordinate system as the MR image and provides the initial position for registration.

Our linear registration is optimized using a downhill simplex algorithm 26 with cross-correlation as the objective function. We found [ ] that it performed effectively and quickly for registering our images of similar image contrast: both processed MR and US images have bright vessels on dark backgrounds. Although other objective functions such as mutual information may be used, we found in our preliminary tests that these methods were significantly slower and did not perform better than cross-corrrelation on our pre-processed vessel images. In fact, observations from our preliminary test are in line with that of 31 , who suggested that registration with mutual [ ] information does not perform well in images with (1) thin structures such as the vessels in retinal images or (2) for cross-modality registration of MR to US images.

We performed a total of three rounds of optimization with sequentially smaller blurring kernels and sampling grids with successively refined linear transforms. Due to the size difference between the US target and the MR source volume, cropping was automatically performed on the MR volume around the location of the US target to improve and speed up the calculation of registration transforms.

Cropping for the first round of registration is set to 50mm, since this is the largest value of brain shift recorded in the literature. Parameters defining the full-width-half-max (FWHM) of the blur (where the FWHM relates to standard deviation by 20 , the sampling grid

[ ]
used, the crop margin, and the type of linear transform used in each round are described in more detail in Table 1.

Non-Linear Registration Phase

Following the linear registration, which provides a robust global alignment, we used the ANIMAL registration algorithm 4 to

[ ]
compute a set of non-linear transformations that corrects for any local vessel deformations that cannot be accounted for by the linear transformations. ANIMAL hierarchically estimates 3D vector field at different scales, mapping deformations progressively from larger to smaller scales and produces a globally non-linear deformation field with point-to-point correspondence between the two volumes.

The affine linear registration transforms recovered from the previous linear registration phase are applied to the preprocessed volumetric intensity source image, which is then used as the source image for this non-linear registration phase of our hybrid registration method. We performed two rounds non-linear registration, once at lower resolution, and once at higher resolution using ANIMAL with cross-correlation objective function, and relatively high stiffness to maintain smoothness in the deformation (See Table 1). The resulting transformation is applied in each round of non-linear registration only if the final value of the objective function is reduced from the starting value.

Through the iterations of non-linear registration, we are able to attain correspondence between the US and the MR vessel images based on their grey scale voxel values. Even though the vessels are not sparse, there exists a significant amount of tissue in between each vessel where the ANIMAL algorithm must interpolate or extrapolate the deformations. This is due to the fact that the algorithm performs registration and regularizes using only the non-zero vessel voxels of the images, and ignores the zero-valued voxels where the non-vascular tissues have been removed. To address the limits of the local elastic regularization used in ANIMAL and to ensure a continuous smooth deformation throughout the tissue, we apply the thin-plate spline (TPS) transform 9 on points extracted from the

[ ]
aligned US and MR vessels to regularize the transformation over the entire image.

One thousand points are randomly sampled from overlaps of the vessels from the two non-linearly registered volumes and inversely transformed using the deformations recovered by ANIMAL. This produces 1000 pairs of homologous points, which are then used as the control points to define the displacements on an interpolating TPS.

Validation Methods

Digital Phantom

As part of the work to validate our hybrid registration algorithm, we developed a non physics-based digital phantom to synthesize images that simulate the appearance of Doppler US vessel images acquired from a brain distorted by brain shift. Qualitative visual observations were made between the vessels acquired using power Doppler US and MR imaging to enable and successfully imitate real power Doppler US images from MR vessel data. We found that the US vessels in comparison to MR vessels tend to be:

This is likely due to the fact that the US operator increases signal gain when acquiring the Doppler US image. This is done to Thick: allow smaller and more distant vessels to be captured in the image.

The surfaces of the acquired US vessels are uneven due, in part, to speckling in US images and artifacts arising from flash Uneven: artifacts caused by sudden probe movements.

Both noise and US image artifacts can be due to the random scattering and reflection of US pulses in Noisy and containing artifacts:

the different tissue being imaged, as well as the noise introduced by the electronics of the US imaging equipment. Some vessels are dimmer or entirely not present in volumes of one modality versus the other (See Fig. 1). This could be due Presence:

to imaging artifacts such as occlusion, shadowing, and the angle of the US probe in respect to the imaged blood flow.

Our digital phantom simulates the appearance of a Doppler US vessel image by taking an MR vessel image as input. The workflow of the digital phantom is shown in Figure 2. The first step to producing the synthetic US vessel involves extracting and processing an MR angiographic image to act as the ground truth image. The second step is to generate a set of random linear and nonlinear transforms to act as the registration ground truth and then applying the transform to an MR vessel image. The third step involves altering the transformed MR vessel image to exhibit the observed characteristic appearance of the US vessels listed above. Finally, noise and speckling is added, and the simulated US volume is masked into the shape of an US sweep. These steps are explained in detail below.

To show that the simulated US vessel volume behaves similarly to a real US vessel volume in MR vessel registration, we also compared the similarity of the cross-correlation objective function curves to each other. The first step in creating the simulated US vessels involves extracting the vessel from a T1 weighted Gadolinium MR image in the manner described in Section 2.1. The extracted MR vessel image is then thresholded at a low intensity to capture anything that may be a vessel and then cropped down to a ROI.

Random Transformations

Parameters for the linear transforms are chosen randomly within translations from 30 to 30mm, rotations from 5 to 5 degrees, -scalings from 0.9 to 1.1, and shears from 0.9 to 1.1, with the fixed transformation point located in the center of the ROI bounding box. The non-linear transforms are generated using a TPS since this method produces realistic deformations that are smooth and regularized. Nine pairs of control points are used to define the TPS transform, with the two points of each pair randomly located in relation to each other at distances between 0 5mm. The first 8 point pairs are located on the corners of the bounding box around the ROI that has been scaled down by 50 to be inside the ROI. The last point pair is positioned in the center of the ROI bounding box (See Fig. 3). The rationale for this % placement scheme is to localize the control points of the thin-plate spline in the center of the ROI where the anatomical structures are located, thereby increasing the local nonlinear deformations there. This set of linear and non-linear transformations are used as the gold-standard for validating our registration method, beyond validation are also applied to the processed MR vessel image to simulate tissue displacement from brain shift.

Simulation of Characteristics

Once the MR vessels have been deformed using the set of random transforms, they are used to simulate the appearance of Doppler images. The simulation occurs in two steps. First, the thick uneven surface and variable presence of the US vessels are simulated through the extensive use of Perlin noise 29 . Afterwards noise and speckle are added onto the simulated vessels.

[ ] Perlin noise is a procedural texturing method commonly used in computer graphics to generate many types of continuously varying pseudo-random textures with specific spatial frequencies. Random gradients are first defined at regularly spaced intervals in a 3D lattice.

These gradients are used to determine the intensity of the field between the nodes of lattice through linear or cubic interpolation. This simple technique allows us to generate random appearing textures of different scales and spatial frequencies that can easily be used as parameters for varying the appearance of generated US vessels. A more detailed explanation of Perlin noise and the technique for generating it can be obtained in 29 .

[ ]

We used the intensity values of a Perlin noise texture with low spatial frequency as the parameter for determining whether the vessel is present in the image (vessel presence) by varying the contrast of a vessel to the background. The noise texture has a value between 0 and 1, which makes the vessel completely invisible and fully visible, respectively. This roughly simulates artifacts from shadowing or from angle dependant modulation of the US signal. Although power Doppler US (signal power) is not as sensitive to probe angle as colour Doppler US (blood flow rate) the probe angle to a vessel may nevertheless alter the intensity of the power signal.

Perlin noise of high spatial frequency is used to determine the radius of a particular voxel s dilation on a visible vessel. The dilation is ' performed using a 3D isotropic Gaussian point spread function with a FWHM of between 3 5mm. This higher spatial frequency Perlin noise, which has a value between 0 and 1, is used to determine extent of the vessel dilation. The value of 0 indicates that the voxel will be dilated to 3mm and 1 indicates that the voxel will be dilated to 5mm. The combined uneven dilation of the original vessels produced a set of resulting vessels that are thick and bumpy. When the resulting dilated vessel volume is thresholded and windowed for intensity, a set of disjoint vessels with uneven surfaces are created, which simulates appearance-wise the increased gain of the power Doppler US imaged vessels and any amplified noise.

Accumulating Noise and Speckling

Various types of artifacts and noise commonly seen in US images and volumes are applied to the generated vessel volume as the final the real US vessels in their thickness, uneven surfaces, noise, artifacts, and their variable presence (see Fig. 1).

Behavior of Simulated Vessels in Registration

To ensure that the simulated US vessels from our digital phantom performed similarly to real Doppler US reconstructed vessels in registration, we measured their cross-correlation values through translations, rotations, and scaling. The simulated US vessels and their corresponding real US vessels, were all registered and aligned to the MR angiography data. The gold-standard transform for aligning the real US image to the MR image was found by first calculating the affine linear transform between manually picked homologous anatomical landmarks between the US and MR image in the manner described in Section 3.3.

All the data sets were blurred with Gaussian kernels of 5 voxel FWHM and the two sets of US vessels were translated, rotated and scaled around their registered position on the MR vessel image. From the test results in Figure 4, we see similar profiles for the cross-correlation objective function for both sets of US data. The correlation of the objective function curves was greater than 0.9 for translation and more than 0.8 for rotation and scaling. This indicates that there is high similarity between the registration performance of our simulated US vessels compared to that of the real US vessels. As well, the generated transforms and distortions from our digital phantom also vary smoothly and are of similar magnitude as that found in clinical cases in the literature. From this, we can conclude that simulated US vessel volumes model the appearance of clinical US vessel images sufficiently well for testing our MR to US registration algorithm. plastic tubing with inner lumens of 1.57, 2.36, and 3.18mm, capable of circulating fluids for simulating blood vessels. The phantom also 7 21 has an inflatable catheter located near its bottom center to simulate brain deformation 32 . The T1 MRI data set from this phantom [ ] consists of 2 image volumes for each catheter inflation of 0, 5, and 10ml of water, giving a total of 6 volumes (See Fig. 5). Three sets of 2 Doppler US images for phantom inflations of 0, 5, and 10ml were also acquired.

PVAc Physical Phantom

The MR images were acquired with a Siemens SonataVision 3T using T1 weighted MR imaging at a voxel resolution of 1mm 1mm × 1mm. The tubes in the phantom were manually segmented from the MR images.

×

The Phantom US images are acquired with an ATL HDI 5000CV US imager using a ATL P7-4 phased-array probe on power Doppler mode. A passive target with markers spheres were attached to the US probe and optically tracked by a Polaris tracker (Northern Digital Incorporated). The ultrasound probe was calibrated with a z-fiducial phantom described in 12 . The US probe is swept smoothly and [ ] steadily by hand over the water-immersed phantom to minimize gaps in the reconstructed volume. The 2D Doppler images are filtered using a colour saturation threshold and reconstructed using the same methods as described in Section 2.1.

Manually labelled homologous landmark points, consisting of 21 fluid filled bubbles, located and distributed evenly throughout the phantom tissue material, were found in each of the 6 physical phantom MR volumes. These landmarks were used to validate the transformation of our algorithm. In 2 separate point picking trials to determine intra-rater variability, the mean difference between the points around each landmark was approximately 0.7mm with a maximum distance of 1.1mm.

Clinical Data Acquisition and Vessel Extraction

Vessel structures in preoperative MR and intraoperative US images must be properly acquired and extracted prior to registration. To acquire the preoperative image, the patient can be scanned using MR angiographic techniques, such as time-of-flight, T1-weighted MPR [ 42 , or contrast enhanced MR angiography.

]

In our institute, open cranial surgery involving such procedures as tumour resection almost always involve Gadolinium (Gd) contrast enhanced MR imaging for enhancement of the pathology, usually within a week prior to their procedure. This also provides us with a contrast enhanced MR vascular image for use in computer assisted neurosurgical guidance. The images were acquired using a GE Signa Excite 1.5T with double-dose Gd enhanced T1 weighted MR imaging at an in-plane resolution of 0.5 0.5mm and a slice thickness of × 1mm. The MR vessels are extracted by filtering the angiographic MR images using Frangi s vesselness method with the same parameters ' as those described in Section 2.1.

We note that although MR angiographic or contrast enhanced MR imaging may not be routinely performed in all institutes for tumour surgery, we nevertheless believe that the development of improved image processing methods may allow blood vessels to be extracted from standard T1-weighted images thus allowing for the use of a technique such as ours.

During the stabilization of the patient s head for operation, the pre-operative MR image was rigidly registered to the patient s head by ' '

the surgeon through manually picking nine facial landmark point pairs. These point pairs were chosen between the skin surfaces of the actual patient and the segmented patient MR image identified using optically tracked pointer and a mouse pointer, respectively. The point pairs consist of the patient s lateral and medial canthus on the left and right eye, the tragus and tragus valley on the left and right ear, as ' well and the most posterior part of the nasal bridge.

The clinical intraoperative US images and their tracking data were acquired and processed in the same manner as the physical phantom data to build the US volume from 2-dimensional slices and extract its vessels (see Section 2.1). The probe was swept in a free-handed manner over the dura of the patient after the craniotomy in a smooth and steady manner to minimize gaps in the reconstructed volume.

Manually labelled landmarks for anatomical features found between the MR volume and the US volume were used as a silver-standard to evaluate registration success. These landmarks consist of blood vessel bifurcations, conjunctions of sulci, and unambiguous distal portions of ventricles located throughout the volume from 4 sets of real clinical data acquired from 3 patients (See Table 2). Depending on the size of the volumes and the availability of landmarks, between 7 and 18 unambiguous points were labeled by an expert. Four separate point picking trials were done to determine intra-rater variability. The mean difference between the points at each landmark was 1.3mm

with the maximum distance of 2.1mm.

Given the clinical nature of these images and their acquisition conditions, it was difficult and took significant amounts of time to accurately identify unambiguous landmark points in both the US and MRI images. This process involves cross-checking between the preoperative and intraoperative images, re-verifying the landmarks, and correcting any inaccuracies, which can take around 15 30 minutes per case or around 2 4 minutes per landmark point. This was especially true for US acquisition 3, where the sweep amounted to a thin slab less than 2 cm in width and only 7 unambiguous landmark points could be indentified.

TRE Calculations

We use the target registration errors (TRE) as a quantitative measure to test and validate our registration methods on our digital phantom trials, physical phantom trials and also our clinical trials.

: The TRE were calculated using an initial set of points, which were evenly sampled from throughout the ROI Digital phantom trials delimited sweep cone of the US images. The initial points were then transformed using both the gold-standard transformations and the recovered transformations from our technique. The root mean square (RMS) errors between the latter two sets of points were then calculated to give us our TRE.

: The TRE were calculated using the labelled homologous fluid filled bubbles in the physical phantom images of Physical phantom each catheter inflation. The transformations recovered through registration of the source MR vessel image to the target US vessel image of a different inflation was applied onto the points picked from the source MR image. The RMS errors between transformed source points and the set of points picked from the target inflation image was then calculated, giving us our TRE.

: The TRE were calculated through a set of silver-standard point pairs, which were manually picked from homologous Clinical trials intra-cranial anatomical landmarks in the source MR image and the target US image. The transformations recovered using our registration technique were applied onto the source points and compared with the locations of the target points through computation of the RMS error as TRE.

Registration Validation

The hybrid registration algorithm was evaluated with the digital phantom (Section 3.1), the PVAc Phantom (Section 3.2), and 4 sets of real clinical data (Section 3.3). Although it is common in the literature to report translation and angular errors, we elected to use TRE for validation of our registration method since it can effectively encompass the displacements incurred from the full affine transformation ( 12degrees-of-freedom) used in our linear registration phase. This includes scaling and shear transformation on top of the translational and rotational transforms used in rigid-body registrations.

The TRE measure is also useful since it allows us effectively compare the quantities of error corrected between our linear and non-linear registration phases. Validation of the registration were done in the following manner for the digital phantom, the physical phantom, and the clinical trials datasets: : Using the digital phantom, we created 2100 registration trials, each with a randomly generated vessel instantiation, Digital phantom transform, and deformations which were used as the gold standard. TREs greater than 2.5mm are considered as failed registration attempts.

: Validation using the physical phantom was performed by registering each of the 6 MR volumes with each of the Physical phantom 3 US volumes to estimate the non-linear transforms. This was done to account for deformations generated from inflating with 0, 5, or 10ml of water.

: Registration validation of the clinical trial was accomplished using manually defined landmark points used as the Clinical trials silver-standards that were chosen from homologous intra-cranial anatomical landmarks in the source MR image and the target US image.

Results

Digital Phantom

We found that our registration method was quite robust and accurate in correcting the simulated brain-shift of our digital phantom, with a high success rate given the 2.5mm threshold criteria. In applying our full registration method (linear and non-linear phases) to the digital phantom we found that our method was successful for 99 , 94 , 92 , and 90 of the digital phantom trials at 5, 10, 15 and 20mm

% % % %
initial misregistration, respectively (see Fig. 6, Non-Linear).

It is interesting to note that the linear step completed in the first phase of registration can significantly reduce the initial TRE. In the same trials as above, when initial misalignment RMS error is less than 5mm, 10mm, 15mm, and 20mm, the linear registration results in TREs less than 2.5mm for 87 , 68 , 65 , and 64 of the trials, respectively (see Fig. 6, Linear). These lower values are expected since % % % % linear transforms cannot completely correct for the local non-linear deformations we generated for our digital phantom. Furthermore, if the linear registration completed with a TRE of less than 6mm in these trials, the non-linear registration was successful more than 96 of the % time, correcting the TRE to below 2.5mm.

Given the initial design criteria of 20mm maximum displacement on the dura, only 10 of the trials failed with TREs greater than % 2.5mm. We found that only 3 of the registration trials attained a greater TRE than the initial alignment (see Fig. 7). The mean reduction % in TRE through non-linear registration from the linear registration results is 0.62mm and quite significant, with 1.73 10 in a p = × 64 two-tailed t-test. This significant reduction in TRE can also be seen in that the mean TRE after linear registration was 2.29 0.62mm with ± the first and third quartile being 1.84mm and 2.66mm, respectively, which reduced after non-linear registration to 1.64 0.47mm with the ± first and third quartile being 1.32mm and 1.82mm, respectively. In general, we found that the regions of the US volume with the least blood vessels had the highest TREs while regions surrounded by vascular structures had consistently low TREs.

We also tried to examine the behavior of the vessel registration algorithm when we push it beyond the 20mm operating range. We see that in this case the success in the non-linear registration was only 84 , 76 , and 70 for 25mm, 30mm, and 35mm initial % % % misregistration.

PVAc Phantom

In our tests with the physical PVAc phantom, our non-linear registration algorithm was able to correct the deformations in all our trials to TREs 2mm or under (See Table 3). This correction corresponds to a reduction of 8 to 45 in TRE remaining from the linear % % registration step after non-linear registration, with the TREs decreasing in all cases in these tests on the physical phantom. This demonstrates that the algorithm is likely converging towards the correct solution.

It should be noted that since the chosen landmarks on the PVAc phantom are spread out evenly over the entire volume of the phantom, the actual amount of non-linear displacement an their corresponding correction through registration has been diluted in calculation of the TRE. To illustrate this, we used the results of the 0ml to 5ml registration trial and calculated the TRE from 5 landmark points located in a region of the phantom closer to the catheter balloon. From this, we noted a reduction of TRE from 2.4mm to 1.3mm, which equals to a 46 reduction in the TRE of this region from our non-linear registration alone. The differences in the magnitude of non-linear deformations % throughout the image can be seen in Figure 8. The degree of non-linear deformation was as high as 3.01 mm, in the region close to the center of the imaged volume where the deformation inducing catheter resides.

Images showing a 5ml inflated phantom corrected with the recovered transformations from registration with a 0ml inflated phantom can be seen in Figure 9. Note that there is almost a complete overlap of the vascularized cerebral region of the phantom images, including the more highly deformed parts of the phantom near the catheter.

Clinical Data

The initial registration of the MR to the patient s head attained from the facial landmarks for the clinical data was achieved with the ' standard homologous landmark point-pair method described in Section 3.3. Only one starting point was evaluated for each clinical case, since experiments using the synthetic phantom data have indicated the algorithm s robustness to the starting location. These initial mean ' TREs from facial registration using the manually identified silver standard landmarks on the four datasets were 7.25mm, 3.51mm, 4.61mm, 3.742mm. It is important to note that although our surgeons sought to keep the root mean square errors of their facial landmark fiducial registration errors (FRE) lower than 1.5mm, the resulting TREs from the initial skin-based registration are nevertheless quite large.

The magnitude of the TREs is likely due to the configuration of skin landmarks in relation to the surgical target. If the geometric center of gravity of the fiducials is not close to the surgical target, the landmark configuration may lead to a lever effect, accentuating the TRE. Such a decorrelation between the measured FRE from the measured TRE has been noted in the past 6,38,

[ ] After the US to MR registration process consisting of the linear and nonlinear registration steps, the TREs were reduced to 2.38mm, 1.23mm, 3.38mm, and 1.70mm on these clinical datasets with improvements of TRE by 67.2 , 64.9 , 26.7 , and 54.5 , respectively.

% % % %
Resulting images from registration of Patient 1 (See Table 2) can be seen in Fig. 10 and 11, which are respectively the registered vessel images and the T1 MR image corrected with the recovered transformation. Resulting registered images for Patient 4 in Fig. 12 show the T1 MR image corrected with the recovered transformation. It can be seen that the amount of non-linear deformations from datasets 1 and 4

are not negligible, with 15.2 and 11.2 of the TRE uncorrected by affine-linear registration (See Table 4). The magnitude of non-linear % % deformations in Patient 4 was as high as 2.5 mm in some local regions of the brain, and can be seen in Fig. 13.Note that this deformation was recovered from our vessel based registration after completion of our linear registration phase, from images acquired prior to dural opening. This indicates that full affine-linear registration alone is insufficient to capture and recover all the brain-shift induced deformations even prior to dural opening.

The average execution time for vessel registration is roughly 6 8 minutes on a 2.4 GHz Intel Q6600 processor. In this time span, the the time taken up by US preprocessing and volume reconstruction, which takes less than a minute.

Discussion

The validation results indicate that the hybrid registration algorithm outlined in this paper can be applied to correct for registration errors due to inaccurate fiducial registration and to update preoperative images to account for brain shift.

Validation 10 21
Our method can successfully correct for brain shifts with displacements as large as 20mm in more than 90 of the cases (See Fig. 6 % and 7) in our digital simulations. As well, we have found in our clinical cases that we can significantly reduce the registration TRE with either only linear registration or both linear and non-linear registration (See Table 4). However, we note that these results need to be interpreted with care given that our digital simulation is not entirely physically realistic and that we have evaluated the technique only on 3 patients in our clinical trials. We plan to apply the method to a larger cohort of patient data in the future to further validate and determine the generalisability of the results. We also intend to improve our methods for simulation of US blood vessels to make them more physically realistic. We also note that while our clinical cases used specified facial landmarks for initial registration, it is also possible to use a random iterative closest point method such as that of 10 for initial rigid body registration, before applying our hybrid method for [ ] non-linear soft tissue registration.

Our digital phantom occasionally removes large regions of vessel segments in its random vessel presence simulation. On some of these simulated US images, our hybrid registration algorithm fails catastrophically, with the resulting TREs significantly higher than 2.5mm even when the initial TREs are less than 20mm. We realized that the removed vessel segments was the cause of the failure when we qualitatively examined several of the generated US digital vessel phantom images that resulted in registration failures. The simulated US volumes were indeed significantly sparser due to vessel removal and as such they aligned themselves to the wrong vessels in the larger sized MR vessel volume image. When the same images were not processed with the random presence portion of the simulation and ran with the same initial misalignments, registration of these same images converge and succeed even with the added noise, artifacts, and large vessel dilations.

The reason for the catestrophic failure in registration is clear. Since the initial vessel images used in digital simulation is already quite sparse, masking and vessel removal can further weaken the global minimum to the point as to cause the algorithm to be trapped at one of the many local minima. This introduced ambiguity in effect increased the possible alignments of the smaller and sparser US image with the large amount of vasculature in the MR vessel image. The problem was further exacerbated when the misalignment distance of the images was greater, which required the algorithm to traverse across more local minima and increased the chances that it would settle into any one of them before reaching the global minimum. Should the vessel removal be extensive enough, the global minimum for registration could be eliminated, causing the registration algorithm to never converge to the correct solution. Although these large failures in alignments are not favorable, there are nevertheless benefits to such a behavior since it allows technicians and surgeons to easily recognize registration failure through visual inspection, which limits the chance that a wrongly registered image is accepted by the surgeon and used for guidance.

We noticed a reduction in the TRE in all the physical phantom registration trials, and although the values were not below 1mm, visual inspection of the registration results show marked improvements in the alignment of the blood vessel tubes in the phantom (See Fig. 9). This is similarly true with the registration results using clinical data, where the corrections showed good alignment not only in vessels (see Fig. 10) but also in anatomical structures such as the ventricles and sulci (see Fig. 11 and12). Although anatomical structures like ventricles can be useful in registration, they are not as topologically distinct as branching vessels and are not well distributed throughout the brain tissue. Nevertheless, structures such as ventricles and sulci can be used as additional registration features in conjunction with vessels in registration. While we realize we cannot generalize with 4 clinical data sets, these tests do demonstrate that our algorithm can be successfully used to register real clinical images acquired intraoperatively during an operation.

Hybrid Vessel Registration

Ideally, we would have liked to compare our techniques to other US to MR registration techniques such as of 16 and 37 , however

[ ] [ ]
these techniques are dissimilar to our own in that (1) they recover linear transformations only and do not recover non-linear deformations in their registration, (2) nor do they take advantage of the vessel information in their images for image registration. As well, these techniques are not publicly available and trying to both implement and run them ourselves might introduce bias. Nevertheless, we have run several tests against Reinertsen et al. [START_REF] Reinertsen | Validation of vessel-based registration for correction of brain shift[END_REF] , which show that our technique, although slower, is more accurate, more robust, and has a higher

[ ]
probability of converging to the correct solution when run on real clinical data and image data from our digital phantom. We intend to continue investigating the use of our technique after and during tumour resection since these situations will likely be more challenging due to the higher tissue distortion and tissue removal. In both cases, the amount of resected tissues will be quite large, but since our registration technique relies on blood vessels for registration and not the removed tissue and can deal with missing vessels segments we believe that it should work well in such situations.

While the total time of 6 8 minutes required for hybrid registration is perhaps pushing the limits of clinical acceptance, this work was done with the intention of evaluating a technique to first determine its functionality and robustness on clinical data prior to further optimization. We believe that optimizing the code for the algorithm as well as adapting it to use hardware acceleration can reduce the registration time required. However, since the present algorithm is highly robust for misalignments up to 20mm (more than 90 success, % See Fig. 6), it should be possible to account for large absolute displacements by progressively estimating smaller relative displacements during the procedure. For example, the deformation estimated using the algorithm before dural opening will be used as the starting transformation for the subsequent deformation estimations during or after a surgical procedure. While this method remains to be tested with additional clinical data, we believe each application of registration algorithm could correct for the continual deformations of 10-15mm that occur throughout a procedure.

In this work, we concern ourselves solely to problem of registering US data acquired on the dura. We believe that by being able to correct brain-shift first on the dura, we can subsequently evaluate our registration algorithm in surgery before proceeding to tackle the potentially more difficult problems of image registration after dural opening and tissue resection. We realize that the characteristics of brain shift and the performance of the algorithm may be altered depending on the size and location of the craniotomy or the presence of tumour tissues, however this is outside the scope of our paper and will be addressed in the future research.

Conclusion

In conclusion, we have presented a new method for correction of preoperative MR images using intraoperative US for non-linear hybrid vessel registration. Since our vessel registration algorithm does not depend heavily on image preprocessing for robust and accurate results, it is able to avoid the pitfalls of feature-based registration methods that relies on their preprocessing methods such as vessel centerline registration algorithms. The minimum processing to extract the vessels allows us to complete registration in less time than purely intensity-based registration techniques.

In this work, we also developed a digital phantom capable of synthetically generating volumes with the appearance of Doppler US vessel images. These images are processed from a MR vessel image and have known ground truths for testing the registration robustness and accuracy of various non-linear MR to US vessel registration algorithms.

Our registration algorithm was then validated for robustness and accuracy through multiple trials using synthetic US vessel images generated from our digital phantom, along with US and MR images scanned from a real physical phantom, and finally 4 sets of real clinical data from 3 patients consisting of preoperative MR vessel images and intraoperative Power Doppler US images. Results from these tests show that the estimated transforms recovered by our novel registration algorithm are not only accurate but also robust when processing images that contain large quantities of noise, artifacts and missing vessels commonly seen in clinical vessel images. PVAc Phantom with 5ml inflation (Left) and its corresponding magnitude of deformation when non-linearly registered to the 0ml inflation phantom (Right). The magnitude of non-linear deformation is highest around the center of the phantom where the inflation catheter resides with the maximum recovered deformation of 3.01mm and the minimum of 0.008mm.

Fig. 9

PVAc Phantom with 5ml inflation image (red) being registered to its 0ml inflation image (green). The image initial alignments shows considerable misalignment. After linear registration, alignment is greatly improved but a large amount of misregistration can be seen in the left side of the image and close to the inflatable catheter. Following non-linear registration, most of the phantom is well aligned. 4.

Fig. 11

Registration results using clinical data, with overlays of the dual-mode US (Red) on the Gadolinium-enhanced T1 MR images (Green). The images show a coronal slice plane of the US volume and the same slice plane through the misaligned MR volume prior to (Left) and after registration (Right). The TREs for this subject (Patient 1 in Table 4) before registration is 7.46mm and 2.88mm afterwards. Note the improved alignment of the sulci and ventricle borders.

Fig. 12

Registration results using clinical data from Patient 4 (see Table 4), with the non-overlayed Gadolinium-enhanced T1 MR images (Left) and overlays of the hot-metal toned dual-mode US on the MR images (Right) in both coronal (Top) and sagittal planes (Bottom) prior to and after registration. The TRE before registration is 3.75mm and 1.7mm afterwards. The arrows on the non-overlayed registered MR image indicate areas of improved alignment at the vessels, sulci, and ventricular borders through the use of our registration method.

Fig. 13

Images from Patient 4 (see Table 4) showing a coronal slice from the reconstructed US image (Left), the Doppler US vessel image (Middle), and the magnitude of non-linear deformation recovered after the linear registration phase (Right). This reveals that significant non-linear deformation occurs even prior to dural opening. The magnitude of non-linear deformation recovered from MR to US registration was very high close to the Doppler vessels with the maximum recovered deformation of 2.5mm and the minimum of 0.006mm. Note that the deformation field can be high in the borders of the image where there are no vessels to lock-down the amount of deformation in these " "

regions.

Table 4

Full registration TREs for our 4 sets of clinical data. The first row shows the TRE from intra-operative patient facial registration, the second row shows the results of our method after the linear registration phase, the third row shows the final results of our method after the linear and non-linear registration phases, the forth and fifth rows show the percent reduction of TRE from the facial registration to the linear and both linear and non-linear registration phases, respectively. The sixth row shows the additional percentage of TRE reduction following the non-linear registration after linear registration. The last row indicate the number of landmarks used to determine TRE on each patient dataset.

Data sets

  [START_REF] Lange | Vessel-based non-rigid registration of MR/CT and 3D ultrasound for navigation in liver surgery[END_REF][START_REF] Nakajima | Use of cortical surface vessel registration for image-guided neurosurgery[END_REF][START_REF] Penney | Registration of freehand 3d ultrasound and magnetic resonance liver images[END_REF][START_REF] Reinertsen | Validation of vessel-based registration for correction of brain shift[END_REF] are used as alignment cues. These anatomical features are extracted and parameterized prior to their use in[ ]registration.

  step to simulate the US vessel volume. Speckle artifacts generated via the methods of Pizurica et al. 30 and grainy Perlin noise are [ ] multiplied to the expanded vessels. Gaussian noise and high intensity voxels from Perlin noise of high frequency are added to the vessel volume. The volume is then masked with a mask generated from the coordinates and orientation of a real US sweep to produce a wedge shape volume similar to a reconstructed US vessel volume. The resulting simulated US vessels qualitatively show a close resemblance to

  MRI and Doppler US vessel data of a polyvinyl alcohol cryogel (PVAc) physical brain phantom described and acquired by Reinertsen et al. 34 were used to validate our algorithm. The phantom is made of 3 types of PVAc of different hardnesses, containing 3 coils of [ ]

-

  multiscale linear registration takes approximately 3 4 minutes and non-linear registration takes approximately 3 5 minutes, with the rest of --

Fig. 1

 1 Fig. 1 Three different cross-section views of the actual T1 Gado-enhanced MR image (MR image), the actual ultrasound volume (Actual US), and the simulated ultrasound volume (Simulated US) imaged in the same brain region.

Fig. 2

 2 Fig. 2 Workflow of the Digital Phantom showing each step of the Doppler US appearance simulation process. Images used and outputted in each processing step are also shown.

Fig. 3

 3 Fig. 3 The placement of TPS control point pairs in relation to a region-of-interest (blue dashed box) for the creation of non-linear distortion. Eight control point pairs (gray) are located at the corners of the 50 scaled (red box) version of the bounding box (black box) surrounding the % region of interest. The ninth point pair (gray) is positioned in the center of the bounding box.

Fig. 4

 4 Fig. 4Cross-correlation curves graphs of the simulated US (red) versus that of the real US (blue) when the vessels are translated (top row), rotated (middle row), or scaled (bottom row). The horizontal axes of the translation graphs, rotational graphs and the scaling graph are labeled in millimeters of translational displacement, radians of rotation, and magnitude of scaling about the center of gravity of the US vessel volume, respectively. All vertical axes of the graphs are cross-correlation values. Cross-correlation between the translation curves in the x, y, and z directions are 0.91, 0.98, and 0.99, respectively. Cross-correlation between the rotation curves in the x, y, and z directions are 0.95, 0.82, and 0.88, respectively. Cross-correlation between the scale curves is 0.86.

Fig. 5 T1

 5 Fig. 5 T1 MR images of the PVAc Phantom inflated with 0, 5, and 10ml of water in the implanted catheter, which can be seen at the bottom center of each image.

Fig. 6

 6 Fig. 6Digital phantom registration percent success on test cases of different initial misalignment RMS errors using the linear and non-linear phases of our registration algorithm. Registration trials with RMS errors of less than 2.5mm are defined as successful.

Fig. 7 64 -Fig. 8

 7648 Fig. 7Digital phantom registration results showing the TRE for all the trials after the nonlinear phase of the registration process. The horizontal axis indicates the initial TRE prior to correction while the vertical axis indicates the final TRE after non-linear registration. The blue circles are the results after non-linear registration and the pink X s are the results after linear registration. The black diagonal line shows where a registration ' trial resulted in the reduction in TRE (point is below the line) or whether it increased the TRE (point is above the line). The green horizontal line shows our cutoff for success (2.5mm). The mean reduction of 0.62mm TRE from linear to non-linear phases of registration is highly significant, with < 2 10 in a two-tailed t-test. P ×

Fig. 10

 10 Fig. 10 Registration results showing improved alignment of the preoperative MR vessels (Green) with the intraoperative US vessels (Red). The images show a transverse slice plane of the US volume and the same slice plane through the misaligned MR volume before (Left) and after (Right) applying the transformations recovered using our hybrid registration method. Notice that the MR vessels previously misaligned and out of plane to the US vessels are now well registered. The images are from Patient 1 in Table4.

  . Previous non-linear feature-based vessel registration algorithms by Reinertsen et al. 33 and Lange et al. 19 preprocess their acquired images by segmenting the vessels and skeletonizing them to a set of discrete vessel centerpoints

	[ ]	[ ]

[ ]

Table 1

 1 Parameters used in multiscale linear registration (Top) and multiscale non-linear registration with ANIMAL (Bottom)

	Round	Blur FWHM	Sampling grid	Crop margin	Transform
	1st	3.3mm	2.5mm	50mm	Rigid-body
	2nd	2.2mm	1.8mm	16mm	Affine w/Isotropic Shear
	3rd	1.5mm	1mm	10mm	Affine w/Anisotropic Shear
	Round		Lattice Diameter	Sampling grid	Rigidity
	1st		15mm	5mm	0.9
	2nd		10mm	2.5mm	0.8

Table 2

 2 Patient information and US scanning parameters for our 4 sets of clinical data

	US Scan

Table 3

 3 TREs resulting from applying the non-linear registration algorithm on 6 linearly registered physical phantom MR images (Final) and the initial TREs (Initial). The MR images were acquired by inflating the physical phantoms twice ( , ) with 3 different volumes of water (0, 5, 10ml). The transforms are computed by registering each physical phantom MR source image with the remaining phantom US 1 2target images of dissimilar inflation volumes. Source images are listed as rows while target images are listed as columns in the table.

	0ml

[START_REF] Letteboer | Brain shift estimation in image-guided neurosurgery using 3-D ultrasound[END_REF]

[START_REF] Letteboer | Brain shift estimation in image-guided neurosurgery using 3-D ultrasound[END_REF]

[START_REF] Letteboer | Brain shift estimation in image-guided neurosurgery using 3-D ultrasound[END_REF]

[START_REF] Letteboer | Brain shift estimation in image-guided neurosurgery using 3-D ultrasound[END_REF]

[START_REF] Letteboer | Brain shift estimation in image-guided neurosurgery using 3-D ultrasound[END_REF]

[START_REF] Letteboer | Brain shift estimation in image-guided neurosurgery using 3-D ultrasound[END_REF]

[START_REF] Letteboer | Brain shift estimation in image-guided neurosurgery using 3-D ultrasound[END_REF]

[START_REF] Letteboer | Brain shift estimation in image-guided neurosurgery using 3-D ultrasound[END_REF]

) 514 -528 1998 ; 16 . Ji S , Wu Z , Hartov A , Roberts D , Paulsen K . Mutual-information-based image to patient re-registration using intraoperative ultrasound in image-guided neurosurgery . Medical physics . 35 :

-4624 2008 ;

[START_REF] Letteboer | Brain shift estimation in image-guided neurosurgery using 3-D ultrasound[END_REF]

[START_REF] Letteboer | Brain shift estimation in image-guided neurosurgery using 3-D ultrasound[END_REF]

[START_REF] Letteboer | Brain shift estimation in image-guided neurosurgery using 3-D ultrasound[END_REF]

[START_REF] Letteboer | Brain shift estimation in image-guided neurosurgery using 3-D ultrasound[END_REF]

[START_REF] Letteboer | Brain shift estimation in image-guided neurosurgery using 3-D ultrasound[END_REF]

[START_REF] Letteboer | Brain shift estimation in image-guided neurosurgery using 3-D ultrasound[END_REF]

[START_REF] Letteboer | Brain shift estimation in image-guided neurosurgery using 3-D ultrasound[END_REF]

Acknowledgements:

We would like to thank Marta Kersten and Lara Bailey for proofreading the manuscript. We are also grateful to the Natural Sciences and Engineering Research Council of Canada (NSERC), the Canadian Institutes of Health Research (CIHR), and Le Fonds Quebec de la Recherche sur la Nature et les Technologies (FQRNT) for funding this project.


S F Eskildsen

J V Manjón

P Coupé

V Fonov

D L Collins

Superresolution improves MRI cortical segmentation with FACE

Introduction

Brain cortical surface extraction from MRI has applications for measurement of gray matter (GM) atrophy, functional mapping, source localization and preoperative neurosurgical planning. Accurate cortex segmentation requires high resolution morphological images and several methods for extracting the cerebral cortex have been developed during the last decade [START_REF] Dale | Cortical surface--based analysis: I. Segmentation and surface reconstruction[END_REF][START_REF] Kim | Automated 3--D extraction and evaluation of the inner and outer cortical surfaces using a Laplacian map and partial volume effect classification[END_REF][START_REF] Eskildsen | Active surface approach for extraction of the human cerebral cortex from MRI[END_REF]). In many studies, the resolution of the morphological image acquisition sequence is chosen to be relatively low (~1mm 3 ) due to time and equipment constraints. To improve segmentation accuracy, such low resolution images can be upsampled using various interpolation techniques. However, many interpolation methods lead to blurred images where high frequency information (e.g., edges) is badly reconstructed. To overcome this issue, superresolution methods have been proposed (Manjón 2010a), that have the ability to effectively increase the image resolution while preserving sharp features of the underlying anatomy. In this study, we investigated the effect of applying superresolution as proposed in (Manjon 2010a) to the accuracy of cerebral cortex segmentation.

Methods

Ten T1--weighted high--resolution (HR) MRI scans were acquired from healthy volunteers on a Philips Achieva 3T system using a 3D fast gradient echo SENSE sequence with TR/TE=9.5/4.6 ms, 8° flip angle, FOV=240 mm, 480×480 matrix, and a slice thickness of 0.5 mm. Fig. 1 shows the processing steps. Images were denoised using the SANLM filter (Manjon 2010b), bias field corrected using N3 [START_REF] Sled | A nonparametric method for automatic correction of intensity nonuniformity in MRI data[END_REF]) and intensity normalized. Linear registration [START_REF] Collins | Automatic 3D intersubject registration of MR volumetric data in standardized Talairach space[END_REF] to MNI space was performed using a HR version of the MNI152 average non--linear template [START_REF] Fonov | Unbiased average age--appropriate atlases for pediatric studies[END_REF] followed by a brain extraction using BEaST [START_REF] Eskildsen | BEaST: Brain Extraction based on nonlocal Segmentation Technique[END_REF]. Low resolution images (LR) were created by downsampling the denoised HR images to 1 mm isotropic voxel sizes. HR images were "reconstructed" by two different methods: i) using B--spline interpolation (BS), and ii) using superresolution (SR) (Manjón 2010a, Fig. 2). LR, BS and SR images were normalized using the bias fields, transformations and brain masks calculated for the original images. Finally, Fast Accurate Cortex Extraction (FACE) [START_REF] Eskildsen | Active surface approach for extraction of the human cerebral cortex from MRI[END_REF] was applied to all images to obtain surfaces delineating the cortex (Fig. 3).

Results

The resulting SR images have sharper features than the BS images and the intensity differences are smaller (Fig. 2). Surfaces from the LR images clearly lack the subtle features found in the surfaces from the ORG images (Fig. 3), while surfaces from BS and SR images regain these features. Though differences between BS and SR surfaces are hard to observe directly, the error maps clearly indicate higher errors in the BS surfaces.

The WM surface errors decrease significantly with the construction of higher resolution images for both BS (p=0.003) and SR (p<0.001) (Fig. 4). The GM surfaces are also improved, however, only significantly with the SR method (p<0.001). Using the SR method improves the reconstruction accuracy significantly compared to the BS method for both WM (p=0.020) and GM surfaces (p<0.001). While increasing the resolution with either improves the Hausdorff distance over the LR image, no difference between BS and SR was detected.

Conclusion

Increasing the resolution of MRI images improves the accuracy of cortical surface based segmentations. Our results show that using the superresolution method provides sharper features and results in higher segmentation accuracy than using upsampling with B--spline interpolation. Superresolution methods can advantageously be used for reconstructing the cerebral cortex and can potentially reveal subtle patterns otherwise hidden where circumstances restrict the acquisition of native high resolution images.