Method for precisely and accurately segmenting a 2d or 3D image

Abstract

The present invention relates to a computer-implemented method for segmenting an image, the method comprising (S1) obtaining a 2D / 3D image I of a region of interest; (S2) initially segmenting the 2D / 3D image I; (S3) processing each image Pj on the basis of one or more margins of error in the initial segmentation of the image I so as to obtain a corresponding image Uj; (S4) processing the pixels / voxels of the 2D / 3D image I corresponding to the uncertain region in the image Uj; (S5) determining an image Ej representing the difference between the pixels / voxels Fjn of the mask Fj and the pixels / voxels Pjn of the image Pj; (S7) assigning each pixel / voxel Ejn representing a difference so as to assign it or not to the object Aj, thereby obtaining a final segmentation mask Mj comprising a refined segmentation of the object Aj.

Description

[0001] DESCRIPTION

[0002] TITLE: A method for the precise and accurate segmentation of a 2D or 3D image

[0003] FIELD OF INVENTION

[0004] The invention relates to a method for segmenting a 2D or 3D image, for example of a subject, a patient, or a structure. In particular, the 2D or 3D image contains an object to be segmented, for example, a bone.

[0005] STATE OF THE ART

[0006] Image segmentation consists of classifying each pixel or voxel of an image (2D or 3D) as belonging or not to an object, for example a bone.

[0007] Regardless of the application, accurate segmentation is desirable, and accurate segmentation should generally follow the most visible contour in the image.

[0008] For example, in orthopedic surgery, bone segmentation in images can be performed manually by an annotator, using interactive image annotation and processing tools, and / or automatically by segmentation algorithms.

[0009] However, manual methods are inherently subject to human bias and are not very precise, while semi-manual methods are complex and time-consuming. With automated methods, the gains in practicality, speed, and accuracy are offset by reduced robustness. These methods are therefore unsatisfactory for obtaining precise and accurate segmentation of objects, for example, in computed tomography (CT) or cone beam computed tomography (CBCT) images.

[0010] DESCRIPTION OF THE INVENTION

[0011] The invention proposes to segment in a precise and accurate manner one or more object(s) contained in an image, preferably a medical image.

[0012] To this end, the invention proposes a computer-implemented method for segmenting an image, comprising

[0013] - Obtaining a 2D / 3D image I of a region of interest of an element comprising at least one object Aj, j = 1 to J, the 2D / 3D image I comprising N pixels / voxels In, n=1 to N;

[0014] - Initial segmentation of the 2D / 3D image I so as to obtain J image Pj of segmented object Aj, each image Pj comprising N pixels / voxels Pjn, each center of which is a probability of belonging to the object Aj; - Processing of each image Pj, from one or more margins of error of the initial segmentation of the image I so as to obtain a corresponding image llj, the image llj comprising a first zone and a second zone, the first zone comprising only pixels / voxels lljn that belong with certainty to the object, the first zone being totally included in the object Aj; the second zone comprising only pixels / voxels lljn that do not belong with certainty to the object Aj; the pixels / voxels not belonging to the first zone nor to the second belonging to an uncertain zone;

[0015] - Processing of the pixels / voxels of the 2D / 3D image I corresponding to the uncertain area of ​​the image Uj so that the values ​​of the pixels / voxels corresponding to the uncertain area of ​​the image Uj are defined with respect to a boundary delimiting the object Aj in the 2D / 3D image I in the uncertain area, the boundary being the most visible in the uncertain area, the pixels / voxels processed in the uncertain area and the pixels / voxels of the first and second areas defining a mask Fj;

[0016] - Determination of a gap image Ej between the pixels / voxels Fjn of the mask Fj and the pixels / voxels Pjn of the image Pj, the non-nuis pixels / voxels Ejn of the gap image Ej being gap pixels / voxels, these pixels / voxels being grouped into M disjoint groups EGjm named "gap", a gap EGjm being a set of gap pixels / voxels Ejn having the same non-zero value, the gap EGjm and therefore the gap pixels / voxels Ejn of the gap EGjm being classified "significant" if the size and / or depth of the set is greater than an EMIN and / or EDEPTH threshold and "not significant" otherwise;

[0017] - Assignment processing of each pixel / voxel Ejn in deviation to assign it or not to the object Aj and thus obtain a final segmentation mask Mj including a refined segmentation of the object Aj, the assignment processing including an assignment rule defined as follows, for each pixel / voxel n of each significant deviation EGjm, if Ejn= 1 then Mjn=0, the pixel / voxel Ejn not belonging to the object Aj, if Ejn=-1 then Mjn=1, the pixel / voxel Ejn belonging to the object Aj, the remaining pixels Mjn being identical to the corresponding pixels Fjn of the mask Fj.

[0018] The invention is advantageously complemented by the following features, taken alone or in any technically feasible combination thereof:

[0019] - The assignment process includes an assignment rule defined as follows:

[0020] VI jn = 1 when the pixel / voxel with significant deviation Ejn belongs to the object Aj and Mjn = 0 when the pixel / voxel with significant deviation Ejn does not belong to the object Aj, the other remaining pixels / voxels Mjn being identical to the corresponding pixels Fjn of the mask Fj. - the process includes a step of classifying each pixel / voxel n of each significant deviation EGjm (class S) according to one of the following classes: Inside the object Aj (IN) or outside the object Aj (OUT) by means of a previously trained classification model.

[0021] - the classification takes into account at least one subclass assigned to each significant difference, this subclass being used according to the use case of the final segmentation Mj to discriminate differences that can be assigned to one or the other of the classes according to the use case.

[0022] - the assignment rule is further defined as follows

[0023] VI jn = 1 or > 0.5 when the pixel / voxel in deviation Ejn belongs to the object Aj and Mjn=0 or < 0.5 when the pixel / voxel in deviation Ejn does not belong to the object Aj, the other remaining pixels / voxels Mjn being identical to the corresponding pixels Fjn of the mask Fj.

[0024] - The processing of each image Pj to obtain the image Uj (S3) consists, for all pixels / voxels n from 1 to N: if PiNj.n > PHIGH WHERE P IN j, n is the minimum value of the points in the image Pj in a convolution kernel K IN j around the pixel / voxel Pjn of Pj, the pixel / voxel Ujn of the image Uj belongs to the first zone; if For j, n < PLOW j where P O uTj,n is the minimum value of the points in the image Pj in a convolution kernel K 0UT j around the pixel / voxel Pjn of Pj or if P IN k , n > PHIGH k for k=1 to J where kj, the pixel / voxel Ujn of the image Uj belongs to the second zone.

[0025] ■ PHIGH j = PLOW j = 1 / 2 and in which the pixels / voxels of the first zone have a probability equal to 1, the pixels / voxels of the second zone have a probability equal to 0, the pixels / voxels of the uncertain zone have a probability equal to 1 / 2.

[0026] - the most visible boundary is obtained by processing the 2D / 3D image I in the uncertain area using a graph cut algorithm, the most visible boundary Fj being the minimum cut of a graph linking neighboring pixels / voxels in the uncertain area by links having a capacity defined by a function f of the image I and possibly of the initial segmentation Pj, the minimum cut being obtained using an algorithm configured to solve a maximum flow problem.

[0027] - the most visible boundary is obtained by processing the 2D / 3D image I in the uncertain zone to obtain a geodesic curve / surface, the geodesic curve / surface being characterized by a function f using the magnitude of the gradient of the image I, the most visible boundary being that which is the solution to the dual problem of the Combinatorial Continuous Max Flow (CCMF) algorithm, solved for example by the Primal-Dual Interior Point (PDIP) algorithm, or an Appleton-Talbot Continuous Max Flow (AT-CMF) type partial differential equation solving algorithm.

[0028] - The function f uses the magnitude of the gradient of the image I and possibly the magnitude of the gradient of the initial segmentation Pj to discriminate an external boundary of Aj from internal boundary(s) of Aj, for example by using the dot product of the normalized gradient vectors of Pjn and In, for example with 0 < a < 1: f (In, Pjn) = Il Vin

[0029] - The determination of the image Ej (S5) includes a binarization of the mask Fj and the image Pj consisting of obtaining a binary image FBj where each pixel FBjn is equal to 1 if Fjn>0.5 otherwise 0 and a binary image PBj where each pixel PBjn is equal to 1 if Pjn>0.5 otherwise 0, the binary mask FBj and the binary image PBj being binary versions of the mask Fj and the image Pj where a pixel / voxel of the binary mask FBj and / or the binary image PBj is equal to 1 for each pixel / voxel greater than 0.5 or 0 for each pixel / voxel less than 0.5, the image Ej comprising pixels / voxels whose levels are 0, +1 or -1, the uncertain pixels / voxels having a value equal to -1 or +1, the gap image Ej being determined from the binary images FBj and PBj.

[0030] - the image Ej of deviation contains pixels / voxels Ejn= 1 if the pixel In of the image I is included in the binary mask FBj and outside the binary image PBj, Ejn=O if the pixel In of the image I is included in the mask FBj and in the mask PBj or outside the mask FBj and the mask PBj, Ejn=-1 if the pixel In of the image I is outside the binary mask FBj and included in the image PBj.

[0031] - the process includes a step of determining a curve / surface Cj delimiting the object Aj in the final segmentation mask Mj, the curve / surface Cj corresponding to a boundary in the final segmentation mask Mj obtained by means of a SurfaceNets, Marching Squares / Cubes, Cubical Marching Squares or Flying Edges algorithm.

[0032] - the image I is a 3D image reconstructed from a number K of 2D images whose projection parameters are obtained, for example K between 1 and 10, and the boundary Fj is then determined on an image which is a linear combination of I and the K 2D images back-projected into the space of I according to their respective projection parameters, for example a weighted average where I has a weight equal to K and the K images each have a weight of 1.

[0033] DESCRIPTION OF THE FIGURES Other features, purposes and advantages of the invention will become apparent from the following description, which is purely illustrative and not limiting, and which should be read in conjunction with the accompanying drawings on which:

[0034] Figure 1 illustrates steps of a segmentation process according to one embodiment of the invention.

[0035] Figure 2 illustrates an initial 2D image I of a region of interest comprising three objects A1, A2, A3, the image I comprising N pixels In, n=1 to N.

[0036] Figure 3 illustrates an initial segmentation of the 2D image I so as to obtain J image Pj of segmented object Aj, each image Pj comprising N pixels Pjn, each value of which is a probability of belonging to the object Aj.

[0037] Figure 4 illustrates a processing of each image Pj, from one or more margins of error of the initial segmentation of the image I for an image comprising a first zone Z1 and a second zone Z2, the first zone Z1 comprising only pixels which are reliably included in the object, the first zone Z1 being totally included in the object; the second zone Z2 comprising only pixels which are not reliably included in the object; the pixels which are neither in the first zone nor in the second belong to an uncertain zone Z3.

[0038] Figure 5 illustrates a processing of the pixels of the initial 2D image I corresponding to the uncertain area of ​​the image such that the values ​​of the pixels corresponding to the uncertain area of ​​the image Uj are defined with respect to the most visible boundary delimiting the object Aj in the 2D / 3D image I in the uncertain area, the pixels processed in the uncertain area and the pixels of the first and second area defining a mask Fj of the most visible boundary and a mask Vj of the relative visibility of this boundary in each pixel / voxel.

[0039] Figure 6 illustrates a determination of an image Ej of difference between the pixels Fjn of the mask Fj and the pixels Pjn of the image Pj.

[0040] Figure 7 illustrates a treatment of significant differences in class(es) or even subclass^).

[0041] Figure 8 illustrates a processing of assigning each pixel in deviation to assign it or not to the object Aj and thus obtain a final segmentation mask Mj including a refined segmentation of the object Aj.

[0042] Across all figures, similar elements relate to similar elements.

[0043] DETAILED DESCRIPTION OF THE INVENTION

[0044] General overview of the process

[0045] The method according to the invention comprises the main steps below enabling the segmentation of objects in 2D or 3D images to be accurate to within one pixel / voxel, or even less, and following the objectively most visible boundary in the initial image if and only if it is determined to be precisely that of the object:

[0046] Step S1: Obtain an initial image I 2D or 3D;

[0047] Step S2: Obtain an initial segmentation of one or more objects, for example semi-automatically or with a machine learning trained model;

[0048] Step S3: Automatically determine an uncertain zone relative to each segmented object, for example according to distance and / or probability errors characteristic of the initial segmentation and determined beforehand;

[0049] The steps are implemented for each segmented object.

[0050] Step S4: Determine a segmentation, following the most visible boundary in image I separating the object from the rest of the image in the uncertain area;

[0051] Step S5: Determine and locate groups of connected pixels / voxels in deviation between the initial segmentation (S2) and the most visible boundary (S4), these deviations being located in the uncertain zone (S3), each deviation being classified according to its dimensions as "significant" or "non-significant";

[0052] Step S6 (optional): Classify each significant deviation according to other characteristics that may be a function of the initial image I;

[0053] Step S7: Decide whether to include or exclude each deviation from the object based on this or these classifications;

[0054] Step S8 (optional): Obtain a curve / surface following a contour of the segmented object in sub-pixel / voxel precision.

[0055] Initial image I (step S1)

[0056] The method includes a step S1 of obtaining an initial 2D or 3D image, denoted I, of a region of interest comprising J object(s) Aj, j=1 to J. An object is preferably bone, cartilage, or ligament. However, the method described here applies to any type of object with visible contours in image I. This image I is acquired using a known type of imaging system, depending on the type of initial image desired, and will not be described further.

[0057] The initial image I partitions the 2D or 3D physical space where the objects Aj are located into pixels or voxels, each representing one or more real numbers. These numbers are influenced by the nature and proportion of the objects Aj present in each pixel or voxel, for example, by X-ray absorption, atomic spin relaxation, wave reverberation, etc., depending on the acquisition modality chosen. The initial image I is, for example, a radiography, fluoroscopy, MRI, CT, or CBCT image.

[0058] In addition, in the case of a 3D image, the initial 3D image can be reconstructed from a small number of smaller dimension images: a 3D image reconstructed with a reduced number of ultrasound slices or; a 3D X-ray absorption image reconstructed with a reduced number of 2D radiographs.

[0059] For example, the number K of images is between 1 and 10.

[0060] This initial image I therefore comprises N pixels / voxels, denoted In = u n with n=1 to N, u n being the value of the pixel / voxel.

[0061] Figure 2 illustrates a 2D image I of N=16x16=256 pixels whose values ​​are influenced mainly by the nature and proportion of each object in each voxel, for example l n e [0, 255] for all ne [1, N]. The proportion is visible in Figure 1 with respect to the gray level. The darker the pixel, the more the object Aj is present proportionally. Conversely, the lighter the pixel, the less the object Aj is present proportionally. When an object fully occupies a pixel / voxel, the pixel / voxel value depends only on the nature of that object and the chosen imaging system.

[0062] Initial segmentation Pi of an object (step 52)

[0063] In a step S2, the initial image I is processed according to at least one initial segmentation so as to obtain, for each initial segmentation, J image(s) Pj of segmented object Aj. At the end of this or these initial segmentation(s), we obtain J*S image(s) Pj, each comprising N pixels / voxels, where each value is a probability of belonging to object Aj, Pj,s denoting the image Pj corresponding to a segmentation s performed, S being the number of different segmentations performed, and s = 1 to S. The set of images Pj of segmented object Aj is called the first set.

[0064] This initial segmentation can indeed be implemented in several ways and must be robust, meaning it must avoid gross errors. However, its accuracy is limited. Such an initial segmentation thus exhibits one or more error margins characteristic of the initial segmentation method in pixel(s) / voxel(s) value and / or in distance; for example, a method might have maximum distance errors of ±2 pixels / voxels, or maximum probability errors of ±0.1, or a combination of both, allowing for a better characterization of the maximum errors and resulting in smaller, average-sized uncertain areas at step S3.

[0065] The initial segmentation can be implemented manually, semi-automatically, or using a neural network trained in image segmentation.

[0066] Performing several initial segmentations also allows us to test several algorithms, different parameters, etc., and to further improve the quality of the result as we will see later.

[0067] Alternatively or in a complementary manner, each initial segmentation can combine several types of segmentations, for example by averaging the annotations of several methods or annotators, possibly weighted by their performance.

[0068] The annotations may be incomplete or even only point-based, and the segmentation obtained by considering for each pixel / voxel In, is for example a probability of belonging to the object of the nearest annotated point, possibly based on distance.

[0069] The initial segmentation can be binary and made non-binary by applying, for example, a convolution operation such as Gaussian blur.

[0070] Figure 3 illustrates an example of an initial segmentation: on the left the initial image I comprising three segmented objects A1, A2, A3 and on the right the different resulting images Pj P1, P2, P3 comprising pixels / voxels at 1 or 0 depending on whether they belong to the object Aj or not.

[0071] Partition Uj of the initial segmentation Pj into 3 zones, one of which is uncertain (step S3)

[0072] During an S3 step, each image Pj is processed, starting from one or more error margins of the initial segmentation of the image I, in order to obtain at least one corresponding image Uj comprising N pixels / voxels.

[0073] We understand that each image Pj of the first set (i.e. obtained for each segmentation) is processed.

[0074] Thus again, in the case of multiple initial segmentations, step S3 is advantageously repeated for each segmentation, so as to process each image Pj,s, from one or more error margins of the s-th initial segmentation of image I in order to obtain one or more corresponding image(s) Uj,s comprising N pixels / voxels and corresponding to an image Pj,s. The set of images Uj obtained is called the second set.

[0075] The objective of this step S3 is to partition each image Pj by defining the following areas, taking into account one or more possible errors in the initial segmentation:

[0076] A first zone comprising only pixels / voxels Ujn of Uj that reliably belong to the object Aj;

[0077] A second area comprising only pixels / voxels Ujn of Uj that do not reliably belong to the object Aj;

[0078] An uncertain area comprising pixels / voxels that belong neither to the first area nor to the second area.

[0079] The pixel / voxel values ​​of the first zone are preferably set to 1, those of the second zone to 0, and those of the uncertain zone to 1 / 2. According to one embodiment, to partition the image Pj automatically, the following procedure is used: For all n in [1 .. N]:

[0080] • lljn = 1 if:

[0081] O PiN jn > PHIGH J (typically > 1 / 2)

[0082] ■ where PIN jn is the minimum value of the points of Pj in a convolution kernel KIN j around Pjn

[0083] • lljn = 0 if: o POUT jn < PLOW j (typically < 1 / )

[0084] ■ where POUT jn is the maximum value of the points of Pj in a convolution kernel KouTj around Pjn o PIN kn > PHIGH k

[0085] ■ for all k in [1 .. J] where kj

[0086] ■ where PIN kn is the minimum value of the points of Pk in a convolution kernel KIN k around Pkn

[0087] • Otherwise Ujn = 1 / 2.

[0088] In other words, the partition Uj of Pj includes a first zone totally included in the object Aj (erosion) and the second zone is a zone which is beyond the object (dilation).

[0089] We thus obtain the image Uj where the pixels Ujn for all n in [1 .. N] have a probability among {1, 1 / 2, 0} for each voxel center In of belonging to Aj, given the characteristic errors of the initial segmentation and the other segmented objects. Note that these characteristic errors are determined before the initial segmentation of the image I by characterizing the errors of the initial segmentation method.

[0090] Alternatively, the uncertain zone can also be obtained by simple thresholding operations and / or discrete mathematical morphology, i.e. in the domain {0,1} or even manually.

[0091] Figure 4 illustrates step S3. By considering the probability of 1D points belonging to an object around x=0 according to a Gaussian distribution, step S3 partitions the points according to different neighborhoods inside (first zone Z1) and outside (second zone Z2), and even according to different directions of the 1D image (KIN and KOUT are anisotropic), as well as according to different probability thresholds to determine an uncertain zone that best accounts for a large number of possible characteristic errors for the initial segmentation method. The uncertain zone Z3 is the complementary zone of the first and second zones Z1 and Z2.

[0092] According to one embodiment, and to allow the next step to determine the most visible boundary of the object in an uncertain area that completely separates the first and second zones Z1, Z2, it is preferable to choose convolution kernels KIN x, KOUT X of size strictly greater than 1 pixel in all dimensions of I and a threshold

[0093] PLOW X - HIGH X-

[0094] In addition, it is possible to process the non-connected components of Z1, Z2 in the image llj by integrating them into the uncertain area to allow the next step to determine whether or not there is a boundary connecting them in the image I.

[0095] In all cases, for each image Pj in the first set, we can optionally obtain several corresponding images llj by implementing step S3 multiple times, notably by varying the processing parameters of the step, known as uncertainty zone parameters, such as error margins or thresholds (in particular PHIGH j and PLOW j), so as to obtain larger or smaller uncertainty zones. Thus, the second set can include more J*S images llj.

[0096] Furthermore, it is quite possible to obtain a variable number of llj images for each initial segmentation, some segmentations allowing or requiring more variability in said uncertain zone parameters than others.

[0097] Mathematically, denoting n s the number of images Uj,s obtained for the image Pj,s (i.e., for the s-th segmentation) at step S3 (with n s ^1), for convenience we set S' = Es=i«s > 5 the total number of images llj for each object Aj, and therefore step S3 sees the obtaining of J*S' images llj (in the second set).

[0098] Most visible boundary Fi in the uncertain zone of llj (stage S4)

[0099] In step S4, the uncertain area determined in step S3 is used to detect the most visible boundary in the initial image I within this uncertain area. The most visible boundary separates the pixels / voxels belonging to object Aj from the others within the uncertain area.

[0100] In particular, the pixels / voxels corresponding to those in the uncertain region identified in the image Uj are associated with a function f measuring the visibility of the object boundary between neighboring pixels / voxels, considering the intensities of the pixels / voxels in the image I and possibly the initial segmentation Pj. This function is used to determine the boundary that is globally the most visible, that is, the one that globally optimizes this function. This boundary is expressed as a segmentation mask denoted Fj.

[0101] We understand that each image Uj of the second set (i.e. obtained for each initial segmentation and / or each different case of uncertain zone parameters) is processed.

[0102] Thus, once again, in the case of multiple images Uj, step S4 is advantageously repeated so as to process each image Uj,s in order to obtain a corresponding mask Fj,s. The set of masks Fj obtained is called the third set.

[0103] In one embodiment, a graph cut algorithm is used, the most visible boundary being the minimum cut of a discrete graph connecting neighboring pixels / voxels in the uncertain region by links with a maximum capacity defined by the function f. The minimum cut is obtained using an algorithm to solve the dual problem of maximum flux passing through the links of the graph (Graph Cut, GC). In this regard, see the document Krcah Marcel, Gâbor Székely, and Rémi Blanc, "Fully automatic and fast segmentation of the femur bone from 3D-CT images with no shape prior," 2011 IEEE International Symposium on Biomedical Imaging: From Nano to Macro, IEEE, 2011.

[0104] In one embodiment, a method for approximating a continuous geodesic curve / surface is used, the geodesic curve / surface being characterized by a metric, the metric being an antimonotonic function of the function f, the most visible boundary being the approximate solution of this geodesic curve / surface, the approximate solution being obtained by means of a known algorithm for solving the continuous maximum flow (CMF) dual problem, for example the Primal-Dual Interior Point (PDIP) algorithm for solving the combinatorial formulation of the problem (Combinatorial Continuous Max Flow, CCMF), or a partial differential equations algorithm for solving a physical formulation of the problem, for example Appleton-Talbot Continuous Max Flow (AT-CMF). For these algorithms, reference can be made to the following documents:

[0105] Couprie, Camille & Grady, Leo & Talbot, Hugues & Najman, Laurent. (2010). Combinatorial Continuous Maximal Flows. CoRR. abs / 1010.2733;

[0106] Appleton, Ben & Talbot, Hugues. (2003). Globally Optimal Surfaces By Continuous Maximal Flows.

[0107] Determining the most visible boundary involves one or more functions f(In, Pjn) that increase as the boundary of Aj becomes more visible in image I. To achieve the greatest possible accuracy in the final segmentation, the function f(In, Pjn) should consider the most local features possible; that is, for each In, it should only consider the values ​​of In and / or its neighbors in image I. For example, a norm of a discretized gradient of image I that does not introduce a shift, such as the central difference of In. When image I contains multiple values, such as a color or, more generally, a multispectral image, the visibility of the boundary may depend on one or more of these values.

[0108] In addition, it may be useful to discriminate between internal and external boundaries of the object Aj. To do this, a directional component is introduced into the metric to ensure that, for example, the external boundary of the cortical bone (decreasing intensities from the inside out) is segmented instead of its internal boundary (increasing intensities from the inside out).

[0109] In this case, the function f uses the vector gradient of the image I and the vector gradient of the initial segmentation Pj to discriminate an external boundary of Aj from internal boundary(s) of Aj, for example by using the dot product of the normalized gradient vectors of Pjn and In, for example with 0 < a < 1: f(In, Pjn) = IlVin

[0110] Such a function f is expressed as a function of I and Pj, to be higher on pixels / voxels where the gradient vector of Pjn is in a direction opposite to the gradient vector of In, thus making these pixels / voxels less "cutting" than pixels with an equally large gradient of In but oriented differently.

[0111] At the end of this step S4, for each image Uj of the second set a mask Fj is obtained which includes the pixels / voxels of the first zone and second zone and in the uncertain zone the values ​​adjusted with respect to the most visible boundary.

[0112] Additionally, when the PDIP algorithm is used to obtain the most visible boundary, it provides a solution v(nu) at each iteration that bounds the desired most visible boundary and generally converges rapidly towards it. The PDIP algorithm also provides a variable A(lambda) quantifying the relative contribution of each pixel / voxel to the visibility of the desired boundary.

[0113] In this case, at the end of step S4, Vj = A local indicator of border visibility Fj that is strictly positive and increasing with the local visibility of the border at the level of In is obtained, Vj can advantageously be normalized between 0 and 1.

[0114] Figure 5 illustrates in 2D the image I with the object Ai in particular, the corresponding image Ui, the resulting mask Fi, and the local visibility indicator Vi. The mask Fi is obtained by approximating a geodesic; the pixels of Fi corresponding to the uncertain region Z3 of U1 receive a value close to 0 (outside Ai) or 1 (inside Ai), except for the pixels on either side of the most visible boundary in the image I, which receive a value on either side of V2. The threshold of V2 allows the pixels on either side of this boundary to be segmented. The mask V1 is obtained by the PDIP algorithm; the pixels have a value close to 0, except for the pixels on either side of the most visible boundary in the image I, which receive a value that is higher the more locally visible the boundary is.

[0115] We observe that from image I and image U1, in the uncertain area Z3, the mask F1 segments the pixels on either side of the most visible border in the initial image I.

[0116] Gap image, grouping and classification by dimension (step S5)

[0117] In step S5, at least one slack image Ej is determined between the pixels / voxels Fjn of a mask Fj and the pixels / voxels Pjn of the image Pj, the unaffected pixels / voxels Ejn of the slack image Ej being slack pixels / voxels. It is understood that, at the end of step S4, for each object Aj, at least one mask Fj has been obtained, and advantageously a plurality of masks Fj corresponding to several initial segmentations and / or various uncertain zone parameters.

[0118] In the case where we have a unique mask Fj for each object Aj (i.e. the third set contains J masks Fj) we can of course only calculate the difference between the pixels / voxels Fjn of the mask Fj and the pixels / voxels Pjn of the image Pj, i.e. determine a difference image Ej for each object Aj.

[0119] If we have at least two masks Fj per object Aj (for example if two different initial segmentations were done in step S2), we can repeat step S5, and for each mask Fj of the third set, determine a gap image Ej between the pixels / voxels Fjn of the mask Fj and the pixels / voxels Pjn of the image Pj (i.e. up to J*S' gap images Ej), so as to exploit the performance of each segmentation / parameterization.

[0120] It is particularly interesting to have several Fj masks corresponding to several initial segmentations and / or various uncertain zone parameters, which makes it possible to be sure of detecting any segmentation error even if it is made both by the initial segmentation and by the search, in different uncertain zones, for the most visible boundary(s) in the image.

[0121] The idea here is to compare the pixels / voxels of the mask(s) Fj with those of the image Pj to determine the differences that should be classified as belonging or not belonging to the object Aj.

[0122] Preferably, this gap image Ej is obtained from binary versions of the mask(s) Fj and the image Pj.

[0123] To do this, the binary version FBj of Fj and the binary version PBj of Pj are obtained according to the following rule:

[0124] For all n in [1, ..., N]:

[0125] - FBjn = 1 if Fjn > 0.5 otherwise 0;

[0126] - PBjn = 1 if Pjn > 0.5 otherwise 0.

[0127] Then each image of the gap between the initial segmentation and the most visible boundary of Aj in the uncertain zone Ej = FBj - PBj is determined where

[0128] Ejn = 1 if In is in FBj (=1) but outside PBj (=0)

[0129] Ejn = 0 if In is in FBj (=1) and PBj (=1), or outside FBj (=0) and PBj (=0)

[0130] Ejn = -1 if In is outside FBj (=0) but in PBj (=1)

[0131] We then obtain, by a known connected component analysis, the set EGj of coordinates, and optionally other features such as bounding boxes, M disjoint groups of connected pixels / voxels EGjm of value 1 and -1 called deviations, the deviation EGjm being classified "significant" (class S) if the size and / or depth of the set of these pixels / voxels is greater than an EMIN and / or EDEPTH threshold and "not significant" (class NS) otherwise.

[0132] The pixels / voxels of the non-significant deviations (class NS) will be assigned at step (S7) so as to follow the most visible boundary to correct the lack of precision of the initial segmentation.

[0133] We note that the number of gap images Ej obtained can be quite high, up to J*S', that is, the maximum number of gaps between the pixels / voxels Fjn of a mask Fj and the pixels / voxels Pjn of the image Pj.

[0134] If more than one gap image Ej is obtained for an object Aj, various methods can be used, preferably combining these gap images Ej into a single overall gap image denoted EZj, notably using a majority system: for a pixel / voxel n, the number of times Ejn = 1 and the number of times Ejn = 0 are counted, and EZjn represents the majority value. Alternatively, one could, for example, average the values ​​and binarize them again, or use any other known technique.

[0135] Next, we can group the non-nuisance pixels / voxels EZjn of the global deviation image EZj into deviations noted EZGjm, classify them as significant or not and continue the process on the basis of the global deviation image EZj and its deviations EZGjm as if it were the only deviation image Ej.

[0136] Alternatively, we can still determine the EGjm gaps for each gap image Ej independently, and combine these gaps directly into EZGjm, again by majority or otherwise, for example by identifying the "groups" of gaps corresponding in practice to the same region of interest as sharing a common non-empty intersection.

[0137] Alternatively, we can keep the plurality of gap images Ej and all gaps EGjm, see further.

[0138] Optional classification of significant differences (step S6)

[0139] When the discrepancies are significant (class S), it may be advantageous to interpret them to determine if the most visible boundary is indeed that of the object Aj and further assign them to one of the following classes: Inside the object Aj (class IN) or outside the object Aj (class OUT).

[0140] Such an interpretation is implemented using a classification model, for example, a neural network trained by deep learning, taking as input a region of interest from I centered on the coordinates or bounding box of a significant deviation (class S) EGjm, and the same region of interest from Ej. In the case of multiple images of deviations Ej, the model can take as input either a global significant deviation (after combination) EZGjm, or a vector of significant deviations EGjm with a non-empty common intersection. In all cases, the classification performed for a deviation EGjm / EZGjm can be applied at step (S7) to all pixels / voxels of said deviation, or just a potentially very small subset of its pixels / voxels located at the center of said deviation, for example, the smallest non-empty set of pixels / voxels after at least one erosion step in at least one dimension of the image.This cleverly avoids small segmentation artifacts at the boundary between 2 differently classified EGjm gaps.

[0141] Furthermore, it can be useful for surgical planning to assign a subclass that can explain the segmentation gap or influence the surgery, for example, for a bone-type object Aj in a CT image:

[0142] • belonging to a neighboring bone (e.g., ribs, ossicles, etc.) or a bone growth not relevant to surgery (e.g., an osteophyte, etc.);

[0143] • an artifact causing values ​​possibly lower than the lower limit of 700 Hll (for example in English, shadowing, beam hardening, ring artifacts, etc);

[0144] • a condition explaining a local under-intensity (e.g. osteoporosis, tumor, geode, etc.);

[0145] • an artifact causing values ​​possibly exceeding the upper limit of 3000 Hll for bones (e.g., streak, ring artifacts, etc.);

[0146] • calcification of another nearby tissue (e.g., ligament, artery, etc.);

[0147] • the presence of hardware (e.g., pedicle screw, etc.);

[0148] Certain subclasses such as "geode" or "belonging to a neighboring bone" can, on their own, justify assigning a deviation to one of the IN or OUT classes.

[0149] Additionally, to assign a significant deviation (class S) to either the IN or OUT classes, it can be helpful to also consider the use case. For example, the subclass "osteophyte" might be assigned to the OUT class for final segmentation intended to calculate the size of an implant, but to the IN class for final segmentation intended to determine the area of ​​bone that a robot should cut.

[0150] Thus, depending on the subclass and optionally the use case, the difference can be integrated in the most appropriate way into the IN or OUT class.

[0151] At the end of this step S6, we obtain ECjm, the set of class(es) and possibly subclass(es). Figure 6 illustrates in 2D a binarized segmentation PBi and a binarized boundary FBi, the difference of which reveals the image of the gaps Ei in three connected components EGjm, one of which is insignificant in size (comprising a single pixel).

[0152] Figure 7 illustrates a 2D image I containing two significant deviations classified by a classification model. The significant deviations correspond to a light artifact and an osteophyte. These must then be classified as belonging to the object or not.

[0153] Assignment of each pixel / voxel to the Ejn gap (step S7)

[0154] In a step S7 each pixel / voxel in gap Ejn (see steps S5 and S6) is then assigned to the final segmentation mask Mj comprising a refined segmentation of the object Aj.

[0155] More specifically, it involves processing each pixel / voxel n of each significant gap EGjm of at least one gap image Ej (including the possible global gap image EZj), to assign it or not to the object Aj and thus obtain a final segmentation mask Mj comprising a refined segmentation of the object Aj (also called final segmentation, as opposed to the initial segmentation(s)).

[0156] If we have a single gap image Ej, this one is used. If we have several gap images Ej, we can either use a global gap image EZj as mentioned before (and therefore process each pixel / voxel n of each significant global gap EZGjm of the global gap image EZj, to assign it or not to the object Aj), or again use all the gaps EGjm of the various Ej.

[0157] We will not be limited in any way to obtain the final segmentation mask Mj from the gap image(s) Ej.

[0158] According to a first embodiment, the assignment is implemented according to the assignment rule defined as follows, for each pixel / voxel n of each significant gap EGjm (class S) of at least a given number of gap images Ej (for example at least 80% of the gap images) or of a combination of gap images Ej (the global gap image EZj), if Ejn=1 then Mjn=0, the pixel / voxel Ejn not belonging to the object Aj, if Ejn=-1 then Mjn=1, the pixel / voxel Ejnp belonging to the object Aj, the remaining pixels Mjn being identical to the corresponding pixels Fjn of the mask Fj to precisely follow the most visible boundary in the absence of gap or when the gaps are not significant (class NS).

[0159] According to a second embodiment, the pixels / voxels with significant deviations (class S) are more finely assigned by taking into account their class determined in step S6 (classes IN or OUT).

[0160] The two modes of implementation can be summarized as follows: either Mj = Fj, then for each gap EGjm of EGj and for each pixel / voxel n of EGjm located in Ej, Fj, Pj and Mj, the assignment is implemented in the following way.

[0161] Thus, according to a first embodiment (only the first two lines of the table above), the pixels / voxels Ejn in deviation are assigned or not to the mask Mj according to their class NS or S and their value +1, -1 only.

[0162] According to a second embodiment, when the pixel / voxel in deviation Ejn is significant (class S) then it is interpreted by a classification model that can take into account the initial image and the use case and its resulting IN or OUT class is taken into account, those that are not significant (class NS) being classified only with respect to their value as in the first embodiment.

[0163] Figure 8 illustrates in 2D an initial segmentation PBi towards the most visible boundary FBi shown in image I to result in the final segmentation Mi.

[0164] As already explained, the classification of an EGjm gap can alternatively be assigned not to the set of pixels / voxels in gap but to a subset of these pixels / voxels, advantageously just a potentially very small subset located at the center of the gap, for example the smallest non-empty set of pixels / voxels after at least one erosion step in at least one dimension of the image.

[0165] In these embodiments, the final segmentation Mj of each object can be obtained similarly to step S4 by determining the best visible boundary of the object Aj, taking into account the reclassification of a reduced number of disjoint pixels / voxels. The technical effect is to provide a semantically more accurate segmentation precisely aligned to the best visible boundary in the image.

[0166] Curve / surface (step S8)

[0167] Additionally, in step S8, a curve / surface Cj delimiting the object Aj in the final image Mj is determined. This is the curve / surface Cj corresponding to a boundary of the image Mj obtained using a known algorithm, for example SurfaceNets, Marching Squares / Cubes, Cubical Marching Squares or Flying Edges.

[0168] Indeed, the geodesic defines the most visible single boundary of Aj, but the accuracy of this geodesic is limited by the resolution of the pixels / voxels where the geodesic metric is calculated. In particular, the visibility indicator given by the CCMF method brackets this geodesic with values ​​that are attracted to the pixel / voxel centers and their midpoints.

[0169] This curve / surface Cj allows us to obtain a contour of the object Aj with an accuracy of the order of ±0.3 pixel / voxel.

Claims

DEMANDS 1. A computer-implemented method for segmenting an image, comprising (S1 ) Obtaining a 2D / 3D image I of a region of interest of an element comprising at least one object Aj, j = 1 to J, the 2D / 3D image I comprising N pixels / voxels In, n=1 to N; (52) At least one initial segmentation of the 2D / 3D image I so as to obtain, for each initial segmentation, J images Pj of segmented object Aj, each image Pj comprising N pixels / voxels Pjn, each center of which is a probability of belonging to the object Aj; (53) Processing each image Pj, from one or more margins of error of the corresponding initial segmentation of the image I so as to obtain at least one corresponding image llj, the image llj comprising a first zone and a second zone, the first zone comprising only pixels / voxels lljn which belong reliably to the object, the first zone being totally included in the object Aj; the second zone comprising only pixels / voxels lljn which do not belong reliably to the object Aj; the pixels / voxels not belonging to the first zone or the second belonging to an uncertain zone; (54) Processing the pixels / voxels of the 2D / 3D image I corresponding to the uncertain area of ​​each image Uj such that the values ​​of the pixels / voxels corresponding to the uncertain area of ​​each image Uj are defined with respect to a boundary delimiting the object Aj in the 2D / 3D image I in the uncertain area, the boundary being the most visible in the uncertain area, the pixels / voxels processed in the uncertain area and the pixels / voxels of the first and second areas defining at least one mask Fj; characterized in that it further comprises: (55) Determination of at least one gap image Ej between the pixels / voxels Fjn of a mask Fj and the pixels / voxels Pjn of the image Pj, the non-nuis pixels / voxels Ejn of the gap image Ej being gap pixels / voxels, these pixels / voxels being grouped into M disjoint groups EGjm called "gap", a gap EGjm being a set of gap pixels / voxels Ejn having the same non-zero value, the gap EGjm and therefore the gap pixels / voxels Ejn of the gap EGjm being classified as "significant" if the size and / or depth of the set is greater than an EMIN and / or EDEPTH threshold and "not significant" otherwise; (S7) Assignment processing of each pixel / voxel n of each significant gap EGjm of at least one gap image Ej to assign it or not to the object Aj and thus obtain a final segmentation mask Mj including a refined segmentation of the object Aj.

2. A method according to claim 1, comprising a step (S6) of classifying each pixel / voxel n of each significant EGjm deviation (class S) according to one of the following classes: Inside object Aj (IN) or outside object Aj (OUT) by means of a previously trained classification model.

3. Method according to claim 2, wherein the classification takes into account at least one subclass assigned to each significant deviation, this subclass being used according to the use case of the final segmentation Mj to discriminate deviations that can be assigned to one or the other of the classes according to the use case.

4. A method according to any one of the preceding claims, wherein the assignment processing comprises an assignment rule defined as follows, for each pixel / voxel Ejn of each significant gap EGjm of at least a given number of gap images Ej or of a combination of gap images Ej, if Ejn=1 then Mjn=0, the pixel / voxel Ejn not belonging to the object Aj, if Ejn=-1 then Mjn=1, the pixel / voxel Ejn belonging to the object Aj, the remaining pixels Mjn being identical to the corresponding pixels Fjn of the mask Fj.

5. A method according to any one of the preceding claims, wherein the processing of each image Pj to obtain an image Uj (S3) consists, for all pixels / voxels n from 1 to N: if PiNj.n > PHIGH WHERE P IN j, n is the minimum value of the points in the image Pj in a convolution kernel K IN j around the pixel / voxel Pjn of Pj, the pixel / voxel Ujn of the image Uj belongs to the first zone; if PouTj.n < PLOW j where PO uTj,n is the minimum value of the points in the image Pj in a convolution kernel K 0UT J around the pixel / voxel Pjn of Pj or if P IN k , n > PHIGH R. P° ur k=1 to J where kj, the pixel / voxel Ujn of the image Uj belongs to the second zone.

6. A method according to claim 5, wherein P H IGH J = l / 2 PLOW J = V 2 and in which the pixels / voxels of the first zone have a probability equal to 1, the pixels / voxels of the second zone have a probability equal to 0, the pixels / voxels of the uncertain zone have a probability equal to 1 / 2.

7. A method according to any one of the preceding claims, wherein the most visible boundary (S4) is obtained by processing the 2D / 3D image I in the uncertain area using a graph cutting algorithm, the most visible boundary Fj being the minimum cut of a graph linking neighboring pixels / voxels in the uncertain area by links having a capacity defined by a function f of the image I and optionally of the initial segmentation Pj, the minimum cut being obtained by means of an algorithm configured to solve a maximum flow problem.

8. A method according to any one of the preceding claims, wherein the most visible boundary (S4) is obtained by processing the 2D / 3D image I in the uncertain area to obtain a geodesic curve / surface, the geodesic curve / surface being characterized by a function f using the magnitude of the gradient of the image I, the most visible boundary being that which is the solution to the dual problem of the Combinatorial Continuous Max Flow (CCMF) algorithm, solved for example by the Primal-Dual Interior Point (PDIP) algorithm, or an Appleton-Talbot Continuous Max Flow (AT-CMF) type partial differential equation solving algorithm.

9. A method according to any one of claims 7 to 8, wherein the function f uses the magnitude of the gradient of the image I and possibly the magnitude of the gradient of the initial segmentation Pj to discriminate an external boundary of Aj from internal boundary(s) of Aj, for example by using the dot product of the normalized gradient vectors of Pjn and In, for example with 0 < a < 1: f(In, Pjn) = IlVin 10. A method according to any one of the preceding claims, wherein the determination of each gap image Ej (S5) comprises a binarization of the mask(s) Fj and the image Pj consisting of obtaining binary images FBj of which each pixel / voxel FBjn is equal to 1 if Fjn>0.5 otherwise 0 and a binary image PBj of which each pixel / voxel PBjn is equal to 1 if Pjn>0.5 otherwise 0, the binary masks being binary versions of the mask(s) Fj and the image Pj, the gap image Ej comprising pixels / voxels whose levels are 0, +1 or -1, the uncertain pixels / voxels having a value equal to -1 or +1, the gap image Ej being determined from the binary images FBj and PBj.

11. Method according to claim 10, wherein the gap image Ej comprises pixels / voxels Ejn=1 if the pixel / voxel In of the image I is included in the binary mask FBj and outside the binary mask PBj, Ejn=0 if the pixel / voxel In of the image I is included in the mask FBj and in the binary mask PBj or outside the binary mask FBj and the binary mask PBj, Ejn=-1 if the pixel In of the image I is outside the binary mask FBj and included in the image PBj.

12. A method according to any one of claims 10 and 11 in combination with claim 4, wherein the assignment rule for step S7 is further defined as follows: VI jn = 1 or > 0.5 when the pixel / voxel in deviation Ejn belongs to the object Aj and Mjn=0 or < 0.5 when the pixel / voxel in deviation Ejn does not belong to the object Aj, the other remaining pixels / voxels Mjn being identical to the corresponding pixels Fjn of the mask Fj.

13. A method according to any one of the preceding claims, comprising a step (S8) of determining a curve / surface Cj delimiting the object Aj in the final segmentation mask Mj, the curve / surface Cj corresponding to a boundary in the final segmentation mask Mj obtained by means of a SurfaceNets, Marching Squares / Cubes, Cubical Marching Squares or Flying Edges algorithm.

14. A method according to any one of the preceding claims, wherein the image I (S1) is a 3D image reconstructed from a number K of 2D images whose projection parameters are obtained, for example K between 1 and 10, and the boundary (S4) Fj is then determined on an image which is a linear combination of I and the K 2D images back-projected into the space of I according to their respective projection parameters, for example a weighted average where I has a weight equal to K and the K images each have a weight of 1.

15. A method according to any one of the preceding claims, wherein a plurality of images Ej is determined in step (S5) for each object Aj, step (S5) further comprising combining the gap images Ej into an overall gap image EZj, step (S7) being implemented on the basis of each significant gap EZGjm of said overall gap image EZj.

Images