WO2014032041A1 - Method and system for image registration - Google Patents

Abstract

A system, method, device, and non-transitory computer readable recording medium storing a program, for processing brain image data using a three-variable solver, the variables being a Point Correspondence, a Deformation Field, and a Resection Region, and resolving the three variables with a Nested Expectation and Maximization (NEM) framework.

Description

METHOD AND SYSTEM FOR IMAGE REGISTRATION CROSS REFERENCE TO RELATED APPLICATIONS

This applicaiion claims the benefit of U.S. Provisional Application No. 61/692,944, filed August 24, 2012, which is hereby incorporated by reference in its entirety. STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH

This invention was made with Government support under contracts CCF- l 139864, CCF-1 1365538, and CSI-1 136536 awarded by National Science Foundation. The Government has certain rights in the invention

FIELD OF THE INVENTION This application relates to the field of image registration and, more particularly, to a method and system for processing images using a nested expectation and maximization framework,

BACKGROUND Image registration is a process of determining positions in a first image to corresponding positions in a second or more other images. Differences between features of the first image and features of the second image, including missing or additional features, produce a correspondence problem. In one example of an application for registration, brain shift can severely compromise the fidelity of Image-Guided Neurosurgery (IGNS). Very few studies in the literature address the difficulties posed by brain deformation during and after tumor resection pose for image registration. Resection creates a cavity, which renders the biomechanical model defined on pre-operative MRI inaccurate, due to the existence of the additional parts of the mode! corresponding to the resection region.

I Most conventional approaches use a biomeclianical model to estimate the brain shift based on sparse intraoperative data after the dura is opened. Very few address brain deformation during and after tumor resection. The difficuity originates from the fact that resection creates a cavity, which renders the biomechanical model defined on pre-operative MRI inaccurate due to the existence of the additional part of the model corresponding to the resection region.

Some have investigated tissue retraction and resection using sparse available OR data and a finite element model, They used a two-step method: (1) remove tissue volume by manual deletion of model elements that coincide with the targeted zone and then (2) apply boundary conditions to the new surfaces created during the excision process. Determining the cavity is challenging because a portion of it will be filled by surrounding tissues.

Based on the bijective Demons algorithm, others have presented a registration framework to handle retraction and resection. They used a level set method to automatically detect resected regions. Also, they have presented an elastic FEM-based registration algorithm and evaluated it on the registration of 2D pre- with intra-operative images, where a superficial tumor has been resected. Some have pursued a semi-automatic method based on postbrain tumor resection and laser range data. Vessels are identified in both pre-operative MRI and laser range image; then a Robust Point Matching (or "RPM") method was used to force the corresponding vessels to exactly match each other under the constraint of the bending energy of the whole image. RPM uses Thin-Plate Splines (or "TPS") as mapping function. The basis function of TPS was a solution of the biharnionie, which does not have compact support and will therefore lead to, in real application, unrealistic deformation in the region far away from the matching points, in another words, RPM is not suitable for estimating deformation using sparse intraoperative data. SUMMARY

Disclosed are embodiments of a system, method, neurosurgery device, and non-transitory recording medium storing a program for processing an MRI resection as a three-variable solver, the variables being a Point Correspondence, a Deformation Field, and a Resection Region, and resolving the three variables with a Nested Expectation and Maximization (NEM) framework. These embodiments may employ a Nested Expectation and

Maximization Non-Rigid Registration (NEMNR ) approach discussed in greater detail below,

The present approach is an improvement, in that it is capable of: (1) removing the tetrahedra in a model corresponding to the resection region; and (2) building a heterogeneous biomechanical model, which is facilitated by a multi-tissue mesh generation method. This approach eliminates the manual removal step by treating the resection region as a variable, which is able to be automatically resolved by the NEM framework, an improvement over traditional EM optimization. In this framework, the present approach uses feature point rather than intensity as the metric. The Point Correspondence is viewed as an additional variable to the Deformation Field and the Resection Region. Unlike the pure geometric transformation characterized by local affinity and global smoothness, a heterogeneous biomechanical model is used to realistically simulate the underlying movement of the brain, which improves accuracy,

BRIEF DESCRIPTION OF DRAWINGS

To the accomplishment of the foregoing and related ends, certain illustrative

embodiments of the invention are described herein in connection with the following description and the annexed drawings. These embodiments are indicative, however, of but a few of the various ways in which the principles of the invention may be employed and the present invention is intended to include all such aspects and their equivalents. Other advantages, embodiments and novel features of the invention may become apparent from the following description of the invention when considered in conjunction with the drawings. The following description is given by way of example, but not intended to limit the invention solely to the specific embodiments described, which can be understood in conjunction with the materials that follow:

Figure 1 is a schematic of a computer based system that may be used in embodiments of the present approach.

Figure 2 is a flowchart of aspects of present method, giving context to the NEMNRR approach, Figure 3 illustrates coarse tissue mesh generation.

Figure 4 illustrates element outliers and point outliers for preoperative MRI and

5 interoperative MRI.

Figure 5 is a schematic of aspects of the nested expectation and maximization framework with Outer EM and Inner EM detailed,

10 Figure 6 is an illustration of aspects of Nested Expectation and Maximization aligned with row for inner EM and column for outer EM,

Figure 7 are results from synthetic data, with 7(a) surgery simulation, 7(b) illustrating central target points with surrounding source points. 7(c) showing central source points I S with surrounding outliers, 7(d) showing target points with noise, 76 showing estimated source points, 7(f) showing estimated target points, 7(g) showing non-resected mesh M and remaining mesh M/M_Rec ,7(h) illustrating a discrepancy between non-resection and true answer, and 7(i) illustrating discrepancy between resection and true answer.

20 Figure 8 are results form NEMNRR and a comparison between NEM RR and BMNRR for low-field MRI.

Figure 9 shows resected elements and mesh quality after deformation. 5 Figure 10 shows a deformation field with deformation magnitude and a portion of the brain, not including ventricles, removed to display the deformation field of ventricles,

Figure 1 1 shows point outlier rejection,. 0 Figure 12 illustrates registration results using NEMNR, with the three colums

corresponding with intraMRI, preMRL and warped preMRL respectively.

Figure 13 shows two rows of to different slices for BOLD deformation, w ith the last row superimposed BOLD on the corrected BOLD to show deformation. Figure 14 shows qualitative evaluation regarding canny edges. Figure 1 5 shows qualitative evaluation regarding resection margin.

Figure 16 is a multi-tissue mesh of brain and ventricle,

DETAILED DESCRIPTION 1. INTRODUCTION

It is noted that in this disclosure and particularly in the claims and/or paragraphs, terms such as "comprises," "comprised," "comprising," and the like can have the meaning attributed to it in U.S. patent law; that is, they can mean "includes," "included,"

"including," "including, but not limited to" and the like, and allow for elements not explicitly recited. Terms such as "consisting essentially of and "consists essentially of have the meaning ascribed to them in U.S. patent law; that is, they allow for elements not explicitly recited, but exclude elements that are found in the prior art or that affect a basic or novel characteristic of the invention. Embodiments of the present invention are disclosed or are apparent from and encompassed by, the following description.

A method and system for the non-rigid registration of images is provided having three beneficial aspects: (1) it does not require point correspondence to be known in advance: (2) it allows the input data to he incomplete; and (3) it simulates the underlying deformation with a heterogeneous biomechanical model. This approach formulates the registration as a three-variable functional minimization problem, with the three variables being Point Correspondence, Deformation Field, and Resection Region, Point Correspondence may be represented by a fuzzy assign matrix. Deformation Field may be represented by a piece- wise linear function regularized by the stress energy of a heterogeneous biomechanical model. Resection Region may be represented by a maximal connected tetrahedrai mesh. A Nested Expectation and Maximization framework has been developed to simultaneously resolve these three-variables. Embodiments of the present approach may include methods, devices, computer readable media, or systems. Method embodiments may include a method for the processing and optionally display of imaging data of the head of a patient, relying on preoperative image data of the head and an intraoperative image data of the head: Preoperative generally means prior to completion of an operation while intraoperative simply means after a first step of an operation. A device embodiment may include a neurosurgery device guided by imaging data of the head of a patient. A neurosurgery device may be considered a device employed for a brain Intervention or operation, such as an ablation, resection, biopsy, injection, implantation, radiosurgery, etc. in other words, a neurosurgery device may be guided by or otherwise employ the present methods or approach. A media embodiment may include non-transitory recording medium storing a program having computer code that configures a computer to implement a method for processing such pre- and intraoperative image data: A non-transitory recording medium is a physical recording medium capable of storing a computer code that can be used to specially configure a computer to implement the present methods or approach. A system embodiment may include a system for performing such a process, the system implementing the process with, among other things, at least one computer

Image registration may be considered as a process in which different sets of image data are brought into a single coordinate system. A mesh is a grid or partition, such as a 3D geometry partitioned into tetrahedrons using a Delaunay triangulation. Image meshing is a process of creating or using computer based models for a desired computational analysis based on image data. The data may be image data drawn from magnetic resonance imaging (MM), computed tomography (CT), or microtomography. For clarity, image data may refer to data from any of the above techniques or other similar approaches. The below description, however, was primarily with image data from MR! for simplification, A mesher is simply a computer based mesh generator or mesh generation tool.

A solver may be considered to include a computer programming tool for solving a computational objective. The present methods or approaches employ a specific form of solver suited to the objects described herein. A brain extraction tool is computer based process for separating or extracting brain image data from image data of a head or skill, which generally includes image data of non-brain tissue. Embodiments of the invention may be implemented by a speciaiiy configured programmable digital computer, or a system using one or more programmable digital computers and computer readable storage media, including those designed for processing CAD-based algorithms as known to those of ordinary skill in the art. In one embodiment, Figure 1 depicts ars example of one such computer system 100, which includes at ieast one processor 1 10, such as, e.g., an Intel or Advanced Micro Devices microprocessor, coupled to a

communications channel or bus 1 12, The computer system 100 further includes at least one input device 1 14 such as, e.g., a keyboard, mouse, touch pad or screen, or other selection or pointing device, at least one output device 1 16 such as, e.g., an electronic display device, at least one communications interface 1 18, at least one computer readable medium or data storage device 120 such as a magnetic disk or an optical disk and memory 122 such as Random-Access Memory (RAM), each coupled to the communications channel 1 12. The communications interface 1 18 may be coupled to a network 142.

One skilled in the art will recognize that many variations of the sy stem 100 are possible, e.g., the system 100 may include multiple channels or buses 1 12, various arrangements of storage devices 120 and memory 122. as different units or combined units, one or more non-transitory computer-readable storage medium (CRSM) readers 136, such as, e.g., a magnetic disk drive, magneto-optical drive, optical disk drive, or flash drive, multiple components of a given type, e.g., processors 1 10, input devices 1 14, communications interfaces 1 18, etc.

In one or more embodiments, computer system 100 communicates over the network 142 with at least one computer 144, which may comprise one or more host computers and/or server computers and/or one or more other computers, e.g. computer system 100, performing host and/or server functions including web server and/or application server functions. In one or more embodiments, a database 146 is accessed by the at least one computer 144. The at least one computer 144 may include components as described for computer system 100, and other components as is well known in the computer arts.

Network 142 may comprise one or more LANS, WANS, intranets, the Internet, and other networks known in the art. In one or more embodiments, computer system 100 is configured as a workstation that communicates with the at Ieast one computer 144 over the network 142, In one or more embodiments, computer system 100 is configured as a client in a client-server system in which the at least one other computer comprises one or more servers, Additional computer systems 100, any of which may be configured as a work station and/or client computer, may communicate with the at least one computer 144 and/or another computer system 100 over the network 142,

II. METHOD

As noted above, brain shift is an example in which a change in images severely

compromises the fidelity of IGNS. The flowchart presented in Fig, 2 describes the context of the present approach, which tbeuses on a framework or three-variable solver referred to in particular as Nested Expectation and Maximization Non-rigid Registration (NEMNRR).

In the brain shift example introduced above, a Brain Extraction Tool (BET) may be used to delete non-brain tissue image data from an image of the head, The BET essentially extracts the brain image data from the skull image, and the ventricle may then be segmented by manually delineating the ventricle boundary. The resulting two-tissue multi-label image (i.e., of the brain and ventricle) may then be fed into a multi-tissue mesher to produce a heterogeneous model in conjunction with specific biomechanical attributes. Edge detection may be performed on both pre- and intra-operative MRI images to produce a source point set and a target point set. Classic Canny edge detection may be facilitated by, for example, the open source Insight Segmentation and Registration Toolkit (IT ) available from the National Library of Medicine, to produce these two point sets.

Among other things, embodiments of the present approach beneficially address a specific feature point-based, non-rigid registration problem, which can be stated as problem 11 (Pi) Given a heterogeneous patient-specific brain model, a source point set in preoperative MRI and a target point set in intra-operative MRI, find Point Correspondence, Deformation Field and Resection Region. The three variables may be simultaneously resolved by a Nested Expectation and

Maximization Non-rigid Registration method with a biomechanical model and two point sets as inputs. The resolved deformation field can be used to warp the segmented preoperative MRI to improve the navigation accuracy, To resolve this problem, the deformation field may be represented by a displacement vector defined on the mesh nodes, The correspondence between two point sets may be represented by a correspondence matrix, and the resection region may be represented by a connected submesh. All three variables may thus be incorporated into one cost function, which is minimized by a Nested Expectation Maximization (NEM or Nested EM) strategy.

Unlike a traditional Point-based Registration (PBR) method, this Nested EM method does not require the correspondence to be known in advance, and allows the input data to be incomplete. To improve accuracy, a heterogeneous biomechanica! model is employed to realistically simulate the underlying movement of the brain. This heterogeneous model is built upon a multi-tissue mesh and specific biomechanica! attributes of each tissue,

A. 3D Multi-tissue Tetrahedral Mesh Generation

A multi-tissue mesher may include two steps: (1) starting from a homogeneous Body- Centered Cubic (BCC) mesh, identifying a coarse multi-tissue mesh by assigning each tetrahedron with a specific tissue label; and (2) deforming the coarse multi-tissue mesh surface to tissue boundaries defined in the multi-label image,

Coarse multi-tissue mesh: BCC mesh is a mesh based on a crystalline structure that is ubiquitous in nature. The nodes of BCC structure are grid points of two interlaced grids, The edges of BCC structure consist of edges of the grid and additional edges between a node and its eight nearest neighbors in the other grid. An advantage of the BCC is that it is highly structured and easily refined during simulation, even after red-green subdivision.

3) Label redistribution may be performed on the homogeneous BCC mesh to produce a coarse multi-tissue mesh, which is deformed subsequently. Given an initial label assignment (Fig. 3(a)), labels are redistributed to produce a surface robust against deformation (see dotted line in Fig. 3(b)). f the surface is not close enough to the tissue boundary (dashed line in Fig. 3(b)), a red-green subdivision may be performed on the tetrahedra across the tissue boundary, as shown in Fig. 3(c), The subdivision is likely to impair the robustness of the surface, and therefore label redistribution is performed again to produce a surface, which is robust and better approximates the tissue boundary (see Fig, 3(d)), The above procedure may be repeated until the multi-tissue surface is weii posed for deformation, and close enough to the tissue boundary. Thus, coarse multi-tissue generation may be as shown in Fig. 3, Li and 12 are tissue labels, the dashed line is the real boundary and the dotted Hne is the submesh interface: 3(a) initiation of labels; 3(b) redistribution of the labels; 3(c) subdivision if the fidelity criterion is not satisfied; and 3(d) redistribution of the labels again.

2) Deform the mesh surface to the tissue boundary; Mesh fidelity is achieved by iteratively deforming the mesh surface toward the tissue boundary. In each iteration, the deformation field, one and the only one variable, is resolved by minimizing the function,

W(U)

+ Xi WHfl - Di jj ²) (1) where 17 is the unknown displacement vector associated with the surface nodes; n is the number of the tissues; K; is the global stiffness matrix assembled by the tetrahedra within i-th tissue. Kj depends on two biomechanical attributes of f^'-th tissue: Young's modulus and Possion's ratio. /¾ is the global linear interpolation matrix assembled by mesh nodes. The assembling of H- is presented below. D, is the distance vector from the i-t surface to the i-Xh tissue boundary. ,, is used to balance the quality (first term) and the fidelity (second term),

W(U) may be minimized by solving:

f - o

UiiiT D_t (2)

B. Nested Expectation and Maximization Non-rigid Registration (NEMNRR)

In this section, a cost function may first be developed step by step from the classic point- based non-rigid registration energy function, and then a Nested Expectation and

Maximization framework may be used to resolve it,

1 ) Cost Function: Given, a Source point set S = (Sij -_! £ R^J and a Target point set T ~ {t₍} __l€ R^a , with known correspondence, i.e., St corresponding to the point-based non- rigid registration problem can be formulated as: =argmin( J_n A'(i ) i il+A∑_Sienl|5j + u(S ) - ti ll²), (3) u

where ihe first term is reguiarization, or smoothing energy, and the second term is similarity energy, is the deformation field and λ controls the trade-off between these two energies. O is the problem domain, namely for this example, the segmented brain.

Brain tissue removal influences both terms in the above equation (3); (i) the reguiarization in terms of the domain on which it is defined; and (ii) the similarity in terms of additional outliers introduced due to tumor resection. Extending equation (3) to equation (4) may be by specifying the reguiarization term with the strain energy of a linear elastic model.

( , c^, IV)=argmin( /_n.ni a( )^t€(u)d(n-n^,)+A∑_Sig_n._n. + u(s,) - ∑tjen_R tj tj \\ ^ (4)

where variable iY' represents the resection region and variable c,y reflects the degree to which point Si corresponds to ί,, The c,y may defined as in RPM with soft assignment. A Classic iterative Closest Point (TCP) method treats the correspondence as a binary variable and assigns the value based on the nearest-neighbor relationship; however, this simple and crude assignment is not valid for non-rigid registration, especially when large deformation and outliers are involved. One may define a range Ω_κ, a sphere centered at the source point with radius R, and only take into account: (1) the target points, which are located in of the source point; and (2) the source points, which have at least one target point in OR. Therefore, with this simple extension of RPM, this method is capable of eliminating outliers existing in both point sets. The homogeneous model employed in the reguiarization term in equation (4) is further extended to the following heterogeneous model:

(u, cT )~zr&MYi(X_a≠it.._{ai t} β¾(ι0*ί<(ι.).ί(-¾) + λ¾_6β__Λ, > + -∑-_{e 2fi}

(5)

where U β,_: ~ Ω - Ω', i ^~ 1 ... n if n = Ι,Ω' - 0 , and <¾· - 1, then equation (5) is reduced to equation (3). This means this approach may be viewed as a generalized point-based NRR method, characterized by: (1) ί 3 employment of a heterogeneous biomechanical model as the regularization; and (2) accommodation of incomplete data; and (3) without a requirement for correspondence.

Equation (5) may be approximated by equation (6) using finite element method:

W(U, C, M_Rec) =∑ U^TK_tU + _t(HY - D C))^TW(HU - D(C)) + λ₂∑ ¾· +

\M_Rec \ (6) where _n , o (ii.)^t _i(u)_£i(/2_i) is approximated by &¾£/. C is a point correspondence matrix with entries c¾-. The entries of the vector D are defined as: dt(cy)— s_t—

∑t_jea_R ify_> sj e M\M_Rec , where is the non-resected mesh that approximates Ω, and MRec is the resected mesh that approximates Ω'. The first term is the strain energy assembled on all elements in M\ M_Rcc , and the second term is similarity energy defined on all source points Si e M\M_Rec . The third and fourth terms prevent too many point and element outliers from being rejected.

W in the second term is the matching stiffness matrix of size 3 I S \ x 3 j 5 j, W is a block- diagonal matrix whose 3 x 3 sub-matrix Wt s defined as, ~~St^V9 where m is the number the vertices of the mesh. ~~ makes the matching term independent of the numbers of the vertices and the registration (source) points. S_k ^s is the average stiffness tensor for k-th registration point. S_k ^va makes the registration point behavior like an elastic node of the finite element model Assume the k-t registration point is located in the tetrahedron defined by vertices C₍, i€ [0: 3] . S ⁸ is calculated by S_k ^vg =∑L₀ , where K_Cl is a 3 x 3 sub-matrix of the global stiffness matrix K, hj is the interpolation factor, which is defined below.

H is the global linear interpolation matrix assembled by registration points. The A th registration point o_t contained in the tetrahedron defined by vertices c_£, £ [0: 3] contributes to four 3 3 sub-matrices: [H]_kc , [H]_kCl , [H]_kC , [H]_kCs . The diagonal matrix [H]_kc.— diag (hi , hi , The linear interpolation factor h_t, is calculated as:

Where v_Cj is the mesh node with index c,-. H is also used in mesh deformation, the second step of multi-tissue mesh generation. In mesh deformation with mesh nodes used as registration points, i.e., <¾ is the same with one of the four tetrahedron nodes, equation (7) is reduced to: f 1 for o_k = v_c.

(8)

( 0 for o_k≠ v_c. Compared to equation (6), equation (1) has only one variable U, and therefore it can be resolved directly by equation (2). The additional two variables in equation (6) result in the difficulty of resolving U.

Finding C and M is equivalent to outlier rejection. In Fig, 4, the identified triangles and rectangles denote point outliers. The original definition of outliers is extended to include elements in addition to points, and therefore the resection region M may be viewed as a collection of element outliers. As a result, the problem for resolving the three unknowns (i.e., Correspondence, Deformation Field and Resection Region) is transformed into the following problem: Given a heterogeneous patien -specific brain model, a source point set in pre-operative MR! and a target point set in intra-operative MRI, reject point outliers from both Source and Target point sets and reject element outliers from the biomechanieal model.

The Nested Expectation and Maxiinization Solver may be used to iterative!)' reject point and element outliers.

2) Nested Expectation and Maximization Solver: The Expectation and Maximization (EM) algorithm is a maximum-likelihood estimation of the model parameter (unknowns) in the presence of missing or hidden data. To estimate the mode! parameters, EM proceeds iteratively. Each iteration of the EM algorithm may be composed of two steps: the E-step and the M-step. in the E step, the missing data is estimated given the observed data and current estimate of the model parameters, in the M step, the likelihood function is maximized under the assumption that the missing data are known. The data of the estimate of the missing data from the E step are used in iieu of the actual missing data. Convergence is assured since the algorithm is guaranteed to increase the likelihood at each iteration.

In considering the problem of registration in the context of EM, the cost function of equation (6), from a probability (Bayesian) point of view, defines the likelihood function, in which the unknown (model parameters) is the deformation vector U. and the missing data are the correspondence C and the resection region M_Rec. Assuming M_REC is known, the more accurate the estimate of C, the more accurate the estimate of £7 is, and vice versa.

Therefore, an EM algorithm is naturally employed to solve U and under a specified ¾_c, To resolve M_ReCs treat {/ and C as one unknown pair < U, C >. The more accurate the estimate ofM _{ec is, the more accurate the estimate of < U, C > is, and therefore another high level EM can resolve the missing data Μ^_ύ,

The basic procedure of EM is: alternate between estimating the unknowns and estimating the missing data. In the Nested EM framework shown in Fig. 5, the inner EM may be used to resolve U and C with M_Rec fixed and the outer EM is used to resolve M_Rec, M_Re is approximated as a collection of tetrahedra located in a region corresponding to the tumor cavity in the intra-operative MR! in the model MR_0C is initialized as empty and updated in each iteration of the outer EM. If ail the tetrhedra contained in the cavity are collected, the outer EM stops, a) inner EM; Inner EM may be used to resolve f/and C given MR_SC. For each source point assume its correspondences are subject to Gaussian distribution, so cy can be estimated (E step) by equation (9),

k-\ ik Ryjln

Once C is estimated, C/may be resolved using equation (2). The resolved C/ is used to warp S closer to 7' and then the correspondence C is estimated again. This inner EM is illustrated in Fig. 6, in which the inner EM iterates along horizontal direction, and the outer EM iterates along vertical direction. In Fig. 6, the subscript ί is used to denote the inner EM and subscript o to denote outer EM. The superscript is used to denote the iteration number. For example, Ef denotes the h iteration of E step in inner EM. In the horizontal direction, inner EM iteratsvely estimates the correspondence and deformation field until no point outliers are detected. Inner EM begins from a search range (circles) with a larger radius R. For each source point, if there are no target points iocated in the circle centered at the source point, this source point will be rejected as outlier. For each target point, if it is outside of the search range, this target point will be rejected as outliers too. Once all the outliers are rejected, C may be estimated by equation (9), and {/can be solved with equation (2). in the next iteration, R is reduced by multiplying a simulated annealing factor 0.93, and the above procedures are repeated. Its pseudo code is presented in Algorithm 1, below, b) Outer EM: Outer EM, illustrated in the vertical direction in Fig. 6, is used to find Af¾_f and < U, C >, In M step, < U. C > is resolved by inner EM. in E step, M_RSC is resolved by an element outlier rejection algorithm. MR_SC is approximated by a collection of tetrahedron outliers, which fall in the cavity in the intra-operative MRL The cavity does not need to be identified in the intra-operative MRI, and it is in fact impossibie to distinguish the cavity from the background. The region Background Image (BGT) including the cavity and the background can be segmented by a simple threshold segmentation method; however, one cannot determine if a tetrahedron is an outlier based only on whether it is located in the BGL because some tetrahedra might happen to fall in the background rather than the cavity, To make the element outlier rejection algorithm robust, one may utilize the fact that the resection region is a collection of tetrahedra, which fall in the BGl of intra-operative MRI, connect with each other, and constitute a maximal connected submesh. This may effectively occlude false outliers. The collection of the outliers proceeds iteratively, and in each iteration, more specifically the E step of outer EM, additional outliers will be added into _¾c if they fail in the BGI md connect with the maximal connected submesh identified in previous iteration.

The element outlier rejeciion algorithm is presented in Aig. 2. Assuming the current deformation field and resection region is ?/and Af_¾c respectively, the new estimated outliers can be obtained by transforming the remaining mesh M \ M_REC using t and then finding all elements, which satisfy the requirements: (1) are completely inside the background of the intra-operative image, and (2) are connected with A/j_¾c. The outer EM iterativelv rejects element outliers using Alg. 2 and computes < U, C > using Alg, 1 until no additional element outliers are detected. Alg. 3 presents the whole pseudo code of the Nested EM algorithm,

Algorithm 1 (Alg, I) Feature Point Outlier Rejection

Input: M non-resected mesh, _¾c: resected mesh, S source points, Γ: target points, e: tolerance, r, annealing rate, R: search range

Output: U displacement vector, C; correspondence matrix

1 : i/ <= /

2: repeat

3: Transform S based on U: S <s= U(S)

4: E step:

// outlier rejection for S

5: S = S si \ if no target points in Ω_κίοτ sj

// outlier rejection for T

6: S = T{ti ] i no source points list it within„¾ }

7: Estimate correspondence C according to equation (9)

8: M step

9: Solve U according to equation (6)

10: change jj i/_j— l/_j-i l! between successive iterations

11 : Decrease R: R <t= R r

12: until change < e

Algorithm 2 (Alg. 2) Element Outlier Rejection Input: M: non-resected mesh, resected mesh, U: displacement vector, BGI; background image

Output: Mn_eC n w resected mesh, S new source points

1 : Obtain deformed remaining mesh M_Dfir ll&^ i M )

2: Find all elements Mj completely contained in the background image BGI and constitute the largest connected mesh with M_Rec

3; Map M_} in ,_¾/to M> in MMn_ec

4: S S{Si \ Si€ M₂}

5: M_Rec ^ M_Rec U M₂

6: Scale Young's modulus for the elements across the boundary

Algorithm 3 (Alg, 3) Feature Point Based NRR (NEMNRR)

Input: MRI; pre-operative MRI, iMRI: intra-operative MRI

Output: U: displacement vector

1 : Segment brain on MRI and mesh generation for M

2: Segment background image BGI on iMRI

3; Canny edge detection on MRI to get S

4: Canny edge detection on iMRI to get T

5: input R, e, and r

6: Initiate M_Rec 0

7: repeat

8: M step: U, C *= PtOutlierRejection(M, M_Rec, S, T, e, r, R)

9; E step: M_Rec, S «= EleO tUerRejection(M, M_Rec, U, BGI)

10: until MR_SC does not change

I II . RESULTS

Experiments have been conducted on synthetic data and clinical data, including low- and high-field MRI. The experiments on low-field MRI represent a typical application, namely using sparse intra-operative information to correct preoperative MRI. The experiments on

37 high-field MRI represent another typical application, namely fusing pre-operatively acquired Blood-Oxygen-Level Dependent scan (or "BOLD"), not available in OR, with intra-operative MRI. in addition to these two applications, experiments were also conducted to compare the present approach and system with a classic point-based NRR and to compare the

homogeneous model with the heterogeneous model,

A. Experiments on Synthetic Data

To generate a synthetic resected brain, a surgery simulation tool was developed to simulate brain resection, as shown in Fig, 1 1(a). The synthetic deformed resected brain was produced by a surface-based registration tool that was capable of deforming the brain based on specific boundary condition: the deformation of the resection surface. The source points S were simulated as the surface nodes of the resection region before deformation, and the target points T' are the surface nodes of the resection region after deformation (Fig. 1 1(b)). All the surface nodes except the green nodes were added into S as the outliers (smaller white points in Fig. 7(c)). The outliers for T, were generated using white Gaussian noise (grey points in Fig. 7(d)).

Fig. 7(e) and Fig. 7(f) show that Alg. 1 correctly detects all the source points and target points. Most outliers are rejected from S and T except three outliers in S (white points in Fig, 7(e)) and corresponding three outliers in 7 (brighter white points in Fig, 7(e)). Fig. 7(g) shows the element outlier removed mesh M \ Mp_wc produced by Aig. 2. The non-resected mesh aud the resected mesh M\ are shown together to illustrate the resection region clearly.

Two experiments were performed to verify that the removal of element outliers from the mode! can improve the accuracy of the registration. Both experiments registered the non- resected brain with the synthetic resected brain, but one experiment rejected element outliers in the model while the other did not. Both experiments used the same source points and target points as shown in Fig. 11(b), therefore the variation of the results was caused by whether outliers are rejected or not. In each experiment, the registration result was compared with synthetic deformed resected brain (true answer) by subtracting one from another to produce a discrepancy image, If a registration result was closer to the true answer, the discrepancy would be smoother. Comparing Fig, 7(h) to Fig. 7(i), the method, involved element outlier rejection, demonstrated a more accurate result. This experiment validated the hypothesis that removal of the resection region from a bsoniechamca! model can improve the accuracy of the registration.

B. Experiments on Clinical MRT

Experiments were conducted on iMRl, including Sow- and high-field iMRL The data include ten public cases from a hospital in the U.S. and four cases from the hospital outside the U.S. Table I lists the patient information including the gender, tumor location, and histopathology of these fourteen cases,

TABLE ί

(Patient information of fourteen cases.)

The first ten cases are drawn from a surgical planning laboratory of the university hospital located within the U.S. while the final four cases are from a hospital withi China. The MRI of the ten U.S. public cases were acquired with a protocol: whole brain sagittal 3D-SPGR (slice thickness 1.3 ram, TE/TR-6/35 ms,

x 256). For the four Chinese eases, the pre-operative MRI and high-field in Ira-operative MRI were acquired in 8 minutes with the same protocol: pre-operative (high-field infra-operative) MRI, 3D ΊΊ - weighied magnetization-prepared rapid gradient echo (MPRAGE) sagittal images with [dimension=256 x 256 x 176, inplane resolution^ .0 x 1.0mm, thickness=1.0 mm, FOV^::::256 x 256].

Low-field MRi, Tl -weighted, three-dimensional, fast spoiled, gradient recalled sequence with [dimension^al 28 x 128 x 35, inplane resolution^ .25 x 0.94, thickness= mm, FOV^™ .160 x 120] was acquired in 7 minutes.

BOLD was acquired in about 20 minutes with the protocol: a gradient-echo EPi sequence [dimension=64 x 64 x 20, inplane reso!ution=4 x 4mm TR=2000ms, slice thickness^ mm, FOV^::::256 256, gap= 1 mm] .

1) Low-Field Intra-operative MRI: Fig. 8 illustrates the results of NEMNRR and a comparison between NEMNRR and BMNRR for low field MRI. The numbered arrows point to the boundaries of the pre-MRI of concern. Fig. 8(a) shows a preoperative MRI (or "preMRI"), and Fig. 8(b) shows an intra-operative MRI. Unlike high-field MR imaging, the low-field MR is incapable of capturing the whole brain due to its limited field of view or FOV (160 x 120). In fact, to save imaging time, only the oval region within the FOV was imaged. As shown in. Fig. 8(b), only the brain near the resection region was displayed. The extracted edge points on pre- (Fig. 8(a)) and intra-operative MRI (Fig. 8(c)) served as the source points and target points, respectively. Both source points and target points were obtained using an implementation of canny edge detection algorithm. In Fig. 8(e), the edges of the iMRI were superimposed on the preMRI to show the discrepancy before registration. The arrows point to the boundaries of the preMRI of concern. After NEMNRR (see Fig. 8(d)), it was observed: (i) the boundary- 1 of the preMRI was not deviated by its surrounding outlier, i.e. the oval boundary, (ii) the boimdary-2 in the vicinity of the tumor in preMRI was deformed to the boundary in the vicinity of the resection region in iMRI, (iii) the boundary-3 still agreed well with the boundary of the iMRI, as it did before registration, even though many outliers surround it. Fig. 8(e) and Fig. 8(f) sho the comparison for another slice (i.e., slice number 92), After registration, both the boundary of the ventricle and the boundary in the vicinity of the tumor in the preMRI matched well with die corresponding boundaries in the iMRI.

Fig. 8(g), Fig. 8(h) and Fig. 8(i) show the results using a block matching or BMNRR method. This method used Block Matching to find the correspondence, and then drove a homogeneous biomechanical model to estimate the deformation. Comparing Fig. 8(g) and Fig. 8(i) to Fig. 8(d) and Fig. 8(f), respectively (i.e., same slice as the number denotes), it can be seen that the larger deformation in the vicinity of the tumor still existed after BMNRR. Two factors account for this: (1) Block Matching cannot find correct correspondence near the 5 resection region due to intensity variations and noise around the resection region; and (2) the resection region was not removed from the biomechanical mode!. BMNRR showed as good results in the deep part of the brain as the present method (e.g., see the ventricle in Fig. 8(h) and Fig. 8(1)), due to fewer outliers and rich texture information, which were helpful for Block Matching.

10

2) High-Field Intra-operative MR J: Fig, 9 shows the result of Alg. 2 for element outlier rejection, with the resected elements and mesh quality after deformation. The quality of the remaining mesh after deformation (i.e., maximal deformation magnitud 18.2mm) is still acceptable, as shown in the bar chart of Fig. 9, in which the minimal dihedral angle was I S 0,212. The reason for this was that the multi-tissue mesh was very robust against larger deformation.

Mesh quality is a measure of how well the elements of a mesh are shaped. Mesh quality can be evaluated using different metrics, such as the minimal dihedral angle, aspect ratio, etc. 0 Mesh quality influences the accuracy and speed of the finite element solver. For example, if the angle between two triangles of a tetrahedron is very small, the assembled stiffness matrix will be ill-posed, characterized by a larger condition number (a measure of the asymptotically worst case of how much the error can be magnified in proportion to small error). This characteristic of "magnifying the error" of an ill-poised matrix severely deteriorates the 5 solution of a linear system of equations because the linear system is usually solved

iteratively.

Fig. 10 shows the deformation field of the heterogeneous model. The grey scale denotes deformation magnitude, A portion of the image of the brain has been cut off to expose the 0 ventricle and its deformation field. The largest deformation reaches 18,2mm, still in the effective range of the linear elastic biomechanical model. The larger deformation occurs in the region near the resection, and the ventricle on the tumor side is squeezed inward as the arrows show. Fig 11 shows the results of Aig, 1 for point outlier rejection. Fig. Ufa) corresponds with case 2, and Fig. 1 1 (b) corresponds to case 6. Comparing to the edges before outlier rejection, most outliers are removed from the preMRl and iMRI after outlier rejection. Table Π lists the numbers of the source/targe points before/after registration.

Fig, 12 shows the registration results of the NEMNRR method for four cases. The three columns correspond with intraMRI, preMRl, and warped preMRl respectively. The four rows correspond with case 2, 7, 12, and 13. Edges detected on iMRI are superimposed onto preMRJ and warped preMRl, respectively, to illustrate the improvement of the boundary matching after registration.

A typical application of the registration method on high-field MRI is to merge pre- operativeiy acquired functional MR! (fMRD such as BOLD with infra-operative MRL BOLD cannot be acquired infra-operatively, which makes the non-rigid registration method the only feasible way to bring BOLD into Operating Room.

BOLD is a type of specialized MRl scan that measures the hemodynamic response (change in blood flow) related to neural activity in the brain or spinal cord of humans or other animals. It is one of the most recently developed forms of neuroimaging. Since the early 1990s, BOLD has come to dominate the brain mapping field due to its relatively low invasiveness, absence of radiation exposure, and relatively wide availability.

BOLD naturally agrees with a reference preoperative MRI, which is used to register with high-field intra-operative MRI to recover the deformation between them. The recovered deformation will be applied on BOLD to produce a corrected BOLD. In Fig. 13 is a merger of the corrected BOLD and iMRI into one image, winch can be used for function-guided neurosurgery navigation.

It Is believed that there is no a feasible way to directly evaluate the accuracy of the corrected BOLD thus far. Resting state fMRi provides a means to intra-operative ly acquire fMRL which will probably become a gold standard in future, but as of now remains under investigation; however, the accuracy of the corrected MRI can be indirectly evaluated by measuring the accuracy of corrected anatomical MRJ based on the fact that function moves with anatomy in the same way, which has been widely accepted in the neurosurgery community.

3) Quantitative evaluation of NEMNRR: A quantitative evaluation of the NEMNRR approach follows, regarding the edges and tumor/resection margin for cases listed in Table Ϊ. The edges come from canny edge detection, which is completely automatic, and therefore the evaluation on canny edges is not subject to the influence of raters. Compared to the evaluation on canny edges, the evaluation on the tumor/resection margin is more specific but involves a manual delineation of the margin.

Hausdorff distance (HD) is employed as the measure of the registration accuracy. This metric has been used for the evaluation of eleven public clinical cases. HD measures the degree of mismatch between two point sets/regions with equation: H(A, B^') ~ ax(h(A, ), k(B, AY) ( 10) where h(A, B) is directed HD defined bv h(A, B)— ^^'ίΐ"" ,, ^ ~ b\\. HD is a very stringent measure because it calculates the maximal mismatch of two point sets in two directions. Effectively using HD for the evaluation requires knowledge of its advantages and disadvantages. Compared to landmark-based evaluation, HD is automatic, repeatable, and not subject to the influence of the raters; however, HD is susceptible to outliers, in the experiment, the outlier rejected source point and target point sets were used to measure the HD. HD relies on the point set pair. For two different pairs, HD could have different values. A comparison of NEMNRR with BMNRR is needed for the next section. One method is to detect edge pair < EIMRI NEMNR E nMRijiE MRR > on iMRI and NEMNRR warped preMRi, and then calculate HDxEmN R* The next step is to detect another edge pair < Ei_mm_BmNFsR

EpreM!u BMNRR> on iMRI and BMNRR warped preMRI, and then calculate HDBMNRR,' however, the comparison between HDNBMNRR and HDBMNRR is meaningless because it is very difficult to guarantee that EIM I NKMNRR is same with EIMRI. s mR," and there are no outliers within each pair. To overcome this challenge, in the experiment, the outlier rejected target point set will be warped using NEMNRR and BMNRR, respectively, and the warped point sets may be taken as Epr mijiEMNRR and EpreMRi SMN Ri respectively. This method applies to the comparison between the homogeneous model and heterogeneous model, as well. To precisely perform the evaluation, surgeons examined each case to determine whether the deformation was induced by brain shift or by resection. The fourteen cases were classified into six brain shift eases and eight resection cases as shown in Table II, in which character 'S' denotes brain shift and character 'K' denotes resection. For the evaluation regarding edge points for fourteen cases, HD was used as the measure. For the evaluation regarding tumor margin for six brain shift cases, both tumor boundaries in (warped) preMRJ and iMRl were delineated, and I-ID was used as the measure. For the evaluation regarding resection margin for eight resection cases, due to the existence of the cut of the resection margin, directed HD (from the resection margin to its superset i.e. the tumor boundary) was used.

Before presenting the quantitative results, the parameter setting for NEM R, including mesh generation, canny edge detection, and EM solver may be considered, Also presented are the parameters for BMNRR., which will be compared to NEMNRR. The BMNRR parameters about the ten public cases are not available. Because the warped pre-operative MRI were already available, those experiments were not repeated.

Table II lists all parameters for fourteen cases. As noted above, the first ten cases were form a U.S. hospital while the next four were from a Chinese hospital,

For each case, both single- and multi-tissue meshes were generated using the algorithm described in ΪΙ-Α. While this example related to a simple two-tissue heterogeneous model (i.e., ventricle plus brain), it may be applied to others, Three parameters are related to this raesher: BCC, Subdivision, and Resolution, BCC specified the size (i.e., physical unit) of the bounding box of the BCC tetrahedron, Subdivision controlled the number of the mesh subdivision for each tissue, and Resolution defined the up-sampling level, For the single mesh, these three parameters were fixed to 8,0.7,0, For the multi-tissue mesh, case 1 required that the original image be up- sampled once (resolution=l) since the ventricle of this patient was very small. For cases 5, 8, 10, 12, 13, and 14, the Subdivision parameter was adjusted to less than 0.7 to avoid more subdivision along the ventricle boundary, Compared to single-tissu mesh, the multi-tissue mesh had more vertices and tetrahedra for precisely describing the geometry of the ventricle.

Canny edge detector was employed to perform the edge detection on both preMRl and inirMRl to obtain source points and target points, Four parameters were related to this detector: Gaussian variance, maximum error, and low and high hysteresis thresholdings. The first two parameters influence the accuracy of the edge location. A non-zero variance will blur the image, thus reducing the outliers at the expense of impairing the accuracy of the edge location. Because the present method is able to reject outliers, the variance was set to zero to avoid impairing the accuracy of the edge location, A zero variance invalidates the parameter: maximum error, so this parameter is not listed in the table. The high and low hysteresis thresholdings divide the potential edge points into three parts: strong edge points (i.e., any points with a gradient magnitude larger than the high thresholding), weak edge points (i.e., any points with a gradient mangitude smaller than the high thresholding but large than the low thresholding), and non-edge points (i.e., any points with a gradient magnitude less than the low thresholding). In the present approach, the strong edge points were only considered by setting the low thresholding equal to the high thresholding due to the following considerations: (1 ) the biomechanical model can work with sparse data; (2) it can avoid too much outliers involved; and (3) users do not need to adjust the low thresholding, enabling ease-of-use. Consequently, the four parameters are reduced to one, i.e., the high hysteresis thresholding, The resulting number of the edge points is listed in Table II, The number listed in the parenthesis is the number of edge points after outlier rejection Alg. 1.

There are four parameters related to NEMNRR. λ; is fixed to 1.0 for all cases because the tensor W makes the registration independent of the number of the vertices and the edge points, e is set to 0,5. The annealing rate is 0,93, Only the search range is listed in the table. A value larger than the maximum displacement magnitude need not be specified because NEMNRR is able to recover the deformation larger than the initial search range. Compared to the cases from the Chinese hospital, the ten public U.S. cases do not involve significant deformation, so the search range was set to 5.0 mm for both resection cases and brain shift cases. The search range was set to 10.0 mm for Chinese cases. Both the homogeneous and heterogeneous models used the same search range, but the heterogeneous model ran faster than the homogeneous model. This was because the ventricle of the heterogeneous model is softer than the homogeneous model, which rendered the heterogeneous model to deform flexibly, and then makes the matching term and regularization term, defined in equation (6), to reach balance quickly.

BMNRR includes three steps: point selection, Block Matching, and solver. The rejection fraction of point selection was set to 0.95, and the half block size of Block Matching was set to 3 x 3 x 3. The parameter listed in the table was the half search window, which functions the same as the search range of NEMNRR, The search window was set equal to the search range for the purpose of fair comparison between NEMNRR and BMNRR. BMNRR was computationally intensive, especially for the Block Matching (e.g., around 10 to 15 minutes), BMNRR was parallelized on a cooperative architecture in previous work. The time listed in the table was the running time of the parallel BMNRR.

Table HI shows the results of quantitative evaluation of NEMNRR on both edges and resection margins (The comparison results are presented in this table as well), To evaluate the accuracy regarding the edges, outlier rejected edge points were used in preMRl and intraMRI to calculate HD before registration in order to alleviate the influence of the outliers on HD, Outlier rejected edge points were used in intraMRI and warped outlier rejected edge points in preMRJ to calculate the HD after NEMNRR registration. To evaluate the accuracy regarding the tumor/resection margin, all fourteen cases were first classified into six. brain shift cases and eight resection cases (see Table I). For the six brain shift cases, the tumor boundaries in preMRl and iMRi were delineated to calculate the HD. One case (case 3) was not included because its tumor boundary was too difficult to be identified in both preMRl and iMRi simultaneously. In each of eight resection cases, the preMRl slice, in which the tumor is completely intra-operatively resected, was chosen so the margin corresponding to the resection margin of iMRi could be identified by using the tumor boundary. Otherwise, there would have been no way to identify the corresponding margin in preMRL The resection margin was delineated in iMRL and directed HD was used for evaluation.

The average error regarding edges was 2.01 mm, and 3,50 mm regarding tumor/resection region. Compared to rigid registration, NEMNRR increased the accuracy by 9,51 mm regarding edges, and 4,01 mm regarding tumor/resection margin on average. Two-tailed t test (significance level 0.5) demonstrated the improvement was statistically significant for both edge points (P-va!ues 6.63E-8) and tumor/resection margin (P-value 0.003).

4) Comparison: This section is a qualitative and quantitative comparison of NEMNRR with BMNRR, and the heterogeneous model with the homogeneous model. Both qualitative and quantitative comparisons were performed on edges and tumor/resection margin. HD was employed to quantitatively measure the degree of the mismatch. The same multi-tissue mesh was employed to build the homogeneous and heterogeneous model, and thereby eliminated the influence of the geometry in comparing these two models. The BMNRR registration results for ten public cases were downloaded from. a) Qualitative comparison: Fig, 14 shows the comparison regarding edge points for four cases. The canny edges detected in intraMRI are superimposed on BMNRR warped preMRi,

NEMNRR warped preMRi, homogeneous model warped preMRi, and heterogeneous modei warped preMRi, respectively in order to illustrate the dismatch. Fig. 15 shows the comparisons regarding tumor/resection margin for the same four cases as Fig. 14. The tumor boundary for brain shift cases and the resection margin for resection cases are identified in intraMRI, and then superimposed on the three warped images to show the dismatch regarding tumor/resection region.

For the evaluation on both edges and resection margin, NEMNR showed better boundary matching than BMNRR. Compared to the homogeneous model, the heterogeneous model did not demonstrate perceptible improvement b) Quantitative comparison: Table 111 shows the results of quantitative comparison between NEMNRR and BMNRR, and between the homogeneous model and heterogeneous modei. To evaluate the accuracy regarding the edges, outlier rejected edge points were used in iMRI and the warped edge points by BMNRR to calculate HD for BMNRR. Outlier rejected points were used in iMRI, and warped points by NEMNRR to calculate the HD for NEMNRR. Due to the unavailability of the deformation field for the ten public cases, the comparison regarding edges was performed on four Chinese cases, To evaluate the accuracy regarding tumor/resection margin, for six brain shift cases, the tumor boundaries in both iMRI and BMNRR warped preMRi were delineated to calculate the HD for BMNRR. The tumor boundaries in both iMRI and NEMNRR warped preMRi were delineated to calculate the HD for NEMNRR, For eight resection cases, the resection margin in iMRI was delineated, and the directed IFD was used for comparison,

As shown in Table III, NEMNRR outperformed BMNRR by increasing the accuracy by 9.92 mm on average regarding edge points, and 2.40 mm on average regarding tumor/resection margin, A two-tailed t test (significance level 0.5) demonstrated that the improvement is statistically significant for both edge points (P-values 0,048) and tumor/resection region (P-vaiue 0.04), The statistical values iisted in the table did not distinguish brain shift cases from resection cases, Performing two-tailed i test on brain shift cases and resection cases yields a P-value 0.491 and 0.033, respectively, in brain shift cases, the NEM RR did not show significant improvement compared to BMNRR, The reason should be that both the intensity variation and the deformation were not large, thus the BMNRR worked well

Both NEMNRR and BMNRR did not perform well in case ten. In examining this case, it was found that few edges were detected, which cannot provide enough information to drive the model. Moreover, the evaluation results on edges and resection margin lacked consistency. ^'The NEMNRR(Ho) reduced the error of rigid registration from 13.30 mm to 1.88 mm regarding the edge points, but the evaluation on resection margin showed the error was reduced only from 13.01 mm to 1 1.09 mm. The reason for this was most likely that the detected edges, although well-aligned, were too far away from the tumor and the resection region, and thus not effective in driving the model to estimate the deformation around the resection margin.

TABLE II

(Parameters for mesh generation. Canny edge detection, NEMNR, and

BMNRR)

For Table Π, the numbers Iisted in the column of Canny is the number of edge points detected by the Canny detector before outlier rejection. The number in parentheses is the number after outlier rejection. Not all parameters of NEMNRR and BMNRR are listed. Only the comparable parameters, i.e., the search range of the NEMNRR and the search window of BMNRR are listed. For BMNRR, the other parameters are block size = 3 3 X 3, block rejection fraction =*·^· .95, block connexity type = 26, A = 1.0, For NEMNRR, c = 0,5, r - 0.93, A = 1.0.

TABLE ill

(Quantitative evaluation and comparison for the fourteen cases,)

"Ho" denotes the homogeneous model, and "Hete" denotes the heterogeneous model.

"Canny" denotes the evaluation was performed on the edge points detected by Canny detector, and "Tumor" denotes the evaluation was performed on the tumor or resection margin and, depending on the case, is brain shift or resection.

To specifically measure the influence of the model to the registration, the multi-tissue mesh, as shown in Fig. 36, was employed in both models. As a result, the influence of the discrepancy of the geometry and topology between single mesh and multi-tissue mesh was eliminated. The only differences between the two models were the biomechanieal attributes of the ventricle. The homogeneous model used Young's modulus E ¹⁼¹ 3000i¾, Poisson's ratio ^:~ 0.45 for all tetrahedra and the heterogeneous model replaced Young's modulus with E - IQP and Poisson's ratio with v = 0.1 for the ventricle. These models were compared for edge points and resection margin. HD was used to measure the dtsmatc of the boundary for fourteen cases regarding edge points and thirteen cases regarding tumor/resection margin. The evaluation results on edge poiiits show the improvement brought by the heterogeneous model is statistically significant (P- value 0.04), but for the resection margin, the improvement is not (P-value 0,777). Because the difference between these two models is so small (0.21 mm on average regarding edge points and 0,05 mm regarding tumor margin), the statistical result for resection margin is more susceptible to the error involved in the manual delineation of the margin. The edge points detected by canny edge detection were not subject to the influence of the raters, so the evaluation result on the edge points was more accurate. This work investigated only a very simple heterogeneous model (brain plus ventricle), which was capable of improving the registration accuracy, but by a small magnitude. IV. DISCUSSION

Tumor resection led to the variation of the intensity and large deformation in the vicinity of the resection region. It not only resulted in pixel/voxel outliers, but also element outliers for finite element-based registration. BMNRR relied on the Block Matching to find the correspondence, Block matching is susceptible to the outliers and does not perform well in the vicinity of the resection region, with which surgeons are concerned most.

Feature point-based method is robust against outliers. Edge points were used specifically due to the following consideration: edges provide valuable information for registration, There is no way to register two homogeneous regions. However, once the intensify changes, (i.e. an appearance of an edge, strong or weak, depending on the magni ude of the gradient), the two regions could be aligned. Intensity- ased methods utilize the edge information implicitly, but edge-based methods utilize it explicitly. The edges are the target, used for comparison, in comparing two images, the matching of the edge is generally of greater interest. Edge-based methods provide a mechanism, as in equation (6), to directly control the degree of the matching. If all the valuable edges (such as the boundaries of the tumor, ventricle, and resection region) can be detected, and all outliers can be removed, the matching degree of the margins can be controlled by the tolerance, if the approximation error induced by the regularization term in equation (6) is not considered. The current edge detector cannot avoid outliers and detect ail valuable edges, which motivates the design of a registration method characterized by utilizing a biomechanical model to accommodate sparse data and reject outliers. In fact, the canny edge detector used in the above examples was not powerful because the functionality of edge-linking is missing. Combined with a powerful edge detector, the NEMNRR could further improve the accuracy of the registration. The use of edges (observation of the data) and a biomechanics! model (priori knowledge) are a trade-off in the registration, if enough edges can be detected, the model is not required any more. Conversely, if the model precisely described the deformation of the soft tissue, the required edges could be reduced to the boundary condition, i.e. brain boundary. Currently, there are no biomechanics 1 models that can precisely describe the deformation of the brain. Thus, in performing edge detection, there is an attempt to detect all brain boundaries, and therefore guarantee that the registration can, at least, reach the accuracy defined by the model. Based on these 14 cases, most brain boundaries, ventricle boundary, and part of the tumor/resection margin can be detected. The additional ventriele/turaor/resection boundaries may be used to reduce the error between the deformation determined by the model and the real deformation. Furthermore, to relax the requirement for the edges, the fidelity of the mode! may be improved by a heterogeneous model, In the above work, only a very simple heterogeneous model was considered, but the presented framework was suitable for any complicated heterogeneous model if corresponding biomechanical properties are available.

For a linear elastic model, it is possible to adjust the parameter .¾ to achieve the same effect of softening or hardening the tissue as adjusting Young's modules; however, this method is a giobai control, making the whole image of the brain softer or harder. A heterogeneous model is able to distinguish different tissues, owing to a multi-tissue mesher, and locally adjust the properties of the brain by specifying appropriate parameter values,

V. CONCLUSION A novel non-rigid registration method to compensate for brain deformation induced by tumor resection has been described. This method does not require the point correspondence to be known in advance and allows the input data to be incomplete, thus producing a more general point-based NRR. This method uses strain energy of the biomechanical model to regularize the solution. To improve the fidelity of the simulation of the underlying deformation field, a heterogeneous model based on a multi-tissue mesher was constructed. To resolve the deformation field with unknown correspondence and resection region, a Nested EM framework was developed, which can effectively resolve these three variables simultaneously. Experiments targeted two typical applications, i.e. the use of low-field MRI to correct brain deformation and the used of high-field MM to bring pre-operative BOLD into OR. The experimental results showed that the registration method, in the vicinity of the tumor resection, was more accurate than the state-of-the-art NRR.

The heterogeneous model was able to improve the registration accuracy, but on an imperceptible scale, it is believed that as more tissues are incorporated into this model, the improvement of the accuracy brought by the heterogeneous model will become noticeable. SAMPLE EMBODIMENTS

Embodiments of the present approach may include a method for the processing, and optionally display, of imaging data of the head of a patient, the imaging data produced by a magnetic resonance imaging scanner to collect pre-operative image data of the head and an intra-operative image data of the head stored in at least one non-transitory data storage device; the processing method may include the steps of (i) accessing the at least one data storage device and extracting pre-operative brain linage data from the pre-operative image data using a brain extraction tool and delineating a ventricle boundary to produce a segmented mask and a segmented pre-operative brain image: (ii) accessing the at least one data storage device and extracting intra-operative brain image data from the intra-operative image data using the brain extraction tool; (iii) producing a heterogeneous biomechanical model having a plurality of mesh nodes using a multi-tissue raesher; (iv) producing a set of source points by edge detection of the segmented pre-operative brain image; (v) producing a set of target points by edge detection of the intra-operative brain image data; (vi) wherein a deformation field may be represented by a deformation vector defined on the plurality of mesh nodes, point correspondence of the set of source points and set of target points may be represented by a correspondence matrix, and a resection region may be represented by a connected submesh; and (v) processing the segmented pre-operative brain image using a three-variable solver, the variables comprising the deformation field, the point

correspondence; and the resection region.

The deformation field may be represented by a deformation vector defined on the plurality of mesh nodes by a piece-wise linear function regularized by a strain energy of a heterogenous biomechanical model; with the point correspondence of the set of source points and set of target points is represented by a correspondence matrix thai is a fuzzy assign matrix, and the resection region is represented by a maximal connected tetrahedrai submesh. The processing of the segmented preoperative brain image may be by using a three-variable solver method comprises resolving the three variables with a Nested

Expectation and Maximization (NEM) framework. This could include iterative ly rejecting element outliers of the resection region and point outliers using the NEM framework; and approaching the remaining non-resection brain domain and resolution correspondence.

In some cases, the mesh may be a multi-tissue tetrahedrai mesh generated by starting from a homogeneous body-eentered-cubic mesh, identifying a coarse multi-tissue mesh by assigning each tetrahedron with a specific tissue label; and deforming the coarse multi-tissue mesh surface to tissue boundaries defined the multi-label image.

Another embodiment may be a neurosurgery device guided by imaging data of the head of a patient, with the neurosurgery device including (i) a magnetic resonant imaging scanner for producing imaging data comprising a pre-operative image of the head and an inira-operative image of the head; (ii) at least one non-transitory data storage device in communication with the magnetic resonant imaging scanner for receiving and storing pre-operative image data and an intra-operative image data; (iii) a computer system having a specially configured computer processor in communication with the non-transitory data storage device, the computer system configured to process the pre-operative image data and an intra-operative image data using a method comprising: (a) accessing the at least one data storage device and extracting pre-operative brain image data using a brain extraction tool and delineating a ventricle boundary to produce a segmented mask and a segmented pre-operative brain image; (b) accessing the at least one data storage device and extracting intra-operative brain image data using the brain extraction tool; (c) producing a heterogeneous biomechanical model having a plurality of mesh nodes using a multi-tissue mesher; (d) producing a set of source points by edge detection of the segmented pre-operative brain image; (e) producing a set of target points by edge detection of the intra-operative brain image data: (f) wherein a deformation field may be represented by a deformation vector defined on the plurality of mesh nodes, point correspondence of the set of source points and set of target points may be represented by a correspondence matrix, and a resection region may be represented by a connected submesh; and (g) processing the segmented pre-operative brain image using a three-variable solver, the variables comprising the deformation field, the point

correspondence; and the resection region,

Optionally, this device may be such thai the deformation field may be represented by a deformation vector defined on the plurality of mesh nodes by a piece-wise linear function regularized by a strain energy of a heterogenous biomecha cal model; the point correspondence of the set of source points and set of target points is represented by a correspondence matrix that is a fuzzy assign matrix, and the resection region is represented by a maximal connected tetrahedral submesh. The processing of the segmented pre- operative brain image may be implemented by using a three-variable solver method comprises resolving the three variables with a Nested Expectation and Maximization (NEM) framework. Further, the method may include iterative ly rejecting element outliers of the resection region and point outliers using the NEM framework; and approaching the remaining non-resection brain domain and resolution correspondence.

The device may generate a multi-tissue tetrahedral mesh by (i) starting from a homogeneous body-centered-cubic mesh; (ii) identifying a coarse multi-tissue mesh by assigning each tetrahedron with a specific tissue label; and (iii) deforming the coarse multi-tissue mesh surface to tissue boundaries defined the multi-label image,

Embodiments may include a system for performing a process, with the process being implemented by a computer system operatively connected to at least one data storage device. The data storage device may store image data of pre-operative image data of a patients head and tra-operative image data of the patient's head, with the system including (i) at least one computer and at least one computer readable medium storing thereon computer code which when executed by the at least one computer performs a method, the method comprising the at least one computer: (a) accessing the at least one data storage device and extracting preoperative brain image data using a brain extraction tool and delineating a ventricle boundary to produce a segmented mask and a segmented pre-operative brain image; (b) accessing the at least one data storage device and extracting mtra-operative brain image data using the brain extraction tool; (c) producing a heterogeneous biomechanical model having a plurality of mesh nodes using a multi-tissue mesher; (d) producing a set of source points by edge detection of the segmented pre-operative brain image; (e) producing a set of target points by edge detection of the intra-operative brain image data; (f) wherein a deformation field may be represented by a deformation vector defined on the plurality of mesh nodes, point correspondence of the set of source points and set of target points may be represented by a correspondence matrix, and a resection region may be represented by a connected submesh; and (g) processing the segmented pre-operative brain image using a three-variable solver, the variables comprising the deformation field, the point correspondence; and the resection region,

Optionally, the system may be such that the deformation field may be represented by a deformation vector defined on the plurality of mesh nodes by a piece-wise linear function regularized by a strain energy of a heterogenous biomechamcai model; the point correspondence of the set of source points and set of target points is represen ted by a correspondence matrix that is a fuzzy assign matrix, and the resection region is represented by a maximal connected tetrahedral submesh. The processing of the segmented preoperative brain image may be implemented by using a three-variable solver method comprises resolving the three variables with a Nested Expectation and Maximization (MEM) framework. Further, the system may include a method that involves iterative!}' rejecting element outliers of the resection region and point outliers using the NEM framework; and approaching the remaining non-resection brain domain and resolution correspondence. The system may generate a multi-tissue tetrahedral mesh generated by starting from a homogeneous hody-eentered-eubic mesh; identifying a coarse multi-tissue mesh by assigning each tetrahedron with a specific tissue label; and deforming the coarse multi- tissue mesh surface to tissue boundaries defined the multi-label image. Embodiments ma include an article of manufacture for configuring computer. For example, embodiments may include a non-transitory recording medium storing a program having computer code that configures a computer to implement a method, (with the computer operatively connected to at least one data storage device in which is stored pre-operative image data of a patients head and inlra-operative image data of the patient's head), such that the computer when configured by the non-transitory recording medium and executing the computer code performs the following method:

(i) accessing the at least one data storage device and extracting pre-operative brain image data from the pre-operative image data using a brain extraction tool and delineating a ventricle boundary to produce a segmented mask and a segmented preoperative brain image;

(is) accessing the at least one data storage device and extracting intraoperative brain image data from the intra-operative image data using the brain extraction tool;

(iii) producing a heterogeneous biomechanics! model having a plurality of mesh nodes using a multi-tissue mesher:

(iv) producing a set of source points by edge detection of the segmented preoperative brain image;

(v) producing a set of target points by edge detection of the intra-operative brain image data;

(vi) wherein a deformation field may be represented by a deformation vector defined on the plurality of mesh nodes, point correspondence of the set of source points and set of target points may he represented by a correspondence matrix, and a resection region may be represented by a connected submesh; and

(v) processing the segmented pre-operative brain image using a three-variable solver, the variables comprising the deformation field, the point correspondence; and the resection region.

The non-transitory recording medium may be such that the deformation field may be represented by a deformation vector defined on the plurality of mesh nodes by a piece-wise linear function regularized by a strain energy of a heterogenous biomechanical model; the point correspondence of the set of source points and set of target points is represented by a correspondence matrix that is a fuzzy assign matrix, and the resection region is represented by a maximal connected tetrahedral submesh. Optionally, the processing of the segmented pre-operative brain image may be implemented by using a three-variable solver method comprises resolving the three variables with a Nested Expectation and Maximization (NEM) framework. In addition, the method may involve iierativeiy rejecting element outliers of the resection region and point outliers using the NEM framework: and approaching the remaining non-resection brain domain and resolution correspondence.

The non-transitory recording medium may be such that the mesh is a multi-tissue tetrahedral mesh generated by starting from a homogeneous body-centered-cubic mesh, identifying a coarse multi-tissue mesh by assigning each tetrahedron with a specific tissue label, and deforming the coarse multi-tissue mesh surface to tissue boundaries defined the multi-label Image,

Non-limiting exemplary embodiments of systems, methods, and recording mediums which can implement the embodiments disclosed herein may be found in U.S. Patent Application Serial No. 13/31 1,335, filed on December 5, 201 1 , which claims priority to Provisional Application No. 61/419,632 filed on December 3, 2010, and to Provisional Application No. 61 /482,797 filed on May 5, 2011, the entirety of each of which is incorporated herein by reference. Non-limiting exemplary embodiments of systems, methods, and recording mediums which can implement the embodiments disclosed herein may also be found in U.S. Patent

Application Serial No. 13/31 1, 289, filed on December 5, 2 1 1, which claims priority to Provisional Application No. 61/41 ,603 filed on December 3, 20.1 , the entirety' of each of which is incorporated herein by reference,

While the invention has been described and illustrated with reference to certain preferred embodiments herein, other embodiments are possible. Additionally, as such, the foregoing illustrative embodiments, examples, features, advantages, and attendant advantages are not meant to be limiting of the present invention, as the invention may be practiced according to various alternative embodiments, as well as without necessarily providing, for example, one or more of the features, advantages, and attendant advantages that may be provided by the foregoing illustrative embodiments.

Systems and modules described herein may comprise software, firmware, hardware, or any combination(s) of software, firmware, or hardware suitable for the purposes described herein. Software and other modules may reside on servers, workstations, personal computers, computerized tablets, PDAs, and other devices suitable for the purposes described herein. Software and other modules may be accessible via local non-transitory memory, via a network or cloud, via a browser or other application in an application service provider or software as a service context, or via other means suitable for the purposes described herein. Data structures described herein may comprise computer files, variables, programming arrays, programming structures, or any electronic information storage schemes or methods, or any combinations thereof, suitable for the purposes described herein. User interface elements described herein may comprise elements from graphical user interfaces, command line interfaces, and other interfaces suitable for the purposes described herein. Except to the extent necessary or inherent in the processes themselves, no particular order to steps or stages of methods or processes described in this disclosure, including the Figures, is implied. In many cases the order of process steps may be varied, and various illustrative steps may be combined, altered, or omitted, without, changing the purpose, effect or import of the methods described.

Accordingly, while the invention has been described and illustrated in connection with preferred embodiments, many variations and modifications as will be evident to those skilled in this art may be made without departing from the scope of the invention, and the invention is thus not to be limited to the precise details of methodology or construction set forth above, as such variations and modification are intended to be included within the scope of the inventiors, Therefore, the scope of the appended claims shoisld not be limited to the description and illustrations of the embodiments contained herein,

Claims

What is claimed is:

1 , A method for the processing of imaging data of the head of a patient the imaging data produced by a magnetic resonance imaging scanner to collect pre-operative image data of the head and an intraoperative image data of the head stored in at least one non-transitory data storage device, die method comprising: accessing the at least one data storage device and extracting pre-operative brain image data from the pre-operative image data using a brain extraction tool and delineating a ventricle boundary to produce a segmented mask and a segmented pre-operative brain image;

accessing the at least one data storage device and extracting intra-operative brain linage data from the intra-operative image data using the brain extraction tool;

producing a heterogeneous biomechanieal model having a plurality of mesh nodes using a multi-tissue rnesher;

producing a set of source points by edge detection of the segmented pre-operative brain image;

producing a set of target points by edge detection of the intra-operative brain image data;

wherein a deformation field may be represented by a deformation vector defined on the plurality of mesh nodes, point correspondence of the set of source points and set of target points may be represented by a correspondence matrix, and a resection region may be represented by a connected submesh; and

processing the segmented pre-operative brain image using a three-variable solver, the variables comprising the deiormation field, the point correspondence; and the resection region.

2. The method of claim I , wherein:

the deformation field may be represented by a deformation vector defined on the plurality of mesh nodes by a piece-wise linear function regularized by a strain energy of a heterogenous biomechanieal model;

the point correspondence of the set of source points and set of target points is represented by a correspondence matrix that is a fuzzy assign matrix, and the resection region is represented by a maximal connected tetrahedral submesh.

3. The method of claim 1 , wherein the mesh is a multi-tissue tetrahedral mesh generated by:

starting from a homogeneous body-eentered-cubic mesh;

identifying a coarse multi-tissue mesh by assigning each tetrahedron with a specific tissue label; and

deforming the coarse multi-tissue mesh surface to tissue boundaries defined the multi-labei image.

4. The method of claim 2, wherein the processing of the segmented preoperative brain image using a three-variable solver method comprises resolving the three variables with a Nested Expectation and Maximization (NEM) framework.

5. The method of claim 4, wherein the method comprises:

iteratively rejecting element outliers of the resection region and point outliers using the NEM framework; and

approaching the remaining non-resection brain domain and resolution

correspondence.

6. A neurosurgery device guided by imaging data of the head of a patient, the neurosurgery device comprising:

a magnetic resonant imaging scanner for producing imaging data comprising a pre- operative linage of the head and an intra-operative image of the head;

at least one non-transitory data storage device in communication with the magnetic resonant imaging scanner for receiving and storing preoperative image data and an intraoperative image data;

a computer system having a specially configured computer processor in

communication with the non-transitory data storage device, the compuier system configured to process the preoperative image data and an intra-operative image data using a method comprising: accessing the at least one data storage device and extracting preoperative brain image data using a brain extraction tool and delineating a veotricie boundary to produce a segmented mask and a segmented pre-operative brain image;

accessing the at least one data storage device and extracting intra-operattve brain image data using the brain extraction tool;

producing a heterogeneous biomechanics! model having a plurality of mesh nodes using a multi-tissue mesher;

producing a set of source points by edge detection of the segmented preoperati e brain image;

producing a set of target points by edge detection of the intra-operative brain image data;

wherein a deformation field may be represented by a deformation vector defined on the plurality of mesh nodes, point correspondence of the set of source points and set of target points may be represented by a correspondence matrix, and a resection region may be represented by a connected submesh; and

processing the segmented pre-operative brain image using a three-variable solver, the variables comprising the deformation field, the point correspondence; and the resection region,

7. The device of claim 6, wherein:

the deformation field may be represented by a deformation vector defined on the plurality of mesh nodes by a piece-wise linear function regularized by a strain energy of a heterogenous biomechanical model;

the point correspondence of the set of source points and set of target points is represented by a correspondence matrix that is a fuzzy assign matrix, and

the resection region is represented by a maximal connected tetrahedral submesh.

8. The device of claim 6, wherein the mesh is a multi-tissue tetrahedral mesh generated by:

starting from a homogeneous body-centered-cubie mesh;

identifying a coarse multi-tissue mesh by assigning each tetrahedron with a specific tissue label; and

deforming the coarse multi-tissue mesh surface to tissue boundaries defined the multi-label image.

9. The device of claim 7, wherein the processing of the segmented pre-operative brain image using a three-variable solver method comprises resolving the three variables with a Nested Expectation and Maximization (NEM) framework.

10. The device of claim 9, wherein the method comprises:

iteratively rejecting element outliers of the resection region and point outliers using the NEM framework; and

approaching th remaining non-resection brain domain and resolution

correspondence.

11. A non-transitory recording medium storing a program having computer code that configures a computer to iinpiement a method, the computer operaiively connected to at least one data storage device in which is stored pre-operative image data of a patients head and intra-operative image data of the patient's head, the computer when configured by the non- transitory recording medium and executing the computer code performs the method comprising;

accessing the at least one data storage device and extracting pre-operative brain image data from the pre-operative image data using a brain extraction tool and delineating a ventricle boundary to produce a segmented mask and a segmented pre-operative brain image;

accessing the at least one data storage device and extracting intra-operative brain image data from the intra-operative image data using the brain extraction tool;

producing a heterogeneous biomechanics! model having a plurality of mesh nodes using a multi-tissue mesher;

producing a set of source points by edge detection of the segmented pre-operative brain image;

producing a set of target points by edge detection of the intra-operative brain image data;

wherein a deformation field may be represented by a deformation vector defined on the plurality of mesh nodes, point correspondence of the set of source points and set of target points may be represented by a correspondence matrix, and a resection region may be represented by a connected subrnesh; and

processing the segmented pre-operative brain image using a three-variable solver, the variables comprising the deformation field, the point correspondence; and the resection region.

12. The non-transitory recording medium of claim 1 1 , wherein:

the deformation field may be represented by a deformation vector defined on the plurality of mesh nodes by a piece-wise linear function regularized by a strain energy of a heterogenous biomechanical model;

the point correspondence of the set of source points and set of target points is represented by a correspondence matrix that is a fuzzy assign matrix, and

the resection region is represented by a maximal connected tetrahedral submesh.

13. The non-transitory recording medium of claim 1 1, wherein the mesh is a multi- tissue tetrahedral mesh generated by:

starting from a homogeneous body-centered-cubic mesh;

identifying a coarse multi-tissue mesh by assigning each tetrahedron with a specific tissue iabei; and

deforming the coarse multi-tissue mesh surface to tissue boundaries defined the multi-label image.

14. The non-transitory recording medium of claim 12, wherein the processing of the segmented pre-operative brain image using a three-variable solver method comprises resolving the three variables with a Nested Expectation and Maximization (NEM) framework.

15. The non -transitory recording medium of claim 14, wherein the method comprises:

iterative!)' rejecting element outliers of the resection region and point outliers using the NEM framework: and

approaching the remaining non-resection brain domain and resolution

correspondence,

16. A system for performing a process, the process being implemented by a computer system operative! y connected to at least one data storage device in which is stored image data pre-operative image data of a patients head and intra-operative image data of the patient's head, the system comprising: at ieast one computer and at least one non-transitory computer readable medium storing thereon computer code which when executed by the at least one computer performs a method, the method comprising the at least one computer:

accessing the at ieast one data storage device and extracting pre-operative brain image data using a brain extraction tool and delineating a ventricle boundary to produce a segmented mask and a segmented pre-operative brain image;

accessing the at least one data storage device and extracting in Ira-operative brain image data using the brain extraction tool;

producing a heterogeneous biomechanics! model having a plurality of mesh nodes using a multi-tissue mesher;

producing a set of source points by edge detection of the segmented pre-operative brain image;

producing a set of target points by edge detection of the intra-operative brain image data;

wherein a deformation field may be represented by a deformation vector defined on the plurality of mesh nodes, point correspondence of the set of source points and set of target points may be represented by a correspondence matrix, and a resection region may be represented by a connected subrnesh; and

processing the segmented pre-operative brain image using a three-variable solver, the variables comprising the deformation field, the point correspondence; and the resection region.

17, The system of claim 16, wherein:

the deformation field may be represented by a deformation vector defined on the plurality of mesh nodes by a piece-wise linear function regularized by a strain energy of a heterogenous biomechanics! model;

the point correspondence of the set of source points and set of target points is represented by a correspondence matrix that is a fuzzy assign matrix, and

the resection region is represented by a maximal connected tetrahedral subrnesh.

18, The system of claim 16, wherein the mesh is a ulti -tissue tetrahedral mesh generated by:

starting from a homogeneous body-centered-cuhic mesh: identifying a coarse multi-tissue mesh by assigning each tetrahedron with a specific tissue label; and

deforming the coarse multi-tissue mesh surface to tissue boundaries defined the multi-label image.

19. The system of claim 17, wherein the processing of the segmented pre-operative brain image using a three-variable solver method comprises resolving the three variables with a Nested Expectation and Maximization (NEM) framework.

20. The device of claim 19, wherein the method comprises:

iterative!}' rejecting element outliers of the resection region and point outliers using the NEM framework; and

approaching the remaining non-resection brain domain and resolution

correspondence.