[go: up one dir, main page]

WO2010037850A2 - Procédé de traitement d'image pour estimer une déformation d'un cerveau d'un patient - Google Patents

Procédé de traitement d'image pour estimer une déformation d'un cerveau d'un patient Download PDF

Info

Publication number
WO2010037850A2
WO2010037850A2 PCT/EP2009/062826 EP2009062826W WO2010037850A2 WO 2010037850 A2 WO2010037850 A2 WO 2010037850A2 EP 2009062826 W EP2009062826 W EP 2009062826W WO 2010037850 A2 WO2010037850 A2 WO 2010037850A2
Authority
WO
WIPO (PCT)
Prior art keywords
patient
brain
points
tree
registration
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Ceased
Application number
PCT/EP2009/062826
Other languages
English (en)
Other versions
WO2010037850A3 (fr
Inventor
Marek Bucki
Yohan Payan
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Centre National de la Recherche Scientifique CNRS
Universite Joseph Fourier Grenoble 1
Original Assignee
Centre National de la Recherche Scientifique CNRS
Universite Joseph Fourier Grenoble 1
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Centre National de la Recherche Scientifique CNRS, Universite Joseph Fourier Grenoble 1 filed Critical Centre National de la Recherche Scientifique CNRS
Priority to US13/122,092 priority Critical patent/US9008381B2/en
Priority to EP09783691A priority patent/EP2353144A2/fr
Publication of WO2010037850A2 publication Critical patent/WO2010037850A2/fr
Publication of WO2010037850A3 publication Critical patent/WO2010037850A3/fr
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/38—Registration of image sequences
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/08—Clinical applications
    • A61B8/0808—Clinical applications for diagnosis of the brain
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30016—Brain

Definitions

  • the present invention relates to the field of computer-assisted neurosurgery, and more particularly to a system and a method for the precise localization of a target.
  • a deformation of the brain tissue may occur after an opening has been made in the skull of the patient.
  • This deformation can be related to physical phenomena (gravity, loss of cerebrospinal fluid, action of the neurosurgeon, etc.) or to physiological phenomena (swelling due to osmotic drugs, anesthetics, etc.), some of which are not known at this time. day.
  • neuro-navigation systems In order to compensate for this deformation of the patient's brain and to retrieve the patient's preoperative data during the surgical procedure, neuro-navigation systems have been developed.
  • a first image of the cortical surface of the patient is acquired before the surgical operation.
  • a second image of the cortical surface of the patient is acquired.
  • System processing means compares these two images and estimates a volume deformation field that aligns the points in the two images to compensate for deformation of the patient's brain.
  • An object of the present invention is to provide a system and a method for compensating for the deformation of the patient's brain making it possible to overcome at least one of the aforementioned drawbacks.
  • the invention provides a passive navigation method for estimating a deformation of a brain of a patient, the method comprising: the treatment of a three-dimensional image of the patient's brain, acquired before a surgical procedure, to obtain a reference cerebral arterial tree of the patient,
  • the initial vascular tree is elastically recalibrated to its deformed, segmented configuration in the US Focused Doppler images, and the resulting displacement field is then extended to the entire organ volume through a biomechanical model specific to the patient.
  • the global deformation thus calculated makes it possible to update the preoperative data and correct the initial surgical planning.
  • the invention provides a fast, robust and accurate response to the surgeon and gives him the opportunity to validate the relevance of the deformation calculated before any modification of the schedule.
  • the method may also include calculating disparities between points of the generated brain image and one or more two-dimensional control images acquired during the operation;
  • the three-dimensional image acquired before the operation may be magnetic resonance angiography
  • the two-dimensional images acquired during the operation may be ultrasound images in Doppler mode
  • the treatment of the three-dimensional image to obtain the reference cerebral arterial tree comprises a step of orienting the three-dimensional image with respect to the patient's head after the installation of the latter on an operating table,
  • the step of orienting the three-dimensional image comprises: locating a set of calibration points on the patient's head after the patient has installed it on the operating table, where the detection of their radiometric counterparts in the three-dimensional image, and o determining rigid transformations for moving from the points of the set of calibration points to the positions of their counterparts in the three-dimensional image, - the calibration points may comprise points of the the surface of the patient's head after installing it on the operating table,
  • the treatment of the two-dimensional images to reconstitute at least partially a current cerebral arterial tree of the patient comprises a segmentation of the two-dimensional images to obtain the median position of the vessels and / or arteries of the brain so as to reconstitute at least partially the current cerebral arterial tree of the patient
  • the step of determining a deformation field comprises the determination of nonlinear transformations, called elastic transformations, allowing the passage of the positions of the points of the cerebral arterial reference tree to the positions of the corresponding points of the current arterial tree.
  • the invention also relates to a system for determining deformation of the brain of a patient, the system comprising means for implementing the method described above.
  • FIG. 1 illustrates an embodiment of the method according to the invention.
  • FIG. 2 illustrates an embodiment of the system for implementing the method according to the invention
  • FIG. 3 illustrates a deformation of a unitary square with a folding of space
  • FIG. 4 illustrates in two dimensions a multi-scale iterative approach described in relation with the step of resetting organic tissue
  • FIG. 5 illustrates cells inserted in the margin of a discretization for two successive refinement levels
  • FIG. 6 illustrates a function of form w
  • FIG. 7 illustrates an example of approximation of a resetting energy curve by a parable
  • FIG. 8 illustrates three steps of the approximation of the inverse of the resetting function
  • FIG. 9 illustrates steps of direct and indirect filtering.
  • the present invention makes it possible to take this deformation into account and to calculate a rectified position of the anatomical structures of the brain in order to locate ancillaries.
  • MRA magnetic resonance angiography
  • the invention proposes a technique for updating the preoperative data according to the deformation of the brain during the procedure. To do this, the invention proposes to follow during the intervention the deformation of the particular anatomical structure that is the cerebral vascular tree.
  • the displacements of the vascular tree are representative of the global deformations of the tissues.
  • the arteries and arterioles irrigating the nerve tissues are present both in the deep parts of the brain, but also on the surface, and more particularly in the vicinity of the tumor whose development passes through an angiogenesis and therefore the formation of new ships.
  • the invention then proposes to use a biomechanical modeling of the patient's brain in order to extrapolate a global deformation from the displacement of the cerebrovascular tree.
  • the global displacement field thus obtained is then passed on to all the preoperative data in order to update them.
  • the image processing method for estimating a displacement of a brain of a patient comprises the following steps:
  • generation 500 from the displacement of the estimated brain, of at least one image of the patient's brain in which the displacement of the brain is compensated.
  • the method comprises a step 100 of processing an image of the brain of the acquired patient prior to the surgical procedure.
  • This treatment aims to locate the cerebrovascular tree in the image acquired before surgery.
  • This cerebral vascular tree in the image acquired before the surgical procedure corresponds to the vascular tree before the deformation of the patient's brain, and is hereinafter referred to as "cerebral vascular tree of reference".
  • the image acquired before the surgical procedure is a three-dimensional image.
  • this three-dimensional image is a magnetic resonance angiography.
  • Magnetic resonance angiography reveals important contrasts between tissues. The analysis of magnetic resonance angiography allows to distinguish with precision the nature of the observed tissues.
  • the data it contains make it possible to locate the reference cerebral vascular tree.
  • the technique allowing the acquisition of magnetic resonance angiography has the advantage of being non-invasive and therefore without danger for the patient, unlike, for example, techniques emitting ionizing radiation such as tomography.
  • the treatment of the three-dimensional image to obtain the reference cerebral arterial tree comprises a step of orienting the three-dimensional image with respect to the patient's head.
  • the spatial position of the reference vascular tree is determined relative to a repository of the physical space of the patient.
  • the reference cerebral vascular tree is transposed into the patient's physical reference system, using a so-called rigid registration technique which relies on the matching of calibration points with their segmented equivalents in the patient. Magnetic resonance angiography acquired before the operation.
  • rigid registration is intended to mean a geometrical transformation composed of a rotation and a translation.
  • the calibration points used may be points or remarkable anatomical surfaces such as the tip of the nose, the eyebrows, the surface of the cheeks or forehead.
  • the calibration points may also be defined by a plurality of adhesive pellets (at least three and preferably ten) adhered to the patient's skin before the acquisition of magnetic resonance angiography. Thus, these pellets are acquired in angiography at the same time as the patient's brain.
  • the radiometric counterparts of the calibration points are identifiable in magnetic resonance angiography. Once the patient is on the operating table, the surgeon tells a tracking system the position of each calibration point using a probe.
  • the location system can be an optical location system (including for example a stereoscopic camera), or a magnetic locating system, or any other location system known to those skilled in the art.
  • the system is able to calculate the spatial position of the patient. reference vascular tree relative to the patient's spatial position.
  • the method also includes a step 200 of processing a plurality of two-dimensional images of the patient's brain acquired during the operation to at least partially reconstitute a cerebral arterial tree of the so-called patient
  • the surgeon performs a two-dimensional image acquisition.
  • the term "plurality of two-dimensional images” refers to 2D images in a determined volume (3D) whose acquisition planes are intersecting or parallel.
  • These 2D images can be obtained from any medical device for acquiring 2D or 3D images such as 3D ultrasound images for example.
  • these two-dimensional images are ultrasound images acquired in Doppler mode, using a probe located by the localization system.
  • Ultrasound images in Doppler mode can visualize blood flow.
  • the Doppler ultrasound technique has the advantage of being benign for the patient. Intraoperative Doppler ultrasound is therefore a good tool for exploring the movements of the brain during the procedure.
  • This visualization mode is based on the Doppler effect and makes it possible to locate, by analyzing the frequency variations of the ultrasounds emitted, displacements within the tissues, such as blood flows. The segmentation of this type of image is made easier because of the color coding used to represent the velocities of the flows. This modality thus makes it possible to visualize in the volume of the encephalon the blood flows that take place along the vascular tree and thus to identify the position of the vessels that irrigate the parenchyma.
  • the acquisition of the two-dimensional images is carried out for example as follows.
  • the localized probe is brought into contact with the patient's brain through the craniotomy performed by the surgeon.
  • a rotational movement is manually printed at the localized probe so that the ultrasound plane scans the widest possible part of the brain around the tumor.
  • the localized probe includes a marker enabling the location system to define its position and orientation relative to the patient's physical reference frame, and thus to spatially reposition the two-dimensional images in the patient's physical reference frame.
  • this After acquiring a series of two-dimensional images, these are processed to locate the patient's current cerebrovascular tree. More precisely, the centers of the vessels contained in each two-dimensional image are selected so as to obtain a cloud of points. This cloud of points at least partially represents the current cerebral vascular tree of the patient.
  • the reference vascular tree and current are matched to determine a deformational field representing the displacement of the current vascular tree relative to the vascular tree. reference.
  • the reference vascular tree is recalibrated with respect to the current vascular tree.
  • elastic registration is meant the fact of determining, from datasets which represent two digital images of the same zone, a transformation making it possible to pass from the reference image to the current image. In other words, it is about superimposing "at best" structures contained in two comparable images.
  • the data from the three-dimensional image and two-dimensional images are of different natures.
  • the reference vascular tree obtained from the preoperative data is almost continuous, the whole of the patient's head is contained in the data volume, and the segmentation performed does not generate, so to speak, artifacts in the reconstruction of the cerebrovascular tree.
  • the current vascular tree obtained from two-dimensional intraoperative images is composed of a cloud of points, these two-dimensional images cover only a fraction of the volume of the brain.
  • two-dimensional images can be more or less evenly distributed in space, some regions with multiple acquisitions while others have no acquisition.
  • step 300 of registration of cerebrovascular trees is therefore to estimate an elastic transformation between:
  • the strain field is then calculated from the correspondence established between the reference vascular tree and the current vascular tree.
  • Each vector of the deformation field originates from a point of the reference vascular tree and indicates a point of the current vascular tree, deformed.
  • the elastic registration is guided by the minimization of a registration energy. This registration energy quantifies the similarity between the current vascular tree and the reference vascular tree.
  • Optimal registration is achieved for zero registration energy.
  • the objective of the elastic registration step is to iteratively minimize the value of this registration energy by changing the points of the current vascular tree towards the points of the reference vascular tree.
  • the elastic registration function is a nonlinear transformation to move from the current vascular tree to the reference vascular tree.
  • conditions related to the nature of the data that are to be recalibrated are determined to ensure the respect of essential physical criteria.
  • a first condition that the elastic registration function must satisfy relates to non-folding of the space. Indeed, the deformation of the patient's brain can not induce the folding of the patient's brain on itself.
  • a second condition that the elastic registration function must satisfy relates to the continuity of the propagation of the deformation. Indeed, reference and current vascular trees are observations of the physical phenomenon of brain deformation. In a similar way to the process of deformation of the patient's brain, without tearing or tissue resection, the registration function must be continuous and even continuously differentiable i.e. differentiable at any point in the field and not having a differential discontinuity.
  • a third condition that the registration function must satisfy is that two distinct points of the current vascular tree can not have the same image in the reference vascular tree, and vice versa.
  • the elastic registration function as proposed in the invention has the advantage of being invertible with the desired degree of accuracy. Once this adjustment function has been calculated, a post-processing of the result produced by the registration makes it possible to eliminate the aberrations due to the presence of noise in the data of the two- or three-dimensional per-operative images.
  • the set of displacements of the vascular tree finally retained forms a deformation field representing the displacement of the current vascular tree relative to the reference vascular tree.
  • the registration step can be applied to the registration of structures other than a cerebrovascular tree.
  • the elastic registration presented above can be applied to the registration of any type of organic tissue.
  • the deformation field of the determined vascular tree is applied to a biomechanical model of the soft tissues of the patient's brain to estimate the overall displacement of the patient's brain.
  • This mathematical model formalizes an a priori knowledge of the behavior of the organ and makes it possible to infer the global behavior under the local strain constraints deduced from the per-operative observation, called boundary conditions.
  • the biomechanical model also known as the deformable model
  • the biomechanical model extrapolates the movements of the vascular tree to the entire volume of the brain to simulate the effect of brain deformation.
  • the biomechanical model makes it possible to regulate the data. Indeed, the estimation of the displacement field of the vascular tree is tainted by errors due to the imaging mode used, the limitations of the segmentation algorithms as well as those of the calculation of the elastic registration carried out to match the current vascular tree and the reference vascular tree.
  • the biomechanical model makes it possible to harmonize the field of displacements by eliminating the errors in the input data of the biomechanical model. This makes the process resistant to the artifacts present in the noisy data collected during the procedure (i.e. two- or three-dimensional ultrasound images).
  • the originality of the process thus concerns the complementarity between a biomechanical model approximating the behaviors of the patient's brain tissues and a follow-up of anatomical structure potentially noisy but representative of the displacement of the brain thanks to the cerebrovascular tree.
  • a final step 600 of the method concerns the generation, from the global displacement of the estimated brain, of at least one image of the brain of the patient in which the movement of the brain is compensated.
  • This image is for example the three-dimensional image (or two-dimensional images in section of this three-dimensional image) acquired before the operation to which the estimated global deformation of the patient's brain is applied so as to update the preoperative data as a function of the estimated global deformation of the patient's brain.
  • the method may comprise a step of calculating the disparities between points of the generated brain image and one or more two-dimensional control images acquired during the operation.
  • the user can for example acquire one or more two-dimensional images of control of the patient's brain to calculate the positional error between the points of the two-dimensional control images and their counterparts in the image generated from the deformation. overall estimate of the patient's brain.
  • the user can scan the brain tissue in the region of interest while projecting in the two-dimensional control images the position of the vascular tree deformed by the model.
  • the superposition of two real-time Doppler data obtained by scanning the organ, on the one hand, and simulated Doppler signal, calculated from the incidence of localized echographic sections relative to the vascular tree deformed by the system, on the other hand, allows the user to assess the coincidence between the model and the reality of the patient.
  • the assessment of the quality of the deformation allows the user to make the decision whether or not to follow the indications of the system on the basis of the preoperative data updated.
  • An elastic registration is an application which minimizes a certain registration error, or rather, a registration energy defined between a set of source points, S and another set of destination points, D.
  • the energy of the registration R is the measure of a "similarity" between the set of transformed points R (S) and the target set D. When the two sets are merged then this energy becomes zero and the registration is completed. If the reshaping energy and the function space in which the R transformation is not defined correctly, then the problem of finding an optimal R is misplaced, ie the existence and uniqueness of R are not guarantees.
  • the registration energy can be freely defined to reflect the nature of the registration problem being addressed. For example, in the case of registration between two images, this measurement can be based on the differences between the values of the pixels, or voxels, of the source and destination images. In this case, the similarity between the S and D datasets is defined by the spatial proximity of the two point clouds. Without further knowledge of the nature of the manipulated data, it is possible, for example, to define a registration energy E reg as the average of the minimal Euclidean distances between R (S) and D.
  • Card (D) is the cardinal of the destination set D
  • d (R (s), D) represents the distance between the set D and the transformed by R of a point s of the set S
  • R (S) represents the distance between a point d of the set D and the transform by R of the set S.
  • the nonlinear transformation R sought is continuous and even continuously differentiable ie differentiable at any point in the domain and does not show any discontinuity. differential. This defines a first restriction on the search space of the registration functions
  • the volume defined by the transformed trihedron, noted dv is then obtained in a similar way.
  • R must be injective. This results in the fact that two distinct points can not have the same image by R and, physically, R must not bring two atoms into the same point of space. If R is an injection, then for each point reached a predecessor can be calculated.
  • the approach presented here is inspired by the mechanical modeling of solids based on finite element discretization.
  • the deformed point cloud is embedded in a resilient virtual solid. Since the displacement of points translates the physical properties of an underlying structure whose geometry is unknown, this virtual solid is arbitrarily assigned the form of a bounding box of points.
  • the set of admissible deformations is then defined as the successive combinations of elementary deformations of the enclosing solid, the effects of which propagate through its volume and modify the position of the internal points.
  • the elementary deformations mentioned here are defined below.
  • the 3D registration technique described here can easily be derived in a 2D registration technique by replacing the notion of "virtual solid" with "virtual plan".
  • the illustrations presented below are based on a 2D implementation of the registration to clarify the concepts.
  • the elastic registration function is assembled iteratively so that the registration energy E r ⁇ g decreases optimally at each step.
  • S 0 : S be the initial source point cloud
  • S 1 be the source point cloud obtained following the iteration i.
  • the domain of the virtual solid is discretized with a certain level of refinement in regular hexahedrons (or squares, in 2D) called "cells".
  • the elementary deformations of the solid are obtained by individually displacing the nodes of the grid thus formed.
  • V (n) The union of the cells surrounding a node n defines its "neighborhood", denoted V (n).
  • Each node of the discretized solid is considered and the gradient of the registration energy, as a function of the displacement of the node, is calculated.
  • the opposite of this vector defines the preferential displacement of the node, that is to say the one which, applied to the node, causes the greatest decrease of the energy of registration.
  • the one that allows the greatest decrease of the gradient is chosen and its preferential displacement is applied to it while the other nodes remain motionless.
  • the local deformation of the solid is calculated by propagating the displacement of the node through its neighborhood and the position of the source points it contains is modified accordingly, thus giving the new configuration of the point cloud Si + 1.
  • the solid returns to its initial configuration, at rest, and the source points Si + 1 are distributed in its volume.
  • a new iteration can then be calculated by taking the new configuration Si + 1 as a starting point.
  • the iterations stop, finally, when no significant decrease in energy can be obtained by moving the nodes of the grid of the solid.
  • the refinement iterations stop when the maximum degree of manually specified subdivision has been reached.
  • the dimension of the smallest cell reached during the refinement process defines the dimensions of the smallest structure present in S and which will not be correctly recalibrated if an equivalent structure is absent from D.
  • Figure 4 illustrates in 2D the multi-scale iterative approach described above.
  • the source points S are represented by the points referenced 1.
  • the set D is not represented, for the sake of clarity.
  • the best nodal displacement is applied to the solid and the source points are moved to S1. (referenced 4).
  • the solid then returns to its initial configuration (a). After several iterations, if no significant improvement in energy can be obtained, the level of refinement of the grid increases (c).
  • the energy gradients are calculated at the nine nodes of the new grid.
  • the best displacement is applied to the solid and the new source point configuration S2 (referenced 6) is calculated.
  • FIG. 5 shows the additional cells 7 inserted in the margin of the domain discretization for two successive refinement levels.
  • G1 ..., GJ be the successive refinement levels, G1 being the highest level
  • a nodal value (displacement or other quantity) is interpolated through the discretized volume by means of shape functions .
  • the value of the shape function associated with a node and evaluated at a point in the domain gives the proportion of the nodal value transferred at that point. So the value of w at its node is 1.
  • the same mechanism for propagating the displacement of a node through its neighborhood will be used.
  • w support is ui is a compact
  • Id3> ⁇ 3 is the identity matrix 3 * 3.
  • the non-folding of r can thus be guaranteed by limiting the amplitude of the nodal displacements.
  • the maximum norm of nodal displacement is given by the inverse of the maximum standard of the gradient of the considered shape function w.
  • V: V (n) be the closed union of the cells surrounding node n.
  • V (n) V (n) be the closed union of the cells surrounding node n.
  • the form function w is defined on the union of neighboring cells around node n by a change of variable and a scaling in order to adapt its canonical definition domain L> J to the dimensions of the cells in the discretization grid of the virtual solid.
  • Figure 6b shows the variable changes needed to assemble from the canonical form function the function of form w over a neighborhood of 2 * 2 cells of the plane.
  • Figure 6c shows the layout of the shape function thus obtained
  • a Euclidean distance map is computed from the points of D.
  • the distance map is computed on a 3D grid discretizing the space around D. At each vertex of the grid gi, the minimum distance to D, d (gi, D) is stored. This data structure is initialized during the preoperative phase, just after the segmentation of the vascular tree D in the ARM images.
  • the contents of the distance card are then saved as a binary file to be read at the beginning of the intervention by the system, before any calculation of elastic registration.
  • the best displacement for each node of the grid that discretizes the virtual solid is calculated using the gradient of Ereg.
  • the calculation of the nodal displacement is carried out in two stages: 1. The gradient of Ereg with respect to the nodal displacements is estimated by finite differences; 2. A search along the direction of greatest descent is performed to estimate the optimal pitch;
  • the minimized function is defined by: or is the elementary registration function associated with the node n, and R is the global registration function assembled before the current iteration.
  • the length of the maximum pitch allowed by the non-folding constraint is 10% of L, which makes it possible to define the maximum vector of displacement of the node n in the direction of descent on f:
  • the best nodal displacement in the sense of minimizing the energy of registration corresponds therefore to the minimum of the function.
  • a suitable estimate of this minimum is obtained by approximating on the interval [0, 1] the function by a parable
  • the three coefficients of parabola A, B and C are determined by the local information on the behavior of f as well as an evaluation of its value at the maximum pitch: . So :
  • FIG. 7 shows a representative example of approximation of the registration energy curve 11 by a parabola P 12 over the interval [0, 1]. Measurements of the difference between the two curves, made on a large number of data sets during the course of the registration procedure shows the qualities of this approximation.
  • the bottom curve 11 is obtained by sampling the resetting energy on the interval [0, 1], that of the top 12 is the plot of the parabola P.
  • the values of t for which these curves reach their minimum are almost identical.
  • the minimization of the energy is thus carried out, step by step, until no significant decrease can be obtained by applying a displacement to one of the nodes. of the mesh.
  • indices i and j will not be specified in the following, and we considered more generally an elementary function r, associated with the displacement of the node n.
  • p r-1 (q) is calculated as follows.
  • p can be written as P ⁇ * ⁇ + * p ⁇
  • tp is the root of the polynomial function f over the interval [ta, tb].
  • f the search interval can be restricted to [max ⁇ -1, ta ⁇ , O].
  • this stopping condition is checked when this inequality is a reliable stopping criterion defined in the undeformed space.
  • the final error can thus be increased by the sum of the errors propagated:
  • each Lipschitz constant ki must be evaluated only once since it depends only on the amplitude of the displacement applied to the node associated with the elementary function ri and the size of the cell at the level of refinement considered. .
  • the role of filtering is to eliminate associations of points that do not satisfy a given accuracy criterion.
  • dmax be the maximum acceptable distance between a source point and its homologue in D.
  • the direct filtering makes it possible to eliminate the points of S that do not fulfill the condition d (R (s), D) ⁇ dmax.
  • the filtered source point set is S f
  • Figure 9a illustrates the direct filtering step. As can be seen in this figure, source points that have not reached D are eliminated.
  • the aim of the registration procedure is to generate a displacement field corresponding to the deformation of the structure of the vascular tree (in the presentation made here, but more generally of the deformed organic tissue) defined from the preoperative images by D , and which bring it into the current configuration S, obtained by analysis of US peroperative Doppler images.
  • Reverse filtering makes it possible to verify that the displacement field obtained by the inversion of the registration function R actually transforms the points from D to S.
  • Figure 9b illustrates the reverse filtering step. As this figure illustrates, displacements of D that do not reach S are eliminated.
  • the elastic registration described in detail above makes it possible to construct an association between two observations of the pre- and perioperative cerebrovascular tree, of very distinct natures.
  • the calculation is guided by the minimization of the registration energy which quantifies the similarity between the two structures.
  • the optimal registration being reached for zero energy, the objective of the registration algorithm is to minimize iteratively. value by changing the current structure, called 'source', to the reference structure, called 'destination'.
  • the speed of the calculation is obtained thanks to the preliminary computation of a distance map on which the evaluation of the energy of registration and the determination of the direction of stronger descent are based, calculated with each iteration of the process of minimization. .
  • the mathematical properties of the registration function family generated by the method described above ensure the respect of essential physical criteria, such as the non-folding of space, the non-overlapping of material and the continuity of the propagation of propagation. the deformation.
  • This elastic registration also has the advantage of being invertible with the desired degree of precision.
  • FIG. 2 illustrates one embodiment of the system for implementing the method described above.
  • the system comprises locating means 10, acquisition means
  • processing means 30 processing means 30 and display means 40.
  • the location means allow the location in space of the instruments used by the user during the surgical procedure. These locating means may comprise, for example, a stereoscopic camera and markers associated with the different objects to be located.
  • the locating means 10 also allow the location in the space of the patient's head during the procedure, in particular for the spatial repositioning of the reference vascular tree calculated from the magnetic resonance angiography.
  • the acquisition means 20 allow the acquisition of two-dimensional images during the operation.
  • These acquisition means 20 comprise for example an ultrasonic probe locatable by the locating means 10.
  • the display means allow for example the display of preoperative data updated.
  • These display means may comprise a screen.
  • the processing means make it possible to implement the various calculation steps of the method, and in particular: the treatment step 100 of the three-dimensional image of the patient's brain, acquired before a surgical procedure, to obtain a reference cerebral arterial tree of the patient,
  • the treatment step 200 of the two-dimensional images of the patient's brain, acquired during the operation, for reconstituting at least partially the current cerebral arterial tree of the patient the determination step 300, from the mapping of the reference and current cerebral arterial tree, the deformation field of the vascular tree representing the displacement of the current vascular tree relative to the reference vascular tree, and
  • These processing means may be a computer and include control means for controlling the various hardware elements of the system.

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)
  • Image Processing (AREA)

Abstract

L'invention concerne un procédé de traitement d'image pour estimer une déformation d'un cerveau d'un patient, le procédé comprenantle traitement (100) d'une image tridimensionnelle du cerveau du patient, acquise avant une intervention chirurgicale, pour obtenir une arborescence artérielle cérébrale de référence du patient, le traitement (200) d'images bidimensionnelles du cerveau du patient, acquises pendant l'opération, pour reconstituer au moins partiellement une arborescence artérielle cérébrale courante du patient, la détermination (300), à partir de la mise en correspondance des arborescences artérielles cérébrales de référence et courante, d'un champ de déformation de l'arbre vasculaire représentant le déplacement de l'arbre vasculaire courant par rapport à l'arbre vasculaire de référence, l'application (400) du champ de déformation de l'arbre vasculaire déterminé à un modèle biomécanique du cerveau du patient pour estimer le déplacement du cerveau du patient et la génération (500), à partir du déplacement du cerveau estimé, d'au moins une image du cerveau du patient dans laquelle le déplacement du cerveau est compensé.

Description

Procédé de traitement d'image pour estimer une déformation d'un cerveau d'un patient
DOMAINE DE L'INVENTION La présente invention concerne le domaine de la neurochirurgie assistée par ordinateur, et plus particulièrement un système et une méthode permettant la localisation précise d'une cible.
PRESENTATION DE L'ART ANTERIEUR La localisation précise d'une cible est essentielle pour réduire le taux de morbidité au cours d'une intervention chirurgicale d'ablation d'une tumeur cérébrale.
Au cours d'une intervention chirurgicale sur le cerveau, une déformation des tissus cérébraux (ou « brain shift » selon la terminologie anglo-saxonne) peut se produire après qu'une ouverture ait été réalisée dans le crâne du patient. Cette déformation peut être liée à des phénomènes physiques (gravité, perte du liquide céphalorachidien, action du neurochirurgien, etc.) ou à des phénomènes physiologiques (tuméfaction due aux médicaments osmotiques, anesthésiques, etc.), dont certains ne sont pas connus à ce jour.
Afin de compenser cette déformation du cerveau du patient et de retrouver les données préopératoires du patient pendant la procédure chirurgicale, des systèmes de neuro-navigation ont été développés.
Ces systèmes de neuro-navigation permettent d'afficher au chirurgien la position des objets chirurgicaux par rapport à la position des caractéristiques anatomiques cérébrales du patient. Dans ce type de systèmes la compensation de la déformation du cerveau du patient est en général réalisée sur le principe suivant.
Une première image de la surface corticale du patient est acquise avant l'opération chirurgicale. Pendant l'opération chirurgicale, une seconde image de la surface corticale du patient est acquise. Des moyens de traitement du système comparent ces deux images et estiment un champ volumique de déformation qui permet d'aligner les points dans les deux images afin de compenser la déformation du cerveau du patient. De tels systèmes sont notamment décrits dans les documents US 7 072 705 et US
2007/0021669.
Toutefois, un inconvénient de cette technique de compensation de la déformation du cerveau du patient à partir de points de la surface corticale est qu'elle laisse une zone importante du cerveau (sous la surface corticale) dans laquelle la déformation réel du cerveau est inconnue.
Un but de la présente invention est de proposer un système et un procédé de compensation de la déformation du cerveau du patient permettant de pallier au moins l'un des inconvénients précités.
RESUME DE L'INVENTION
A cet effet l'invention prévoit un procédé de navigation passive pour estimer une déformation d'un cerveau d'un patient, le procédé comprenant : - le traitement d'une image tridimensionnelle du cerveau du patient, acquise avant une intervention chirurgicale, pour obtenir une arborescence artérielle cérébrale de référence du patient,
- le traitement d'images bidimensionnelles du cerveau du patient, acquises pendant l'opération, pour reconstituer au moins partiellement une arborescence artérielle cérébrale courante du patient,
- la détermination, à partir de la mise en correspondance des arborescences artérielles cérébrales de référence et courante, d'un champ de déformation de l'arbre vasculaire représentant le déplacement de l'arbre vasculaire courant par rapport à l'arbre vasculaire de référence, - l'application du champ de déformation de l'arbre vasculaire déterminé à un modèle biomécanique du cerveau du patient pour estimer le déplacement du cerveau du patient .
- la génération, à partir du déplacement du cerveau estimé, d'au moins une image du cerveau du patient dans laquelle le déplacement du cerveau est compensé. L'arbre vasculaire initial est recalé élastiquement vers sa configuration déformée, segmentée dans les images US Doppler focalisées, et le champ de déplacements résultant est ensuite étendu à l'ensemble du volume de l'organe par le biais d'un modèle biomécanique spécifique au patient.
La déformation globale ainsi calculée permet de mettre à jour les données préopératoires et de rectifier le planning chirurgical initial.
L'invention fournit une réponse rapide, robuste et précise au chirurgien et lui donne la possibilité de valider la pertinence de la déformation calculée avant toute modification du planning.
Des aspects préférés mais non limitatifs de la présente invention sont les suivants : le procédé peut également comprendre le calcul des disparités entre des points de l'image du cerveau générée et une ou plusieurs images bidimensionnelles de contrôle acquise pendant l'opération ;
- l'image tridimensionnelle acquise avant l'opération peut être une angiographie par résonance magnétique, les images bidimensionnelles acquises pendant l'opération peuvent être des images ultrasonores en mode Doppler,
- le traitement de l'image tridimensionnelle pour obtenir l'arborescence artérielle cérébrale de référence comprend une étape d'orientation de l'image tridimensionnelle par rapport à la tête du patient après l'installation de celui-ci sur une table d'opération,
- l'étape d'orientation de l'image tridimensionnelle comprend : o la localisation d'un jeu de points de calibrage sur la tête du patient après l'installation de celui-ci sur la table d'opération, o la détection de leurs homologues radiométriques dans l'image tridimensionnelle, et o la détermination de transformations rigides permettant de passer des positions des points du jeu de points de calibrage aux positions de leurs homologues dans l'image tridimensionnelle, - les points de calibrage peuvent comprendre des points de la surface de la tête du patient après l'installation de celui-ci sur la table d'opération,
- le traitement des images bidimensionnelles pour reconstituer au moins partiellement une arborescence artérielle cérébrale courante du patient comprend une segmentation des images bidimensionnelles pour obtenir la position médiane des vaisseaux et/ou artères du cerveau de sorte à reconstituer au moins partiellement l'arborescence artérielle cérébrale courante du patient, l'étape de détermination d'un champ de déformation comprend la détermination de transformations non-linéaires, dites transformations élastiques, permettant le passage des positions des points de l'arborescence artérielle cérébrale de référence aux positions des points correspondants de l'arborescence artérielle courante.
L'invention concerne également un système pour la détermination d'une déformation du cerveau d'un patient, le système comprenant des moyens pour la mise en œuvre du procédé décrit ci-dessus. BREVE DESCRIPTION DES DESSINS
D'autres avantages et caractéristiques ressortiront mieux de la description qui va suivre de plusieurs variantes d'exécution, données à titre d'exemples non limitatifs, à partir des dessins annexés sur lesquels - la figure 1 illustre un mode de réalisation du procédé selon l'invention, la figure 2 illustre un mode de réalisation du système pour la mise en œuvre du procédé selon l'invention, la figure 3 illustre une déformation d'un carré unitaire avec un repliement d'espace,
- la figure 4 illustre en deux dimensions une approche itérative multi-échelle décrite en relation avec l'étape de recalage de tissu organique,
- la figure 5 illustre des cellules insérées en marge d'une discrétisation pour deux niveaux de raffinement successifs, la figure 6 illustre une fonction de forme w, la figure 7 illustre un exemple d'approximation d'une courbe d'énergie de recalage par une parabole,
- la figure 8 illustre trois étapes de l'approximation de l'inverse de la fonction de recalage,
- la figure 9 illustre des étapes de filtrage direct et indirect.
DESCRIPTION DETAILLEE DE L'INVENTION
On va maintenant décrire plus en détail un mode de réalisation du procédé selon l'invention.
Principe général de l'invention Comme expliqué précédemment, la localisation précise de la cible chirurgicale est essentielle pour réduire la morbidité au cours de l'exérèse chirurgicale d'une tumeur cérébrale. Lorsque les dimensions de la craniotomie sont importantes, une déformation des tissus mous du cerveau peut survenir au cours de l'intervention.
Du fait de cette déformation du cerveau, les données préopératoires ne correspondent plus à la réalité, et la neuronavigation s'en trouve fortement compromise.
La présente invention permet de prendre en compte cette déformation et de calculer une position rectifiée des structures anatomiques du cerveau afin de localiser des ancillaires.
Avant une intervention chirurgicale sur le cerveau d'un patient, une angiographie par résonance magnétique (ARM) du cerveau du patient peut être acquise. Au cours de l'intervention, suite à une déformation importante du cerveau, le chirurgien effectue un balayage échographique Doppler de la région d'intérêt.
L'invention propose une technique permettant la mise à jour des données préopératoire en fonction de la déformation du cerveau au cours de l'intervention. Pour ce faire, l'invention propose de suivre au cours de l'intervention la déformation de la structure anatomique particulière qu'est l'arbre vasculaire cérébral.
Du fait de sa répartition dans le volume du parenchyme, les déplacements de l'arbre vasculaire sont représentatifs des déformations globales des tissus.
En effet, les artères et artérioles irriguant les tissus nerveux sont présentes à la fois dans les parties profondes de l'encéphale, mais aussi en surface, et plus particulièrement au voisinage de la tumeur dont le développement passe par une angiogenèse et donc la formation de nouveaux vaisseaux.
De la localisation de l'arbre vasculaire par rapport à sa position initiale en début d'intervention il est possible de déduire les déplacements de ses vaisseaux, provoqués par la déformation du cerveau.
Cette information, loin d'être dense dans le volume global, reste localisée à certaines régions.
L'invention propose alors d'utiliser une modélisation biomécanique du cerveau du patient afin d'extrapoler une déformation globale à partir du déplacement de l'arbre vasculaire cérébral.
Le champ de déplacements global ainsi obtenu est ensuite répercuté à l'ensemble des données préopératoires afin de les mettre à jour.
Description d'un mode de réalisation du procédé selon l'invention En référence à la figure 1, on a illustré différentes étapes du procédé selon l'invention.
Le procédé de traitement d'image pour estimer un déplacement d'un cerveau d'un patient comprend les étapes suivantes :
- le traitement 100 d'une image tridimensionnelle du cerveau du patient, acquise avant une intervention chirurgicale, pour obtenir une arborescence artérielle cérébrale de référence du patient, le traitement 200 d'images bidimensionnelles du cerveau du patient, acquises pendant l'opération, pour reconstituer au moins partiellement une arborescence artérielle cérébrale courante du patient, - la détermination 300, à partir de la mise en correspondance des arborescences artérielles cérébrales de référence et courante, d'un champ de déformation de l'arbre vasculaire représentant le déplacement de l'arbre vasculaire courant par rapport à l'arbre vasculaire de référence, l'application 400 du champ de déformation de l'arbre vasculaire déterminé à un modèle biomécanique du cerveau du patient pour estimer le déplacement global du cerveau du patient ; ce modèle sert d'interpolateur mécanique afin de calculer au cours de l'intervention l'ensemble des déplacements volumiques à l'intérieur du cerveau
- la génération 500, à partir du déplacement du cerveau estimé, d'au moins une image du cerveau du patient dans laquelle le déplacement du cerveau est compensé.
Détermination de l'arbre vasculaire cérébral de référence
Comme décrit précédemment, le procédé comprend une étape 100 consistant à traiter une image du cerveau du patient acquise avant l'intervention chirurgicale. Ce traitement a pour but de localiser l'arbre vasculaire cérébral dans l'image acquise avant l'intervention chirurgicale.
Cet arbre vasculaire cérébral dans l'image acquise avant l'intervention chirurgicale correspond à l'arbre vasculaire avant la déformation du cerveau du patient, et est appelé dans la suite « arbre vasculaire cérébral de référence ». Préférentiellement, l'image acquise avant l'intervention chirurgicale est une image tridimensionnelle.
Préférentiellement encore, cette image tridimensionnelle est une angiographie par résonnance magnétique.
L'angiographie par résonance magnétique permet de révéler des contrastes importants entre les tissus. L'analyse de l'angiographie par résonance magnétique permet de distinguer avec précision la nature des tissus observés.
En particulier, les données qu'elle contient permettent de localiser l'arbre vasculaire cérébral de référence. Par ailleurs, la technique permettant l'acquisition d'une angiographie par résonance magnétique présente l'avantage d'être non invasive et donc sans dangers pour le patient, contrairement par exemple à des techniques émettant un rayonnement ionisant comme la tomographie.
Une fois le patient installé sur la table opératoire, le traitement de l'image tridimensionnelle pour obtenir l'arborescence artérielle cérébrale de référence comprend une étape d'orientation de l'image tridimensionnelle par rapport à la tête du patient. En d'autres termes, la position spatiale de l'arbre vasculaire de référence est déterminée par rapport à un référentiel de l'espace physique du patient. Plus précisément, l'arbre vasculaire cérébral de référence est transposé dans le référentiel physique du patient, à l'aide d'une technique dite de recalage rigide qui s'appuie sur la mise en correspondance de points de calibrages avec leurs équivalents segmentés dans l'angiographie par résonance magnétique acquise avant l'opération. On entend, dans le cadre de la présente invention, par « recalage rigide », une transformation géométrique composée d'une rotation et d'une translation.
Un calcul permet alors de trouver la transformation entre les deux systèmes de coordonnées.
Les points de calibrage utilisés peuvent être des points ou des surfaces anatomiques remarquables comme la pointe du nez, les arcades sourcilières, la surface des joues ou du front.
Les points de calibrage peuvent aussi être définis par une pluralité de pastilles adhésives (au moins trois et préférentiellement dix) collées sur la peau du patient avant l'acquisition de l'angiographie par résonance magnétique. Ainsi, ces pastilles sont acquises dans l'angiographie en même temps que le cerveau du patient.
Les homologues radiométriques des points de calibrage sont identifiables dans l'angiographie par résonance magnétique. Une fois le patient installé sur la table d'opération, le chirurgien indique à un système de localisation la position de chaque point de calibrage au moyen d'un palpeur. Le système de localisation peut être un système optique de localisation (comprenant par exemple une caméra stéréoscopique), ou un système magnétique de localisation, ou tout autre système de localisation connu de l'homme du métier.
Connaissant la position de l'arbre vasculaire cérébral de référence par rapport aux points de calibrage, et connaissant la position des points de calibrage une fois le patient installé sur la table d'opération, le système est apte à calculer la position spatiale de l'arbre vasculaire de référence par rapport à la position spatiale du patient.
Le recalage entre les images et le patient, réalisé en début d'intervention, permet par la suite de localiser dans le champ opératoire tous les éléments identifiés dans l'image par résonance magnétique et, réciproquement, de transposer la position des ancillaires dans le volume de données afin de suivre et contrôler la bonne réalisation de l'intervention.
Détermination de l'arbre vasculaire cérébral courant
Le procédé comprend également une étape 200 consistant à traiter une pluralité d'images bidimensionnelles du cerveau du patient, acquises pendant l'opération, pour reconstituer au moins partiellement une arborescence artérielle cérébrale du patient dite
« arborescence artérielle cérébrale courante ».
Durant l'opération, suite à une déformation importante du cerveau du patient, le chirurgien effectue une acquisition d'images bidimensionnelles. On entend, dans le cadre de la présente invention, par « pluralité d'images bidimensionnelles », des images 2D dans un volume (3D) déterminé dont les plans d'acquisition sont sécants ou parallèles. Ces images 2D peuvent être obtenues à partir de tout dispositif médical d'acquisition d'images 2D ou 3D telles que des images échographiques 3D par exemple. Préférentiellement, ces images bidimensionnelles sont des images échographiques acquises en mode Doppler, à l'aide d'une sonde localisée par le système de localisation.
Les images échographiques en mode Doppler permettent de visualiser les flux sanguins. Par ailleurs, la technique d'échographie Doppler présente l'avantage d'être bénigne pour le patient. L'échographie Doppler per-opératoire est donc un bon outil pour l'exploration des déplacements du cerveau au cours de l'intervention. Ce mode de visualisation repose sur l'effet Doppler et permet de localiser, en analysant les variations de fréquence des ultrasons émis, des déplacements au sein des tissus, tels que les écoulements de sang. La segmentation de ce type d'images est rendu plus aisée du fait du codage en couleur utilisé pour représenter les vélocités des écoulements. Cette modalité permet donc de visualiser dans le volume de l'encéphale les flux sanguins qui s'effectuent le long de l'arbre vasculaire et d'identifier ainsi la position des vaisseaux qui irriguent le parenchyme.
L'acquisition des images bidimensionnelles est réalisée par exemple comme suit. La sonde localisée est mise en contact avec le cerveau du patient à travers la craniotomie effectuée par le chirurgien. Un mouvement de rotation est manuellement imprimé à la sonde localisée de sorte que le plan échographique balaye la partie du cerveau la plus étendue possible autour de la tumeur.
La sonde localisée comprend un marqueur permettant au système de localisation de définir sa position et son orientation par rapport au référentiel physique du patient, et donc de repositionner spatialement les images bidimensionnelles dans le référentiel physique du patient.
Après l'acquisition d'une série d'images bidimensionnelles, celles-ci sont traitées pour localiser l'arbre vasculaire cérébral courant du patient. Plus précisément, les centres des vaisseaux contenus dans chaque image bidimensionnelle sont sélectionnés de sorte à obtenir un nuage de points. Ce nuage de points représente au moins partiellement l'arbre vasculaire cérébral courant du patient.
Recalaae des arbres vasculaires de référence et courant Dans une autre étape 300 du procédé, l'arbre vasculaire de référence et courant sont mis en correspondance pour déterminer un champ de déformation représentant le déplacement de l'arbre vasculaire courant par rapport à l'arbre vasculaire de référence.
Plus spécifiquement, l'arbre vasculaire de référence est recalé par rapport à l'arbre vasculaire courant. On entend par « recalage élastique », le fait de déterminer, à partir de jeux de données qui représentent deux images numériques d'une même zone, une transformation permettant de passer de l'image de référence à l'image courante. En d'autres termes, il s'agit de superposer "au mieux" des structures contenues dans deux images comparables. Les données issues de l'image tridimensionnelle et des images bidimensionnelles sont de natures distinctes.
D'un côté, l'arbre vasculaire de référence obtenu à partir des données préopératoires est presque continu, l'ensemble de la tête du patient est contenu dans le volume de données, et la segmentation réalisée n'engendre pour ainsi dire pas d'artefacts dans la reconstruction de l'arbre vasculaire cérébral.
D'un autre coté, l'arbre vasculaire courant obtenu à partir des images bidimensionnelles per-opératoires est composé d'un nuage de points, ces images bidimensionnelles ne couvrent qu'une fraction du volume du cerveau. Selon l'habileté de l'opérateur, les images bidimensionnelles peuvent être plus ou moins régulièrement réparties dans l'espace, certaines régions comportant des acquisitions multiples tandis que d'autres ne comportent aucune acquisition.
La tâche qui incombe à l'étape 300 de recalage des arbres vasculaires cérébraux est donc d'estimer une transformation élastique entre :
- la configuration courante, per-opératoire, de l'arbre vasculaire définie par un nuage de points épars, redondants, incomplet et bruités, et
- la configuration de référence définie par les données complètes et propres obtenues à partir de l'imagerie préopératoire.
Le champ de déformation est alors calculé à partir de la correspondance établie entre l'arbre vasculaire de référence et l'arbre vasculaire courant. Chaque vecteur du champ de déformation a pour origine un point de l'arbre vasculaire de référence et désigne un point de l'arbre vasculaire courant, déformé. Le recalage élastique est guidé par la minimisation d'une énergie de recalage. Cette énergie de recalage quantifie la similarité entre l'arbre vasculaire courant et l'arbre vasculaire de référence.
Le recalage optimal est atteint pour une énergie de recalage nulle. L'objectif de l'étape de recalage élastique est de minimiser itérativement la valeur de cette énergie de recalage en faisant évoluer les points de l'arbre vasculaire courant vers les points de l'arbre vasculaire de référence.
La rapidité du calcul est obtenue grâce :
- au calcul préalable d'une carte de distances - entre chaque point de l'arbre vasculaire courant et son homologue au sens radiométrique (i.e. le point qui désigne la même structure) dans l'arbre vasculaire de référence - sur laquelle s'appuie l'évaluation de l'énergie de recalage et
- la détermination de la direction de plus forte descente, calculée à chaque itération du processus de minimisation. Etant donnée la nature éparse et irrégulière des données, la recherche d'une fonction de recalage élastique doit être contrainte afin de ne pas aboutir à un recalage aberrant. La fonction de recalage élastique est une transformation non-linéaire permettant de passer de l'arbre vasculaire courant à l'arbre vasculaire de référence.
Plus précisément, des conditions liées à la nature des données que l'on cherche à recaler sont déterminées pour garantir le respect de critères physiques essentiels.
S'agissant de données organiques, une première condition que doit satisfaire la fonction de recalage élastique concerne le non-repliement de l'espace. En effet, la déformation du cerveau du patient ne peut induire le repliement du cerveau du patient sur lui-même. Une deuxième condition que doit satisfaire la fonction de recalage élastique concerne la continuité de la propagation de la déformation. En effet, les arbres vasculaires de référence et courant sont des observations du phénomène physique de déformation du cerveau. De façon similaire au processus de déformation du cerveau du patient, sans déchirement ou résection de tissu, la fonction de recalage doit être continue et même continûment différentiable i.e. différentiable en tout point du domaine et ne présentant pas de discontinuité de différentielle.
Une troisième condition que doit satisfaire la fonction de recalage est que deux points distincts de l'arbre vasculaire courant ne peuvent avoir la même image dans l'arbre vasculaire de référence, et vice versa. La fonction de recalage élastique telle qu'elle est proposée dans l'invention présente l'avantage d'être inversible avec le degré de précision souhaité. Une fois cette fonction de recalage calculée, un post-traitement du résultat produit par le recalage permet d'éliminer les aberrations dues à la présence de bruit dans les données des images bi- ou tridimensionnelles per-opératoires.
L'ensemble de déplacements de l'arbre vasculaire finalement retenu forme un champ de déformation représentant le déplacement de l'arbre vasculaire courant par rapport à l'arbre vasculaire de référence.
Le lecteur appréciera que l'étape de recalage peut être appliquée au recalage de structures autre qu'un arbre vasculaire cérébral. Notamment, le recalage élastique présenté ci-dessus peut s'appliquer au recalage de tout type de tissu organique.
Extrapolation du champ de déplacements de l'arbre vasculaire cérébral à l'ensemble du volume de l'organe par le biais d'un modèle biomécanique
Dans une autre étape 500 du procédé, le champ de déformation de l'arbre vasculaire déterminé est appliqué à un modèle biomécanique des tissus mous du cerveau du patient pour estimer le déplacement global du cerveau du patient.
Ce modèle mathématique formalise une connaissance a priori du comportement de l'organe et permet d'en inférer le comportement global sous les contraintes de déformation locales déduites de l'observation per-opératoire, appelées conditions aux limites. Le modèle biomécanique (aussi connu sous le nom de modèle déformable) joue un double rôle.
En premier lieu, le modèle biomécanique permet d'extrapoler les déplacements de l'arbre vasculaire à l'ensemble du volume du cerveau afin de simuler l'effet de la déformation du cerveau. En deuxième lieu, le modèle biomécanique permet de réguler les données. En effet, l'estimation du champ de déplacements de l'arbre vasculaire est entachée d'erreurs dues au mode d'imagerie utilisé, aux limitations des algorithmes de segmentation ainsi qu'à celles du calcul du recalage élastique effectué pour mettre en correspondance l'arbre vasculaire courant et l'arbre vasculaire de référence. Le modèle biomécanique permet d'harmoniser le champ de déplacements en supprimant les erreurs dans les données d'entrée du modèle biomécanique. Ceci permet de rendre le procédé résistant aux artefacts présents dans les données bruitées recueillies au cours de l'intervention (i.e. les images échographiques bi- ou tridimensionnelles).
L'originalité du procédé concerne ainsi la complémentarité entre un modèle biomécanique approximant les comportements des tissus cérébraux du patient et un suivi de structure anatomique potentiellement bruité mais représentatif du déplacement du cerveau grâce à l'arbre vasculaire cérébral.
A l'issu de l'étape 500 du procédé, le déplacement global du cerveau du patient est estimé. Une dernière étape 600 du procédé concerne la génération, à partir du déplacement global du cerveau estimé, d'au moins une image du cerveau du patient dans laquelle le déplacement du cerveau est compensé. Cette image est par exemple l'image tridimensionnelle (ou des images bidimensionnelles en coupe de cette image tridimensionnelle) acquise avant l'opération à laquelle est appliquée la déformation globale estimée du cerveau du patient de sorte à mettre à jour les données préopératoire en fonction de la déformation globale estimée du cerveau du patient.
Vérification de la cohérence de la déformation calculée
Avantageusement, le procédé peut comprendre une étape de calcul des disparités entre des points de l'image du cerveau générée et une ou plusieurs images bidimensionnelles de contrôle acquises pendant l'opération.
Ceci permet à l'utilisateur de vérifier la cohérence et la précision du calcul réalisé par le procédé décrit ci-dessus. Pour ce faire, l'utilisateur peut par exemple acquérir une ou plusieurs images bidimensionnelles de contrôle du cerveau du patient pour calculer l'erreur de position entre les points des images bidimensionnelles de contrôle et leurs homologues dans l'image générée à partir de la déformation globale estimée du cerveau du patient.
Par exemple, l'utilisateur peut scruter les tissus cérébraux dans la région d'intérêt tout en projetant dans les images bidimensionnelles de contrôle la position de l'arbre vasculaire déformé par le modèle. La superposition des deux données i.e. Doppler en temps réel obtenu par balayage de l'organe, d'une part et signal Doppler simulé, calculé à partir de l'incidence des coupes échographiques localisées par rapport à l'arbre vasculaire déformé par le système, d'autre part, permet à l'utilisateur d'évaluer la coïncidence entre le modèle et la réalité du patient. L'appréciation de la qualité de la déformation permet à l'utilisateur de prendre la décision de suivre ou non les indications du système sur la base des données préopératoires mises à jour.
Le déplacement du cerveau étant un processus évolutif, cette validation est réalisée aussitôt après l'acquisition des données sur lesquelles se base l'algorithme de déformation afin d'en garantir la pertinence. Le procédé décrit ci-dessus est parfaitement adapté à un suivi continu de la déformation des tissus mous du cerveau, tant par ses temps de réponse qui sont de l'ordre de quelques secondes, que par la répétabilité de la procédure d'acquisition de données qui ne nécessite pas une interruption majeure de l'intervention ou la mise en place d'équipements encombrants.
Théorie associée au recalape élastique
On va maintenant décrire plus en détail la théorie associée à l'étape de recalage élastique. Les explications qui vont suivre sont faites en relation avec le recalage élastique de l'arbre vasculaire cérébral, mais il est bien entendu que ce recalage peut s'appliquer à d'autres types de tissus organiques. Un recalage élastique est une application
Figure imgf000015_0002
qui minimise une certaine erreur de recalage, ou plutôt, une énergie de recalage définie entre un ensemble de points sources, S et un autre ensemble de points destination, D. L'énergie du recalage R est la mesure d'une « similarité » entre l'ensemble de points transformés R(S) et l'ensemble cible D. Lorsque les deux ensembles se confondent alors cette énergie devient nulle et le recalage est accompli. Si l'énergie de recalage et l'espace de fonctions dans lequel évolue la transformation R ne sont pas définis correctement, alors le problème consistant à trouver un R optimal est mal posé, i.e. l'existence et l'unicité de R ne sont pas garanties.
Nous présentons dans cette partie, le cadre de calcul de la déformation élastique qui permet d'établir de façon robuste une correspondance élastique entre les données considérées. Les spécificités de cette approche résident dans les propriétés mathématiques de la transformation R, démontrées plus loin, qui en garantissent le réalisme physique.
Energie de recalaαe
L'énergie de recalage peut être librement définie afin de refléter la nature du problème de recalage traité. Par exemple, dans le cas de recalage entre deux images, cette mesure peut être basée sur les différences entre les valeurs des pixels, ou voxels, des images source et destination. Dans le cas présent, la similarité entre les jeux de données S et D est définie par la proximité spatiale des deux nuages de points. Sans autre connaissance de la nature des données manipulées, il est possible, par exemple, définir une énergie de recalage Ereg comme la moyenne des distances euclidiennes minimales entre R(S) et D.
Figure imgf000015_0001
Où : - R représente la fonction de recalage élastique,
- Card(S) est le cardinal de l'ensemble source S,
- Card(D) est le cardinal de l'ensemble destination D,
- d(R(s),D) représente la distance entre l'ensemble D et le transformé par R d'un point s de l'ensemble S,
- d(d, R(S) représente la distance entre un point d de l'ensemble D et la transformée par R de l'ensemble S.
Dans ce cas un recalage optimal est réalisé lorsque la valeur de Erβg s'annule. L'énergie symétrique présentée ci-dessus suppose que chaque structure présente dans S peut trouver, au moyen d'une transformation R appropriée, une correspondance parmi D, et réciproquement. Les désignations « S » et « D » des deux jeux de données sont interchangeables.
Dans le cas présent, toutefois, les structures présentes dans S et D peuvent ne pas être équivalentes et l'hypothèse de symétrie ne tient pas puisque les structures dans S ne représentent qu'un sous-ensemble de D. Une définition alternative doit être trouvée pour en tenir compte. Dorénavant, on désignera par « S » l'ensemble de points définissant la configuration courante, per-opératoire, de l'arbre vasculaire, et D ceux qui en définissent la configuration de référence, préopératoire. La fonction de recalage élastique R est donc celle qui amène les points de S sur D et une définition asymétrique de l'énergie de ce recalage peut être formulée en tronquant l'expression symétrique ci- dessus :
Figure imgf000016_0001
La signification des appellations « source » et « destination » peut être mieux saisie en voyant le recalage élastique comme un processus itératif où les points source S sont déplacés vers l'ensemble de destination selon une trajectoire optimale sur une surface d'énergie potentielle dont le relief est défini par la configuration spatiale des points de D. L'évaluation d'une configuration candidate R(S) est faite sur l'ensemble contenant moins de données tandis que la fonction d'énergie est plus strictement définie par D, lequel contient l'information de référence sur l'ensemble de la structure de l'arbre vasculaire. Des solutions alternatives au problème de l'asymétrie figurent dans la littérature. La définition ci-dessus de Ereg peut aboutir à des situations où un groupe de points sources est recalé vers une partie confinée de D, en d'autres termes, cette énergie ne favorise pas la répartition des points sources sur la structure destination. Pour éviter cela, il est possible d'en étendre l'expression en rajoutant un terme qui quantifie cette répartition mais qui opère uniquement dans un voisinage restreint de points de D qui contiennent des points de R(S).
Finalement, étant donné que le cardinal de l'ensemble S n'évolue pas au cours du recalage, il est encore possible de simplifier l'expression précédente de l'énergie Erβg en éliminant la constante multiplicative qui la précède, ce qui nous amène à la définition suivante de Ereg :
Figure imgf000017_0001
Espace des fonctions de déformation
Etant donné la nature éparse et irrégulière des données, la recherche d'une fonction de recalage élastique doit être contrainte afin de ne pas aboutir à un recalage aberrant. Cet écueil est principalement dû aux nombreux minima locaux qui jalonnent le produit de la fonction d'énergie pour une configuration D donnée. Les deux ensembles S et D sont des observations du phénomène physique étudié i.e. la déformation du cerveau.
On suppose donc que de façon similaire au processus du « brain-shift », sans déchirement ou résection de tissu, la transformation non-linéaire R cherchée est continue et même continuement différentiable i.e. différentiable en tout point du domaine et ne présentant pas de discontinuité de différentielle. Ceci définit une première restriction sur l'espace de recherche des fonctions de recalage
Figure imgf000017_0003
Une autre contrainte forte sur la nature de la transformation est qu'elle ne doit pas replier l'espace. La contrainte de non-repliement peut s'exprimer mathématiquement comme suit. Considérons un trièdre orienté positivement constitué de trois vecteurs infinitésimaux, dXi, dX2 and dX3, placés en un point p de l'espace. Soit dV le volume, positif, défini par ces trois vecteurs et dont la valeur est donnée par le déterminant de la matrice formée par les trois vecteurs.
Figure imgf000017_0002
Après application de la transformation R au point p les trois vecteurs se transforment en dx1 , dx2 et dx3. La nature infinitésimale de ces vecteurs nous autorise l'approximation suivante :
Figure imgf000018_0001
Le volume défini par le trièdre transformé, noté dv est alors obtenu de façon similaire.
\ JR(P) I W
Figure imgf000018_0002
La contrainte de non-repliement s'énonce alors tout simplement : dv > 0. Ce qui signifie que le volume infinitésimal déformé n'a pas été « retourné » (dv < 0) par R ou totalement écrasé (dv = 0). Ainsi l'application R ne replie pas localement l'espace si et seulement si son jacobien en p est strictement positif. Dans notre cas le recalage élastique ne doit replier l'espace en aucun point, ce qui nous amène à l'expression de la contrainte de non-repliement :
Figure imgf000018_0003
Si en un point donné, le jacobien est plus grand que 1 cela signifie que le recalage élastique étire localement le volume, si au contraire il est plus petit que 1 , alors le volume est contracté. Enfin le produit de la matrice jacobienne de R calculée en un point par un vecteur unitaire nous donne le taux d'allongement ou de contraction de l'espace dans une direction donnée au voisinage de ce point. La figure 3 montre une déformation d'un carré unitaire avec un repliement d'espace.
Finalement pour les mêmes motivations physiques, l'application R doit être injective. Ce qui se traduit par le fait que deux points distincts ne peuvent pas avoir la même image par R et, physiquement, R ne doit pas amener deux atomes en un même point de l'espace. Si R est une injection, alors pour chaque « point atteint
Figure imgf000018_0006
un prédécesseur peut être calculé.
Dans le cas présent où S représente une fraction de la structure présente dans D, la transformation cherchée in fine est celle qui amène D sur S, soit R"1, afin de définir le champ de déplacements de l'arbre vascuiaire de la configuration de référence D vers la configuration courante S, comme décrit ci-dessus. Pour que l'expression R"1 ait un sens, il faut que
Figure imgf000018_0004
Afin de simplifier ces considérations, on considère les fonctions R surjectives, Le. telles que
Figure imgf000018_0005
Il est donc possible de conclure que la fonction de recalage R cherchée est bijective. La contrainte de non-repliement garantit que le jacobien est strictement positif en tout point de m ceci est cependant insuffisant pour garantir la bijectivité globale.
L'espace de fonctions doit alors être convenablement défini afin d'assurer la bijectivité de la fonction R. Les sections suivantes décrivent la stratégie de recalage adoptée afin d'aboutir à une minimisation de l'énergie de recalage tout en respectant les contraintes énoncées plus haut de différentiabilité, non-repliement et bijection.
Calcul du recalaαe
L'approche présentée ici est inspirée de la modélisation mécanique des solides basée sur une discrétisation de type éléments finis. Afin d'approximer le comportement mécanique des données non-structurées, on suppose que le nuage de points soumis à la déformation est englobé dans un solide virtuel élastique. Puisque le déplacement des points traduit les propriétés physiques d'une structure sous-jacente dont la géométrie est inconnue, on attribue arbitrairement à ce solide virtuel la forme d'une boîte englobante des points. On définit alors l'ensemble des déformations admissibles comme les combinaisons successives de déformations élémentaires du solide englobant, dont les effets se propagent à travers son volume et modifient la position des points internes. Les déformations élémentaires mentionnées ici sont définies plus loin. La technique de recalage 3D décrite ici peut aisément être dérivée en une technique de recalage 2D en remplaçant la notion de « solide virtuel » par « plan virtuelle ». Les illustrations présentées ci-après se basent sur une implémentation en 2D du recalage pour en clarifier les concepts.
La fonction de recalage élastique est assemblée itérativement de manière à ce que l'énergie de recalage Erβg diminue de façon optimale à chaque pas. Soit S0 := S, le nuage de points source initial, et soit S1 le nuage de points source obtenu suite à l'itération i.
Le domaine du solide virtuel est discrétisé avec un certain niveau de raffinement en hexaèdres réguliers (ou carrés, en 2D) appelés « cellules ». Les déformations élémentaires du solide sont obtenues en déplaçant individuellement les nœuds de la grille ainsi formée.
L'union des cellules entourant un nœud n définit son « voisinage », noté V (n). Une fonction de recalage élémentaire r est définie en propageant le déplacement du nœud n à l'ensemble des points de sous les contraintes suivantes : La frontière du voisinage de n, i)V(n) est inchangée i.e. Vp E dV(n), r(p) = p ; Le voisinage du n reste globalement inchangé i.e. r(V(n)) = V(n) ;
Tous les points à l'extérieur de ve voisinage restent inchangés i.e. Vp φ V(n), r(p) = p.
Chaque nœud du solide discrétisé est considéré et le gradient de l'énergie de recalage, en tant que fonction du déplacement du nœud, est calculé. L'opposé de ce vecteur définit le déplacement préférentiel du nœud, c'est-à-dire celui qui, appliqué au nœud, entraîne la diminution la plus importante de l'énergie de recalage.
De tous les nœuds, celui qui permet la diminution la plus importante du gradient est choisi et son déplacement préférentiel lui est appliqué tandis que les autres nœuds demeurent immobiles. La déformation locale du solide est calculée en propageant le déplacement du nœud à travers son voisinage et la position des points source qu'il contient est modifiée en conséquence, donnant ainsi la nouvelle configuration du nuage de points Si+1.
Une fois la déformation élémentaire appliquée, le solide retourne dans sa configuration initiale, au repos, et les points source Si+1 sont distribués dans son volume. Une nouvelle itération peut alors être calculée en prenant la nouvelle configuration Si+1 comme point de départ. Les itérations s'arrêtent, enfin, lorsqu'aucune diminution significative de l'énergie ne peut être obtenue en déplaçant les nœuds de la grille du solide.
La boucle itérative décrite ci-dessus suppose un degré de raffinement du domaine
* - fixe. Afin de maintenir la consistance spatiale de S au cours de l'assemblage de R et en éviter le morcellement, une approche multi-échelle a été adoptée. L'assemblage itératif de R commence au niveau de raffinement le plus grossier et lorsque l'énergie de recalage ne peut être améliorée davantage à un niveau donné, les cellules sont subdivisées en 8 sous-cellules (ou 4, en 2D) à la manière d'un octree.
Les itérations de raffinement s'arrêtent lorsque le degré maximal de subdivision spécifié manuellement a été atteint. La dimension de la plus petite cellule atteinte au cours du processus de raffinement définit les dimensions de la plus petite structure présente dans S et qui ne sera pas recalée correctement si une structure équivalente est absente de D.
La figure 4 illustre en 2D l'approche itérative multi-échelle décrite plus haut. Les points source S sont représentés par les points référencés 1. L'ensemble D n'est pas représenté, pour plus de clarté. En (a) l'ensemble initial SO := S est placé au sein du solide virtuel (référencé 2) discrétisé par une cellule (niveau 1) et les gradients de l'énergie de recalage sont calculés en ses quatre nœuds 3. En (b) le meilleur déplacement nodal est appliqué au solide et les points source sont déplacés en S1. (référencé 4). Le solide retourne alors à sa configuration initiale (a). Après plusieurs itérations, si aucune amélioration significative de l'énergie ne peut être obtenue, le niveau de raffinement de la grille augmente (c). Les gradients de l'énergie sont calculés aux neuf nœuds 5 de la nouvelle grille. En (d) le meilleur déplacement est appliqué au solide et la nouvelle configuration de points source S2 (référencé 6) est calculée.
Grâce à la technique de « remise à zéro » de la déformation du solide après chaque itération, de grandes déformations de S peuvent être envisagées sans avoir à maintenir une grille irrégulière, soumise à une somme de déformations importantes ce qui aurait été le cas si le solide suivait strictement les déformations de S. La régularité de la grille avant chaque itération permet également un gain substantiel d'effort calculatoire, étant donné que l'interpolation des déplacements nodaux peut être calculée sur la base d'une cellule type cubique. Avant de détailler davantage le calcul de R, nous devons étendre le concept intuitif du solide virtuel de manière à ce que chaque déformation élémentaire soit définie comme une application ™ ~~ * ™" , afin de répondre à la restriction énoncée plus haut.
Considérons pour cela un point source placé à la frontière du solide virtuel , par exemple à proximité du bord supérieur du carré sur la figure 4. Si le nœud est déplacé comme cela est illustré en (d), ce point va tomber en dehors du solide dans sa configuration au repos (carré) et deviendra un « point mort », dont la position ne pourra être modifiée par les itérations suivantes, lesquelles n'ont d'effet que sur les cellules de Ω.
Pour pallier cet inconvénient, on ajoute autour du domaine, une marge de cellules de dimensions identiques à celles utilisées pour la discrétisation. Cet ajout est réalisé à chaque niveau de raffinement.
Les nœuds extérieurs de ces cellules sont fixes et assurent ainsi un rôle d'interface entre la zone déformable et le reste de ιx , qui demeure inchangé. La figure 5 montre les cellules supplémentaires 7 insérées en marge de la discrétisation du domaine pour deux niveaux de raffinement successifs.
Propriétés du recalage élastique
Soient G1, ...,GJ les niveaux de raffinement successifs, G1 étant le niveau le plus
grossier et GJ le plus fin. So l'ensemble des
Figure imgf000021_0001
déformations élémentaires calculées à chaque niveau de raffinement j (où Ij est le nombre total de déformations élémentaires au niveau j). L'expression de la fonction de recalage élastique R assemblée itérativement est :
Figure imgf000022_0001
II est clair que les propriétés de différentiabilité, non-repliement et bijection de R dépendent des propriétés des fonctions de recalage élémentaires ή. Si chaque ή est une τπ>3 j, ira 3 bijection alors la fonction de recalage R est également une bijection et son inverse est donnée par :
Figure imgf000022_0002
Si chaque rj est une application continuement différentiable
Figure imgf000022_0005
alors la
fonction de recalage R l'est également. Soit
Figure imgf000022_0006
un point quelconque et pk le point transformé après application de k déformations élémentaires. Soit K = 11 + 12 + ... +
IJ , alo
Figure imgf000022_0007
st l'ensemble de toutes les positions successives de pO et R(p) = pK. La différentielle de R en pO est donnée par la règle de composition suivante :
Figure imgf000022_0003
Finalement, si tous les Hj vérifient la condition de non-repliement, alors cela est vrai également de l'application R. Ceci est la conséquence de l'expression du jacobien de R :
Figure imgf000022_0004
Ainsi pour prouver que toute fonction R de la famille de fonctions de recalage élastique générée à l'aide de la procédure décrite plus haut vérifie les contraintes de différentiabilité, bijection et non-repliement, il suffit de prouver que ces propriétés sont bien vérifiées pour chacune des fonctions de recalage élémentaire ή qui composent R. Etude des fonctions de recalaae élémentaire
Expression générale de ή
A ce point, la famille de fonctions de recalage élémentaires doit être spécifiée. Dans la méthode des Eléments Finis, une valeur nodale (déplacement ou autre quantité) est interpolée à travers le volume discrétisé au moyen de fonctions de forme
Figure imgf000023_0003
. La valeur de la fonction de forme associée à un nœud et évaluée en un point du domaine donne la proportion de la valeur nodale transférée en ce point. Ainsi la valeur de w en son nœud est 1. Le même mécanisme pour propager le déplacement d'un nœud à travers son voisinage va être utilisé.
On considère la ième itération de recalage effectuée au niveau de raffinement j et mettant en œuvre la fonction de recalage élémentaire r'j. Les indices d'ordre et de niveau de raffinement i et j seront omis dans la suite de ce paragraphe. Soit u le déplacement et w le fonction de forme associés au nœud dont le déplacement est réalisé lors de cette itération. L'expression de r est donc :
Figure imgf000023_0002
Les conditions sur la zone d'influence de r se traduisent par des contraintes sur le support de w de la façon suivante :
Figure imgf000023_0004
_. . . . . . . x Remarquons que le support de w est
Figure imgf000023_0005
ui est un compact
Figure imgf000023_0006
Différentiabilité La différentiabilité de r découle de celle de
Figure imgf000023_0007
alors
Figure imgf000023_0008
^ ' et la matrice jacobienne de r s'écrit :
Figure imgf000023_0001
où Id3><3 est la matrice identité 3 * 3.
Donc pour que r so
Figure imgf000023_0009
suffit que w le soit aussi. Non-repliement
En réécrivant la contrainte de non-repliement à l'aide de l'expression de la matrice jacobienne de r et en développant le déterminant de la matrice, une nouvelle expression de la contrainte de non-repliement est obtenue :
Figure imgf000024_0003
Puisque P
Figure imgf000024_0005
est une application continue à support compact, soit
Figure imgf000024_0004
te point où
Figure imgf000024_0001
Alors l'équation précédente est équivalente à :
Figure imgf000024_0006
Le non-repliement de r peut donc être garanti en limitant l'amplitude des déplacements nodaux. La norme maximale du déplacement nodal est donnée par l'inverse de la norme maximale du gradient de la fonction de forme considérée w.
Bijection Soit V := V (n) l'union fermée des cellules entourant le nœud n. Dans cette partie nous allons montrer que si une application élémentaire r ne replie pas l'espace, alors elle est également bijective de V dans lui-même mais également de R3 ! R3, puisque les points à l'extérieur de V ne sont pas affectés par r.
Implémentation
Choix d'une expression de w Afin de simplifier le calcul du recalage élastique, une expression polynomiale de w est choisie. Soit TT un polynôme de degré 3, défini par :
Figure imgf000024_0002
Son expression es ,t
Figure imgf000024_0007
Une fonction de forme w définie sur la cellule « canonique » i"1 1J et prenant la valeur 1 au noeud peut ainsi être obtenue comme :
Figure imgf000024_0008
Figure imgf000025_0001
La figure 6a montre la fonction de forme restreinte à une dimension, i.e. w = 7S sur l'intervalle [0, 1]. La fonction de forme w est définie sur l'union des cellules voisines autour du nœud n par un changement de variable et une mise à l'échelle afin d'adapter son domaine de définition canonique
Figure imgf000025_0006
L > J aux dimensions des cellules dans la grille de discrétisation du solide virtuel. La figure 6b montre les changements de variable nécessaires pour assembler à partir de la fonction de forme canonique, la fonction de forme w sur un voisinage de 2*2 cellules du plan. La figure 6c montre le tracé de la fonction de forme ainsi obtenue
Figure imgf000025_0005
Contrainte de non-repliement et de bijection
L'expression de w étant choisie, il est possible de calculer l'amplitude du déplacement nodal maximal autorisé pour une fonction de recaiage élémentaire r. Pour définir cette amplitude, une borne supérieure M est calculée pour la quantité
Figure imgf000025_0007
ur le domaine
Figure imgf000025_0008
La dérivée de
Figure imgf000025_0009
_ Elle atteint son maximum sur [0, 1] en tmax
Figure imgf000025_0010
On obtient l'expression de
Figure imgf000025_0011
Figure imgf000025_0002
Su
Figure imgf000025_0012
supérieure est valide pour chaque coordonnée, ainsi
Figure imgf000025_0003
Figure imgf000025_0004
Etant donnée cette borne supérieure pour la norme du gradient, une condition suffisante qui garantit la bijection et par là même, le non repliement de l'espace pour des fonctions de déformation élémentaires définies pour une grille composée de cellules
canoniques [0, 1] est donc
Figure imgf000025_0013
Ce résultat peut être étendu à n'importe quelle dimension de cellule de la façon su .ivant .e. S eoi .t.
Figure imgf000026_0003
une cellule cubique
d'arête L. Soit
Figure imgf000026_0008
la fonction
Figure imgf000026_0004
qui transforme cette cellule en cellule canonique
Figure imgf000026_0007
I L J , définie par :
Figure imgf000026_0006
La fonction de forme w définie sur CeII est
Figure imgf000026_0005
et :
Figure imgf000026_0001
De la même façon que pour la cellule canonique, il est possible de conclure qu'une borne supérieure pour la norme du gradient est :
Figure imgf000026_0002
Ce qui conduit à la condition suffisante garantissant pour une cellule de dimension L la bijection et le non-repliement de la déformation élastique élémentaire r :
Figure imgf000026_0009
Etant donné que l'approche présentée ici est basée sur la modélisation mécanique des solides, on choisit de restreindre l'amplitude des déplacements individuels des noeuds au domaine dit des « petites déformations ». Pour ce faire on limite les amplitudes des déformations imposées aux cellules du solide à 10% de leur dimension, soit :
Figure imgf000026_0010
Une conséquence de cette restriction est que la compression et l'étirement de l'espace appliqués au solide à chaque itération sont bornés. Puisque
Figure imgf000026_0011
et
Figure imgf000026_0012
), alors :
Figure imgf000026_0013
L'équation ci-dessus fournit une approximation
Figure imgf000027_0002
des bornes inférieure et supérieure du jacobien de r valable pour toutes les dimensions de cellules :
Figure imgf000027_0001
Minimisation de l'énergie
On considère deux ensembles S et D ayant les cardinaux nS = Card(S) et nD = Card(D). La minimisation de l'énergie de recalage Erβg nécessite le calcul de la distance minimale entre chaque point de S et l'ensemble D. La complexité du calcul brutal, considérant pour chaque point de S sa distance à chaque point de D est s est nS * nD. De nombreuses techniques ont été proposées pour accélérer le calcul de la distance la plus proche entre deux ensembles de points, principalement dans le cadre du recalage rigide basé sur le principe de I1ICP pour « Itérative Closest Point ».
Afin d'accélérer le calcul de Ereg à chaque itération, une carte de distance euclidienne est calculée à partir des points de D. La carte de distance est calculée sur une grille 3D discrétisant l'espace autour de D. A chaque sommet de la grille gi, la distance minimale à D, d(gi,D) est mémorisée. Cette structure de données est initialisée lors de la phase pré-opératoire, juste après la segmentation de l'arbre vasculaire D dans les images ARM.
Le contenu de la carte de distance est ensuite enregistré sous forme de fichier binaire pour être lu au début de l'intervention par le système, avant tout calcul de recalage élastique.
Pour un point source donné s, la valeur de d(s,D) est calculée par interpolation trilinéaire parmi les 8 valeurs lues dans la grille de la carte de distance, associées aux sommets qui entourent s, {gi}i=s1,...,s8. Lors du recalage, à chaque itération, le meilleur déplacement pour chaque nœud de la grille qui discrétise le solide virtuel est calculé à l'aide du gradient de Ereg.
Le calcul du déplacement nodal est réalisé en deux étapes : 1. Le gradient de Ereg par rapport aux déplacements nodaux est estimé par différences finies ; 2. Une recherche le long de la direction de plus forte descente est effectuée afin d'estimer le pas optimal ;
3. L'amplitude du déplacement ainsi obtenu est limitée de manière à satisfaire une contrainte. Etant donné que l'on ne considère que les déplacements d'un nœud à la fois, n, seules les trois coordonnées du vecteur de déplacement nodal
Figure imgf000028_0004
sont impliquées dans
le calcul de la minimisation. La fonction minimisée est
Figure imgf000028_0003
définie par :
Figure imgf000028_0005
où
Figure imgf000028_0006
est la fonction de recalage élémentaire associée au nœud n, et R est la fonction de recalage global assemblée avant l'itération courante.
La stratégie consistant à revenir à la forme initiale du solide virtuel à l'issue de chaque itération permet de supposer que l'ensemble des points source au début de l'itération a déjà été transformé par R selon l'algorithme décrit ci-dessus Ceci nous permet de laisser de coté l'historique du recalage que constitue R. So
Figure imgf000028_0007
ensemble de points source dont la position appartient au support de la fonction de forme w associée au nœud du maillage n. Puisque ce sont les seuls points concernés par le déplacement de n à l'itération de recalage courante, le calcul de l'énergie de recalage peut être restreint à cet ensemble. Ceci nous permet d'exprimer la fonction à minimiser f comme :
Figure imgf000028_0001
La fonction f n'est pas systématiquement différentiable par rapport à
Figure imgf000028_0008
, du fait de la racine carrée portée par la distance, et de l'opérateur min(x, y) qui opère sur les distances entre les nuages de points. La direction de plus grande descente doit donc être estimée localement par différences finies centrées. On obtient ainsi :
Figure imgf000028_0002
où £ est une valeur très petite par rapport à l'amplitude maximale du déplacement possible de n, qui est inférieure à 10% de la taille de la cellule L, d Les
Figure imgf000028_0012
vecteur
Figure imgf000028_0010
quant à eux, forment la base canonique de
Figure imgf000028_0011
Le vecteur de plus grande descente sur f peut s'écrire :
Figure imgf000028_0009
La longueur du pas maximal autorisé par la contrainte de non-repliement est 10% de L, ce qui permet de définir le vecteur maximal de déplacement du nœud n dans la direction de descente sur f :
Figure imgf000029_0001
Le meilleur déplacement nodal au sens de la minimisation de l'énergie de recalage correspond donc au minimum de la fonctio
Figure imgf000029_0005
Une estimation convenable de ce minimum est obtenue en approximant sur l'intervalle [0, 1] la fonction
Figure imgf000029_0003
par une parabole
Figure imgf000029_0004
Les trois coefficients de la parabole A, B et C sont déterminés par les informations locales sur le comportement de f ainsi qu'une évaluation de sa valeur au pas maximal :
Figure imgf000029_0006
. Ainsi :
Figure imgf000029_0002
Deux situations sont alors possibles. Si A > 0 alors P est une parabole convexe et son minimum est atteint en tmin = -B/(2A). Dans ce cas, si tmin > 1 alors la valeur du pas est limité à 1. Si au contraire, A < 0 alors P est une parabole concave et son minimum sur est atteint en tmin = 1.
La figure 7 montre un exemple représentatif d'approximation de la courbe d'énergie de recalage 11 par une parabole P 12 sur l'intervalle [0, 1]. Des mesures d'écart entre les deux courbes, effectuées sur un grand nombre de jeux de données lors du déroulement de la procédure de recalage montre les qualités de cette approximation. La courbe du bas 11 est obtenue en échantillonnant l'énergie de recalage sur l'intervalle [0, 1], celle du haut 12 est le tracé de la parabole P. Malgré un léger écart entre les valeurs des deux courbes, on remarque que les valeurs de t pour lesquelles ces courbes atteignent leur minimum sont quasiment identiques. Pour un niveau de raffinement de la discrétisation du solide virtuel donné, la minimisation de l'énergie est ainsi menée, pas à pas, jusqu'à ce qu'aucune diminution significative ne puisse être obtenue par application d'un déplacement à un des nœuds du maillage.
Dans ce cas, si le niveau de raffinement maximal n'a pas été atteint, la grille est raffinée et la même procédure est répétée. Inversion
Le calcul de l'inverse de la fonction de recalage élastique R va maintenant être décrit. Ce calcul nécessite l'inversion successive de toutes les fonctions de recalage élémentaires : (^)"1. Etant donné qu'aucune expression analytique ne peut être donnée pour (ή)"1, le calcul de (^)"1 (q) est réalisé individuellement pour chaque point considéré q.
Pour simplifier les écritures, les indices i et j ne seront pas précisés dans la suite, et on considéra plus généralement une fonction élémentaire r, associée au déplacement du nœud n. Pour un point donn
Figure imgf000030_0007
é Q e W, son antécédent p = r-1(q) est calculé comme suit. On a vu précédemment que p peut s'écrire comme
Figure imgf000030_0008
P ^ *ï + *p ^ où tp est la racine de la fonction polynomiale f sur l'intervalle [ta, tb]. En fait, puisque
Figure imgf000030_0006
— ' ^ "f — l'intervalle de recherche peut être restreint à [max{-1 , ta}, O].
Le choix de l'expression de la fonction de forme
Figure imgf000030_0009
amène à ^pression suivante de f :
Figure imgf000030_0010
π est un polynôme de degré 3 en t, donc f est de degré 9. Il n'y a pas de formule explicite permettant de trouver les racines d'un polynôme de degré 9 aussi doit on faire appel à une technique itérative de type Newton. En supposant que la valeur de ta est déjà initialisée à max{-1, ta}, on commence avec l'estimation de la solution suivante :
Figure imgf000030_0001
Les itérations sont ensuite menées de manière classique :
μ+i = fi _ J ^p)
Figure imgf000030_0002
p ) L'approximation de
Figure imgf000030_0005
x^-' à chaque itération est donnée par
Figure imgf000030_0003
P *~ ^ ' p u et l'erreur correspondante dans l'espace déformé est ^
Figure imgf000030_0004
La précision de l'approximation de l'inverse est mesurée dans l'espace initial, non déformé, en évaluant la quantité
Figure imgf000031_0002
p n'est pas connu, ce qui rend difficile l'évaluation cette erreur. Afin d'évalue
Figure imgf000031_0003
on s'appuie sur le fait que l'application inverse r"1 est k-lipschitzienne, en d'autres termes :
Figure imgf000031_0001
Cette propriété permet alors de majorer l'erreu par :
Figure imgf000031_0004
On désire poursuivre les itérations de l'algorithme de Newton jusqu'à ce que la précision voulue soit atteinte i.e. lorsque
Figure imgf000031_0005
: Au vu de ce qui précède, cette condition d'arrêt est vérifiée lorsque Cette inégalité
Figure imgf000031_0006
constitue un critère d'arrêt fiable défini dans l'espace non-déformé.
Considérons maintenant le cas d'une fonction de recalage R composée de plusieurs fonctions élémentaires rj . Dans ce cas, la précision du calcul de l'inverse doit être contrôlée à chaque étape d'inversion
Figure imgf000031_0008
d'une fonction élémentaire. Voyons comment nous pouvons évaluer la précision du calcul de R- 1 à l'aide des constantes de Lipschitz des fonctions
Figure imgf000031_0007
Pour plus de clarté, on suppose que
Figure imgf000031_0009
On veut évalue
Figure imgf000031_0010
en ayant, pour chacune des , une
Figure imgf000031_0012
constante de Lipschitz kj .
Soit l'antécédent cherché de p3 = R(p) ayant pour images intermédiaires a figure 8 illustre, de droite à gauche,
Figure imgf000031_0011
les trois étapes de l'approximation de ainsi que l'accumulation successive des
Figure imgf000031_0013
erreurs commises à chaque pas : 1. q2 approxime
Figure imgf000031_0014
L'erreur dans l'espace déformé est _3 := kq3-p3k.
Cette erreur se propage vers l'espace non-déformé comme
Figure imgf000031_0016
et finalement
Figure imgf000031_0015
2. ii approxime
Figure imgf000032_0008
avec une erreur dans l'espace déformé
Figure imgf000032_0007
Cette erreur se propage vers la gauche, i.e. vers l'espace non-déformé, comme
Figure imgf000032_0003
Figure imgf000032_0004
1 3. Finalement s approxime avec une erreur dans l'espace déformé
d
Figure imgf000032_0006
e . Cette erreur se propage comme
Figure imgf000032_0005
L'erreur finale peut ainsi être majorée par la somme des erreurs propagées :
Figure imgf000032_0001
L'expression ci-dessus peut être immédiatement étendue à toute fonction de recalage élastique R. Si R est composée de N fonctions élémentaires et la précision voulue dans l'espace non-déformé est £• alors l'effort d'approximation peut être réparti parmi les N fonctions élémentaires en ajustant le seuil de précision de l'inversion de Newton à :
€
Figure imgf000032_0002
où chaque constante de Lipschitz ki ne doit être évaluée qu'une seule fois étant donné qu'elle ne dépend que de l'amplitude du déplacement appliqué au nœud associé à la fonction élémentaire ri et de la taille de la cellule au niveau de raffinement considéré.
Filtrage du recalaαe élastique et calcul du champ de déplacements Comme décrit ci-dessus, les images per-opératoires sont bruitées et le recalage élastique peut produire occasionnellement des artefacts de recalage. Il est donc important de mettre en place une stratégie permettant de réduire le nombre d'association entre points source et points destination erronées, chacune de ses associations se traduisant par la suite par un déplacement, leur effet risque de se propager jusqu'au calcul de la déformation biomécanique du modèle avec pour conséquence, une perte de précision importante.
Une technique robuste et simple a donc été implémenté pour éliminer les erreurs de recalage grossières, dues à la présence de faux vaisseaux dans les données issues de la segmentation. Cette méthode est divisée en deux phases : un filtrage direct, qui s'appuie sur le calcul de la fonction de recalage R, suivi par un filtrage inverse, lequel fait intervenir l'inverse du recalage : R"1.
Le rôle du filtrage est d'éliminer les associations de points qui ne vérifient pas un critère de précision donné. Soit dmax la distance maximale acceptable entre un point source et son homologue dans D. Le filtrage direct permet d'éliminer les points de S qui ne remplissent pas la condition d(R(s),D) < dmax. L'ensemble de points source filtré est Sf
, défini par :
Figure imgf000033_0002
La figure 9a illustre l'étape de filtrage direct. Comme on peut le constater dans cette figure, les points sources n'ayant pas atteint D sont éliminés.
L'objectif de la procédure de recalage est de générer un champ de déplacements correspondant à la déformation de la structure de l'arbre vasculaire (dans la présentation faite ici, mais plus généralement du tissus organique déformé) définie à partir des images préopératoires par D, et qui l'amènent dans la configuration courante S, obtenue par analyse des images US Doppler per-opératoires. Le filtrage inverse permet de vérifier que le champ de déplacements obtenu par l'inversion de la fonction de recalage R transforme effectivement les points de D en S. Afin de ne pas tenir compte des points source éliminés précédemment, nous restreignons l'ensemble S à Sf et définissons l'ensemble de points destination filtré Db par :
Figure imgf000033_0001
La figure 9b illustre l'étape de filtrage inverse. Comme l'illustre cette figure, les déplacements de D qui n'atteignent pas S sont éliminés.
Finalement, l'ensemble de déplacements de l'arbre vasculaire, de sa configuration de référence vers sa configuration courante, est assemblé en s'appuyant sur ce dernier ensemble de points. Le champ de déplacements destiné au modèle biomécanique est ainsi donné par :
Figure imgf000033_0003
Le recalage élastique décrit en détail ci-dessus permet de construire une association entre deux observations de l'arbre vasculaire cérébral pré- et per-opératoire, de natures très distinctes. Le calcul est guidé par la minimisation de l'énergie de recalage qui quantifie la similarité entre les deux structures. Le recalage optimal étant atteint pour une énergie nulle, l'objectif de l'algorithme de recalage est d'en minimiser itérativement la valeur en faisant évoluer la structure courante, dite 'source', vers la structure de référence, dite 'de destination'.
La rapidité du calcul est obtenue grâce au calcule préalable d'une carte de distance sur laquelle s'appuient l'évaluation de l'énergie de recalage ainsi que la détermination de la direction de plus forte descente, calculées à chaque itération du processus de minimisation.
Les propriétés mathématiques de la famille de fonction de recalage engendrée par la méthode décrite ci-dessus garantissent le respect de critères physiques essentiels, tels que le non-repliement de l'espace, le non-recouvrement de matière et la continuité de la propagation de la déformation. Ce recalage élastique présente également l'avantage d'être inversible avec le degré de précision souhaité.
Sur la base de cette dernière caractéristique, un post-traitement du résultat produit par le recalage est proposé. Ce post-traitement permet d'éliminer les aberrations dues à la présence de bruit dans les données d'imagerie per-opératoire. L'ensemble de déplacements de l'arbre vasculaire finalement retenu constitue la donnée d'entrée du modèle biomécanique dont le rôle est d'en étendre la portée à l'ensemble du volume du cerveau du patient afin de simuler l'effet du brain-shift.
On a illustré à la figure 2 un mode de réalisation du système pour la mise en œuvre du procédé décrit ci-dessus. Le système comprend des moyens de localisation 10, des moyens d'acquisition
20, des moyens de traitement 30 et des moyens d'affichage 40.
Les moyens de localisation permettent la localisation dans l'espace des instruments utilisés par l'utilisateur lors de l'intervention chirurgicale. Ces moyens de localisation peuvent comprend par exemple une caméra stéréoscopique et des marqueurs associés aux différents objets à localiser. Les moyens de localisation 10 permettent également la localisation dans l'espace de la tête du patient lors de l'intervention, notamment pour le repositionnement spatial de l'arbre vasculaire de référence calculé à partir de l'angiographie par résonance magnétique.
Les moyens d'acquisition 20 permettent l'acquisition des images bidimensionnelles pendant l'opération. Ces moyens d'acquisition 20 comprennent par exemple une sonde à ultrasons localisable par les moyens de localisation 10.
Les moyens d'affichage permettent par exemple l'affichage des données préopératoire mise à jour. Ces moyens d'affichage peuvent comprend un écran.
Les moyens de traitement permettent de mettre en œuvre les différentes étapes de calcul du procédé, et notamment : l'étape de traitement 100 de l'image tridimensionnelle du cerveau du patient, acquise avant une intervention chirurgicale, pour obtenir une arborescence artérielle cérébrale de référence du patient,
- l'étape de traitement 200 des images bidimensionnelles du cerveau du patient, acquises pendant l'opération, pour reconstituer au moins partiellement l'arborescence artérielle cérébrale courante du patient, l'étape de détermination 300, à partir de la mise en correspondance des arborescences artérielles cérébrales de référence et courante, du champ de déformation de l'arbre vasculaire représentant le déplacement de l'arbre vasculaire courant par rapport à l'arbre vasculaire de référence, et
- l'application 400 du champ de déformation de l'arbre vasculaire déterminé à un modèle biomécanique du cerveau du patient pour estimer le déplacement du cerveau du patient.
Ces moyens de traitement peuvent être un ordinateur et comprendre des moyens de commande pour le contrôle des différents éléments matériels du système.
L'homme du métier aura compris que de nombreuses modifications peuvent être apportées au procédé décrit sans sortir matériellement des nouveaux enseignements présentés ici. Notamment, le procédé de recalage décrit précédemment en relation avec le recalage de l'arbre vasculaire cérébral peut être mis en œuvre pour permettre le recalage de tout type de tissu organique ayant subit une déformation. Par ailleurs, même si cette étape de recalage a été présentée en relation avec des ensembles de points sources S et destinations D de natures distinctes du fait de la différence entre les zones de recouvrement des données ou de la présence de bruit, l'étape de recalage peut être mise en œuvre sur des jeux de données S et D homogènes. Il est donc bien évident que les exemples qui viennent d'être donnés ne sont que des illustrations particulières en aucun cas limitatives.

Claims

REVENDICATIONS
1. Procédé de traitement d'image pour estimer une déformation d'un cerveau d'un patient, le procédé comprenant : - le traitement (100) d'une image tridimensionnelle du cerveau du patient, acquise avant une intervention chirurgicale, pour obtenir une arborescence artérielle cérébrale de référence du patient, le traitement (200) d'images bidimensionnelles du cerveau du patient, acquises pendant l'opération, pour reconstituer au moins partiellement une arborescence artérielle cérébrale courante du patient, la détermination (300), à partir de la mise en correspondance des arborescences artérielles cérébrales de référence et courante, d'un champ de déformation de l'arbre vasculaire représentant le déplacement de l'arbre vasculaire courant par rapport à l'arbre vasculaire de référence, - l'application (400) du champ de déformation de l'arbre vasculaire déterminé à un modèle biomécanique du cerveau du patient pour estimer le déplacement du cerveau du patient . la génération (500), à partir du déplacement du cerveau estimé, d'au moins une image du cerveau du patient dans laquelle le déplacement du cerveau est compensé.
2. Procédé selon la revendication 1 , caractérisé en ce qu'il comprend en outre le calcul (600) des disparités entre des points de l'image du cerveau générée et une ou plusieurs images bidimensionnelles de contrôle acquise pendant l'opération.
3. Procédé selon la revendication 1 ou la revendication 2, caractérisé en ce que l'image tridimensionnelle acquise avant l'opération est une angiographie par résonance magnétique.
4. Procédé selon l'une des revendications 1 à 3, caractérisé en ce que les images bidimensionnelles acquises pendant l'opération sont des images ultrasonores en mode Doppler.
5. Procédé selon l'une des revendications 1 à 4, caractérisé en ce que le traitement de l'image tridimensionnelle pour obtenir l'arborescence artérielle cérébrale de référence comprend une étape d'orientation de l'image tridimensionnelle par rapport à la tête du patient après l'installation de celui-ci sur une table d'opération.
6. Procédé selon la revendication 5, caractérisé en ce que l'étape d'orientation de l'image tridimensionnelle comprend : la localisation d'un jeu de points de calibrage sur la tête du patient après l'installation de celui-ci sur la table d'opération, la détection de leurs homologues radiométriques dans l'image tridimensionnelle, et - la détermination de transformations rigides permettant de passer des positions des points du jeu de points de calibrage aux positions de leurs homologues dans l'image tridimensionnelle.
7. Procédé selon la revendication 6, caractérisé en ce que les points de calibrage comprennent des points de la surface de la tête du patient après l'installation de celui- ci sur la table d'opération.
8. Procédé selon l'une des revendications 1 à 7, caractérisé en ce que le traitement des images bidimensionnelles pour reconstituer au moins partiellement une arborescence artérielle cérébrale courante du patient comprend une segmentation des images bidimensionnelles pour obtenir la position médiane des vaisseaux et/ou artères du cerveau de sorte à reconstituer au moins partiellement l'arborescence artérielle cérébrale courante du patient.
9. Procédé selon l'une des revendications 1 à 8, caractérisé en ce que l'étape de détermination d'un champ de déformation comprend la détermination de transformations non-linéaires, dites transformations élastiques, permettant le passage des positions des points de l'arborescence artérielle cérébrale de référence aux positions des points correspondants de l'arborescence artérielle courante.
10. Système pour la détermination d'une déformation du cerveau d'un patient, caractérisé en ce qu'il comprend des moyens pour la mise en œuvre du procédé selon l'une des revendications 1 à 9.
PCT/EP2009/062826 2008-10-03 2009-10-02 Procédé de traitement d'image pour estimer une déformation d'un cerveau d'un patient Ceased WO2010037850A2 (fr)

Priority Applications (2)

Application Number Priority Date Filing Date Title
US13/122,092 US9008381B2 (en) 2008-10-03 2009-10-02 Image processing method for estimating a brain shift
EP09783691A EP2353144A2 (fr) 2008-10-03 2009-10-02 Procédé de traitement d'image pour estimer une déformation d'un cerveau d'un patient

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
FR0856712 2008-10-03
FR0856712A FR2936890B1 (fr) 2008-10-03 2008-10-03 Procede de traitement d'image pour estimer une deformation d'un cerveau d'un patient

Publications (2)

Publication Number Publication Date
WO2010037850A2 true WO2010037850A2 (fr) 2010-04-08
WO2010037850A3 WO2010037850A3 (fr) 2010-06-10

Family

ID=40561817

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/EP2009/062826 Ceased WO2010037850A2 (fr) 2008-10-03 2009-10-02 Procédé de traitement d'image pour estimer une déformation d'un cerveau d'un patient

Country Status (4)

Country Link
US (1) US9008381B2 (fr)
EP (1) EP2353144A2 (fr)
FR (1) FR2936890B1 (fr)
WO (1) WO2010037850A2 (fr)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2468207A1 (fr) * 2010-12-21 2012-06-27 Renishaw (Ireland) Limited Procédé et appareil pour analyser des images
US9668710B2 (en) 2011-03-18 2017-06-06 Koninklijke Philips N.V. Tracking brain deformation during neurosurgery

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2706372A1 (fr) * 2012-09-10 2014-03-12 Esaote S.p.A. Procédé et appareil d'acquisition d'images à ultra-son
WO2014165787A1 (fr) * 2013-04-05 2014-10-09 Old Dominion University Research Foundation Système et procédé de détection de structures
CN106890025B (zh) * 2017-03-03 2020-02-28 浙江大学 一种微创手术导航系统和导航方法
FR3073135B1 (fr) * 2017-11-09 2019-11-15 Quantum Surgical Dispositif robotise pour une intervention medicale mini-invasive sur des tissus mous
US10292769B1 (en) 2018-08-07 2019-05-21 Sony Corporation Surgical assistive device and method for providing assistance in surgery of anatomical portions of internal organ affected by intraoperative shift
EP4102513A1 (fr) * 2021-06-10 2022-12-14 Inteneural Networks Inc. Procédé et appareil permettant d'enregistrer un patient en neurochirurgie et déterminer le déplacement du cerveau pendant la chirurgie en utilisant l'apprentissage machine et une caméra de profondeur tridimensionnelle stéréoptique avec un système de cartographie de surface

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7103399B2 (en) 2003-09-08 2006-09-05 Vanderbilt University Apparatus and methods of cortical surface registration and deformation tracking for patient-to-image alignment in relation to image-guided surgery
WO2005025404A2 (fr) * 2003-09-08 2005-03-24 Vanderbilt University Appareil et procedes de reperage de surface corticale et de depistage de deformation pour alignement patient-sur-image en rapport avec une chirurgie guidee par l'image
JP4630564B2 (ja) * 2004-03-30 2011-02-09 国立大学法人浜松医科大学 手術支援装置、方法及びプログラム
US7620223B2 (en) * 2004-04-22 2009-11-17 Siemens Medical Solutions Usa, Inc. Method and system for registering pre-procedural images with intra-procedural images using a pre-computed knowledge base
US7711405B2 (en) * 2004-04-28 2010-05-04 Siemens Corporation Method of registering pre-operative high field closed magnetic resonance images with intra-operative low field open interventional magnetic resonance images
US20060036167A1 (en) * 2004-07-03 2006-02-16 Shina Systems Ltd. Vascular image processing
WO2006036842A2 (fr) * 2004-09-24 2006-04-06 The University Of North Carolina At Chapel Hill Procedes, systemes et produits-programmes informatiques pour superposition hierarchique entre un modele de vaisseau sanguin et de surface de tissu pour un sujet et une image de vaisseau sanguin et de surface de tissu pour ce sujet
DE102005036322A1 (de) * 2005-07-29 2007-02-15 Siemens Ag Registrieren intraoperativer Bilddatensätze mit präoperativen 3D-Bilddatensätzen auf Basis optischer Oberflächenextraktion
DE602005009370D1 (de) * 2005-10-06 2008-10-09 Medcom Ges Fuer Medizinische B Registrierung von 2D Ultraschall-Daten und einem 3D Bilddatensatz

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2468207A1 (fr) * 2010-12-21 2012-06-27 Renishaw (Ireland) Limited Procédé et appareil pour analyser des images
WO2012085511A1 (fr) 2010-12-21 2012-06-28 Renishaw (Ireland) Limited Procédé et appareil pour l'analyse d'images
US9463073B2 (en) 2010-12-21 2016-10-11 Renishaw (Ireland) Limited Method and apparatus for analysing images
US9668710B2 (en) 2011-03-18 2017-06-06 Koninklijke Philips N.V. Tracking brain deformation during neurosurgery

Also Published As

Publication number Publication date
FR2936890B1 (fr) 2010-12-03
FR2936890A1 (fr) 2010-04-09
US9008381B2 (en) 2015-04-14
WO2010037850A3 (fr) 2010-06-10
US20110257514A1 (en) 2011-10-20
EP2353144A2 (fr) 2011-08-10

Similar Documents

Publication Publication Date Title
WO2010037852A1 (fr) Procédé pour le recalage d&#39;un ensemble de points dans des images
WO2010037850A2 (fr) Procédé de traitement d&#39;image pour estimer une déformation d&#39;un cerveau d&#39;un patient
JP5832523B2 (ja) 光コヒーレンストモグラフィのための動き補正および画像改善の方法および装置
Ralston et al. Deconvolution methods for mitigation of transverse blurring in optical coherence tomography
EP2078222B1 (fr) Microscopie à ouverture synthétique interférométrique
US20150371400A1 (en) Segmentation and identification of closed-contour features in images using graph theory and quasi-polar transform
FR2862791A1 (fr) Procede de segmentation d&#39;une structure et systeme d&#39;imagerie associe
FR3095332A1 (fr) Méthode de planification d’une ablation de tissus basée sur l’apprentissage profond
FR3035531A1 (fr)
Giovannelli et al. Regularization and Bayesian Methods for Inverse Problems in Signal and Image Processing
Ma et al. Large area kidney imaging for pre-transplant evaluation using real-time robotic optical coherence tomography
WO2018104609A1 (fr) Procede de suivi de cible dans une sequence d&#39;images medicales, dispositif, equipement terminal et programmes d&#39;ordinateurs associes
Pennec Statistical computing on manifolds for computational anatomy
Saha Non-rigid registration of De-noised ultrasound breast tumors in image guided breast-conserving surgery
Al Moussawi et al. 3D Reconstruction of blood vessels
EP3509483B1 (fr) Systeme et procede pour reconstruire un signal physiologique d&#39;un systeme dynamique artere/tissu/veine d&#39;un organe dans un espace surfacique
FR3143314A1 (fr) Dispositif d’assistance à la planification d’une intervention mini-invasive
WO2010046573A1 (fr) Procédé d&#39;estimation de la concentration d&#39;un traceur dans un ensemble de structures tissulaires, support d&#39;information et dispositif correspondants.
FR2838628A1 (fr) Procede d&#39;assistance et de guidage de navigation d&#39;un outil dans des structures anatomiques
Carton Image segmentation and registration using machine learning for brain shift compensation in image-guided neurosurgery
Nair et al. Convolutional Neural Networks Enable Direct Strain Estimation in Quasistatic Optical Coherence Elastography
Zhang et al. 3D Fourier ptychographic microscopy based on the beam propagation method and time-reversal scheme
CA2277788C (fr) Methode de detection d&#39;une discontinuite geologique presente dans un milieu par utilisation du flux optique
FR3143313A1 (fr) Dispositif d’assistance à la planification d’une intervention mini-invasive sur un os
FR3145860A1 (fr) Procédé de génération d’un modèle d’un ensemble osseux

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 09783691

Country of ref document: EP

Kind code of ref document: A2

NENP Non-entry into the national phase

Ref country code: DE

REEP Request for entry into the european phase

Ref document number: 2009783691

Country of ref document: EP

WWE Wipo information: entry into national phase

Ref document number: 2009783691

Country of ref document: EP

WWE Wipo information: entry into national phase

Ref document number: 13122092

Country of ref document: US