WO2006048419A1 - Procede d'estimation du mouvement et/ou de la deformation d'au moins un objet contenu dans un ensemble d'images 3d, et procede de recalage tridimensionnel correspondant - Google Patents

Procede d'estimation du mouvement et/ou de la deformation d'au moins un objet contenu dans un ensemble d'images 3d, et procede de recalage tridimensionnel correspondant Download PDF

Info

Publication number
WO2006048419A1
WO2006048419A1 PCT/EP2005/055690 EP2005055690W WO2006048419A1 WO 2006048419 A1 WO2006048419 A1 WO 2006048419A1 EP 2005055690 W EP2005055690 W EP 2005055690W WO 2006048419 A1 WO2006048419 A1 WO 2006048419A1
Authority
WO
WIPO (PCT)
Prior art keywords
image
current
images
resolution
entities
Prior art date
Application number
PCT/EP2005/055690
Other languages
English (en)
Inventor
Antoine Simon
Mireille Garreau
Original Assignee
Universite De Rennes 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 Universite De Rennes 1 filed Critical Universite De Rennes 1
Publication of WO2006048419A1 publication Critical patent/WO2006048419A1/fr

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/20—Analysis of motion
    • G06T7/246—Analysis of motion using feature-based methods, e.g. the tracking of corners or segments
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR 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

Definitions

  • the field of the invention is that of 3D images. More specifically, the invention relates to a technique for estimating the movement and / or deformation of at least one object contained in each 3D image of a set of 3D images acquired at different times and / or from several imaging modalities and / or with an imaging modality and an atlas (pre-existing model).
  • the technique of the invention is for example implemented in multi-bar scanner imaging (MSCT), a new technological advance in scanner imaging and non-invasive imaging, to acquire images of high spatial and temporal resolutions.
  • MSCT multi-bar scanner imaging
  • Such an implementation of the invention in multi-strip scanner imaging has great advantages over other modes of cardiac imaging, in the sense that the scanner modality is a non-invasive imaging technique, allowing lower cost, in a limited time of examination while leaving a comfort to the patient, to visualize all the cardiac structures of interest in a single acquisition.
  • the invention can also be implemented with other imaging modalities, such as MRI (Magnetic Resonance Imaging) and echocardiography, since these techniques (still evolving) can provide a sufficient spatial and temporal resolution.
  • imaging modalities such as MRI (Magnetic Resonance Imaging) and echocardiography, since these techniques (still evolving) can provide a sufficient spatial and temporal resolution.
  • the set of 3D images it is for example: either of a set comprising two 3D images acquired at two different instants and / or with two different imaging modalities; a temporal sequence of 3D images comprising at least three 3D images acquired at successive instants within the same time range or of several non-adjacent time ranges, and with the same imaging modality or several modalities of imaging.
  • the invention applies in particular, but not exclusively, to the extraction and characterization of cardiac motion in multi-bar scanner imaging.
  • cardiac motion For example, it is a question of estimating the cardiac motion of Myocardial structures, for the moment limited to cavities, but can extend to the muscular wall with detection of both the endocardium and the epicardium.
  • This application is based on the estimation of three-dimensional movements (3D) of complex and deformable shapes (cardiac cavities), subjected to movements of a non-rigid nature with the objective of extracting quantitative information of cardiac dynamics (ejection fraction , volume variations, movement amplitudes by region).
  • the 3D motion estimation is carried out here at the beginning of a sequence of CT volumes reconstructed by post-synchronization on the ECG at different instants of the cardiac cycle, and a semi-interactive segmentation of the left cavities.
  • the analysis of these scanner images is subject to new constraints such as, in particular, the presence of movement artifacts in the images, the presence of fine structures accessible by the good spatial resolution of the images and large volumes of data. treat.
  • Atherosclerosis the main arterial disease characterized by arterial parietal lesions or atheroma plaques
  • Atherosclerosis the main arterial disease characterized by arterial parietal lesions or atheroma plaques
  • ischemic phenomena that result in alterations of myocardial dynamics.
  • the diagnosis is based on the state of the coronary circulation but also on the analysis of global and local ventricular kinetics.
  • MSCT multi-bar systems
  • the technique of the invention can be applied to the estimation of deformations of other structures.
  • the registration three-dimensional requires a prior estimation of the movement and / or deformation of one or more objects between the two 3D images considered.
  • the latter are acquired at two different instants and / or with two different imaging modalities and / or with an imaging modality and an atlas (pre-existing model).
  • atlas pre-existing model
  • the technique of the invention because it makes it possible to obtain correspondences between two objects inscribed in three-dimensional volumes, can be generalized to other applications.
  • Three-dimensional acquisition while particularly relevant for medical imaging, is also used in other areas such as marine imagery (sonar), satellite imagery.
  • Other fields such as the manufacture of materials can also benefit from the present invention, for the measurement of deformations of flexible materials or for the evaluation of resistance of materials under the effect of pressure for example.
  • the deformable models In methods based on deformable models, the deformable models, geometrical, parametric [Bardinet '96, Chen '94] or non-parametric models [Choi '01, Lin '00, Benayoun '95], after having been initialized, undergo deformations until an optimal match between model and surface is achieved.
  • parametric models the evolution of the parameters over time makes it possible to have a condensed description of the movement.
  • non-parametric models most often based on mass-spring models and / or finite elements with differential stress resolution, the temporal tracking of surface points makes it possible to obtain the estimation of their movement. The main defect of these geometric models remains the difficulty to model complex surfaces.
  • the optical flow methods are based on three constraints: the conservation of the density of the points during the movement, the absence of divergence and the smoothing of the movement. Once these three constraints are respected, the minimization of a penalty function makes it possible to obtain a dense field of displacement defined over the entire volume [Song '91, Coatrieux '94, Gorce '97].
  • the application of this type of method in this case encounters two difficulties: on the one hand, the size of the databases used makes this type of calculation extremely expensive in terms of machine resources.
  • the conservation constraint of the density of points is not a priori respected.
  • the duration of acquisition (the data can be acquired over several cardiac cycles) as well as the acquisition method itself (dynamic visualization of the concentration of contrast medium, subject to diffusion) and that the principle of reconstruction of the images by post-synchronization, makes that one point can change density during acquisition.
  • this intensity conservation constraint may be somewhat relaxed in this type of method, the computation of a dense displacement field over the entire volume studied is difficult to envisage on a large volume of data to be processed.
  • Shape tracking, or mapping methods [Amini '92, Kambhamettu '03, Ruan '94, Shi '0O], are based on the selection of entities and their tracking from a matching measure established on several descriptive parameters.
  • the technique of the present invention lies in this third family of motion estimation methods, namely that of the methods for mapping descriptive entities, between the different 3D images taken at different moments of observation of the object. and / or with different imaging modalities.
  • the matching methods offer advantages over the methods of the other two families mentioned above.
  • the matching methods allow, by the selection of fine entities, an accurate description of a complex surface.
  • the larger the number of entities the longer the computation time to find the optimal matches.
  • the matching methods can be limited to significant data (the myocardial wall in the case of the heart). However, they require having access to these data beforehand by segmentation tools.
  • mapping has the major disadvantage of requiring the segmentation of all the 3D images studied.
  • the segmentation tools are not effective in all cases. They do not guarantee the temporal coherence of the extracted surfaces, since the 3D images are segmented independently of each other, with a real risk of changing from one image to another of the acquisition conditions, the structures This is therefore a source of error for the estimation of the movement and / or the deformation.
  • the segmentation of a 3D image remains a difficult and tedious step of treatment, especially as it is often the case if it is conducted in an interactive or semi-interactive way (that is to say, if it requires intervention of the user to fix certain parameters of the algorithm for example).
  • the invention particularly aims to overcome this major drawback of known methods of matching.
  • one of the objectives of the present invention in at least one embodiment, is to provide a technique for estimating the movement and / or deformation of at least one object contained in each 3D image of a set of 3D images, this technique requiring the segmentation only of a 3D image of this set.
  • the invention also aims, in at least one embodiment, to provide such a technique to reduce the calculation time.
  • Another object of the invention in at least one embodiment, is to provide such a technique which, not only is not penalized by a great topological complexity of the objects studied, but on the contrary takes advantage of this complexity.
  • 3D of a set of 3D images comprising at least two 3D images acquired at different times and / or from several imaging modalities and / or with an imaging modality and an atlas.
  • the method according to the invention comprises the following steps: a) obtaining a 3D surface of said at least one object contained in one of the 3D images of said set, said initial image, by processing the image of initially, said 3D surface being represented in the form of a mesh comprising nodes; b) selection of the starting image as the current image; c) determining optimal matches between entities of the 3D surface obtained for the current 3D image, and entities of a next 3D image of said set, so as to estimate a field of displacement vectors and / or deformation associated with entities of the 3D surface and representing the movement and / or deformation of said at least one object between the current and next images.
  • the invention holds its original originality because it is based on a mapping of an object surface (obtained for the current image, acquired for example at time t) and a volume (which is the following 3D image, acquired for example at a time following t + 1).
  • This surface / volume mapping gives access to calculated displacement and / or deformation vectors for each descriptive entity of the object surface.
  • the invention thus differs from the known methods of matching by the fact that a single 3D image (the starting 3D image) requires a segmentation-type processing, making it possible to extract the 3D surface of the object or objects. objects of interest.
  • the method of the invention only requires the segmentation of one of the 3D images (volumes) of all the 3D images considered.
  • the fact of having only one segmentation to realize for a single 3D image is here an important asset.
  • the known methods of matching, as well as many known methods based on a deformable model require the pre-segmentation of all the processed 3D images (volumes).
  • the developed method makes it possible to take advantage of the topological complexity of the objects studied.
  • the complexity of the studied surfaces will always increase. Unlike methods based on deformable models that are problematic in the case of too complex shapes, the method developed relies on this complexity to extract the matches accurately.
  • the main constraint of the method of the invention is related to the amplitude of the extracted matches. This should not be too important, as this could jeopardize the validity of the results. However, in the case of the applications concerned, this limit has a very moderate impact. In the case of motion extraction, advances in acquisition methods make it possible to obtain higher and higher temporal resolutions, thereby reducing the amplitude of the movements between each instant considered.
  • said set of 3D images comprises two 3D images acquired at two different instants and / or with two different imaging modalities and / or with an imaging modality and an atlas.
  • said set of 3D images is a sequence of 3D images comprising at least three 3D images acquired at successive instants within the same time range or several non-time ranges. adjacent, and with the same imaging modality or multiple imaging modalities.
  • said starting image is the first 3D image of said 3D image sequence
  • the method further comprises the following steps: d) obtaining a 3D surface of said at least one object contained in the next 3D image, from the 3D surface obtained for the current 3D image and from the displacement and / or deformation vector field; e) iteration of steps c) and d) until all the 3D images of the sequence have been processed, the current image used during each new iteration being the next image used during the previous iteration.
  • the invention leads to the coupled detection of displacement and / or deformation vectors associated with the surface of the object between the two instants t and t + 1 (result of step c) and the surface of the object at time t + 1 (result of step d).
  • the starting 3D image requires segmentation type processing.
  • Other images 3D of the analyzed sequence do not undergo this treatment, the search for matches then being directly between the 3D surface of the reconstructed object at a time and the entire image volume at the next instant.
  • said starting image is different from the first 3D image of said 3D image sequence
  • the method further comprises the following steps: d) obtaining a 3D surface of said at least one object contained in the next 3D image, from the 3D surface obtained for the current 3D image and from the displacement and / or deformation vector field; e ') iterating steps c) and d), by traversing said sequence of 3D images in a first direction and until all the 3D images from the initial image to the last image of the sequence have been used, the current image used during each new iteration being the following image, in the first sense, used during the previous iteration; e ") iterating steps c) and d), by traversing said sequence of 3D images in a second direction and until all the 3D images from the starting image to the first image of the sequence have been used, the current image used during each new iteration being the following image, according to the second direction, used during the preceding iteration.
  • the entities of a given 3D image are voxels of said given 3D image.
  • said step c) comprises a step of minimizing a global energy consisting of a sum of local energies, each local energy to compare an entity of the 3D surface, obtained for the current 3D image, and an entity of the next 3D image.
  • each local energy is composed of a weighted sum of: a local correspondence term E 1 (Cd), expressing a coherence measurement of a pairing between an entity of the 3D surface, obtained for the current 3D image, and an entity of the next 3D image; a regularization term E 2 (T), expressing spatio ⁇ temporal smoothing constraints of the estimated movement and / or deformation, said constraints being such that displacement vectors associated with entities close to the 3D surface must be close to each other; amplitude and / or direction.
  • said local correspondence term E 1 (Cd) is a function of a term expressing a topological correspondence between the entity of the 3D surface, obtained for the current 3D image, and the entity of the following 3D image.
  • said local correspondence term E 1 (Cd) is composed of a weighted sum of said term expressing a topological correspondence and at least one of the terms belonging to the group comprising: a term expressing a probability that the entity of the following 3D image belongs to a contour element, of said at least one object contained in said next 3D image; a term expressing a distance between the entity of the 3D surface, obtained for the current 3D image, and the entity of the next 3D image; a term related to gray level densities.
  • contour element can be either previously detected or searched in a specific manner (if necessary) when calculating the energy term expressing a probability and involving the knowledge of this contour element.
  • said global energy minimization step is based on a statistical formulation of an estimation process according to a Markov model, such as the movement and / or the deformation to be estimated.
  • a Markov model such as the movement and / or the deformation to be estimated.
  • the 3D surface obtained for the current 3D image, during step a) or d), is considered as a 3D maximum resolution surface.
  • the next 3D image included in said set of 3D images is considered a 3D image of maximum resolution.
  • said step c) is performed in a multi-resolution mode and comprises the following successive steps: cl) obtaining on the one hand at least one lower resolution 3D surface for the current 3D image, starting from the surface 3D maximum resolution for the current 3D image, and secondly at least one subsequent 3D image of lower resolution, from the next maximum resolution 3D image; c2) choosing the lowest resolution as the current resolution, and initializing matches to the current resolution; c3) for the current resolution, determining optimal matches between entities of the 3D surface having said current resolution for the current 3D image, and entities of the next 3D image also having said current resolution; c4) if the current resolution is not equal to the maximum resolution:
  • step c3 * propagation of the results of the execution of step c3), the current resolution to a next higher resolution, so as to allow a resetting of the matches;
  • step c3 * return to step c3), taking as new current resolution the next higher resolution; otherwise stop of step c).
  • said step a) comprises the following steps: segmentation of the initial image, so as to extract contours and points inside said at least one object contained in the initial image; reconstruction of a 3D surface, represented in the form of a mesh comprising nodes, from the contours and interior points obtained during the segmentation step.
  • step c2) the initialization of the correspondences at the lowest resolution is carried out by null movements, each entity of the surface
  • step c4) the resetting of the correspondences is obtained by interpolation, over the entire mesh of the 3D surface having said next higher resolution, optimal matches determined in step c3).
  • said step a) further comprises the following step: pre-treatment of the reconstructed 3D surface, making it possible to modify, if necessary, said mesh so that the length of each of the edges constituting a mesh is equal to one distances between voxels of 3D images.
  • the imaging modalities belong to the group comprising: scanner imaging; magnetic resonance imaging; ultrasound imaging; nuclear imaging. This list is not exhaustive.
  • the method according to the invention is advantageously applied in an imaging field belonging to the group comprising: medical imaging, marine imaging, satellite imagery and imaging of materials. This list is not exhaustive.
  • the invention also relates to a computer program product comprising program code instructions for executing the steps of the method according to the invention, when said program is executed on a computer.
  • the invention also relates to a storage medium, possibly completely or partially removable, readable by a computer, storing a set of instructions executable by said computer to implement the method according to the invention.
  • the invention also relates to a method for three-dimensional registration of at least one object contained in a 3D image, comprising the following steps: estimation (41) of the movement and / or of the deformation undergone (s), between a first and a second 3D images acquired at two different instants and / or from two different imaging modalities and / or with an imaging modality and an atlas, by at least one object contained in said first and second 3D images; application (42) to said at least one object contained in the second image of a function correcting the movement and / or the deformation estimated during the estimation step.
  • said estimation step (41) is implemented according to the estimation method of the invention.
  • FIG. 1 presents a flowchart of a particular embodiment of the motion estimation and / or deformation method according to the invention
  • Figure 2 details an embodiment of the step referenced 11 in Figure 1
  • Figure 3 details an embodiment of the step referenced 13 in Figure 1
  • FIG. 4 presents a flowchart of a particular embodiment of the three-dimensional registration method according to the invention. 6. Description of an embodiment of the invention
  • the set of 3D images considered is a temporal sequence of 3D images acquired within the same time range and with the same imaging modality (for example in multi-strip scanner imaging (MSCT)). It is clear, however, that the invention can also be implemented with 3D images acquired in non-adjacent time slots and / or with several imaging modalities. Moreover, the invention can also be applied if the set of 3D images comprises only two images (see hereinafter the description of the three-dimensional registration method, in relation with FIG. 4).
  • MSCT multi-strip scanner imaging
  • the method according to the invention allows the spatio / temporal detection and the follow-up of a movement and / or a deformation of object, by a mapping principle of a surface represented as a mesh and a 3D image (also called 3D volume).
  • the starting data consists of a sequence of 3D images representing an object (or set of objects) in its (or their) environment, this object being subjected, between each 3D image, to a movement and / or to non-rigid deformation and moderate amplitudes.
  • the general principle of the invention consists in extracting the local correspondences between the initial mesh (3D surface) describing the considered object and the next 3D image, corresponding to the first deformation of this object.
  • the segmentation step which is crucial and difficult in image analysis, is therefore only necessary here for one (preferably the first) 3D image of the analyzed sequence.
  • the 3D surface of the object contained in the following 3D image is obtained, not by segmenting this next image, but by applying displacement and / or deformation vectors, resulting from the extraction of the correspondences, to the 3D surface of the object contained in the current 3D image.
  • This process is applied iteratively in order to extract the correspondences between the different 3D images of the sequence, so as to describe all the movements and / or deformations of the object. It allows to extract, simultaneously, the position of the deformed object in the analyzed volume (next 3D image) on the one hand, the amplitude and the direction of displacements and / or deformations of the object of other on the other hand, estimated locally. As illustrated in FIG.
  • the method according to the invention for estimating motion and / or deformation comprises the following steps: in a first step 11, the initial image is chosen (the first 3D image of the sequence in this example) as the current image; in a second step 12, this current image (starting image) is subjected to a treatment making it possible to obtain a 3D suriàce of the object that it contains.
  • This 3D surface is represented in the form of a mesh comprising nodes; in a third step 13, optimum matches are established between entities of the 3D surface obtained for the current 3D image, and entities of a next 3D image within the sequence, so as to estimate a vector field displacement and / or deformation associated with the entities of the 3D surface and representing the movement and / or deformation of the object between the current and next images; in a fourth step 14, a 3D surface of the object contained in the next 3D image is obtained from the 3D surface obtained for the current 3D image and from the displacement and / or deformation vector field; in a fifth step 15, it is determined whether the aforementioned next image is the last image of the sequence, that is to say if all the 3D images of the sequence have been processed; if so, we move on to a sixth end step 16; in the negative, a seventh step 17 consisting in taking the next image mentioned above as a new current image (with the 3D surface resulting from the fourth step 14) and the following image of the following one as the
  • the second step 12 is performed with a starting image that is not the first 3D image of the sequence. In this case, it is advisable to treat: all the images situated before the starting image, within the sequence (the current image used during each new iteration being the following image, according to a first direction of travel of the sequence, used during the previous iteration ); and - all the images located after the initial image, within the sequence (the current image used during each new iteration being the following image, according to a second direction of travel of the sequence, used during the previous iteration).
  • the second step 12 comprises, for example, the following steps: segmentation (121) of the initial image, so as to obtain the contours and interior points of the object contained in the initial image. This segmentation step is for example performed by combining a region growth approach and contour information extracted by a gradient operator (see [Guillaume'03]). It allows to provide as a result a set of points
  • reconstruction (122) of a 3D surface represented in the form of a mesh comprising nodes, from the points identified as belonging to the object during the segmentation step.
  • This reconstruction step is for example carried out by submitting the set of 3D points retained to the classic algorithm of Marching Cubes (creation of local topologies, represented as nodes and links between these nodes, thanks to the application of masks characteristics) ; preprocessing (123) of the reconstructed 3D surface, making it possible to modify, if necessary, the mesh so that the mesh size of the 3D surface is regular and close to the dimension of each of the voxels of the 3D images.
  • the preprocessing step (123) makes it possible to ensure that the edges constituting a mesh (a mesh being delimited by three nodes) are of length equal to the inter-voxel distances (they take their values in the set ⁇ 1 , V2, V3 ⁇
  • This pretreatment is ensured by a passage in integer values of the coordinates of all the nodes (approximation to the nearest integer), a merger of nodes which can be rendered confused by the preceding stage, a re-triangularization of the mesh (decomposition or fusion of meshes not having three nodes)
  • the third step 13 for determining optimal matches between entities of the 3D surface obtained for the current 3D image, and entities of the next 3D image is now described in more detail.
  • mapping method In order to establish optimal matches between a surface and a complete volume, the present invention provides a mapping method. This type of method requires considering three steps:
  • mapping method the minimization of a global energy consisting of the sum of the local energies.
  • One of the original features of the mapping method according to the invention is that it establishes correspondences between entities not belonging to the same domain. Indeed, the mapping is carried out between entities of a pre-extracted mesh (defined by nodes and cells (or meshes)) and entities of a complete volume (composed of voxels, three-dimensional elementary volumes). It should be noted that these two sets are very dissimilar: the mesh brings more knowledge
  • the entities used for the description of the mesh are the nodes that constitute it.
  • the entities to match with these nodes are the voxels of the volume.
  • the local energy must make it possible to compare two entities, each belonging to its respective set (nodes on one side, voxels on the other) while taking into account the constraints established on the neighborhood correspondences of each entity.
  • the probability that the voxel considered belongs to the studied boundary is evaluated by the application of a contour detector (which may be, for example, Canny filtering) on the analyzed image volume.
  • a contour detector which may be, for example, Canny filtering
  • a voxel detected as belonging to the contour is privileged to a voxel detected as not belonging to a boundary.
  • Other types of contour detectors may be retained in this phase, applied globally on the volume or locally.
  • a topological correspondence term is defined and allows an evaluation of the configuration compatibility between the node and the voxel considered and their respective neighborhood (the neighborhood of a node of the mesh is defined by the nodes having a mesh in common with this central node, and the neighborhood of a voxel, by the voxels verifying a first connection with the latter).
  • This topological correspondence can be based on various descriptive criteria of the local topology.
  • the mean and Gaussian curvatures are good descriptive characteristics.
  • the search for the local topology of the voxel studied is carried out.
  • Another technique for describing the local topology proposed in this invention consists in transposing the topology of the neighborhood of the studied node to the level of the candidate voxel.
  • This topology is described in terms of differences in the coordinates of the neighbors relative to the coordinates of the studied node.
  • the values of the voxels corresponding to this topology are summed.
  • the higher the sum the closer the local topologies of the two entities are.
  • the defined neighborhood can be extended to a larger number of meshes.
  • the Euclidean distance between node considered and candidate voxel is calculated. Short-distance connections are preferred.
  • a term of energy could be defined according to the densitometric information of the images from the moment when a correspondence between intensities of level of gray can be established between the two volumes considered ( 1st volume from which is defined the 3D surface and 2 e " 16 volume in which is sought the correspondence.)
  • the regularization term E 2 [J 1 ) makes it possible to avoid obtaining aberrant correspondences (for example the movements of two neighboring nodes that intersect) . She is based on the fact that the correspondences of neighboring nodes must generate close movements, whether in amplitude or in direction. Different implementations of this energy can be envisaged (knowing that the energy finally retained can also be a combination of different formulations).
  • a first formulation may correspond to the difference between the evaluated displacement and the average displacement of the neighbors of the considered node.
  • Other formulations may also be suitable.
  • a second formulation is proposed in the context of the Markovian formulation (see detailed discussion below), making it possible to maintain the homogeneity of the pairings.
  • a Markovian formulation is used, for example.
  • the Markov model makes it possible to formulate an estimation problem conditional on observations, in the form of a global model described by local specifications (introducing the definition of neighborhood).
  • the sites of the field considered are given by all the nodes representative of the mesh of the surface.
  • the labels of these sites or labels searched for (estimates f t ) are the corresponding voxels of these nodes.
  • the neighborhood of a site (or node) is composed of nodes having a mesh in common with the studied node.
  • T is a constant value called temperature and U (f) is a global energy function
  • the first term is expressed as a quadratic error so that the minimization of this energy component favors the estimated configurations close to the observations. It is calculated from the local matching measure (or distance) previously developed, on all the sites of the field:
  • ⁇ 2 (/) ⁇ V ⁇ f, > f j ) > or ⁇ , is the first connectivity neighborhood associated with m, site i (7)
  • v (f t , f j ) is the associated interaction potential to a neighborhood clique (a clique being represented here by two sites) with:
  • ⁇ (.) is a delta Kronecker function that takes the value 1 if the labels are close and the value 0 otherwise.
  • a click defined on a neighborhood consists of a subset of neighborhood elements, containing from 1 to N elements (the site included), N being the number of neighborhood elements. These elements constituting the clique must check a connectivity (point or mesh edge in common within the scope of this invention).
  • the minimization of the global energy U (f, d) makes it possible to find the optimal correspondences and thus the estimate of deformation sought for the whole object.
  • the different embodiments / can be obtained according to several optimization algorithms.
  • One method to solve this problem, in the context of the Markov fields, is the random search combined with a simulated annealing technique.
  • the stochastic relaxation algorithm of Metropolis is for example chosen, combined with a simulated annealing, to better ensure convergence to a global minimum, compared to a deterministic relaxation algorithm.
  • Other optimization algorithms are possible here.
  • An example of an optimization algorithm is written as follows:
  • Number_to_T Number_to_T + 1 Calculate AU (f, d)
  • Proba (exp (- AU (f, d) / T)) If ⁇ Proba> random (0, l) ⁇
  • T: T-step_T (decrease in T)
  • nodes For each iteration (at a given temperature), several nodes are chosen (randomly) and processed, until a sufficient number of nodes have been considered (number of tests at T) or until a given number of successes
  • T 0 the starting temperature
  • Not T the temperature decrease step
  • T max the maximum number of successes, at a given temperature, with AU (f, d) (change in overall energy) positive (that is, causing an increase in overall energy) .
  • the stopping criterion can be defined in different ways: either the temperature T must reach a threshold (which can be the temperature 0), or the number of correspondence modifications at a given temperature must be less than a threshold.
  • the selection of the studied node may either be completely random or guided by exhaustivity criteria (for example preventing the same node from being selected twice for a given temperature).
  • the selection of the candidate voxel is carried out by choosing (either randomly or according to criteria of exhaustivity or non-repetition) a voxel neighbor of the current corresponding voxel.
  • This neighborhood can be defined as being a first-order neighborhood (the two voxels then have a face or a stop in common), or of higher order (for example for order two, the two voxels have a face or a stop or a vertex in common) .
  • the third step 13 of determining optimal matches is applied in a multi-resolution mode.
  • the 3D surface obtained for the current 3D image, during the second step 12 or the fourth step 14 is considered as a 3D maximum resolution surface.
  • the next 3D image included in the set of 3D images is considered a 3D image of maximum resolution.
  • Multi-resolution mode requires the ability to go from lower resolution to lower resolution, whether for original data volumes (ie original 3D images of the source movie) or for 3D meshes (that is to say the 3D surfaces).
  • original data volumes ie original 3D images of the source movie
  • 3D meshes that is to say the 3D surfaces.
  • the transition from a given resolution representation to a lower resolution representation is for example achieved by a smoothing (which may be a Gaussian filtering) followed by a subsampling applied to the image volume.
  • the transition from a given resolution representation to a lower resolution representation is performed for example by smoothing (which may be Gaussian filtering) followed by subsampling of the coordinates of the nodes, then an operation of merge of the confused nodes and of re-triangulation of the mesh (so that all the meshes consist of three nodes).
  • smoothing which may be Gaussian filtering
  • the change of resolution at t is applied to the surface mesh, while the change of resolution at t + 1 is performed on the original volume (image of levels of Grey).
  • the change of resolution at t + k is applied to the surface mesh (resulting from the correspondence estimate made at the previous instants and having undergone the pre- already mentioned treatment (regularization and re-triangulation)), while the change of resolution at t + k + 1 is performed on the original volume.
  • the aforementioned third step 13 comprises the following successive steps: in a step referenced 131, obtaining on the one hand at least one 3D surface of lower resolution (for the 3D image current, and on the other hand at least one subsequent 3D image of lower resolution (see the detailed discussion above giving an example of a transition from a given resolution to a lower resolution, whether for the original 3D images of the source sequence or for 3D surfaces) ...
  • This step 133 is for example performed according to the example of optimization algorithm described above (Metropolis algorithm combined with simulated annealing), in the context of a Markovian formulation of the estimation problem; in a step referenced 134, it is determined whether the current resolution is equal to the maximum resolution; if so, the third step 13 is complete (and proceed to the next step 14, see Figure 1); in the negative, we return to the step referenced 133, after execution of a step referenced 135 consisting of: * propagate the results of the execution of the referenced step 133 (that is to say the vector fields displacement and / or deformation estimated for the nodes of the 3D surface), from the current resolution to a next higher resolution, so as to allow a resetting of the correspondences; * Take as new current resolution the next higher resolution.
  • a step referenced 135 consisting of: * propagate the results of the execution of the referenced step 133 (that is to say the vector fields displacement and / or deformation estimated for the nodes of the 3D surface), from the current resolution
  • the initialization of the correspondences at the lowest resolution is for example performed by null motions, each node (of the 3D surface obtained for the current image) being matched with the most voxel. near (within the next 3D image).
  • the propagation of the results, from a given resolution to the higher resolution, that is to say the resetting of the correspondences is for example ensured by interpolation, over the entire mesh of the 3D surface having this higher resolution, matches from the lower resolution (i.e. determined during the step referenced 133).
  • the node or nodes corresponding to the lower resolution are selected and their average displacement.
  • This average displacement is transposed to the higher resolution (by multiplication by the subsampling factor) and applied after passing in integer values (a corresponding voxel is estimated).
  • a field is thus available. displacement vectors and / or deformation estimated at each point of the 3D surface obtained for the current image.
  • This vector field represents the movement and / or the deformation of the object between the instant t of acquisition of the current 3D image and the time t + 1 of acquisition of the next 3D image.
  • These vectors located in 3D space, are characterized by an orientation and an amplitude. As already indicated above, these estimated vectors also make it possible to extract the contour points of the object at the instant t + 1 and thus to calculate the new area corresponding to the instant t + 1 (step referenced 14 on Figure 1). Then, this matching process is applied in the same way on all the successive instants of the sequence (between t + 1 and t + 2, then t + 2 and t + 3, etc.), starting from the new calculated surface and the image volume of the next moment. It thus provides a displacement vector field at each instant of the sequence, associated with the mesh nodes of the successively estimated surfaces.
  • such a registration method comprises: a step (41) for estimating the movement and / or the deformation undergone (s), between a first and a second 3D image acquired at two different times and / or from two different imaging modalities and / or from an imaging modality and an atlas (pre-existing model), by (at least) an object contained in these first and second 3D images; a step (42) for applying to the object contained in the second image a function correcting the movement and / or the deformation estimated during the estimation step.
  • step (41) itself comprises the steps referenced 11 to 13 of the estimation method discussed above in relation to FIGS. 1 to 3. These steps 11 to 13 are not so not described again.

Landscapes

  • Engineering & Computer Science (AREA)
  • Multimedia (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

L'invention concerne un procédé d'estimation du mouvement et/ou de la déformation d'au moins un objet contenu dans chaque image 3D d'un ensemble d'images 3D comprenant au moins deux images 3D acquises à différents instants et/ou à partir de plusieurs modalités d'imagerie et/ou avec une modalité d'imagerie et un atlas. Selon l'invention, ce procédé comprend les étapes suivantes : a) obtention (12) d'une surface 3D dudit au moins un objet contenu dans l'une des images 3D dudit ensemble, dite image de départ, par traitement de l'image de départ, ladite surface 3D étant représentée sous la forme d'un maillage comprenant des noeds ; b) choix (11) de l'image de départ comme image courante ; c) détermination (13) de correspondances optimales entre des entités de la surface 3D obtenue pour l'image 3D courante, et des entités d'une image 3D suivante dudit ensemble, de façon à estimer un champ de vecteurs de déplacement et/ou de déformation associés aux entités de la surface 3D et représentant le mouvement et/ou la déformation dudit au moins un objet entre les images courante et suivante.

Description

Procédé d'estimation du mouvement et/ou de la déformation d'au moins un objet contenu dans un ensemble d'images 3D, et procédé de recalage tridimensionnel correspondant. 1. Domaine de l'invention Le domaine de l'invention est celui des images 3D. Plus précisément, l'invention concerne une technique d'estimation du mouvement et/ou de la déformation d'au moins un objet contenu dans chaque image 3D d'un ensemble d'images 3D acquises à différents instants et/ou à partir de plusieurs modalités d'imagerie et/ou avec une modalité d'imagerie et un atlas (modèle préexistant). La technique de l'invention est par exemple mise en œuvre en imagerie scanner multi-barrette (MSCT), nouvelle avancée technologique en imagerie scanner et imagerie non invasive, permettant d'acquérir des images de hautes résolutions spatiales et temporelles. Une telle mise en œuvre de l'invention en imagerie scanner multi-barrette présente de gros avantages par rapport à d'autres modalités de l'imagerie cardiaque, dans le sens où la modalité scanner est une technique d'imagerie non invasive, permettant à moindre coût, dans un temps limité d'examen tout en laissant un confort au patient, de visualiser l'ensemble des structures cardiaques d'intérêt en une seule acquisition.
Il est clair cependant que l'invention peut également être mise en œuvre avec d'autres modalités d'imagerie, telles que l'IRM (Imagerie par Résonance Magnétique) et l'échocardiographie, dès lors que ces techniques (toujours en évolution) peuvent fournir une résolution spatio-temporelle suffisante.
En ce qui concerne l'ensemble d'images 3D, il s'agit par exemple : soit d'un ensemble comprenant deux images 3D acquises à deux instants différents et/ou avec deux modalités d'imagerie différentes ; soit d'une séquence temporelle d'images 3D comprenant au moins trois images 3D acquises à des instants successifs au sein d'une même plage temporelle ou de plusieurs plages temporelles non adjacentes, et avec une même modalité d'imagerie ou plusieurs modalités d'imagerie.
Dans le domaine médical, l'invention s'applique notamment, mais non exclusivement, à l'extraction et la caractérisation du mouvement cardiaque en imagerie scanner multi-barrette. Il s'agit par exemple d'estimer le mouvement cardiaque des structures du myocarde, pour l'instant limitées aux cavités, mais pouvant s'étendre à la paroi musculaire avec détection à la fois de l'endocarde et de l'épicarde. Cette application relève de l'estimation de mouvements tridimensionnels (3D) de formes complexes et déformables (cavités cardiaques), soumises à des mouvements de nature non-rigide avec pour objectif d'extraire des informations quantitatives de la dynamique cardiaque (fraction d'éjection, variations de volumes, amplitudes de mouvements par région). L'estimation de mouvement 3D est ici conduite en disposant au départ d'une séquence de volumes scanner reconstruits par post-synchronisation sur l'ECG à différents instants du cycle cardiaque, et d'une segmentation semi-interactive des cavités gauches. L'analyse de ces images scanner est soumise à de nouvelles contraintes telles que, en particulier, la présence d'artefacts de mouvement dans les images, la présence de structures fines accessibles de par la bonne résolution spatiale des images et de gros volumes de données à traiter.
Cette application médicale particulière présente un intérêt important. En effet, en terme de santé publique, les maladies cardio-vasculaires représentent aujourd'hui 40% des examens pathologiques réalisés et constituent la première cause de mortalité. L'athérosclérose, principale maladie artérielle caractérisée par des lésions pariétales artérielles ou plaques d'athérome, est à l'origine de phénomènes ischémiques qui se traduisent par des altérations de la dynamique du myocarde. En imagerie cardio- vasculaire, le diagnostic repose sur l'état de la circulation coronaire mais aussi sur l'analyse de la cinétique ventriculaire globale et locale. Au niveau des modalités d'observation et plus particulièrement, en imagerie scanner, l'apparition récente des systèmes multi-barrette (MSCT) (multi-détection, amélioration des résolutions spatiale et temporelle), permet aujourd'hui l'acquisition de structures en mouvement et ouvre une nouvelle voie en ce qui concerne l'imagerie non invasive cardiaque.
Hors du champ médical, la technique de l'invention peut être appliquée à l'estimation de déformations d'autres structures.
Elle peut aussi répondre aux problèmes de recalage tridimensionnel d'objets qui se posent en imagerie médicale, mais également dans d'autres domaines liés à la vision par ordinateur, quand les observations de scènes sont acquises à différents instants et/ou à partir de plusieurs modalités d'imagerie. En effet, dans ce cas, le recalage tridimensionnel nécessite une estimation préalable du mouvement et/ou de la déformation d'un ou plusieurs objets entre les deux images 3D considérées. Ces dernières sont acquises à deux instants différents et/ou avec deux modalités d'imagerie différentes et/ou avec une modalité d'imagerie et un atlas (modèle préexistant). Dans le dernier cas, on cherche à superposer une image réelle avec des données de type atlas
(atlas anatomique par exemple), qui permettent de partir d'une image modèle. Ces objets sont par exemple des formes tridimensionnelles déformables soumises à des mouvements non-rigides. C'est bien sûr le cas de nombreux organes du corps humain, pour lesquels le diagnostic ou la conduite de gestes interventionnels en thérapie ou chirurgie assistée par ordinateur nécessite l'acquisition et le traitement d'images 3D.
Mais le domaine de la chirurgie esthétique ou réparation tissulaire est aussi concerné.
La technique de l'invention, du fait qu'elle permet d'obtenir des correspondances entre deux objets inscrits dans des volumes tridimensionnels, peut être généralisée à d'autres applications. L'acquisition tridimensionnelle, si elle est particulièrement présente dans le cadre de l'imagerie médicale, est aussi utilisée dans d'autres domaines comme l'imagerie marine (sonar), l'imagerie satellitaire. D'autres domaines tels que la fabrication de matériaux peuvent aussi bénéficier de la présente invention, pour la mesure de déformations de matières souples ou pour l'évaluation de résistances de matériaux sous l'effet de pression par exemple. 2. Art antérieur
II existe de nombreuses techniques connues d'estimation de mouvement. Dans un souci de simplification, on limite ci-après la discussion aux techniques connues d'estimation de mouvement non rigide, type de mouvement dont sont dotées les structures cardiaques. En effet, la plupart des travaux méthodologiques connus relatifs à l'extraction de mouvement sont appliqués à l'analyse du mouvement cardiaque, du fait que le cœur est soumis à des déformations complexes (translation, torsion, compression, expansion, élongation, raccourcissement) qui relèvent de mouvements non-rigides et en particulier élastiques.
On résume ci-dessous ces travaux, toutes modalités d'acquisition d'images confondues (IRM taggué, IRM à contraste de phase, imagerie par ultrasons, ...), en retenant les avantages et inconvénients de chaque type de méthode, et en les situant par rapport à la technique de l'invention. Il est également fait mention de diverses publications, dont les références précises sont fournies ci-après dans une annexe intitulée « bibliographie ».
En vision par ordinateur, trois principales méthodes d'estimation de mouvement élastique existent et ont été développées dans la littérature :
- les méthodes basées sur les modèles déformables ;
- les méthodes de flot optique ; et
- les méthodes de suivi de forme, ou de mise en correspondance.
Dans les méthodes basées sur les modèles déformables, les modèles déformables, modèles géométriques, paramétriques [Bardinet '96, Chen '94] ou non paramétriques [Choi '01, Lin '00, Benayoun '95], après avoir été initialisés, subissent des déformations jusqu'à obtention d'une correspondance optimale entre modèle et surface. Dans le cas des modèles paramétriques, l'évolution des paramètres au cours du temps permet d'avoir une description condensée du mouvement. Dans le cas des modèles non paramétriques, le plus souvent basés sur des modèles de type masse-ressort et/ou éléments finis avec résolution de contraintes différentielles, le suivi temporel de points de surface permet d'obtenir l'estimation de leur mouvement. Le principal défaut de ces modèles géométriques reste la difficulté à modéliser des surfaces complexes. Dans le cas des volumes cardiaques acquis par scanner multi-barrette, la très haute qualité de la résolution spatiale fait que la surface peut être par endroit très découpée. Ces aspérités locales seront difficilement restituées par ce type de modélisation alors qu'elles peuvent être significatives du point de vue anatomique.
Les méthodes de flot optique se basent sur trois contraintes : la conservation de la densité des points lors du mouvement, l'absence de divergence et le lissage du mouvement. Une fois ces trois contraintes respectées, la minimisation d'une fonction de pénalité permet d'obtenir un champ dense de déplacement défini sur tout le volume [Song '91, Coatrieux '94, Gorce '97]. L'application de ce type de méthodes dans le présent cas se heurte à deux difficultés : d'une part, la taille des bases de données utilisées fait que ce type de calcul s'avérerait extrêmement coûteux en ressources machine. D'autre part, et surtout, la contrainte de conservation de la densité des points n'est pas a priori respectée. En effet, la durée d'acquisition (les données pouvant être acquises sur plusieurs cycles cardiaques) ainsi que la méthode d'acquisition en elle- même (visualisation dynamique de la concentration en produit de contraste, soumis à la diffusion) et que le principe de reconstruction des images par post-synchronisation, font qu'un point peut changer de densité au cours de l'acquisition. Bien que cette contrainte de conservation d'intensité puisse être quelque peu relâchée dans ce type de méthodes, le calcul d'un champ de déplacement dense sur tout le volume étudié est difficilement envisageable sur un gros volume de données à traiter.
Les méthodes de suivi de forme, ou de mise en correspondance [Amini '92, Kambhamettu '03, Ruan '94, Shi '0O], sont basées sur la sélection d'entités et leur suivi à partir d'une mesure de correspondance établie sur plusieurs paramètres descriptifs.
La technique de la présente invention se situe dans cette troisième famille de méthodes d'estimation de mouvement, à savoir celle des méthodes de mise en correspondance d'entités descriptives, entre les différentes images 3D prises à différents instants d'observation de l'objet et/ou avec différentes modalités d'imagerie. En effet, les méthodes de mise en correspondance offrent des avantages par rapport aux méthodes des deux autres familles précitées.
Tout d'abord, contrairement aux méthodes basées sur les modèles déformables, les méthodes de mise en correspondance permettent, par la sélection d'entités fines, une description précise d'une surface complexe. Cependant, plus le nombre d'entités est important, plus le temps de calcul pour trouver les correspondances optimales est important.
En outre, contrairement aux méthodes de flot optique, les méthodes de mise en correspondance permettent de se limiter aux données significatives (la paroi myocardique dans le cas du cœur). Cependant, elles imposent d'avoir préalablement accès à ces données par des outils de segmentation.
De plus, les méthodes connues de mise en correspondance présentent l'inconvénient majeur de nécessiter la segmentation de toutes les images 3D étudiées.
Or, les outils de segmentation, très nombreux, ne sont pas efficaces dans tous les cas. Ils ne garantissent pas la cohérence temporelle des surfaces extraites, puisque les images 3D sont segmentées indépendamment les unes des autres, avec un risque réel de modification d'une image à l'autre des conditions d'acquisition, des structures environnantes, de la qualité des images ,... Ceci est donc source d'erreur pour l'estimation du mouvement et/ou de la déformation.
D'autre part, il est pénalisant en termes de performances de répéter la phase de segmentation pour chaque image 3D. En effet, la segmentation d'une image 3D demeure une étape de traitement difficile et fastidieuse, surtout comme c'est souvent le cas si elle est menée de manière interactive ou semi-interactive (c'est-à-dire si elle nécessite l'intervention de l'utilisateur pour fixer certains paramètres de l'algorithme par exemple).
3. Objectifs de l'invention L'invention a notamment pour objectif de pallier cet inconvénient majeur des méthodes connues de mise en correspondance.
Plus précisément, l'un des objectifs de la présente invention, dans au moins un mode de réalisation, est de fournir une technique d'estimation du mouvement et/ou de la déformation d'au moins un objet contenu dans chaque image 3D d'un ensemble d'images 3D, cette technique ne nécessitant la segmentation que d'une image 3D de cet ensemble.
L'invention a également pour objectif, dans au moins un mode de réalisation, de fournir une telle technique permettant de réduire les temps de calcul.
Un autre objectif de l'invention, dans au moins un mode de réalisation, est de fournir une telle technique qui, non seulement n'est pas pénalisée par une grande complexité topologique des objets étudiés, mais au contraire tire parti de cette complexité.
4. Caractéristiques essentielles de l'invention
Ces différents objectifs, ainsi que d'autres qui apparaîtront par la suite, sont atteints selon l'invention à l'aide d'un procédé d'estimation du mouvement et/ou de la déformation d'au moins un objet contenu dans chaque image 3D d'un ensemble d'images 3D comprenant au moins deux images 3D acquises à différents instants et/ou à partir de plusieurs modalités d'imagerie et/ou avec une modalité d'imagerie et un atlas. Le procédé selon l'invention comprend les étapes suivantes : a) obtention d'une surface 3D dudit au moins un objet contenu dans l'une des images 3D dudit ensemble, dite image de départ, par traitement de l'image de départ, ladite surface 3D étant représentée sous la forme d'un maillage comprenant des nœuds ; b) choix de l'image de départ comme image courante ; c) détermination de correspondances optimales entre des entités de la surface 3D obtenue pour l'image 3D courante, et des entités d'une image 3D suivante dudit ensemble, de façon à estimer un champ de vecteurs de déplacement et/ou de déformation associés aux entités de la surface 3D et représentant le mouvement et/ou la déformation dudit au moins un objet entre les images courante et suivante. L'invention tient son originalité première du fait qu'elle est basée sur une mise en correspondance d'une surface d'objet (obtenue pour l'image courante, acquise par exemple à l'instant t) et d'un volume (qui est l'image 3D suivante, acquise par exemple à un instant suivant t+1). Cette mise en correspondance surface/volume donne accès à des vecteurs de déplacement et/ou de déformation calculés pour chaque entité descriptive de la surface d'objet.
L'invention se démarque donc des méthodes connues de mise en correspondance par le fait qu'une seule image 3D (l'image 3D de départ) nécessite un traitement de type segmentation, permettant d'extraire la surface 3D de l'objet ou des objets d'intérêt. En d'autres termes, le procédé de l'invention ne nécessite la segmentation que d'une des images 3D (volumes) de l'ensemble des images 3D considérées. Le fait de n'avoir qu'une segmentation à réaliser pour une seule image 3D est ici un atout important. En effet, les méthodes connues de mise en correspondance, ainsi que bon nombre de méthodes connues basées sur un modèle déformable, nécessitent la pré-segmentation de toutes les images 3D (volumes) traitées. Le procédé développé permet de tirer parti de la complexité topologique des objets étudiés. Dans un contexte où les progrès de l'acquisition tridimensionnelle permettent d'obtenir une résolution spatiale de plus en plus élevée, la complexité des surfaces étudiées ira toujours croissante. Contrairement aux méthodes basées sur les modèles déformables qui posent problème dans le cas de formes trop complexes, la méthode développée s'appuie sur cette complexité afin d'extraire les correspondances avec précision. II est à noter que la principale contrainte du procédé de l'invention est liée à l'amplitude des correspondances extraites. Celle-ci ne doit pas être trop importante au risque de compromettre la validité des résultats. Cependant, dans le cas des applications concernées, cette limite a un impact très modéré. Dans le cas de l'extraction de mouvement, les progrès des méthodes d'acquisition permettent d'obtenir des résolutions temporelles de plus en plus élevées, entraînant par là la diminution de l'amplitude des mouvements entre chaque instant considéré. Par exemple, pour l'extraction du mouvement cardiaque, les scanners MSCT permettent d'obtenir une résolution temporelle nettement suffisante pour l'utilisation du procédé selon l'invention. Dans un premier mode de réalisation avantageux de l'invention, ledit ensemble d'images 3D comprend deux images 3D acquises à deux instants différents et/ou avec deux modalités d'imagerie différentes et/ou avec une modalité d'imagerie et un atlas.
Dans un second mode de réalisation avantageux de l'invention, ledit ensemble d'images 3D est une séquence d'images 3D comprenant au moins trois images 3D acquises à des instants successifs au sein d'une même plage temporelle ou de plusieurs plages temporelles non adjacentes, et avec une même modalité d'imagerie ou plusieurs modalités d'imagerie.
De façon avantageuse, dans le cas du second mode de réalisation précité, ladite image de départ est la première image 3D de ladite séquence d'images 3D, et en ce que le procédé comprend en outre les étapes suivantes : d) obtention d'une surface 3D dudit au moins un objet contenu dans l'image 3D suivante, à partir de la surface 3D obtenue pour l'image 3D courante et du champ de vecteurs de déplacement et/ou de déformation ; e) itération des étapes c) et d) jusqu'à ce que toutes les images 3D de la séquence aient été traitées, l'image courante utilisée lors de chaque nouvelle itération étant l'image suivante utilisée lors de l'itération précédente.
Ainsi, dans le cas d'une séquence d'images 3D, l'invention conduit à la détection couplée des vecteurs de déplacement et/ou déformation associés à la surface de l'objet entre les deux instants t et t+1 (résultat de l'étape c) et de la surface de l'objet à l'instant t+1 (résultat de l'étape d). De cette façon, dans le cas d'une séquence d'images 3D, seule l'image 3D de départ nécessite un traitement de type segmentation. Les autres images 3D de la séquence analysée ne subissent pas ce traitement, la recherche de correspondances se faisant ensuite directement entre la surface 3D de l'objet reconstruite à un instant et le volume image entier à l'instant suivant.
Selon une variante avantageuse, ladite image de départ est différente de la première image 3D de ladite séquence d'images 3D, et en ce que le procédé comprend en outre les étapes suivantes : d') obtention d'une surface 3D dudit au moins un objet contenu dans l'image 3D suivante, à partir de la surface 3D obtenue pour l'image 3D courante et du champ de vecteurs de déplacement et/ou de déformation ; e') itération des étapes c) et d'), en parcourant ladite séquence d'images 3D selon un premier sens et jusqu'à ce que toutes les images 3D depuis l'image de départ jusqu'à la dernière image de la séquence aient été utilisées, l'image courante utilisée lors de chaque nouvelle itération étant l'image suivante, selon le premier sens, utilisée lors de l'itération précédente ; e") itération des étapes c) et d'), en parcourant ladite séquence d'images 3D selon un second sens et jusqu'à ce que toutes les images 3D depuis l'image de départ jusqu'à la première image de la séquence aient été utilisées, l'image courante utilisée lors de chaque nouvelle itération étant l'image suivante, selon le second sens, utilisée lors de l'itération précédente. Préférentiellement, les entités d'une surface 3D donnée sont des nœuds du maillage de ladite surface 3D donnée.
De façon préférentielle, les entités d'une image 3D donnée sont des voxels de ladite image 3D donnée.
Au lieu des noeuds et des voxels, des entités d'autre nature peuvent également être considérées, telle que « la maille » ou encore « un patch de surface (ou ensemble de mailles)» pour la surface à mettre en correspondance avec « un ensemble de voxels » sur le volume. Les voisinages de ces entités et les énergies locales associées à ces entités sont bien sûr dans ce cas à re-définir, tout en pouvant les construire sur le même type de descripteurs de base (position, courbure, ...). Avantageusement, ladite étape c) comprend une étape de minimisation d'une énergie globale constituée d'une somme d'énergies locales, chaque énergie locale permettant de comparer une entité de la surface 3D, obtenue pour l'image 3D courante, et une entité de l'image 3D suivante.
Selon une caractéristique avantageuse, chaque énergie locale est composée d'une somme pondérée de : - un terme de correspondance locale E1(Cd), exprimant une mesure de cohérence d'un appariement entre une entité de la surface 3D, obtenue pour l'image 3D courante, et une entité de l'image 3D suivante ; un terme de régularisation E2(T), exprimant des contraintes de lissage spatio¬ temporel du mouvement et/ou de la déformation estimée, lesdites contraintes étant telles que des vecteurs de déplacement associés à des entités voisines de la surface 3D doivent être proches en amplitude et/ou en direction. Avantageusement, ledit terme de correspondance locale E1(Cd) est fonction d'un terme exprimant une correspondance topologique entre l'entité de la surface 3D, obtenue pour l'image 3D courante, et l'entité de l'image 3D suivante. De façon avantageuse, ledit terme de correspondance locale E1(Cd) est composé d'une somme pondérée dudit terme exprimant une correspondance topologique et d'au moins un des termes appartenant au groupe comprenant : un terme exprimant une probabilité que l'entité de l'image 3D suivante appartienne à un élément de contour, dudit au moins un objet contenu dans ladite image 3D suivante ; un terme exprimant une distance entre l'entité de la surface 3D, obtenue pour l'image 3D courante, et l'entité de l'image 3D suivante ; un terme lié aux densités de niveaux de gris.
Il est à noter que l'élément de contour peut être soit préalablement détecté, soit recherché de manière ponctuelle (si besoin) au moment du calcul du terme d'énergie exprimant une probabilité et impliquant la connaissance de cet élément de contour.
Dans un mode de réalisation particulier de l'invention, ladite étape de minimisation d'une énergie globale est basée sur une formulation statistique d'un processus d'estimation selon un modèle de Markov, telle que le mouvement et/ou la déformation à estimer est considéré comme une réalisation f = {f, , i = \,...,NS} d'un champ aléatoire F = {Ft , i = l,..., NS} regroupant toutes les réalisations possibles, les NS sites du champ aléatoire F étant donnés par toutes les entités (par ex. nœuds) de la surface 3D, obtenue pour l'image 3D courante, représentatives du maillage de ladite surface 3D, des labels desdits sites, aussi appelés étiquettes recherchées, étant les entités (par ex. voxels) de l'image 3D suivante jugées correspondantes desdites entités de la surface 3D.
Cette formalisation du problème d'estimation dans le cadre statistique des champs markoviens 3D permet d'optimiser les calculs.
De façon préférentielle, la surface 3D obtenue pour l'image 3D courante, lors de l'étape a) ou d), est considérée comme une surface 3D de résolution maximale. En outre, l'image 3D suivante comprise dans ledit ensemble d'images 3D est considérée comme une image 3D de résolution maximale. Enfin, ladite étape c) est effectuée dans un mode multi-résolution et comprend les étapes successives suivantes : cl) obtention d'une part d'au moins une surface 3D de résolution inférieure pour l'image 3D courante, à partir de la surface 3D de résolution maximale pour l'image 3D courante, et d'autre part d'au moins une image 3D suivante de résolution inférieure, à partir de l'image 3D suivante de résolution maximale ; c2) choix de la plus basse des résolutions comme résolution courante, et initialisation des correspondances à la résolution courante ; c3) pour la résolution courante, détermination de correspondances optimales entre des entités de la surface 3D possédant ladite résolution courante pour l'image 3D courante, et des entités de l'image 3D suivante possédant également ladite résolution courante ; c4) si la résolution courante n'est pas égale à la résolution maximale :
* propagation de résultats de l'exécution de l'étape c3), de la résolution courante à une résolution supérieure suivante, de façon à permettre une réinitialisation des correspondances ;
* retour à l'étape c3), en prenant comme nouvelle résolution courante la résolution supérieure suivante ; sinon arrêt de l'étape c). Ainsi, l'intégration d'un processus de multi-résolution permet d'affiner l'estimation de manière guidée et progressive. On tient compte et on exploite les différents niveaux de détails topologiques des surfaces complexes, et on réduit le temps de calcul.
Avantageusement, ladite étape a) comprend les étapes suivantes : segmentation de l'image de départ, de façon à extraire des contours et points intérieurs audit au moins un objet contenu dans l'image de départ ; reconstruction d'une surface 3D, représentée sous la forme d'un maillage comprenant des nœuds, à partir des contours et points intérieurs obtenus lors de l'étape de segmentation.
De façon avantageuse, dans l'étape c2), l'initialisation des correspondances à la résolution la plus basse est réalisée par des mouvements nuls, chaque entité de la surface
3D obtenue pour l'image 3D courante étant mise en correspondance avec l'entité la plus proche dans l'image 3D suivante.
Avantageusement, dans l'étape c4), la réinitialisation des correspondances est obtenue par interpolation, sur tout le maillage de la surface 3D possédant ladite résolution supérieure suivante, des correspondances optimales déterminées lors de l'étape c3).
De façon avantageuse, ladite étape a) comprend en outre l'étape suivante : pré¬ traitement de la surface 3D reconstruite, permettant de modifier, si nécessaire, ledit maillage de façon que la longueur de chacune des arêtes constituant une maille soit égale à une des distances entre voxels des images 3D.
Selon une caractéristique avantageuse, les modalités d'imagerie appartiennent au groupe comprenant : l'imagerie scanner ; l'imagerie par résonance magnétique ; l'imagerie par ultrasons ; l'imagerie nucléaire. Cette liste n'est pas exhaustive.
Le procédé selon l'invention s'applique avantageusement dans un domaine d'imagerie appartenant au groupe comprenant : l'imagerie médicale, l'imagerie marine, l'imagerie satellitaire et l'imagerie des matériaux. Cette liste n'est pas exhaustive.
L'invention concerne également un produit programme d'ordinateur comprenant des instructions de code de programme pour l'exécution des étapes du procédé selon l'invention, lorsque ledit programme est exécuté sur un ordinateur. L'invention concerne aussi un moyen de stockage, éventuellement totalement ou partiellement amovible, lisible par un ordinateur, stockant un jeu d'instructions exécutables par ledit ordinateur pour mettre en œuvre le procédé selon l'invention.
L'invention concerne encore un procédé de recalage tridimensionnel d'au moins un objet contenu dans une image 3D, comprenant les étapes suivantes : estimation (41) du mouvement et/ou de la déformation subi(s), entre une première et une seconde images 3D acquises à deux instants différents et/ou à partir de deux modalités d'imagerie différentes et/ou avec une modalité d'imagerie et un atlas, par au moins un objet contenu dans lesdites première et seconde images 3D ; application (42) audit au moins un objet contenu dans la seconde image d'une fonction corrigeant le mouvement et/ou la déformation estimé lors de l'étape d'estimation. En outre, ladite étape d'estimation (41) est mise en œuvre selon le procédé d'estimation de l'invention.
5. Liste des figures
D'autres caractéristiques et avantages de l'invention apparaîtront à la lecture de la description suivante d'un mode de réalisation préférentiel de l'invention, donné à titre d'exemple indicatif et non limitatif, et des dessins annexés, dans lesquels : la figure 1 présente un organigramme d'un mode de réalisation particulier du procédé d'estimation de mouvement et/ou de déformation selon l'invention ; la figure 2 détaille un exemple de réalisation de l'étape référencée 11 sur la figure 1 ; la figure 3 détaille un exemple de réalisation de l'étape référencée 13 sur la figure 1 ; et la figure 4 présente un organigramme d'un mode de réalisation particulier du procédé de recalage tridimensionnel selon l'invention. 6. Description d'un mode de réalisation de l'invention
On présente ci-dessous, en relation avec les figures 1 à 3, un mode de réalisation particulier du procédé selon l'invention d'estimation du mouvement et/ou de la déformation d'un objet contenu dans chaque image 3D d'un ensemble d'images 3D.
Dans ce mode de réalisation particulier, l'ensemble d'images 3D considéré est une séquence temporelle d'images 3D acquises au sein d'une même plage temporelle et avec une même modalité d'imagerie (par exemple en imagerie scanner multi-barrette (MSCT)). Il est clair cependant que l'invention peut également être mise en œuvre avec des images 3D acquises dans des plages temporelles non adjacentes et/ou avec plusieurs modalités d'imagerie. Par ailleurs, l'invention peut aussi s'appliquer si l'ensemble d'images 3D ne comprend que deux images (voir ci-après la description du procédé de recalage tridimensionnel, en relation avec la figure 4).
Dans ce mode de réalisation particulier, à partir d'images 3D acquises à différents instants (séquence temporelle), le procédé selon l'invention permet la détection spatio/temporelle et le suivi d'un mouvement et/ou d'une déformation d'objet, par un principe de mise en correspondance d'une surface représentée sous forme de maillage et d'une image 3D (aussi appelée volume 3D). Les données de départ sont constituées d'une séquence d'images 3D représentant un objet (ou un ensemble d'objets) dans son (ou leur) environnement, cet objet étant soumis, entre chaque image 3D, à un mouvement et/ou à une déformation de type non-rigide et d'amplitudes modérées.
On rappelle que le principe général de l'invention consiste à extraire les correspondances locales entre le maillage initial (surface 3D) décrivant l'objet considéré et l'image 3D suivante, correspondant à la première déformation de cet objet. L'étape de segmentation, qui est cruciale et difficile en analyse d'images, n'est donc ici nécessaire que pour une (préférentiellement la première) image 3D de la séquence analysée. Une fois extraite, ces correspondances permettent alors de modifier le maillage surfacique initial afin d'obtenir un maillage décrivant de façon satisfaisante l'état de l'objet détecté dans le volume considéré après cette première déformation. En d'autres termes, la surface 3D de l'objet contenu dans l'image 3D suivante est obtenue, non pas en effectuant une segmentation de cette image suivante, mais grâce à l'application de vecteurs de déplacement et/ou de déformation, résultant de l'extraction des correspondances, à la surface 3D de l'objet contenu dans l'image 3D courante. Ce processus est appliqué de manière itérative afin d'extraire les correspondances entre les différentes images 3D de la séquence, de façon à décrire tous les mouvements et/ou déformations de l'objet. Il permet d'extraire, de façon simultanée, la position de l'objet déformé dans le volume analysé (image 3D suivante) d'une part, l'amplitude et la direction des déplacements et/ou déformations de l'objet d'autre part, estimés de façon locale. Comme illustré sur la figure 1, dans le mode de réalisation particulier précité, le procédé selon l'invention d'estimation de mouvement et/ou de déformation comprend les étapes suivantes : dans une première étape 11, on choisit l'image de départ (la première image 3D de la séquence dans cet exemple) comme image courante ; dans une seconde étape 12, cette image courante (image de départ) est soumise à un traitement permettant d'obtenir une suriàce 3D de l'objet qu'elle contient. Cette surface 3D est représentée sous la forme d'un maillage comprenant des nœuds ; - dans une troisième étape 13, on détermine des correspondances optimales entre des entités de la surface 3D obtenue pour l'image 3D courante, et des entités d'une image 3D suivante au sein de la séquence, de façon à estimer un champ de vecteurs de déplacement et/ou de déformation associés aux entités de la surface 3D et représentant le mouvement et/ou la déformation de l'objet entre les images courante et suivante ; dans une quatrième étape 14, on obtient une surface 3D de l'objet contenu dans l'image 3D suivante, à partir de la surface 3D obtenue pour l'image 3D courante et du champ de vecteurs de déplacement et/ou de déformation; dans une cinquième étape 15, on détermine si l'image suivante précitée est la dernière image de la séquence, c'est-à-dire si toutes les images 3D de la séquence ont été traitées ; dans l'affirmative, on passe à une sixième étape de fin 16 ; dans la négative, on effectue une septième étape 17 consistant à prendre l'image suivante précitée comme nouvelle image courante (avec la surface 3D issue de la quatrième étape 14) et l'image suivante de la suivante précitée comme nouvelle image suivante, avant de revenir à la troisième étape 13.
Dans une variante de réalisation, la deuxième étape 12 est effectuée avec une image de départ qui n'est pas la première image 3D de la séquence. Il convient dans ce cas de traiter : toutes les images situées avant l'image de départ, au sein de la séquence (l'image courante utilisée lors de chaque nouvelle itération étant l'image suivante, selon un premier sens de parcours de la séquence, utilisée lors de l'itération précédente) ; et - toutes les images situées après l'image de départ, au sein de la séquence (l'image courante utilisée lors de chaque nouvelle itération étant l'image suivante, selon un second sens de parcours de la séquence, utilisée lors de l'itération précédente). Comme illustré sur la figure 2, la deuxième étape 12 comprend par exemple les étapes suivantes : segmentation (121) de l'image de départ, de façon à obtenir les contours et points intérieurs de l'objet contenu dans l'image de départ. Cette étape de segmentation est par exemple réalisée en combinant une approche de croissance de régions et des informations contour extraites par un opérateur gradient (voir [Guillaume'03]). Elle permet de fournir comme résultat un ensemble de points
3D supposés appartenir à l'objet étudié ; reconstruction (122) d'une surface 3D, représentée sous la forme d'un maillage comprenant des nœuds, à partir des points identifiés comme appartenant à l'objet lors de l'étape de segmentation. Cette étape de reconstruction est par exemple réalisée en soumettant l'ensemble des points 3D retenus à l'algorithme classique des Marching Cubes (création de topologies locales, représentées sous forme de noeuds et de liens entre ces noeuds, grâce à l'application de masques caractéristiques) ; pré-traitement (123) de la surface 3D reconstruite, permettant de modifier, si nécessaire, le maillage de façon que la dimension des mailles de la surface 3D soit régulière et proche de la dimension de chacun des voxels des images 3D. L'étape (123) de pré-traitement permet de s'assurer que les arêtes constituant une maille (une maille étant délimitée par trois nœuds) sont de longueur égale aux distances inter- voxels (elles prennent leur valeurs dans l'ensemble {1, V2, V3}. Ce pré-traitement est assuré par un passage en valeurs entières des coordonnées de tous les nœuds (approximation à l'entier le plus proche), une fusion des nœuds qui peuvent être rendus confondus par l'étape précédente, une re- triangularisation du maillage (décomposition ou fusion des mailles ne comportant pas trois nœuds)
II est à noter que cette contrainte sur la longueur des arêtes des mailles devant être égale à une distance inter-voxels peut être relâchée tout en restant dans le cadre de la présente invention, mais implique de re-défïnir le voisinage des entités de base choisies et de re-spécifier également la formulation des énergies.
On décrit maintenant plus en détail la troisième étape 13 de détermination de correspondances optimales entre des entités de la surface 3D obtenue pour l'image 3D courante, et des entités de l'image 3D suivante.
On rappelle que dans la présente description, on utilise indifféremment les termes « image 3D » et « volume » .
Afin d'établir des correspondances optimales entre une surface et un volume complet, la présente invention propose une méthode de mise en correspondance. Ce type de méthode nécessite de considérer trois étapes :
- la sélection des entités à mettre en correspondance,
- la définition d'une énergie locale permettant de comparer deux entités,
- la minimisation d'une énergie globale constituée de la somme des énergies locales. Une des originalités de la méthode de mise en correspondance selon l'invention est qu'elle établit des correspondances entre des entités n'appartenant pas au même domaine. En effet, la mise en correspondance est réalisée entre des entités d'un maillage pré-extrait (défini par des noeuds et des cellules (ou mailles)) et des entités d'un volume complet (composés de voxels, volumes tridimensionnels élémentaires). Il est à noter que ces deux ensembles sont très dissemblables : le maillage apporte plus de connaissance
(notamment sur l'emplacement de la frontière à étudier, potentiellement difficile à segmenter) que le volume, "données brutes" issues de l'acquisition. Les entités utilisées pour la description du maillage sont les noeuds qui le constituent. Les entités à faire correspondre avec ces noeuds sont les voxels du volume. L'énergie locale doit permettre de comparer deux entités, chacune appartenant à son ensemble respectif (noeuds d'un côté, voxels de l'autre) tout en tenant compte des contraintes établies sur les correspondances du voisinage de chaque entité. Elle est par exemple composée de la somme pondérée de deux termes principaux (équation 1) : d'une part un terme de correspondance locale ou terme d'attache aux données E1[J1, d) (exprimant une distance entre estimation du mouvement (fj associée à une entité i (i=l, ...,NE) (NE étant le nombre d'entités considérées) et observation (d)),, d'autre part un terme de régularisation E2 (ft ) (exprimant des contraintes de lissage spatio-temporel des déformations).
E (fι,d) = aιEι(fι, d)+a2E2(fι ) (1)
Le terme d'attache aux données E1[J1, d) est directement lié aux données du problème, c'est-à-dire à l'estimation du déplacement d'un noeud (ft ) entre les instants t et t+1, conditionnellement aux observations (volumes images à t et t+1) (d).
Il est donné par une mesure de cohérence d'un appariement et a pour but de quantifier :
- la probabilité que le voxel considéré à t+1 appartienne à la frontière étudiée, - la correspondance topologique entre le noeud à t et le voxel à t+1 considérés,
- la distance entre le noeud et le voxel.
La probabilité que le voxel considéré appartienne à la frontière étudiée est évaluée grâce à l'application d'un détecteur de contour (qui peut être par exemple le filtrage de Canny) sur le volume image analysé. Ainsi, un voxel détecté comme appartenant au contour est privilégié à un voxel détecté comme n'appartenant pas à une frontière. D'autres types de détecteurs de contour peuvent être retenus dans cette phase, appliqués globalement sur le volume ou localement.
Un terme de correspondance topologique est défini et permet d'effectuer une évaluation de la compatibilité de configuration entre le noeud et le voxel considérés et leur voisinage respectif (le voisinage d'un nœud du maillage est défini par les nœuds ayant une maille en commun avec ce nœud central, et le voisinage d'un voxel, par les voxels vérifiant une connexité première avec ce dernier). Cette correspondance topologique peut se baser sur différents critères descriptifs de la topologie locale. Ainsi, dans le cas de la surface, les courbures moyenne et gaussienne sont de bonnes caractéristiques descriptives. Afin d'implémenter le calcul de ces courbures dans le cas du volume 3D et d'obtenir des descripteurs qui soient comparables à ceux calculés sur la surface, la recherche de la topologie locale du voxel étudié est effectuée. Celle-ci est basée sur l'information de contour fournie par le détecteur de contour et correspond à la recherche de la surface d'iso-intensité (points appartenant au contour) au voisinage du voxel. Une fois cette surface locale trouvée (modélisée par un maillage), le calcul de courbure peut s'effectuer de la même façon que pour la surface, mais est appliqué de manière très locale.
Une autre technique de description de la topologie locale proposée dans cette invention consiste à transposer la topologie du voisinage du nœud étudié au niveau du voxel candidat. Cette topologie est décrite en termes de différences de coordonnées des voisins par rapport aux coordonnées du nœud étudié. Les valeurs des voxels correspondants à cette topologie (dans le volume filtré avec le détecteur de contours) sont sommées. Ainsi, plus cette somme est élevée, plus les topologies locales des deux entités sont considérées comme proches. A ce niveau, le voisinage défini peut être étendu à un nombre plus important de mailles. Enfin, afin de privilégier les correspondances représentant des mouvements faibles, la distance euclidienne entre noeud considéré et voxel candidat est calculée. Les correspondances à faible distance sont privilégiées.
Le terme d'attache aux données E1[J1, d) peut donc être exprimé pour chaque couple Nœud (N1 , i = \, ...,NS, NS : nombre de nœuds sur la surface à t) - Voxel (Vk , Jc = I, ...,NV, NV : nombre de voxels du volume à t+1), à partir de trois termes représentatifs des potentiels d'énergie précédemment décrits.
EU,d) [N1 , Vk ) = ac Econ1our [N1 , Vk)+at E1opol [N1 , Vk )+ad Eώs1 [N1 , Vk ) A ces trois termes composant la fonction d'énergie locale, d'autres termes peuvent également être ajoutés, introduisant par exemple une contrainte sur l'intensité de niveaux de gris de l'image. En effet, un terme d'énergie pourrait être défini en fonction de l'information densitométrique des images à partir du moment où une correspondance entre intensités de niveau de gris peut être établie entre les deux volumes considérés (1er volume à partir duquel est défini la surface 3D et 2e"16 volume dans lequel est recherché la correspondance). Le terme de régularisation E2[J1 ) permet d'éviter l'obtention de correspondances aberrantes (par exemple les mouvements de deux noeuds voisins qui se croisent). Elle est basée sur le fait que les correspondances de noeuds voisins doivent engendrer des mouvements proches, que ce soit en amplitude ou en direction. Différentes implémentations de cette énergie peuvent être envisagées (sachant que l'énergie finalement retenue peut être également une combinaison de différentes formulations). Par exemple, une première formulation peut correspondre à l'écart entre le déplacement évalué et la moyenne des déplacements des voisins du noeud considéré. D'autres formulations peuvent également convenir. Par exemple, une deuxième formulation est proposée dans le cadre de la formulation Markovienne (voir discussion détaillée ci- après), permettant de maintenir l'homogénéité des appariements. Afin d'estimer globalement les correspondances à partir des énergies locales formulées précédemment, une formulation markovienne est par exemple utilisée.
Le modèle de Markov permet de formuler un problème d'estimation conditionnellement à des observations, sous la forme d'un modèle global décrit par des spécifications locales (introduisant la définition de voisinage). Dans sa définition, le mouvement 3D recherché est considéré comme une réalisation f = {ft , i = \,...,NS} d'un champ aléatoire F = {Ft , i = l,..., NS} regroupant toutes les réalisations possibles. Les sites du champ considéré sont donnés par tous les nœuds représentatifs du maillage de la surface. Les labels de ces sites ou étiquettes recherchées (estimations ft) sont les voxels jugés correspondants de ces noeuds. Le voisinage d'un site (ou d'un nœud) est composé des noeuds ayant une maille en commun avec le noeud étudié.
La réalisation la plus probable du champ des étiquettes est donnée par la maximisation d'une probabilité P {f) qui s'écrit :
Figure imgf000022_0001
où T est une valeur constante appelée température et U (f) est une fonction d'énergie globale,
Z = ∑ exp \ -ψ U (/) L Ω étant l'ensemble des étiquettes possibles.
La réalisation la plus probable revient alors minimiser la fonction d'énergie globale U (f). La fonction U [f) ou plus généralement U (f, d) s'exprime par l'équation suivante :
U [f, d)=aι Ud [f,d)+a2 U, [f) (4)
Où Ud [f, d) modélise l'erreur de l'estimation par rapport aux observations et le terme U1 (/) modélise l'énergie interne du champ correspondant au mouvement à restituer. Ce deuxième terme joue un rôle de stabilisation dans le problème de régularisation posé.
Le premier terme s'exprime comme une erreur quadratique de manière que la minimisation de cette composante d'énergie favorise les configurations estimées proches des observations. Il est calculé à partir de la mesure d'appariement locale (ou distance) développée précédemment, sur tous les sites du champ :
Ud (f, d) = ∑ (E1 (f,, d)f , L étant l'ensemble des sites du champ (nœuds du
maillage) (5)
Le second terme est exprimé de la façon suivante :
U, [f) = ∑ E2 [f,) (6) v≡L
Une formulation de E2 (/) est également proposée par l'équation suivante :
^2(/) = ∑ V \f, > fj)> ou Α, est le voisinage de première connexité associé au m, site i (7) où v (ft , fj) est le potentiel d'interaction associé à une clique du voisinage (une clique étant représentée ici par deux sites) avec :
Figure imgf000023_0001
où δ(.) est une fonction delta Kronecker qui prend la valeur 1 si les étiquettes sont proches et la valeur 0 sinon.
Plusieurs formulations de E2 (/) sont ici possibles en fonction des cliques définies sur le voisinage et de leur potentiel d'interaction associé.
On rappelle qu'une clique définie sur un voisinage est constituée d'un sous- ensemble d'éléments du voisinage, contenant de 1 à N éléments (le site compris), N étant le nombre d'éléments du voisinage. Ces éléments constituant la clique doivent vérifier une connexité (point ou arête de maille en commun dans le cadre de cette invention). La minimisation de l'énergie globale U (f,d) permet de trouver les correspondances optimales et donc l'estimation de déformation recherchée pour tout l'objet. Les différentes réalisations /peuvent être obtenues selon plusieurs algorithmes d'optimisation. Une méthode pour résoudre ce problème, dans le contexte des champs de Markov, est la recherche aléatoire combinée avec une technique de recuit simulé.
L'algorithme de relaxation stochastique de Metropolis est par exemple choisi, combiné avec un recuit simulé, permettant de mieux assurer la convergence vers un minimum global, comparé à un algorithme de relaxation déterministe. D'autres algorithmes d'optimisation sont ici possibles. Un exemple d'algorithme d'optimisation s'écrit comme suit :
Initialisation des correspondances T := T0
Tant que {Critère d'arrêt non vérifié} Nbre essais à T := 0
Nbre succès := 0 ;
Tant que {Nbre essais à T < Nbre essais à T max} ou {Nbre succès à T < Nbre succès à T max}
Sélection du nœud étudié (qui peut être choisi aléatoirement) Sélection du voxel candidat
Nbre_succès_à_T : = Nbre_succès_à_T + 1 Calcul de AU(f,d)
Si { AU(f,d) <=0}
Validation de la correspondance candidate Sinon
Proba = (exp(- AU(f,d) /T)) Si {Proba > random(0,l)}
Validation de la correspondance candidate Nbre_succès := Nbre_succès+1 ; Fm
Fin Fin
T :=T-pas_T (diminution de T) Fin
La valeur initiale des paramètres impliqués dans cette procédure d'optimisation
(température T0, Nbre essais à T max (nombre d'essais maximum à une température donnée) et Nbre succès max (nombre maximum de succès à AU(f,d) positif) a été déterminée empiriquement.
La démarche développée dans ce schéma est la suivante. C'est un procédé itératif dans lequel une température donnée initialement va décroître jusqu'à la rencontre d'un critère d'arrêt, pour lequel la solution de correspondance globale est proposée. La loi de décroissance de la température à chaque itération peut prendre plusieurs formes, par exemple résulter de la décrémentation de la température d'un pas.
Pour chaque itération (à une température donnée), plusieurs nœuds sont choisis (soit de manière aléatoire) et traités, ceci jusqu'à ce qu'un nombre suffisant de nœuds ait été considéré (Nbre essais à T ) ou jusqu'à ce qu'un nombre donné de succès
(Nbre succès à T) ait été atteint. Pour chacun des nœuds sélectionnés, la possibilité de modifier la correspondance en cours est étudiée et acceptée en fonction de la variation d'énergie globale obtenue par ce changement. Un succès correspond à un changement de correspondance pour un nœud, qui a été accepté avec une variation d'énergie globale positive.
La valeur initiale des paramètres impliqués dans cette procédure d'optimisation a été déterminée empiriquement. Ces paramètres sont :
T0 : la température de départ, Pas T : le pas de décroissance de la température,
Nbre essais à T max : le nombre maximum d'essais à une température donnée,
Nbre succès à T max : le nombre maximum de succès, à une température donnée, avec AU(f,d) (variation de l'énergie globale) positive (c'est-à-dire entraînant une augmentation de l'énergie globale). Le critère d'arrêt peut être défini de différentes façons : soit la température T doit atteindre un seuil (qui peut être la température 0), soit le nombre de modifications de correspondance à une température donnée doit être inférieur à un seuil.
La sélection du nœud étudié peut-être soit complètement aléatoire, soit guidée par des critères d'exhaustivité (par exemple empêcher de sélectionner deux fois le même nœud pour une température donnée).
La sélection du voxel candidat s'effectue en choisissant (soit aléatoirement, soit suivant des critères d'exhausitivité ou de non répétition) un voxel voisin du voxel correspondant courant. Ce voisinage peut être défini comme étant un voisinage de premier ordre (les deux voxels ont alors une face ou une arrête en commun), ou d'ordre supérieur (par exemple pour l'ordre deux, les deux voxels ont une face ou une arrête ou un sommet en commun).Optionnellement, afin de tenir compte et d'exploiter les différents niveaux de détails topologiques des surfaces complexes (surface 3D de l'objet), de converger plus vite vers la solution optimale et dans le but de réduire le temps de calcul, la troisième étape 13 de détermination de correspondances optimales est appliquée dans un mode multi-résolution.
Dans ce cas, la surface 3D obtenue pour l'image 3D courante, lors de la deuxième étape 12 ou la quatrième étape 14, est considérée comme une surface 3D de résolution maximale. En outre, l'image 3D suivante comprise dans l'ensemble d'images 3D est considérée comme une image 3D de résolution maximale.
Le mode multi-résolution nécessite de pouvoir passer d'une résolution à une résolution inférieure, que ce soit pour les volumes de données originaux (c'est-à-dire les images 3D originales de la séquence source) ou pour les maillages 3D (c'est-à-dire les surfaces 3D). Pour les volumes originaux, le passage d'une représentation de résolution donnée à une représentation de résolution inférieure est par exemple réalisé par un lissage (qui peut être un filtrage gaussien) suivi d'un sous-échantillonnage appliqués au volume image.
Pour les maillages 3D, le passage d'une représentation de résolution donnée à une représentation de résolution inférieure est par exemple réalisé par un lissage (qui peut être un filtrage gaussien) suivi d'un sous-échantillonnage des coordonnées des nœuds, puis une opération de fusion des nœuds confondus et de re-triangularisation du maillage (afin que toutes les mailles soient constituées de trois noeuds).
Pour traiter le premier intervalle temporel, t - 1+1, le changement de résolution à t est appliqué sur le maillage de surface, tandis que le changement de résolution à t+1 est effectué sur le volume d'origine (image de niveaux de gris). Pour traiter chaque intervalle temporel suivant, t+k - t+k+1, le changement de résolution à t+k est appliqué sur le maillage de surface (issu de l'estimation de correspondance réalisée aux instants précédents et ayant subi le pré-traitement déjà cité (régularisation et re- triangularisation)), tandis que le changement de résolution à t+k+1 est effectué sur le volume d'origine.
Dans ce cas, et comme illustré sur la figure 3, la troisième étape 13 précitée comprend les étapes successives suivantes : dans une étape référencée 131, obtention d'une part d'au moins une surface 3D de résolution inférieure (pour l'image 3D courante, et d'autre part d'au moins une image 3D suivante de résolution inférieure (voir la discussion détaillée ci- dessus donnant un exemple de passage d'une résolution donnée à une résolution inférieure, que ce soit pour les images 3D originales de la séquence source ou pour les surfaces 3D)... On travaille par exemple avec trois niveaux de résolution, y inclus le niveau de résolution maximal ; - dans une étape référencée 132, choix de la plus basse des résolutions comme résolution courante, et initialisation des correspondances à la résolution courante ; dans une étape référencée 133, pour la résolution courante, détermination de correspondances optimales entre des noeuds de la surface 3D (possédant cette résolution courante) contenue dans l'image 3D courante, et des voxels de l'image 3D suivante (possédant cette même résolution courante) de la séquence.
Cette étape 133 est par exemple réalisée selon l'exemple d'algorithme d'optimisation décrit ci-dessus (algorithme de Metropolis combiné au recuit simulé), dans le cadre d'une formulation markovienne du problème d'estimation; - dans une étape référencée 134, on détermine si la résolution courante est égale à la résolution maximale ; dans l'affirmative, la troisième étape 13 est terminée (et on passe à l'étape 14 suivante, voir la figure 1) ; dans la négative, on retourne à l'étape référencée 133, après exécution d'une étape référencée 135 consistant à : * propager les résultats de l'exécution de l'étape référencée 133 (c'est-à- dire les champs de vecteurs de déplacement et/ou de déformation estimés pour les nœuds de la surface 3D), de la résolution courante à une résolution supérieure suivante, de façon à permettre une réinitialisation des correspondances ; * prendre comme nouvelle résolution courante la résolution supérieure suivante.
Dans l'étape référencée 132, l'initialisation des correspondances à la résolution la plus basse est par exemple réalisée par des mouvements nuls, chaque nœud (de la surface 3D obtenue pour l'image courante) étant mis en correspondance avec le voxel le plus proche (au sein de l'image 3D suivante).
Dans l'étape référencée 135, la propagation des résultats, d'une résolution donnée à la résolution supérieure, c'est-à-dire la réinitialisation des correspondances est par exemple assurée par interpolation, sur tout le maillage de la surface 3D possédant cette résolution supérieure, des correspondances issues de la résolution inférieure (c'est- à-dire déterminées au cours de l'étape référencée 133).
Pour ce faire, pour chaque nœud à la résolution supérieure, le ou les nœuds correspondants à la résolution inférieure (suivant que le sous-échantillonnage donne des coordonnées entières ou non) sont sélectionnés et leur déplacement moyenne. Ce déplacement moyen est transposé à la résolution supérieure (par multiplication par le facteur de sous-échantillonnage) et appliqué après passage en valeurs entières (un voxel correspondant est estimé).Quand la troisième étape 13 est terminée, on dispose ainsi d'un champ de vecteurs de déplacement et/ou de déformation estimés en chaque point de la surface 3D obtenue pour l'image courante. Ce champ de vecteurs représente le mouvement et/ou la déformation de l'objet entre l'instant t d'acquisition de l'image 3D courante et l'instant t+1 d'acquisition de l'image 3D suivante. Ces vecteurs, localisés dans l'espace 3D, sont caractérisés par une orientation et une amplitude. Comme déjà indiqué ci-dessus, ces vecteurs estimés permettent également d'extraire les points de contours de l'objet à l'instant t+1 et donc de calculer la nouvelle surface correspondant à l'instant t+1 (étape référencée 14 sur la figure 1). Puis, ce processus de mise en correspondance est appliqué de la même manière sur tous les instants successifs de la séquence (entre t+1 et t+2, puis t+2 et t+3 etc), à partir de la nouvelle surface calculée et du volume image de l'instant suivant. Il fournit ainsi un champ de vecteurs de déplacement à chaque instant de la séquence, associé aux nœuds des mailles des surfaces successivement estimées.
On présente maintenant, en relation avec la figure 4, un mode de réalisation particulier du procédé de recalage tridimensionnel selon l'invention.
De façon classique, un tel procédé de recalage comprend : une étape (41) d'estimation du mouvement et/ou de la déformation subi(s), entre une première et une seconde images 3D acquises à deux instants différents et/ou à partir de deux modalités d'imagerie différentes et/ou à partir d'une modalité d'imagerie et un atlas (modèle préexistant), par (au moins) un objet contenu dans ces première et seconde images 3D ; une étape (42) d'application à l'objet contenu dans la seconde image d'une fonction corrigeant le mouvement et/ou la déformation estimé lors de l'étape d'estimation. L'originalité réside ici dans le fait que l'étape (41) comprend elle-même les étapes référencées 11 à 13 du procédé d'estimation discuté ci-dessus en relation avec les figures 1 à 3. Ces étapes 11 à 13 ne sont donc pas décrites à nouveau. On notera simplement que l'ensemble d'images 3D considéré ne comprend que deux images 3D (il ne s'agit donc pas d'une séquence temporelle comprenant plus de deux images). Bien que l'invention ait été décrite ci-dessus en relation avec un nombre limité de modes de réalisation, l'homme du métier, à la lecture de la présente description, comprendra que d'autres modes de réalisation peuvent être imaginés sans sortir du cadre de la présente invention. En conséquence, la portée de l'invention n'est limitée que par les revendications ci-jointes. Annexe : bibliographie
[Amini '92] Amini, A. A. and Duncan, J. S. (1992). "Bending and stretching models for LV wall motion analysis from curves and surfaces." Image and Vision Computing 10(6): 418-430. [Bardinet '96] Bardinet, E., Cohen, L. D. et al. (1996). "Tracking and motion analysis of the left ventricle with deformable superquadrics." Médical Image Analysis 1(2): 129-49.
[Benayoun '95] Benayoun S., Nastar C, Ayache N. (1995), "Dense Non-Rigid Motion Estimation in Séquences of 3D images Using Differential Constraints", Computer Vision, Virtual Reality and Robotics in Medicine, pp. 309-318.
[Chen '94] Chen, C. W., Huang, T. S. et al. (1994). "Modeling, analysis and visualization of left ventricle shape and motion by hierarchical décomposition." IEEE Transactions on Pattern Analysis and Machine Intelligence 16(4): 342-356.
[Choi 1Ol] Choi, S.-M. and Kim, M. -H. (2001). "Motion visualization of human left ventricle with a time-varying deformable model for cardiac diagnosis." The journal of visualization and computer animation 12: 55-66.
[Coatrieux '94] Coatrieux, J.L., Garreau M., Ruan S., Mao F. (1994). "On motion analysis in médical imaging". Innov. Technol. Biol. Med., 3: 253-267.
[Gorce '97] Gorce J.M., Friboulet D., Magnin LE. (1997), "Estimation of three- dimensional cardiac velocity fields : assessment of a differential method and application to three-dimensional et data", Médical Image Analysis, 1(3) : 245-61. [Guillaume '03] Guillaume H., Garreau M., Boldak C, Larralde A. (2003). "Segmentation de cavités cardiaques en imagerie scanner multbarettes", 12ème Forum des Jeunes Chercheurs en Génie biologique et Médical, Nantes, 92-93. [Kambhamettu '03] Kambhamettu, C, Goldgof, D., et al. (2003). "3D nonrigid motion analysis under small déformations." Image and Vision Computing 21(3): 229- 245.
[Lin '0O] Lin, W. and Robb, R. A. (2000). "Visualization of cardiac dynamics using physics-based deformable mode"l. SPIE Médical Imaging 2000: Image Display and Visualization. [Ruan '94] Ruan S., Bruno A., Coatrieux J. L. (1994). "Three-dimensional motion and reconstruction of coronary arteries from biplane cineangiography", Image and Vision Computing, 12,10: 683-689.
[Shi 1OO] Shi, P., Sinusas, A. J., et al. (2000). "Point-tracked quantitative analysis of left ventricular motion from 3D image séquences." IEEE Transactions on Médical Imaging 19(1): 36-50.
[Song '91] Song, S. M. and Leahy, R. M. (1991). "Computation of 3-D velocity fields from 3-D cine et images of the human heart." IEEE Transactions on Médical Imaging 10(3): 295-306.

Claims

REVENDICATIONS
1. Procédé d'estimation du mouvement et/ou de la déformation d'au moins un objet contenu dans chaque image 3D d'un ensemble d'images 3D comprenant au moins deux images 3D acquises à différents instants et/ou à partir de plusieurs modalités d'imagerie et/ou avec une modalité d'imagerie et un atlas, caractérisé en ce qu'il comprend les étapes suivantes : a) obtention (12) d'une surface 3D dudit au moins un objet contenu dans l'une des images 3D dudit ensemble, dite image de départ, par traitement de l'image de départ, ladite surface 3D étant représentée sous la forme d'un maillage comprenant des nœuds ; b) choix (11) de l'image de départ comme image courante ; c) détermination (13) de correspondances optimales entre des entités de la surface 3D obtenue pour l'image 3D courante, et des entités d'une image 3D suivante dudit ensemble, de façon à estimer un champ de vecteurs de déplacement et/ou de déformation associés aux entités de la surface 3D et représentant le mouvement et/ou la déformation dudit au moins un objet entre les images courante et suivante.
2. Procédé selon la revendication 1, caractérisé en ce que ledit ensemble d'images 3D comprend deux images 3D acquises à deux instants différents et/ou avec deux modalités d'imagerie différentes et/ou avec une modalité d'imagerie et un atlas.
3. Procédé selon la revendication 1, caractérisé en ce que ledit ensemble d'images 3D est une séquence d'images 3D comprenant au moins trois images 3D acquises à des instants successifs au sein d'une même plage temporelle ou de plusieurs plages temporelles non adjacentes, et avec une même modalité d'imagerie ou plusieurs modalités d'imagerie.
4. Procédé selon la revendication 3, caractérisé en ce que ladite image de départ est la première image 3D de ladite séquence d'images 3D, et en ce que le procédé comprend en outre les étapes suivantes : d) obtention (14) d'une surface 3D dudit au moins un objet contenu dans l'image 3D suivante, à partir de la surface 3D obtenue pour l'image 3D courante et du champ de vecteurs de déplacement et/ou de déformation ; e) itération des étapes c) et d) jusqu'à ce que toutes les images 3D de la séquence aient été traitées, l'image courante utilisée lors de chaque nouvelle itération étant l'image suivante utilisée lors de l'itération précédente.
5. Procédé selon la revendication 3, caractérisé en ce que ladite image de départ est différente de la première image 3D de ladite séquence d'images 3D, et en ce que le procédé comprend en outre les étapes suivantes : d') obtention d'une surface 3D dudit au moins un objet contenu dans l'image 3D suivante, à partir de la surface 3D obtenue pour l'image 3D courante et du champ de vecteurs de déplacement et/ou de déformation ; e') itération des étapes c) et d'), en parcourant ladite séquence d'images 3D selon un premier sens et jusqu'à ce que toutes les images 3D depuis l'image de départ jusqu'à la dernière image de la séquence aient été utilisées, l'image courante utilisée lors de chaque nouvelle itération étant l'image suivante, selon le premier sens, utilisée lors de l'itération précédente ; e") itération des étapes c) et d'), en parcourant ladite séquence d'images 3D selon un second sens et jusqu'à ce que toutes les images 3D depuis l'image de départ jusqu'à la première image de la séquence aient été utilisées, l'image courante utilisée lors de chaque nouvelle itération étant l'image suivante, selon le second sens, utilisée lors de l'itération précédente.
6. Procédé selon l'une quelconque des revendications 1 à 5, caractérisé en ce que les entités d'une surface 3D donnée sont des nœuds du maillage de ladite surface 3D donnée.
7. Procédé selon l'une quelconque des revendications 1 à 6, caractérisé en ce que les entités d'une image 3D donnée sont des voxels de ladite image 3D donnée.
8. Procédé selon l'une quelconque des revendications 1 à 7, caractérisé en ce que ladite étape c) comprend une étape de minimisation d'une énergie globale constituée d'une somme d'énergies locales, chaque énergie locale permettant de comparer une entité de la surface 3D, obtenue pour l'image 3D courante, et une entité de l'image 3D suivante.
9. Procédé selon la revendication 8, caractérisé en ce que chaque énergie locale est composée d'une somme pondérée de : un terme de correspondance locale E1(Cd), exprimant une mesure de cohérence d'un appariement entre une entité de la surface 3D, obtenue pour l'image 3D courante, et une entité de l'image 3D suivante ; un terme de régularisation E2(T), exprimant des contraintes de lissage spatio- temporel du mouvement et/ou de la déformation estimée, lesdites contraintes étant telles que des vecteurs de déplacement associés à des entités voisines de la surface 3D doivent être proches en amplitude et/ou en direction.
10. Procédé selon la revendication 9, caractérisé en ce que ledit terme de correspondance locale E1(Cd) est fonction d'un terme exprimant une correspondance topologique entre l'entité de la surface 3D, obtenue pour l'image 3D courante, et l'entité de l'image 3D suivante.
11. Procédé selon la revendication 10, caractérisé en ce que ledit terme de correspondance locale E1(Cd) est composé d'une somme pondérée dudit terme exprimant une correspondance topologique et d'au moins un des termes appartenant au groupe comprenant : un terme exprimant une probabilité que l'entité de l'image 3D suivante appartienne à un élément de contour, dudit au moins un objet contenu dans ladite image 3D suivante ; un terme exprimant une distance entre l'entité de la surface 3D, obtenue pour l'image 3D courante, et l'entité de l'image 3D suivante ; un terme lié aux densités de niveaux de gris.
12. Procédé selon l'une quelconque des revendications 8 à 11, caractérisé en ce que ladite étape de minimisation d'une énergie globale est basée sur une formulation statistique d'un processus d'estimation selon un modèle de Markov, telle que le mouvement et/ou la déformation à estimer est considéré comme une réalisation f = {f, , i = l,...,NS} d'un champ aléatoire F = [F1 , i = 1,..., NS} regroupant toutes les réalisations possibles, les NS sites du champ aléatoire F étant donnés par toutes les entités (par ex. nœuds) de la surface 3D, obtenue pour l'image 3D courante, représentatives du maillage de ladite surface 3D, des labels desdits sites, aussi appelés étiquettes recherchées, étant les entités (par ex. voxels) de l'image 3D suivante jugées correspondantes desdites entités de la surface 3D.
13. Procédé selon l'une quelconque des revendications 1 à 12, caractérisé en ce que la surface 3D obtenue pour l'image 3D courante, lors de l'étape a) ou d), est considérée comme une surface 3D de résolution maximale, en ce que l'image 3D suivante comprise dans ledit ensemble d'images 3D est considérée comme une image 3D de résolution maximale, et en ce que ladite étape c) est effectuée dans un mode multi-résolution et comprend les étapes successives suivantes : cl) obtention (131) d'une part d'au moins une surface 3D de résolution inférieure pour l'image 3D courante, à partir de la surface 3D de résolution maximale pour l'image 3D courante, et d'autre part d'au moins une image 3D suivante de résolution inférieure, à partir de l'image 3D suivante de résolution maximale ; c2) choix (132) de la plus basse des résolutions comme résolution courante, et initialisation des correspondances à la résolution courante ; c3) pour la résolution courante, détermination (133) de correspondances optimales entre des entités de la surface 3D possédant ladite résolution courante pour l'image 3D courante, et des entités de l'image 3D suivante possédant également ladite résolution courante ; c4) si la résolution courante n'est pas égale à la résolution maximale :
* propagation (135) de résultats de l'exécution de l'étape c3), de la résolution courante à une résolution supérieure suivante, de façon à permettre une réinitialisation des correspondances ;
* retour à l'étape c3), en prenant comme nouvelle résolution courante la résolution supérieure suivante ; sinon arrêt de l'étape c).
14. Procédé selon l'une quelconque des revendications 1 à 13, caractérisé en ce que ladite étape a) (12) comprend les étapes suivantes : segmentation (121) de l'image de départ, de façon à extraire des contours et points intérieurs audit au moins un objet contenu dans l'image de départ ; reconstruction (122) d'une surface 3D, représentée sous la forme d'un maillage comprenant des nœuds, à partir des contours et points intérieurs obtenus lors de l'étape de segmentation.
15. Procédé selon l'une quelconque des revendications 13 et 14, caractérisé en ce que, dans l'étape c2), l'initialisation des correspondances à la résolution la plus basse est réalisée par des mouvements nuls, chaque entité de la surface 3D obtenue pour l'image 3D courante étant mise en correspondance avec l'entité la plus proche dans l'image 3D suivante.
16. Procédé selon l'une quelconque des revendications 13 à 15, caractérisé en ce que, dans l'étape c4), la réinitialisation des correspondances est obtenue par interpolation, sur tout le maillage de la surface 3D possédant ladite résolution supérieure suivante, des correspondances optimales déterminées lors de l'étape c3).
17. Procédé selon l'une quelconque des revendications 14 à 16, caractérisé en ce que ladite étape a) comprend en outre l'étape suivante : pré-traitement (123) de la surface 3D reconstruite, permettant de modifier, si nécessaire, ledit maillage de façon que la longueur de chacune des arêtes constituant une maille soit égale à une des distances entre voxels des images 3D.
18. Procédé selon l'une quelconque des revendications 1 à 17, caractérisé en ce que lesdites modalités d'imagerie appartiennent au groupe comprenant : l'imagerie scanner ; l'imagerie par résonance magnétique ; l'imagerie par ultrasons ; - l'imagerie nucléaire.
19. Application du procédé selon l'une quelconque des revendications 1 à 18 dans un domaine d'imagerie appartenant au groupe comprenant : l'imagerie médicale ; l'imagerie marine ; - l'imagerie satellitaire ; l'imagerie des matériaux.
20. Produit programme d'ordinateur, caractérisé en ce qu'il comprend des instructions de code de programme pour l'exécution des étapes du procédé selon l'une quelconque des revendications 1 à 18, lorsque ledit programme est exécuté sur un ordinateur.
21. Moyen de stockage, éventuellement totalement ou partiellement amovible, lisible par un ordinateur, stockant un jeu d'instructions exécutables par ledit ordinateur pour mettre en œuvre le procédé selon l'une quelconque des revendications 1 à 18.
22. Procédé de recalage tridimensionnel d'au moins un objet contenu dans une image 3D, caractérisé en ce qu'il comprend les étapes suivantes : estimation (41) du mouvement et/ou de la déformation subi(s), entre une première et une seconde images 3D acquises à deux instants différents et/ou à partir de deux modalités d'imagerie différentes et/ou avec une modalité d'imagerie et un atlas, par au moins un objet contenu dans lesdites première et seconde images 3D ; application (42) audit au moins un objet contenu dans la seconde image d'une fonction corrigeant le mouvement et/ou la déformation estimé lors de l'étape d'estimation, caractérisé en ce que ladite étape d'estimation (41) est mise en œuvre selon le procédé d'estimation de l'une quelconque des revendications 1, 2 et 6 à 18, sauf quand elles dépendent de l'une des revendications 3 à 5.
PCT/EP2005/055690 2004-11-02 2005-11-02 Procede d'estimation du mouvement et/ou de la deformation d'au moins un objet contenu dans un ensemble d'images 3d, et procede de recalage tridimensionnel correspondant WO2006048419A1 (fr)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
FR0411672A FR2877465B1 (fr) 2004-11-02 2004-11-02 Procede d'estimation du mouvement et/ou de la deformation d'au moins un objet dans un ensemble d'images 3d, et procede de recalage tridimensionnel correspondant
FR0411672 2004-11-02

Publications (1)

Publication Number Publication Date
WO2006048419A1 true WO2006048419A1 (fr) 2006-05-11

Family

ID=34951935

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/EP2005/055690 WO2006048419A1 (fr) 2004-11-02 2005-11-02 Procede d'estimation du mouvement et/ou de la deformation d'au moins un objet contenu dans un ensemble d'images 3d, et procede de recalage tridimensionnel correspondant

Country Status (2)

Country Link
FR (1) FR2877465B1 (fr)
WO (1) WO2006048419A1 (fr)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110008959A (zh) * 2019-03-26 2019-07-12 北京博瑞彤芸文化传播股份有限公司 一种医学数据处理方法及系统

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6289135B1 (en) * 1997-03-20 2001-09-11 Inria Institut National De Recherche En Informatique Et En Antomatique Electronic image processing device for the detection of motions
US20010043722A1 (en) * 2000-03-10 2001-11-22 Wildes Richard Patrick Method and apparatus for qualitative spatiotemporal data processing

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6289135B1 (en) * 1997-03-20 2001-09-11 Inria Institut National De Recherche En Informatique Et En Antomatique Electronic image processing device for the detection of motions
US20010043722A1 (en) * 2000-03-10 2001-11-22 Wildes Richard Patrick Method and apparatus for qualitative spatiotemporal data processing

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
PENNEC X ET AL: "Tracking brain deformations in time sequences of 3D US images", PATTERN RECOGNITION LETTERS, NORTH-HOLLAND PUBL. AMSTERDAM, NL, vol. 24, no. 4-5, February 2003 (2003-02-01), pages 801 - 813, XP004391218, ISSN: 0167-8655 *
POSITANO V ET AL: "Automatic time sequence alignment in contrast enhanced MRI by maximization of mutual information", PROCEEDINGS OF THE 23RD. ANNUAL INTERNATIONAL CONFERENCE OF THE IEEE ENGINEERING IN MEDICINE AND BIOLOGY SOCIETY. 2001 CONFERENCE PROCEEDINGS. (EMBS). INSTANBUL, TURKEY, OCT. 25 - 28, 2001, ANNUAL INTERNATIONAL CONFERENCE OF THE IEEE ENGINEERING IN M, vol. VOL. 1 OF 4. CONF. 23, 25 October 2001 (2001-10-25), pages 2407 - 2410, XP010592138, ISBN: 0-7803-7211-5 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110008959A (zh) * 2019-03-26 2019-07-12 北京博瑞彤芸文化传播股份有限公司 一种医学数据处理方法及系统

Also Published As

Publication number Publication date
FR2877465B1 (fr) 2006-12-29
FR2877465A1 (fr) 2006-05-05

Similar Documents

Publication Publication Date Title
Saha et al. A survey on skeletonization algorithms and their applications
EP1851722B1 (fr) Dispositif et procede de traitement d&#39;image
Xia et al. Super-resolution of cardiac MR cine imaging using conditional GANs and unsupervised transfer learning
US8111885B2 (en) Image processing device and method
RU2522038C2 (ru) Сегментация при мр-визуализации сердца в проекции по длинной оси с поздним усилением контраста
Ouzir et al. Motion estimation in echocardiography using sparse representation and dictionary learning
Dakua LV segmentation using stochastic resonance and evolutionary cellular automata
Gao et al. Simultaneous multi-object segmentation using local robust statistics and contour interaction
Mirunalini et al. Segmentation of Coronary Arteries from CTA axial slices using Deep Learning techniques
Jaffar et al. Anisotropic diffusion based brain MRI segmentation and 3D reconstruction
Gao et al. 4D cardiac reconstruction using high resolution CT images
Zhang et al. 3D anatomical shape atlas construction using mesh quality preserved deformable models
Li et al. 3D coronary artery reconstruction by 2D motion compensation based on mutual information
Khaled et al. Automatic fuzzy-based hybrid approach for segmentation and centerline extraction of main coronary arteries
Archip et al. Ultrasound image segmentation using spectral clustering
WO2006048419A1 (fr) Procede d&#39;estimation du mouvement et/ou de la deformation d&#39;au moins un objet contenu dans un ensemble d&#39;images 3d, et procede de recalage tridimensionnel correspondant
Duan et al. Surface function actives
Arega et al. Using Polynomial Loss and Uncertainty Information for Robust Left Atrial and Scar Quantification and Segmentation
van Walsum et al. Averaging centerlines: mean shift on paths
Behbahani et al. Analysing optical flow based methods
Fritz et al. Automatic 4D segmentation of the left ventricle in cardiac-CT-data
Metaxas et al. A segmentation and tracking system for 4D cardiac tagged MR images
Upendra CNN-Based Cardiac Motion Extraction to Generate Deformable Geometric Ventricular Models from Cine MRI
Beitone et al. Mutual cineMR/RT3DUS cardiac segmentation
CN114882050A (zh) 一种心肌灌注spect图像分割方法

Legal Events

Date Code Title Description
AK Designated states

Kind code of ref document: A1

Designated state(s): AE AG AL AM AT AU AZ BA BB BG BR BW BY BZ CA CH CN CO CR CU CZ DE DK DM DZ EC EE EG ES FI GB GD GE GH GM HR HU ID IL IN IS JP KE KG KM KN KP KR KZ LC LK LR LS LT LU LV LY MA MD MG MK MN MW MX MZ NA NG NI NO NZ OM PG PH PL PT RO RU SC SD SE SG SK SL SM SY TJ TM TN TR TT TZ UA UG US UZ VC VN YU ZA ZM ZW

AL Designated countries for regional patents

Kind code of ref document: A1

Designated state(s): BW GH GM KE LS MW MZ NA SD SL SZ TZ UG ZM ZW AM AZ BY KG KZ MD RU TJ TM AT BE BG CH CY CZ DE DK EE ES FI FR GB GR HU IE IS IT LT LU LV MC NL PL PT RO SE SI SK TR BF BJ CF CG CI CM GA GN GQ GW ML MR NE SN TD TG

DPEN Request for preliminary examination filed prior to expiration of 19th month from priority date (pct application filed from 20040101)
121 Ep: the epo has been informed by wipo that ep was designated in this application
NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 05810972

Country of ref document: EP

Kind code of ref document: A1