[go: up one dir, main page]

WO2010037852A1 - Procédé pour le recalage d'un ensemble de points dans des images - Google Patents

Procédé pour le recalage d'un ensemble de points dans des images Download PDF

Info

Publication number
WO2010037852A1
WO2010037852A1 PCT/EP2009/062829 EP2009062829W WO2010037852A1 WO 2010037852 A1 WO2010037852 A1 WO 2010037852A1 EP 2009062829 W EP2009062829 W EP 2009062829W WO 2010037852 A1 WO2010037852 A1 WO 2010037852A1
Authority
WO
WIPO (PCT)
Prior art keywords
registration
points
function
point
brain
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/062829
Other languages
English (en)
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 EP09783693A priority Critical patent/EP2353145A1/fr
Priority to US13/122,116 priority patent/US8724926B2/en
Publication of WO2010037852A1 publication Critical patent/WO2010037852A1/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
    • 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 surgery, and more particularly to a system and method for resetting an organic tissue in two images.
  • the precise location of a target is essential during a surgical procedure to reduce the morbidity rate.
  • this deformation can 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.
  • compensation for the deformation of the organic tissues can be performed on the following principle.
  • a first image of a region of the organic tissues is acquired before the surgical operation.
  • a second image of the same region is acquired.
  • System processing means process these two images to align the dots in the two images to compensate for the deformation of the organic tissues.
  • the processing means implement an image registration method.
  • the registration methods do not make it possible to precisely align the points in two images when the characteristics of the data present in the two images are heterogeneous.
  • An object of the present invention is to propose a method of resetting a set of source points in a first image and a set of destination points in a second image, the set of source points being of a nature different from the set of destination points because of the difference between the areas of data overlap or the presence of noise.
  • the invention provides a method for resetting a set of source points S in a reference image representing an organic tissue and a set of destination points D in a current image representing the deformed organic tissue, the method comprising determining a reset function R by iterative minimizing of a registration energy E reg , the registration function corresponding to a transformation for passing source points S of the reference picture to the corresponding destination points D in the frame.
  • current image characterized in that the registration function R is a resilient registration function satisfying conditions to ensure compliance with physical criteria of the organic tissue, the conditions comprising:
  • the Jacobian of the elastic registration function is greater than 0, the registration function is continuously differentiable, the registration function is bijective.
  • the elastic registration function is determined so as to ensure the respect of physical criteria associated with the organic tissue that is to be recalibrated. These physical criteria to be respected are the following: non-folding of the organic tissue on itself (which results in the fact that the Jacobian of the function of registration must be greater than 0), - regularity of the propagation of the deformation of the organic tissue (which results in the fact that the registration function is continuously differentiable), non-superposition of two distinct points of the initial space at a single point in the deformed space and range of the deformation function extended to the whole space (which results in the fact that the function is bijective).
  • the set of source points is of a nature distinct from the set of destination points
  • the registration energy E i of the registration function (R) satisfies the relationship next : with: - s a point of the set of source points S, - R (s) the transform of the point s by the resetting function R 1 - d (R (s), D) the Euclidean distance between the transform of the point s by the registration function and the set of destination points
  • the method comprises a direct filtering step in which each point s of the set of source points S whose association with a point d of the set of destination points D checks a rejection criterion is deleted so as to obtain a filtered set of source points S ', the rejection criterion is that the distance between: o transforming the point s of the set of source points S by the resetting function and o the corresponding point d of the set of destination points D is greater than a predetermined maximum value
  • the invention also relates to a system for resetting a set of source points in a reference image representing an organic tissue and a set of destination points in a current image representing the deformed organic tissue, characterized in that comprises means for implementing the method described above.
  • the invention also relates to a computer program product comprising program code instructions recorded on a medium that can be used in a computer, characterized in that it comprises instructions for implementing the method described above.
  • 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 recalibration energy curve by a parabola
  • FIG. 8 illustrates three steps of the approximation of the inverse of the resetting function
  • FIG. 9 illustrates steps of direct and indirect filtering.
  • This registration method will now be described in more detail. This registration method will be presented as part of an image processing method for estimating a deformation of a brain.
  • this registration method is completely independent of the treatment method for estimating a deformation of the brain.
  • - can be implemented for the registration of other organic tissues than the brain, and / or
  • the precise location of the surgical target is essential for reducing morbidity during surgical excision of a brain tumor.
  • deformity of the soft tissues of the brain may occur during the procedure. Because of this deformation of the brain, preoperative data no longer correspond to reality, and neuronavigation is greatly compromised.
  • 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 treatment method for estimating a deformation of the brain allows the preoperative data to be updated according to the deformation of the brain during the procedure.
  • the deformation of the particular anatomical structure that is the cerebrovascular tree Because of its distribution in the volume of the parenchyma, 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 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 deformation of the brain comprises the following steps:
  • 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.
  • the term "rigid registration” means a geometric transformation composed of a rotation and a translation. A calculation then makes it possible to find the transformation between the two coordinate systems.
  • 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.
  • 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. Knowing the position of the reference cerebral vascular tree with respect to the calibration points, and knowing the position of the calibration points once the patient is installed on the operating table, 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 patient known as the "common cerebral arterial tree".
  • 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. After acquiring a series of two-dimensional images, these are processed to locate the patient's current cerebrovascular tree.
  • 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 strain field representing the displacement of the current vascular tree relative to the reference vascular tree.
  • the reference vascular tree is recalibrated with respect to the current vascular tree.
  • 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 cerebral vascular 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.
  • three-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 resection of tissue, the registration function must be continuous and even continuously differentiable ie differentiable at any point in the domain and having no 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.
  • 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 tissue and a follow-up of potentially noisy anatomical structure representative of the displacement of the brain thanks to the cerebrovascular tree.
  • step 500 of the method the overall displacement of the patient's brain is estimated.
  • 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 homologues in the image generated from the estimated global deformation 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.
  • this validation is carried out immediately after the acquisition of the data on which the deformation algorithm is based in order to guarantee its relevance.
  • the method described above is perfectly adapted to a continuous monitoring of the deformation of the soft tissues of the brain, both by its response times which are of the order of a few seconds, and by the repeatability of the data acquisition procedure. which does not require a major interruption of the intervention or the installation of bulky equipment.
  • An elastic registration is an application i L • ⁇ m that minimizes some 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 resetting 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.
  • we present the elastic strain calculation framework which allows to robustly establish an elastic correspondence between the data. considered.
  • the specificities of this approach lie in the mathematical properties of the transformation R, demonstrated below, which guarantee its physical realism.
  • 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
  • - d (d, R (S) represents the distance between a point d of the set D and the transform by R of the set S.
  • 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 P e ⁇ m ⁇ # e W, a predecessor can be computed.
  • 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 reg 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).
  • An elementary registration function r is defined by propagating the displacement of the node n to
  • 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 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.
  • a nodal value (displacement or other quantity) is interpolated through the discretized volume by means of functions of form w.
  • a function function associated with a node and evaluated at a point of the domain gives the proportion of the nodal value transferred at this 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.
  • 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: V (n) be the closed union of the cells surrounding node n.
  • the form function w is defined on the union of the neighboring cells around the node n by a change of variable and a scaling in order to adapt its canonical definition domain I 'J to the dimensions of the cells in the grid of discretization 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 plot of the shape function thus obtained w: M. ⁇ M.
  • Tt ⁇ ⁇ 0.38 ⁇ -A canonical [0, 1] is therefore: 3 ⁇ 3
  • ⁇ be the function that transforms this cell into a canonical cell
  • 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;
  • a search along the direction of greatest descent is performed to estimate the optimal pitch
  • f (H) E ⁇ fig (ro R) where r 'P * ⁇ * P ⁇ * ⁇ is the elementary registration function associated with the node n, and R is the global registration function assembled before the current iteration.
  • Our f £ [U, Ij A suitable estimate of the minimum is obtained by approximating on the interval [0, 1] the function th ⁇ f ( tu max) by a parabola P ' *> ⁇ At 2 + Bt + C
  • 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:
  • 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. Despite a slight difference between the values of the two curves, we note that 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. In this case, if the maximum refinement level has not been reached, the grid is refined and the same procedure is repeated.
  • indices i and j will not be specified in the following, and one considered more generally an elementary function r, associated with the displacement u of the node n.
  • the intraoperative images are noisy and the elastic registration may occasionally produce registration artifacts. It is therefore important to put in place a strategy that makes it possible to reduce the number of associations between source points and erroneous destination points, each of its associations being subsequently translated by a displacement, their effect is likely to propagate up to the calculation of the biomechanical deformation of the model resulting in a significant loss of precision.
  • 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 objective 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 O , 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. Finally, the set of displacements of the vascular tree, from its reference configuration to its current configuration, is assembled based on this last set of points.
  • the displacement field for the biomechanical model is thus given by:
  • 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 20, 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 setting in correspondence of the reference and current cerebral arterial trees, the deformation field of the vascular tree representing the displacement of the current vascular tree relative to the reference vascular tree, and the application 400 of the deformation field of the vascular tree.
  • vascular tree determined to a biomechanical model of the patient's brain to estimate the displacement of the brain of the patient.
  • 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é pour le recalage d'un ensemble de points sources (S) dans une image de référence représentant un tissu organique et d'un ensemble de points destinations (D) dans une image courante représentant le tissu organique déformé, le procédé comprenant la détermination d'une fonction de recalage (R) par minimisation itérative d'une énergie de recalage Ereg, la fonction de recalage correspondant à une transformation permettant de passer des points sources (S) de l'image de référence aux points destinations (D) correspondants dans l'image courante, caractérisé en ce que la fonction de recalage (R) est une fonction de recalage élastique satisfaisant des conditions pour garantir le respect de critères physiques du tissu organique, les conditions comprenant : - le Jacobien de la fonction de recalage élastique est supérieur à 0, la fonction de recalage est continûment différentiable, la fonction de recalage est bijective.

Description

Procédé pour le recalage d'un ensemble de points dans des images
DOMAINE DE L'INVENTION
La présente invention concerne le domaine de la chirurgie assistée par ordinateur, et plus particulièrement un système et une méthode permettant le recalage d'un tissu organique dans deux images.
PRESENTATION DE L'ART ANTERIEUR
La localisation précise d'une cible est essentielle au cours d'une intervention chirurgicale pour réduire le taux de morbidité.
Au cours d'une intervention chirurgicale, une déformation des tissus organiques peut se produire.
Par exemple, au cours d'une intervention chirurgicale d'ablation d'une tumeur cérébrale, cette déformation peut intervenir 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 des tissus organiques et de retrouver les données préopératoires du patient pendant la procédure chirurgicale, une compensation de la déformation des tissus organiques peut être réalisée sur le principe suivant.
Une première image d'une région des tissus organiques est acquise avant l'opération chirurgicale. Pendant l'opération chirurgicale, une seconde image de la même région est acquise. Des moyens de traitement du système traitent ces deux images pour aligner les points dans les deux images afin de compenser la déformation des tissus organiques. Afin d'aligner les points dans les deux images, les moyens de traitement mettent en œuvre un procédé de recalage d'images.
Toutefois, les procédés de recalage ne permettent pas d'aligner avec précision les points dans deux images lorsque les caractéristiques des données présentes dans les deux images sont hétérogènes.
Un but de la présente invention est de proposer un procédé de recalage d'un ensemble de points sources dans une première image et d'un ensemble de points destinations dans une deuxième image, l'ensemble de points sources étant de nature distincte de l'ensemble de points destinations, du fait de la différence entre les zones de recouvrement des données ou de la présence de bruit. RESUME DE L'INVENTION
A cet effet l'invention prévoit un procédé pour le recalage d'un ensemble de points sources S dans une image de référence représentant un tissu organique et d'un ensemble de points destinations D dans une image courante représentant le tissu organique déformé, le procédé comprenant la détermination d'une fonction de recalage R par minimisation itérative d'une énergie de recalage Ereg, la fonction de recalage correspondant à une transformation permettant de passer des points sources S de l'image de référence aux points destinations D correspondants dans l'image courante, caractérisé en ce que la fonction de recalage R est une fonction de recalage élastique satisfaisant des conditions pour garantir le respect de critères physiques du tissu organique, les conditions comprenant :
- le Jacobien de la fonction de recalage élastique est supérieur à 0, la fonction de recalage est continûment différentiable, la fonction de recalage est bijective. Ainsi, la fonction de recalage élastique est déterminée de sorte à garantir le respect de critères physiques associés au tissu organique que l'on souhaite recaler. Ces critères physiques à respecter sont les suivants : non repliement du tissu organique sur lui-même (ce qui se traduit par le fait que le jacobien de la fonction de recalage doit être supérieur à 0), - régularité de la propagation de la déformation du tissu organique (ce qui se traduit par le fait que la fonction de recalage est continûment différentiable), non-superposition de deux points distincts de l'espace initial en un point unique de l'espace déformé et portée de la fonction de déformation étendue à l'espace entier (ce qui se traduit par le fait que la fonction est bijective).
Des aspects préférés mais non limitatifs de la présente invention sont les suivants : l'ensemble de points sources est de nature distincte de l'ensemble de points destinations, l'énergie de recalage E^ de la fonction de recalage (R) satisfait la relation suivante :
Figure imgf000004_0001
avec : - s un point de l'ensemble de points sources S, - R(s) le transformé du point s par la fonction de recalage R1 - d(R(s), D) la distance euclidienne entre le transformé du point s par la fonction de recalage et l'ensemble de points destinations, noté D, le procédé comprend une étape de filtrage direct dans laquelle chaque point s de l'ensemble de points sources S dont l'association avec un point d de l'ensemble de points destinations D vérifie un critère de rejet est supprimé de sorte à obtenir un ensemble filtré de points source S', le critère de rejet est que la distance entre : o le transformé du point s de l'ensemble de points sources S par la fonction de recalage et o le point correspondant d de l'ensemble de points destinations D est supérieure à une valeur maximale prédéterminée, Le procédé comprend en outre une étape de filtrage inverse dans laquelle chaque point d de l'ensemble de points destinations D dont l'association avec un point s de l'ensemble de points sources S vérifie un autre critère de rejet est supprimé de sorte à obtenir un ensemble filtré de points destinations D', l'autre critère de rejet est que la distance entre : o le transformé inverse du point d de l'ensemble de points destinations D par la fonction de recalage inverse et o le point correspondant s' de l'ensemble filtré de points sources S1 est supérieure à la valeur maximale prédéterminée. L'invention concerne également un système pour le recalage d'un ensemble de points source dans une image de référence représentant un tissu organique et d'un ensemble de points destinations dans une image courante représentant le tissu organique déformé, caractérisé en ce qu'il comprend des moyens pour la mise en œuvre du procédé décrit ci-dessus. L'invention concerne également un produit programme d'ordinateur comprenant des instructions de code programme enregistré sur un support utilisable dans un ordinateur, caractérisé en ce qu'il comprend des instructions 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 oeuvre 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é de recalage selon l'invention. Ce procédé de recalage sera présenté dans le cadre d'un procédé de traitement d'image pour estimer une déformation d'un cerveau.
Toutefois, il est bien entendu que ce procédé de recalage est totalement indépendant du procédé de traitement pour estimer une déformation du cerveau.
Plus précisément, le procédé de recalage selon l'invention :
- peut être mis en œuvre pour le recalage d'autres tissus organiques que le cerveau, et/ou
- peut être mis en œuvre dans d'autres types de procédés de traitement d'image.
Principe général du procédé de traitement pour estimer une déformation du cerveau
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. Le procédé de traitement pour estimer une déformation du cerveau permet 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, on 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.
On 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é pour estimer une déformation du cerveau
En référence à la figure 1 , on a illustré différentes étapes du procédé pour estimer la déformation du cerveau. Le procédé de traitement d'image pour estimer une déformation du cerveau 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 vascυlaires 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 vascυlaire 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 tridimensionnelles 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 recalage é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.
D . jrp3 πτ)3
Un recalage élastique est une application i L • ^ m 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 recalaae
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.
v ' sÇS v ' d€D 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 Ereg 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 imgf000017_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 Ereg en éliminant la constante multiplicative qui la précède, ce qui nous amène à la définition suivante de Ereg :
Figure imgf000017_0002
ses
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 : ^ ' .
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, dX^ 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.
dV = (dXi x dX2) • dX3 = dX.% dJζ.2 cDζ.3
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.
|JΛ(p) | dV
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 :
Vp e R\ |Jβ(p)| > 0 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 P e ^mα#e W, 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 vasculaire 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 ^ ^- * 1^CLgG[H) Af1n ^6 simplifier ces considérations, on considère les fonctions R surjectives, i.e. telles que lma9e{R) = "*- . || 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 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 Ereg 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 noeuds 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 à
'ensemble des points de m sous les contraintes suivantes :
La frontière du voisinage de n, dV(n) est inchangée i.e. Vp G dV(ή), r(p) = p ;
LP voisinage du n reste globalement inchangé i.e. r(V(n)) = V(n) ;
Tous les points à l'extérieur de ce 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
^ I 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 noeuds 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 τπ>3 τrp3 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 Λ , 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. Soit ^ •? " ' II ' J Jl J l'ensemble des 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 :
/? : p t→ [ r1/ o . . . o r j o . . . o rj1 o . . . o r J | (p)
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 Hj. Si chaque ή est une τπ>3 y τπ>3 bijection alors la fonction de recalage R est également une bijection et son inverse est donnée par :
Z?"1 : p ι→ (p)
Figure imgf000022_0002
]πι3 ^ τπ)3
Si chaque ή est une application continuement différentiable alors la
fonction de recalage R l'est également. Soit ^0 ^ un point quelconque et pk le point transformé après application de k déformations élémentaires. Soit K = 11 + 12 + ... +
IJ , alors ^ „ ' "• ' * ' est 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 imgf000023_0001
Finalement, si tous les r1, 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 imgf000023_0002
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 r1, qui composent R.
Etude des fonctions de recalage élémentaire Expression générale de f,
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 w . K -→ [U, IJ j_a va|eur fje |a 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 r1,. 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 la 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 :
r : , p \→ p + w(p)!î
Les conditions sur la zone d'influence de r se traduisent par des contraintes sur le support de w de la façon suivante :
Vp e dV{μ), r(p) = p =* w(p) = 0 Vp i V(n), r(p) = p => w{p) = 0 Remarquons que le support de w est W^ ' ^ ', qui est un compact ^ .
Différentiabilité
La différentiabilité de r découle de celle de w. Si w • ' ^3 → M e C1O^3), alors r G ^ v^ ) et la matrice jacobienne de r s'écrit :
d %w" uvτ dwi_
Figure imgf000024_0001
où Id3><3 est la matrice identité 3 >< 3.
Donc pour que r soit ^ V^ ) W 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 :
Vp 6 Ë3, l + Vzufp) - 11 > O
Puisque P .
Figure imgf000024_0002
II est une application continue à support compact, soit
P«- e 13Ie point où II ^ Hatteint son maximum Af := || Vw(PfIMJ||, M > O- Alors l'équation précédente est équivalente à :
I. → M 1
M
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 imgf000025_0001
Son expression es *t τr( vt ') = t2( vS - 2t) \
Une fonction de forme w définie sur la cellule « canonique » I ' J et prenant la valeur 1 au nœud n = \ ^ΛΛ) peut ainsi être obtenue comme : ιυ(p) = ?t;(pi, p2, p3) '= π(pi)7r(p2)π(p3)
La figure 6a montre la fonction de forme restreinte à une dimension, i.e. w = n, 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 I ' 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 w : M. → M.
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 recalage élémentaire r. Pour définir cette amplitude, une borne supérieure M est calculée pour la quantité H w H sur le domaine 1^' M .
La dérivée de π est π'(f) = 6^1 ~ 0. Elle atteint son maximum sur [0, 1] en tmax = 1/2 où ^(U) = 3/2 Qn obtient l'expression de V W : V u'(p) =
Figure imgf000026_0001
SuJM]3, on .(^(PO^M*))' ≤ (§ * 1 * tf = (f)2. Cette borne
supérieure est valide pour chaque coordonnée, ainsi " w" — V ^2/ t et donc
Figure imgf000026_0002
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
Tt\\ < 0.38 < -A canoniques [0, 1] est donc : 3^3
Ce résultat peut être étendu à n'importe quelle dimension de cellule de la façon su .ivant .e. S oo .it. CeIl „ := [ [α λ, a + L J] x \ ιβy,' βH + Ll J x [ L7/.î 7 ' + Ll J une cel „lul .e cubique
d'arête L. Soit Φ la fonction qui transforme cette cellule en cellule canonique
[O. lf , définie par : φ(Pi> P2i P3) Ψ2{P2), ΨsiPs)) = ((Pi - a) /L, (p2 - β)/L, (p3 - l)jL) La fonction de forme w définie sur CeII est lϋ(p) = W(P1, p2, p^) = ττ(φ1{p1))τî(φ2(p2))τr{ψ2(P2))eV
Vw(p)
Figure imgf000026_0004
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 : 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 :
||1?|| < 0.38 L
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 imgf000027_0001
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 II u II — -W I^ et lfv »(p)|| < (3V3)/(2I) alore :
L'équation ci-dessus fournit une approximation ( (3 v 3)/20 ~ 0.25980 • • • ) des bornes inférieure et supérieure du jacobien de r valable pour toutes les dimensions de cellules :
Vp € CeIl, 0.74 < |Jr(p)| < 1.26
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 Ereg 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 u sont impliquées dans
le calcul de la minimisation. La fonction minimisée est ~ définie par : f(H) = Eτfig(r o R ) où r ' P *~* P ~*~
Figure imgf000028_0001
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. rsn (- 9 Soit •- ' *& L l'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 : f(H) = Ere_g(r) = J2d(φi), D) iel La fonction f n'est pas systématiquement différentiable par rapport à u , 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 imgf000029_0001
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'où e ^ •"/ *-". Les vecteurs ι * J ^=I.2, 3_ quant à eux, forment la base canonique de IR . Le vecteur de plus grande descente sur f peut s'écrire :
Figure imgf000029_0002
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_0003
Le meilleur déplacement nodal au sens de la minimisation de l'énergie de recalage correspond donc au minimum de la fonction f ' 0 + * « max) = /(* u max) pOur f £ [U, Ij Une estimation convenable de ce minimum est obtenue en approximant sur l'intervalle [0, 1] la fonction t h→ f(t u max) par une parabole P ' * >→ At2 + Bt + C
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 :
J \ U mΑx) Ajnsj •
P(O) = C = f(t)
P'(0) = £ = v7 - l?max = -L/10 0 f II
P(I) = A + B + C = /( t?max) => A = /(!?„,„) + L/10 || V/ || - f(~0 ) 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 (Hj)"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 u du nœud n.
Pour un point donné ^ e \n>, son antécédent p = r-1(q) est calculé comme suit. On a vu précédemment que p peut s'écrire comme P = Q + ^p w ou tp est la racine de la fonction polynomiale f sur l'intervalle [ta, tb]. En fait, puisque
— v / — l'intervalle de recherche peut être restreint à [max{-1, ta}, O]. Le choix de l'expression de la fonction de forme w{puP2,p3) = π(Pi)*{pMP3)amène à reXpression suivante de f :
/ : t H→ t + π(qι + tu1)π(q2 4- tιi2)π(q3 + tu2) π 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 imgf000031_0001
Les itérations sont ensuite menées de manière classique :
Figure imgf000031_0002
L'approximation de VMy à chaque itération est donnée par P '~ Q ' τp u et l'erreur correspondante dans l'espace déformé est
Hr(PO - QlI = LWi)IIKII.
La précision de l'approximation de l'inverse est mesurée dans l'espace initial, non déformé, en évaluant la quantité IIP ~~ PII. p n'est pas connu, ce qui rend difficile l'évaluation cette erreur. Afin d'évaluer IIP ~ PII1 on s'appuie sur le fait que l'application inverse r"1 est k-lipschitzienne, en d'autres termes :
3k > 0, Vp,q G R2, Wr-1M - r"1^ < fc||p - q||
Cette propriété permet alors de majorer l'erreur HP ~ P 11 par :
IIP' - PlI = Hr- )) - T-1Cq)H < λφ-(p') - q|| = λ;|/(Q|||^||
On désire poursuivre les itérations de l'algorithme de Newton jusqu'à ce que la précision voulue soit atteinte i.e. lorsque IIP ~ P lI ^ e: Au vu de ce qui précède, cette condition d'arrêt est vérifiée lorsque U VS) I
Figure imgf000031_0003
H u II ). Cette inégalité 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
11 doit être contrôlée à chaque étape d'inversion i 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 5 .
Pour plus de clarté, on suppose que — 1 3 ° r2 ° π
On veut évaluer /» = l 2 3 en ayant, pour chacune des i , une constante de Lipschitz kj .
Soit p € M l'antécédent cherché de p3 = R(p) ayant pour images intermédiaires p1 = r1(p) et p2 = r2(r1(p)) = r2(p1). La figure 8 illustre, de droite à gauche, les trois étapes de l'approximation de ' ^*3' ainsi que l'accumulation successive des erreurs commises à chaque pas : 1. q2 approxime P2 TT r3 'Ps . L'erreur dans l'espace déformé est _3 := kq3-p3k.
Cette erreur se propage vers l'espace non-déformé comme 3 '3' 3 3Λ/2 e{ finalement N " P» ≤ ^ i
2. r1 approxime Qi = '2 (<Ï2}avec une erreur dans l'espace déformé
Figure imgf000032_0001
Cette erreur se propage vers la gauche, i.e. vers l'espace non-déformé, comme
Figure imgf000032_0002
r = rr1 (ri)
3. Finalement s approxime avec une erreur dans l'espace déformé
€1 := Ilsi — 1*1 II Ik - r II < ^-, Jh de ' ' . Cette erreur se propage comme H a 1 II -^ tIM.
L'erreur finale peut ainsi être majorée par la somme des erreurs propagées : e = ||s - p|| < ||s - r|| + ||r - q|| + ||q - p|| < e^i + e2k2k1 + e^^h
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_0003
i=l 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ée 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 :
S* := {s G S | d(R(s), D) < dmαx}
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 O, 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 : Db := {d € D I d{R-χ{d), S1) < dmαx}
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 :
Disp := {u = âb, a e Dh, b = FT1 (a)} 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 oeuvre 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 présentant des nombres de points différents, l'étape de recalage peut être mise en œuvre sur des jeux de données S et D comportant le même nombre de points. 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é pour le recalage d'un ensemble de points sources (S) dans une image de référence représentant un tissu organique et d'un ensemble de points destinations (D) dans une image courante représentant le tissu organique déformé, le procédé comprenant la détermination d'une fonction de recalage (R) par minimisation itérative d'une énergie de recalage Ereg, la fonction de recalage correspondant à une transformation permettant de passer des points sources (S) de l'image de référence aux points destinations (D) correspondants dans l'image courante, caractérisé en ce que la fonction de recalage (R) est une fonction de recalage élastique satisfaisant des conditions pour garantir le respect de critères physiques du tissu organique, les conditions comprenant : le Jacobien de la fonction de recalage élastique est supérieur à 0, la fonction de recalage est continûment différentiable, - la fonction de recalage est bijective.
2. Procédé selon la revendication 1 , caractérisé en ce que l'ensemble de points sources est de nature distincte de l'ensemble de points destinations.
3. Procédé selon l'une des revendications 1 ou 2, caractérisé en ce que l'énergie de recalage Ereg de la fonction de recalage (R) satisfait la relation suivante :
Figure imgf000037_0001
avec : - s un point de l'ensemble de points sources S,
- R(s) le transformé du point s par la fonction de recalage R, - d(R(s), D) la distance euclidienne entre le transformé du point s par la fonction de recalage et l'ensemble de points destinations, noté D.
4. Procédé selon l'une des revendications 1 à 3, caractérisé en ce qu'il comprend une étape de filtrage direct dans laquelle chaque point s de l'ensemble de points sources S dont l'association avec un point d de l'ensemble de points destinations
D vérifie un critère de rejet est supprimé de sorte à obtenir un ensemble filtré de points source S'.
5. Procédé selon la revendication 4, caractérisé en ce que le critère de rejet est que la distance entre : o le transformé du point s de l'ensemble de points sources S par la fonction de recalage et o le point correspondant d de l'ensemble de points destinations D est supérieure à une valeur maximale prédéterminée.
6. Procédé selon l'une des revendications 4 ou 5, caractérisé en ce qu'il comprend en outre une étape de filtrage inverse dans laquelle chaque point d de l'ensemble de points destinations D dont l'association avec un point s de l'ensemble de points sources S vérifie un autre critère de rejet est supprimé de sorte à obtenir un ensemble filtré de points destinations D'..
7. Procédé selon la revendication 6, caractérisé en ce que l'autre critère de rejet est que la distance entre : o le transformé inverse du point d de l'ensemble de points destinations D par la fonction de recalage inverse et o le point correspondant s' de l'ensemble filtré de points sources S' est supérieure à la valeur maximale prédéterminée.
8. Système pour le recalage d'un ensemble de points source dans une image de référence représentant un tissu organique et d'un ensemble de points destinations dans une image courante représentant le tissu organique déformé, caractérisé en ce qu'il comprend des moyens pour la mise en œuvre du procédé selon l'une des revendications 1 à 7.
9. Produit programme d'ordinateur comprenant des instructions de code programme enregistré sur un support utilisable dans un ordinateur, caractérisé en ce qu'il comprend des instructions pour la mise en œuvre du procédé selon l'une des revendications 1 à 7.
PCT/EP2009/062829 2008-10-03 2009-10-02 Procédé pour le recalage d'un ensemble de points dans des images Ceased WO2010037852A1 (fr)

Priority Applications (2)

Application Number Priority Date Filing Date Title
EP09783693A EP2353145A1 (fr) 2008-10-03 2009-10-02 Procédé pour le recalage d'un ensemble de points dans des images
US13/122,116 US8724926B2 (en) 2008-10-03 2009-10-02 Method for registering a set of points in images

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
FR0856711 2008-10-03
FR0856711A FR2936889B1 (fr) 2008-10-03 2008-10-03 Procede pour le recalage d'un ensemble de points dans des images

Publications (1)

Publication Number Publication Date
WO2010037852A1 true WO2010037852A1 (fr) 2010-04-08

Family

ID=40552662

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/EP2009/062829 Ceased WO2010037852A1 (fr) 2008-10-03 2009-10-02 Procédé pour le recalage d'un ensemble de points dans des images

Country Status (4)

Country Link
US (1) US8724926B2 (fr)
EP (1) EP2353145A1 (fr)
FR (1) FR2936889B1 (fr)
WO (1) WO2010037852A1 (fr)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8724926B2 (en) 2008-10-03 2014-05-13 Universite Joseph Fourier—Grenoble 1 Method for registering a set of points in images

Families Citing this family (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
BR112013015181A2 (pt) * 2010-12-15 2016-09-13 Koninkl Philips Electronics Nv método e sistema
DE102011120937B4 (de) * 2011-12-14 2019-03-14 Carl Zeiss Meditec Ag Anordnung und Verfahren zur Registrierung von Gewebeverschiebungen
WO2013108182A1 (fr) * 2012-01-16 2013-07-25 Koninklijke Philips N.V. Appareil d'imagerie
PT2999392T (pt) * 2013-05-21 2019-09-12 Neurovision Imaging Inc Sistema e método de imagiologia de uma parte do corpo de um utilizador paciente
US20150356704A1 (en) * 2014-06-08 2015-12-10 Yeda Research And Development Co. Ltd. Good planar mappings and controlling singular values with semidefinite programming
US9582882B2 (en) 2015-03-02 2017-02-28 Nokia Technologies Oy Method and apparatus for image registration in the gradient domain
US9466108B1 (en) * 2015-03-31 2016-10-11 Nokia Technologies Oy Method and apparatus for multiple image registration in the gradient domain
CN107771055B (zh) 2015-06-19 2021-02-26 圣犹达医疗用品心脏病学部门有限公司 用于装置导航的电磁动态配准
CN107750148B (zh) * 2015-06-19 2021-01-01 圣犹达医疗用品心脏病学部门有限公司 阻抗位移及漂移检测和校正
US10467497B2 (en) * 2016-02-22 2019-11-05 Sony Corporation System and method for providing assistance in surgery in presence of tissue deformation
US11010630B2 (en) * 2017-04-27 2021-05-18 Washington University Systems and methods for detecting landmark pairs in images
US11918334B2 (en) 2018-11-07 2024-03-05 St Jude Medical International Holding, Sa.R.L. Impedance transformation model for estimating catheter locations
US12051198B2 (en) 2019-03-29 2024-07-30 Howmedica Osteonics Corp. Pre-morbid characterization of anatomical object using statistical shape modeling (SSM)
CN115334518B (zh) * 2022-06-30 2024-08-06 中国船舶重工集团公司第七0七研究所 一种面向特定海域的水下导航网络节点区域布局方法
CN117351052B (zh) * 2023-10-16 2024-09-20 北京科技大学顺德创新学院 一种基于特征一致性和空间一致性的点云精配准方法

Family Cites Families (3)

* 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
US20060036167A1 (en) * 2004-07-03 2006-02-16 Shina Systems Ltd. Vascular image processing
FR2936889B1 (fr) 2008-10-03 2012-09-28 Univ Grenoble 1 Procede pour le recalage d'un ensemble de points dans des images

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
ARMSPACH J-P ET AL: "3-D Deformable Image Registration: A Topology Preservation Scheme Based on Hierarchical Deformation Models and Interval Analysis Optimization", IEEE TRANSACTIONS ON IMAGE PROCESSING, IEEE SERVICE CENTER, PISCATAWAY, NJ, US, vol. 14, no. 5, 1 May 2005 (2005-05-01), pages 553 - 566, XP011130121, ISSN: 1057-7149 *
M. I. MILLER; L. YOUNES: "Group Actions, Homeomorphisms, and Matching: A General Framework", INTERNATIONAL JOURNAL OF COMPUTER VISION, vol. 41, no. 1-2, February 2001 (2001-02-01), Kluwer Academic Publishers Hingham, MA, USA, pages 61 - 84, XP002525069, ISSN: 0920-5691 *
MAREK BUCKI; YOHAN PAYAN: "Framework and Bio-Mechanical Model for a Per-Operative Image-Guided Neuronavigation Including 'Brain-Shift' Compensation", 23 October 2006 (2006-10-23), XP002525068, Retrieved from the Internet <URL:http://hal.archives-ouvertes.fr/docs/00/10/85/09/PDF/Brainshift_Santiago.pdf> [retrieved on 20090423] *
NOBLET V ET AL: "Retrospective evaluation of a topology preserving non-rigid registration method", MEDICAL IMAGE ANALYSIS, OXFORD UNIVERSITY PRESS, OXOFRD, GB, vol. 10, no. 3, 1 June 2006 (2006-06-01), pages 366 - 384, XP025154092, ISSN: 1361-8415, [retrieved on 20060601] *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8724926B2 (en) 2008-10-03 2014-05-13 Universite Joseph Fourier—Grenoble 1 Method for registering a set of points in images

Also Published As

Publication number Publication date
FR2936889B1 (fr) 2012-09-28
FR2936889A1 (fr) 2010-04-09
EP2353145A1 (fr) 2011-08-10
US8724926B2 (en) 2014-05-13
US20110176746A1 (en) 2011-07-21

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) 光コヒーレンストモグラフィのための動き補正および画像改善の方法および装置
EP2078222B1 (fr) Microscopie à ouverture synthétique interférométrique
Rajendran et al. Photoacoustic imaging aided with deep learning: a review
US9135682B2 (en) Image recovery from single shot digital hologram
JP6375619B2 (ja) オクルージョンを意識した明視野像からの３次元シーンの再構成
FR2918205A1 (fr) Procede et systeme de rendu volumique multivue
US20150371400A1 (en) Segmentation and identification of closed-contour features in images using graph theory and quasi-polar transform
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
EP2961316B1 (fr) Procede automatique de determination predictive de la position de la peau
Shenoy et al. R-SLAM: Optimizing eye tracking from rolling shutter video of the retina
WO2018104609A1 (fr) Procede de suivi de cible dans une sequence d&#39;images medicales, dispositif, equipement terminal et programmes d&#39;ordinateurs associes
EP0814430B1 (fr) Procédé de production d&#39;une séquence restaurée d&#39;images d&#39;un objet en mouvement à partir de mesures bruitées
WO2024126656A1 (fr) Dispositif d&#39;assistance à la planification d&#39;une intervention mini-invasive
WO2024213641A1 (fr) Methode mise en œuvre par ordinateur et dispositif associe pour modeliser une articulation d&#39;un patient
FR3055998A1 (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
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
Song et al. A multi-grid method for joint reconstruction of initial pressure and speed of sound in photoacoustic tomography
Zhang et al. 3D Fourier ptychographic microscopy based on the beam propagation method and time-reversal scheme
Nair et al. Convolutional Neural Networks Enable Direct Strain Estimation in Quasistatic Optical Coherence Elastography
FR3143313A1 (fr) Dispositif d’assistance à la planification d’une intervention mini-invasive sur un os

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: 09783693

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 13122116

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

WWE Wipo information: entry into national phase

Ref document number: 2009783693

Country of ref document: EP